Ctrl+K

Mixed elements

6.4. Mixed elements#

Flat low-order elements are prone to shear locking.

A simple visualization: on a two-element mesh of a cantilever beam, prescribe the vertical displacement at the tip. Then, visualize the shear strain for high and low order approximations.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
# from netgen.webgui import Draw as DrawGeo

SetNumThreads(4)

from ngsolve.meshes import MakeStructured2DMesh

length = 10
height = 1
def mapping(x, y): return (length*x, height*(y-0.5))

mesh = MakeStructured2DMesh(nx = 2, ny = 1, mapping=mapping)
scene = Draw(mesh)

E, nu = 21000, 0.3
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def psi(u):
    strain = Sym(Grad(u))
    return mu*InnerProduct(strain,strain) + lam/2*Trace(strain)**2


fes_u = VectorH1(mesh, order=1, dirichlet="left")
fes_lag = FESpace("number", mesh)
fes = fes_u * fes_lag
u, f = fes.TrialFunction()
q = GridFunction(fes)


a = BilinearForm(fes, symmetric=True)
a += Variation(psi(u)*dx)
a += Variation(-f*(u[1]+2)*ds("right"))

solvers.Newton(a, q, printing=False)

disp = q.components[0]
strain = Sym(Grad(disp))
scene = Draw(strain[1], mesh, "shear", deformation=disp)

print(f"force: {q.components[-1].vec[0]}")

force: -114.09779811266799

Shear strain is around 5%, force necessary to achieve this deformation is \(-114\text{ N/mm}^2\).

Try the same for higher order - shear strain is much lower (below 0.1%), necessary force only \(-11.5\text{ N/mm}^2\).

fes_u = VectorH1(mesh, order=5, dirichlet="left")
fes_lag = FESpace("number", mesh)
fes = fes_u * fes_lag
u, f = fes.TrialFunction()
q = GridFunction(fes)


a = BilinearForm(fes, symmetric=True)
a += Variation(psi(u)*dx)
a += Variation(-f*(u[1]+2)*ds("right"))

solvers.Newton(a, q, printing=False, maxit=2)

disp = q.components[0]
strain = Sym(Grad(disp))
scene = Draw(strain[1], mesh, "shear", deformation=disp)
print(f"force: {q.components[-1].vec[0]}")

force: -11.542769435564512

This indicates that flat low order elements seem much too stiff. Shear strain (and stress) are over-estimated. The TDNNS method is an approach in relaxing the shear locking by relaxing continuity assumptions.

6.4.1. The TDNNS idea#

Relax normal continuity of the displacement at element interfaces - allow for gaps between elements.

To control the gaps, introduce normal stresses as independent unknowns (Lagrangian multiplier concept)

6.4.1.1. Implementation in a nutshell#

  • relax normal continuity of \(\vec u\): use tangential continuous Nedelec elements well-known from Maxwell’s equations (and implemented in ngsolve)

  • introduce the stress tensor as independent unknown via Legendre transform

  • ensure normal-normal continuity of the stress tensor, such that normal stresses can act as Lagrangian multipliers:

  • approach 1: normal-normal continuous symmetric elements \(H(\operatorname{div}\operatorname{div})\)

  • approach 2: hybridization to generate this continuity through another set of Lagrangian multipliers

No details, just visualize the shear strain for the lowest order elements.

Discontinuity (i.e. the gaps) is not visible. Strain (as computed from stress) is small (below \(0.1\)%), necessary force is under-estimated at \(-11.45 \text{ N/mm}^2\)

fes_u = HCurl(mesh, order=1, dirichlet="left")
fes_sigma = HDivDiv(mesh, order=1, dirichlet="bottom|top|right")
fes_lag = FESpace("number", mesh)
fes = fes_u * fes_sigma * fes_lag
u, sigma, f = fes.TrialFunction()
q = GridFunction(fes)

n = specialcf.normal(2)

def psi_s(sigma): 
    D = mesh.dim
    trsigma = Trace(sigma)
    devsigma = sigma - 1/D*trsigma*Id(D)
    return -0.5*(1/2/mu*InnerProduct(devsigma,devsigma) + 1/D/(D*lam+2*mu)*trsigma**2)

a = BilinearForm(fes, symmetric=True)
a += Variation(psi_s(sigma)*dx)
a += Variation(InnerProduct(sigma, Sym(Grad(u)))*dx)
a += Variation(- InnerProduct(sigma*n, n) * InnerProduct(u,n) * dx(element_boundary=True))
a += Variation(-f*(u.Trace()[1]+2)*ds("right"))

solvers.Newton(a, q, maxit=2, printing=False)

disp = q.components[0]
stress = q.components[1]
strain = psi_s(stress).Diff(stress)
scene = Draw(strain[1], mesh, "shear", deformation=disp)

