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.749999999999973
it 1 energy =  8.811178009630002
it 2 energy =  8.748116754115088
it 3 energy =  8.747829199068358
it 4 energy =  8.747829118196048
it 5 energy =  8.747829118181393
it 6 energy =  8.747829118181397
it 7 energy =  8.747829118181397
it 8 energy =  8.747829118181397
it 9 energy =  8.747829118181397

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

it 0 energy =  8.743591365303297
it 1 energy =  8.781886726588903
it 2 energy =  8.741979242436797
it 3 energy =  8.741888204568422
it 4 energy =  8.741861486112727
it 5 energy =  8.741850797186476
it 6 energy =  8.741850632571305
it 7 energy =  8.741850632185969
it 8 energy =  8.74185063218597
it 9 energy =  8.741850632185972

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

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

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