6.5. TDNNS method for nonlinear elasticity#
For nonlinear, invertible material laws the TDNNS method for linear elasticity can not be directly extended to the nonlinear regime. The main problem is that the gradient of the displacement field
We use the Hu-Washizu principle following [Neunteufel, Pechstein, Schöberl. Three-field mixed finite element methods for nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, (2021).] introducing a new independent strain field. We “lift” the distribution
For a hyperelastic potential
with the first Piola-Kirchhoff stress tensor
As
where
Therefore, the TDNNS method for nonlinear elasticity reads: Find
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
import ipywidgets as widgets
def GenerateMesh(t=0.1):
bar = MoveTo(0, -t / 2).Rectangle(1, t).Face()
bar.edges.Min(X).name = "left"
bar.edges.Max(X).name = "right"
bar.edges.Min(Y).name = "bottom"
bar.edges.Max(Y).name = "top"
bar.edges.Min(X).maxh = t / 8
mesh = Mesh(OCCGeometry(bar, dim=2).GenerateMesh(maxh=t / 2.5))
return mesh
mesh = GenerateMesh()
Draw(mesh)
E, nu = 210, 0.2
mu = E / 2 / (1 + nu)
lam = E * nu / ((1 + nu) * (1 - 2 * nu))
def NeoHooke(F):
C = F.trans * F
return mu / 2 * (Trace(C - Id(2)) - log(Det(C))) + lam / 2 * (Det(F) - 1) ** 2
def SkewT(v):
if v.dim == 3:
return 0.5 * CoefficientFunction(
(0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0), dims=(3, 3)
)
else:
return 0.5 * CoefficientFunction((0, -v, v, 0), dims=(2, 2))
order = 3
par = Parameter(0)
force = CF((0, -5 * par))
# clamped boundary at left, free elsewhere
V = HCurl(mesh, order=order, dirichlet="left")
P_sym = HDivDiv(mesh, order=order, dirichlet="bottom|right|top")
P_skw = L2(mesh, order=order)
F = MatrixValued(L2(mesh, order=order))
fes = V * P_sym * P_skw * F
gf_solution = GridFunction(fes)
u, Psym, Pskw, F = fes.TrialFunction()
P = Psym + CoefficientFunction((0, Pskw, -Pskw, 0), dims=(2, 2))
n = specialcf.normal(2)
a = BilinearForm(fes, symmetric=True, condense=True)
a += Variation((NeoHooke(F)) * dx).Compile()
a += Variation(InnerProduct(F - Grad(u) - Id(2), P) * dx).Compile()
a += Variation(u * n * Psym[n, n] * dx(element_boundary=True)).Compile()
a += Variation(-(force * u.Trace()) * ds("right")).Compile()
gf_u, gf_Psym, gf_Pskw, gf_F = gf_solution.components
# F=I
gf_F.Set(Id(2))
# set boundary data
gf_Psym.Set((force * n) * OuterProduct(n, n), definedon=mesh.Boundaries("right"))
scene = Draw(
gf_u,
mesh,
deformation=True,
settings={
"camera": {
"transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -2.5}]
}
},
)
tw = widgets.Text(value="step = 0")
display(tw)
gf_history = GridFunction(fes, multidim=0)
num_steps = 20
with TaskManager():
for step in range(num_steps):
par.Set((step + 1) / num_steps)
solvers.Newton(
a,
gf_solution,
maxerr=1e-8,
printing=False,
inverse="",
)
scene.Redraw()
tw.value = f"step = {step+1}/{num_steps}"
gf_history.AddMultiDimComponent(gf_solution.vec)
Draw(
gf_history.components[0],
mesh,
animate=True,
min=0,
max=1,
autoscale=True,
deformation=True,
settings={
"camera": {
"transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -2.5}]
}
},
);
Noting that
Find
from ngsolve import *
from ngsolve.webgui import Draw
import ipywidgets as widgets
from ngsolve.meshes import MakeQuadMesh
L_1 = 48 # 48m
L_2 = 44 # 44m
L_3 = 16 # 16m
def GenerateMesh(num_el=2):
mesh = MakeQuadMesh(
nx=num_el,
ny=num_el,
mapping=lambda x, y: (L_1 * x, L_2 * x + L_2 * y - (L_2 - L_3) * x * y),
)
return mesh
mesh = GenerateMesh()
mu = Parameter(80.194) # N/mm^2
lam = Parameter(400889.8) # N/mm^2
par = Parameter(0)
force = CoefficientFunction((0, par * 32))
def NeoHooke(F):
C = F.trans * F
return mu / 2 * (Trace(C - Id(2)) - log(Det(C))) + lam / 2 * (Det(F) - 1) ** 2
def SkewT(v):
if v.dim == 3:
return 0.5 * CoefficientFunction(
(0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0), dims=(3, 3)
)
else:
return 0.5 * CoefficientFunction((0, -v, v, 0), dims=(2, 2))
order = 3
# clamped at left, free elsewhere
V = HCurl(mesh, order=order, dirichlet="left")
Sigma = Discontinuous(HDivDiv(mesh, order=order))
Hyb = NormalFacetFESpace(mesh, order=order, dirichlet="left")
Gamma = MatrixValued(L2(mesh, order=order + 1), symmetric=True)
fes = V * Sigma * Gamma * Hyb
gf_solution = GridFunction(fes)
u, Psym, Fsym, uh = fes.TrialFunction()
Fmat = Fsym + SkewT(curl(u))
n = specialcf.normal(2)
a = BilinearForm(fes, symmetric=True, condense=True)
a += Variation((NeoHooke(Fmat)) * dx).Compile()
a += Variation(InnerProduct(Fsym - Grad(u) - Id(2), Psym) * dx).Compile()
a += Variation((u - uh) * n * Psym[n, n] * dx(element_boundary=True)).Compile()
a += Variation(-(force * uh.Trace() + force * u.Trace()) * ds("right")).Compile()
gf_u, gf_Psym, gf_Fsym, gf_uh = gf_solution.components
gf_F = gf_Fsym + SkewT(curl(gf_u))
# F=I
gf_Fsym.Set(Id(2))
scene = Draw(
gf_u,
mesh,
deformation=True,
settings={
"camera": {
"transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -1.5}]
}
},
)
tw = widgets.Text(value="step = 0")
display(tw)
gf_history = GridFunction(fes, multidim=0)
num_steps = 40
with TaskManager():
for step in range(num_steps):
par.Set((step + 1) / num_steps)
solvers.Newton(
a,
gf_solution,
maxerr=1e-8,
printing=False,
inverse="sparsecholesky",
)
scene.Redraw()
tw.value = f"step = {step+1}/{num_steps}"
gf_history.AddMultiDimComponent(gf_solution.vec)
Draw(
gf_history.components[0],
mesh,
animate=True,
min=0,
max=25,
autoscale=True,
deformation=True,
settings={
"camera": {
"transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -1.5}]
}
},
);
Alternative TDNNS formulations for nonlinear elasticity:
Use Cauchy-Green strain tensor
instead of as additional strain field. Advantage: symmetric field independent of rigid body motions. Challenge: Handling nonlinear distribution [Neunteufel, Pechstein, Schöberl. Three-field mixed finite element methods for nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, (2021).]Use Updated Lagrangian scheme for linear TDNNS formulation [Pechstein. Large deformation mixed finite elements for smart structures. Mechanics of Advanced Materials and Structures, (2020).]