This page was generated from wta/elasticity.ipynb.

Elasticity

import netgen.geom2d as geom2d
from ngsolve import *
from ngsolve.webgui import Draw

a rectangle with refinement at corners:

geo = geom2d.SplineGeometry()
pnums = [geo.AddPoint (x,y,maxh=0.01) for x,y in [(0,0), (1,0), (1,0.1), (0,0.1)] ]
for p1,p2,bc in [(0,1,"bot"), (1,2,"right"), (2,3,"top"), (3,0,"left")]:
     geo.Append(["line", pnums[p1], pnums[p2]], bc=bc)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))

Cauchy-Green tensor and hyperelastic energy density:

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def C(u):
    F = Id(2) + Grad(u)
    return F.trans * F

def NeoHooke (C):
    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)

stationary point of total energy:

\[\delta \int W(C(u)) - f u = 0\]

factor = Parameter(0)
force = CoefficientFunction( (0,factor) )

fes = H1(mesh, order=4, dirichlet="left", dim=mesh.dim)
u  = fes.TrialFunction()

a = BilinearForm(fes, symmetric=True)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(force,u)).Compile()*dx)

u = GridFunction(fes)
u.vec[:] = 0

a simple Newton solver, using automatic differentiation for residual and tangential stiffness:

def SolveNewton():
    res = u.vec.CreateVector()

    for it in range(10):
        print ("it", it, "energy = ", a.Energy(u.vec))
        a.Apply(u.vec, res)
        a.AssembleLinearization(u.vec)
        inv = a.mat.Inverse(fes.FreeDofs() )
        u.vec.data -= inv*res

factor.Set(0.4)
SolveNewton()
scene = Draw (C(u)[0,0], mesh, deformation=u)

it 0 energy =  8.749999999999964
it 1 energy =  8.81117801242968
it 2 energy =  8.748116754101272
it 3 energy =  8.747829199029733
it 4 energy =  8.747829118157389
it 5 energy =  8.747829118142732
it 6 energy =  8.747829118142734
it 7 energy =  8.747829118142734
it 8 energy =  8.747829118142734
it 9 energy =  8.747829118142734

factor.Set(factor.Get()+0.4)
SolveNewton()
scene.Redraw()

it 0 energy =  8.743591365193794
it 1 energy =  8.7818867275068
it 2 energy =  8.741979242316216
it 3 energy =  8.741888204442981
it 4 energy =  8.74186148598637
it 5 energy =  8.741850797059199
it 6 energy =  8.741850632444025
it 7 energy =  8.741850632058682
it 8 energy =  8.74185063205868
it 9 energy =  8.74185063205868

Compute \(2^{nd}\) Piola Kirchhoff stress tensor by symbolic differentiation:

C_=C(u).MakeVariable()
sigma = NeoHooke(C_).Diff(C_)

Draw (sigma[0,0], mesh, "Sxx", deformation=u, min=-10.001, max=10.001);