print(f"force: {q.components[-1].vec[0]}")

Warning: Newton might not converge! Error =  8.55632255931183e-11

force: -11.452567283826077

6.4.2. A Hellinger-Reissner formulation#

First step, introduce stress tensor \(\boldsymbol{\sigma}\) as independent unknown. Do so via a Legendre transformation, aka co-energy, explicitely computeable for linear isotropic elasticity

\[ \psi_s(\boldsymbol\sigma) = \min_{\boldsymbol \varepsilon} \psi - \boldsymbol \sigma : \boldsymbol \varepsilon = -\frac{1}{2}\left(\frac{1}{2\mu} \operatorname{dev}\boldsymbol\sigma : \operatorname{dev}\boldsymbol\sigma + \frac{1}{d(d\lambda + 2 \mu)} (\operatorname{tr}\boldsymbol\sigma)^2\right) =: \frac{1}{2} \boldsymbol\sigma : \mathbf A : \boldsymbol \sigma. \]

Stress-strain relation, using flexibility tensor \(\mathbf A = \mathbf C^{-1}\)

\[ \boldsymbol\varepsilon = -\frac{\partial \psi_s}{\partial \boldsymbol \sigma} = \mathbf A : \boldsymbol \sigma \]

Then, use this in a Hellinger-Reissner type formulation

\[ \int_\Omega \delta (\psi_s(\boldsymbol\sigma) + \boldsymbol \sigma : \boldsymbol \varepsilon)\, dV - \delta W_{ext} = 0. \]

Explicitely put, using kinematically admissible virtual displacement \(\delta \vec u\), according strain \(\delta \boldsymbol \varepsilon\), and admissible virtual stress \(\delta \boldsymbol \sigma\),

\[ \int_\Omega (\boldsymbol\sigma : \mathbf A : \delta \boldsymbol \sigma + \delta \boldsymbol \sigma : \boldsymbol \varepsilon + \boldsymbol \sigma : \delta \boldsymbol \varepsilon)\, dV - \delta W_{ext} = 0. \]

Kinematically admissible usually means that \(\vec u\) satisfies kinematic (displacement) boundary conditions. TDNNS choice is different. On a finite element mesh (Sobolev space analysis see papers)

  • \(\vec u\) and \(\delta \vec u\) tangentially continuous, i.e. \(\vec u_t = \vec u - u_n \vec n\) with \(u_n = \vec u \cdot \vec n\) is equivalent on all elements attached to one edge.

  • \(\boldsymbol \sigma\) and \(\delta \boldsymbol \sigma\) normal-normal continuous, i.e. \(\vec n \cdot \boldsymbol \sigma \cdot \vec n =: \sigma_{nn}\) is equivalent on both elements attached to one face/edge.

  • \(\vec u_t = 0\) on \(\Gamma_D\) and \(\delta \vec u_t = 0\) on \(\Gamma_D\)

  • \(\sigma_{nn} = t_n\) on \(\Gamma_N\) and \(\delta \sigma_{nn} = 0\) on \(\Gamma_N\)

Tangential displacements and normal stresses are degrees of freedom and constitute essential boundary conditions.

These spaces are implemented in ngsolve. Consider a plate with hole example.

from netgen.meshing import IdentificationType 

ea = { "euler_angles" : [-50,3,-30] }

rad=15/2 #radius of hole

#substrate geometry
len_a = 25 #width
he = 0.5 #height of substrate

BCDisp="fix"

rect = WorkPlane().Rectangle(len_a, len_a).Face()
circ = WorkPlane().Circle(len_a/2, len_a/2,rad).Face()
circ.maxh=5

plate = (rect-circ).Extrude(he)

plate.faces.name = "free"
plate.faces.Min(X).name = "fix"
plate.faces.Max(X).name = "force"

plate.faces.Min(Z).Identify(plate.faces.Max(Z), name="plate_identification", type=IdentificationType.CLOSESURFACES)

geo = OCCGeometry(plate)

ngmesh = geo.GenerateMesh(maxh=2)


mesh = Mesh(ngmesh)


scene = Draw(mesh, **ea)

order = 1
mesh.Curve(order+1)

fes_u = HCurl(mesh, order=order, dirichlet="fix")
fes_sigma = HDivDiv(mesh, order=order, dirichlet="force|free")

fespace = fes_u * fes_sigma

q = GridFunction(fespace)

n = specialcf.normal(3)

For \(\vec u \in H(\operatorname{curl})\), the strain \(\boldsymbol \varepsilon\) is not defined in weak sense. However, for the variational formulation, it is necessary to evaluate work pairs of the form

\[ \int_\Omega \boldsymbol \sigma : \boldsymbol \varepsilon\, dV \]

