56. Elastodynamics with Newmark time-stepping

\(u\) … displacement, \(v = \dot u \) … velocity, \(a = \dot v\) .. acceleration

Newmark scheme:

\[\begin{align*} \frac{u^{n+1}-u^n}{\tau} = \frac{v^n+v^{n+1}}{2} \\ \frac{v^{n+1}-v^n}{\tau} = \frac{a^n+a^{n+1}}{2} \\ \end{align*}\]

with new acceleration, with elasticity operator \(K\):

\[ a^{n+1} = f - K(u^{n+1}) \]

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *
shape = Rectangle(1,0.1).Face()
shape.edges.Max(X).name="right"
shape.edges.Min(X).name="left"
shape.edges.Max(Y).name="top"
shape.edges.Min(Y).name="bot"
shape.vertices.Min(X+Y).maxh=0.01
shape.vertices.Min(X-Y).maxh=0.01
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))
Draw (mesh);

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)

tau = 0.02
tend = 20
force = CF( (0, -1) )

fes = VectorH1(mesh, order=3, dirichlet="left")
u,v = fes.TnT()

gfu = GridFunction(fes)
gfv = GridFunction(fes)
gfa = GridFunction(fes)
gfuold = GridFunction(fes)
gfvold = GridFunction(fes)
gfaold = GridFunction(fes)

bfa = BilinearForm(fes)
bfa += Variation(NeoHooke(C(u))*dx)

vel_new = 2/tau * (u-gfuold) - gfvold
acc_new = 2/tau * (vel_new-gfvold) - gfaold

bfa += acc_new*v*dx
bfa += -force*v*dx

def MyNewton(a, u, printing=False):
    for it in range(10):
        res = a.Apply(gfu.vec)
        if printing:
            print ("it=", it, "err=", Norm(res))
        a.AssembleLinearization(gfu.vec)
        inv = a.mat.Inverse(fes.FreeDofs() ) 
        gfu.vec.data -= inv*res

from ngsolve.solvers import Newton
sceneu = Draw (gfu, deformation=True)
scenev = Draw (gfv)
scenea = Draw (gfa)
gfu.vec[:] = 0
t = 0
while t < tend:
    t += tau
    # MyNewton(a=bfa, u=gfu, printing=False)
    Newton(a=bfa, u=gfu, printing=False, inverse="sparsecholesky")

    gfv.vec[:] = 2/tau * (gfu.vec-gfuold.vec) - gfvold.vec
    gfa.vec[:] = 2/tau * (gfv.vec-gfvold.vec) - gfaold.vec

    sceneu.Redraw()
    scenev.Redraw()
    scenea.Redraw()
    
    gfuold.vec[:] = gfu.vec
    gfvold.vec[:] = gfv.vec
    gfaold.vec[:] = gfa.vec
    # input("step")