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.749999999999972
it 1 energy =  8.811178274176932
it 2 energy =  8.74811675282931
it 3 energy =  8.747829195439456
it 4 energy =  8.747829114564311
it 5 energy =  8.747829114549656
it 6 energy =  8.74782911454966
it 7 energy =  8.74782911454966
it 8 energy =  8.74782911454966
it 9 energy =  8.74782911454966

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

it 0 energy =  8.743591354967782
it 1 energy =  8.781886826184756
it 2 energy =  8.741979230980744
it 3 energy =  8.74188819248704
it 4 energy =  8.741861474021357
it 5 energy =  8.74185078504465
it 6 energy =  8.741850620429458
it 7 energy =  8.741850620044117
it 8 energy =  8.74185062004412
it 9 energy =  8.74185062004412

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)

BaseWebGuiScene