As we have additional continuity for the stress, such an evaluation is possible in the sense of distributions,

\[ \int_\Omega \boldsymbol \sigma : \boldsymbol \varepsilon\, dV \qquad \text{evaluated via} \qquad \langle \boldsymbol{\varepsilon} ,\boldsymbol{\sigma} \rangle = \sum_{T} \left( \int_T \boldsymbol{\varepsilon} : \boldsymbol{\sigma}\, dV - \int_{\partial T} u_n \sigma_{nn}\, dS\right) \]

As expected, \(\sigma_{nn}\) acts as Lagrangian multiplier for the displacement gap \(u_n\)!

An alternative, equivalent representation needs the tangetial components of stress and displacement vectors, \(\vec u_t = \vec u - \vec u \cdot \vec n \, \vec n\), \(\vec \sigma_{nt} = (\boldsymbol{\sigma} \cdot \vec n) - \sigma_{nn} \vec n\),

\[ \langle \boldsymbol{\varepsilon} ,\boldsymbol{\sigma} \rangle = \sum_{T} \left( -\int_T \operatorname{div} \boldsymbol{\sigma}\cdot \vec u \, dV + \int_{\partial T} \vec u_t \vec \sigma_{nt}\, dS\right) \]
force = -0.1
u_, sigma_ = fespace.TrialFunction()

def tang(u): return u - InnerProduct(u,n)*n
def normal(u): return InnerProduct(u,n)*n

def psi_s(sigma): 
    D = mesh.dim
    trsigma = Trace(sigma)
    devsigma = sigma - 1/D*trsigma*Id(D)
    return -0.5*(1/2/mu*InnerProduct(devsigma,devsigma) + 1/D/(D*lam+2*mu)*trsigma**2)


a = BilinearForm(fespace, eliminate_internal=True, symmetric=True)
a += Variation ((psi_s(sigma_))*dx).Compile()

a += Variation ( -InnerProduct(div(sigma_),u_)*dx).Compile()
a += Variation ( InnerProduct(sigma_*n, tang(u_))*dx(element_boundary=True)).Compile()

a += Variation ( - force*u_.Trace()[2]*ds(definedon=mesh.Boundaries("force")))

solvers.Newton(a, q, maxit=2, printing=False)

disp = q.components[0]
stress = q.components[1]
strain = psi_s(stress).Diff(stress)
scene = Draw(BoundaryFromVolumeCF(stress[0]), mesh, "shear", deformation=BoundaryFromVolumeCF(disp), **ea)

Warning: Newton might not converge! Error =  4.706126909417359e-11

6.4.3. Hybridization#

Break nn-continuity of \(\boldsymbol \sigma\), add another Lagrangian multiplier to reinforce it (equivalently).

Interpretation: Lagrangian multiplier corresponds to normal displacement

Benefit: static condensation leads to positive definite system matrix (better condition number, solvers, preconditioning etc). Also, more standard choice of essential boundary conditions.

order = 2
mesh.Curve(order+1)

fes_u = HCurl(mesh, order=order, dirichlet="fix")
fes_sigma = Discontinuous(HDivDiv(mesh, order=order))
fes_un = HDiv(mesh, order=order, orderinner=0, dirichlet="fix")

fespace = fes_u * fes_un * fes_sigma

q = GridFunction(fespace)

n = specialcf.normal(3)

force = -0.1
u_, un_, sigma_ = fespace.TrialFunction()

def tang(u): return u - InnerProduct(u,n)*n
def normal(u): return InnerProduct(u,n)*n

def psi_s(sigma): 
    D = mesh.dim
    trsigma = Trace(sigma)
    devsigma = sigma - 1/D*trsigma*Id(D)
    return -0.5*(1/2/mu*InnerProduct(devsigma,devsigma) + 1/D/(D*lam+2*mu)*trsigma**2)


a = BilinearForm(fespace, eliminate_internal=True, symmetric=True)
a += Variation ((psi_s(sigma_))*dx).Compile()

a += Variation ( -InnerProduct(div(sigma_),u_)*dx).Compile()
a += Variation ( InnerProduct(sigma_*n, normal(un_) + tang(u_))*dx(element_boundary=True)).Compile()

a += Variation ( - force*u_.Trace()[2]*ds(definedon=mesh.Boundaries("force")))

solvers.Newton(a, q, maxit=2, printing=False)

disp = q.components[0]
stress = q.components[2]
strain = psi_s(stress).Diff(stress)
scene = Draw(BoundaryFromVolumeCF(stress[0]), mesh, "shear", deformation=BoundaryFromVolumeCF(disp), **ea)

Warning: Newton might not converge! Error =  1.2807284671622876e-08

6.4.3.1. Application#

see piezo tutorial example

Contents