12.7. Elastodynamics

Idea: Discretize velocity by push-forward from deformation

12.7.1. Scalar wave equation:

\[ \ddot u = \Delta u \]

reduction to first order in time:

\[\begin{eqnarray*} \dot u & = & v \\ \dot v & = & -\Delta u \end{eqnarray*}\]

Solution space (Bochner spaces):

\[ u \in L_2([0,T], H^1(\Omega)), \qquad v \in L_2([0,T], L_2(\Omega)) \]

Do our discretizations respect these spaces?

12.7.2. Elastic wave equation:

\[ u \in L_2([0,T], H^1(\Omega)^d), \qquad v \in L_2([0,T], H(\operatorname{curl})) \]

See Michi Neunteufel’s master thesis, Chapter 6.4

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

geo = WorkPlane().RectangleC(10,1).Face()
mesh = geo.GenerateMesh(maxh=0.4, dim=2)

tend = 10
taumax = 0.2
tau = Parameter(taumax)

order = 1

dirichlet = "xxx"

Vu = VectorH1(mesh, order=order, dirichlet=dirichlet)
Vv = Discontinuous(HCurl(mesh, order=order, dirichlet=dirichlet)) 
Vp = Discontinuous(HCurl(mesh, order=order, dirichlet=dirichlet))
V = Vu*Vv*Vp

u, v, p = V.TrialFunction()
ut, vt, pt = V.TestFunction()

p = p.Operator("dual")
pt = pt.Operator("dual")

gfu = GridFunction(V)
gfuold = GridFunction(V)

uold = gfuold.components[0]
vold = gfuold.components[1]
pold = gfuold.components[2].Operator("dual")

um = 0.5*(u+uold)
vm = 0.5*(v+vold)
# vm = v
pm = 0.5*(p+pold)

# Id = CoefficientFunction( (1,0,0,1), dims=(2,2) )
F = Grad(u) + Id(2)

Fold = Grad(uold) + Id(2)
Ft = grad(ut)
# Fm = 0.5 * (F+Fold)
Fm = 2*Inv (Inv(F)+Inv(Fold))
# Cm = 0.5 * (F.trans*F + Fold.trans*Fold)
C = F.trans*F
Cold = Fold.trans*Fold
Cm = 0.5*(C+Cold)

def Rot(u): return CF( ( u[1], -u[0]) )
E = 0.5 * (F.trans * F - Id(2))
def Pow(a, b): return a**b 

truecompile = True
nonlinear_material = True

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

def NeoHook (C):
    if nonlinear_material:
        return 0.5 * mu * (Trace(C-Id(2)) + 2*mu/lam * (Pow(Det(C), -lam/2/mu) - 1))
    else:
        return 0.5 * mu * InnerProduct(C-Id(2),C-Id(2))

def Stress (C):
    if nonlinear_material:    
        return 0.5 * mu * (Id(2) - Pow(Det(C), -lam/2/mu) * Inv(C))
    else:
        return mu * (C-Id(2))

n = specialcf.normal(2)
def Ptau (v):
    return v - (v*n)*n
    

a = BilinearForm(V, condense=True)
a += (1/tau*Fm.trans*(u-uold)-vm) * pt * dx(element_boundary=False)
a += (1/tau*Fm.trans*(u-uold)-vm) * pt * dx(element_boundary=True)

a += (p-pold)*vt * dx(element_boundary=False)
a += (p-pold)*vt * dx(element_boundary=True)

a += Inv(Cm)*(-v+vold)*vt*dx
a += -0.5 * (Inv(Cm)*(Cold-C)*Inv(Cm)*vm)*vt *dx 
a +=  (((tau*0.5/Det(Cm)*(0.5*curl(v)+0.5*curl(vold)))*Rot(vm))) * vt*dx 
    
a +=  (1/tau * (Fm*(p-pold)) * ut)*dx(element_boundary=True)
a +=  (1/tau * (Fm*(p-pold)) * ut)*dx(element_boundary=False)

a += 2*InnerProduct(Stress(Cm), Fm.trans*Ft)*dx
# a += SymbolicBFI ( 5 * ut[1] )



akin = BilinearForm(V)
akin += Variation ((0.5*(Inv(F.trans)*v) * (Inv(F.trans)*v)) * dx)

apot = BilinearForm(V)
apot += Variation ( NeoHook(C)*dx )
# apot += Variation ( 5*u[1]*dx)

u = gfu.components[0]
F = grad(gfu.components[0]).trans + Id(2)

print ("Velocity y-component:")
scene_velhat = Draw (gfu.components[1][1], mesh, "vel-hat", deformation=u)

Velocity y-component:

vold.Set( CF( (-y,x) ))

res = gfu.vec.CreateVector()
w = gfu.vec.CreateVector()

t = 0
gfu.vec.data = gfuold.vec

errE = []
from time import sleep
with TaskManager():
  while t < tend+1e-8:
    # print (t)

    newton_conv = True
    lasterr = 1
    for step in range(15):
        a.AssembleLinearization(gfu.vec)
        a.Apply (gfu.vec, res)
        inv = a.mat.Inverse(V.FreeDofs(True), inverse="umfpack")

        res.data += a.harmonic_extension_trans * res
        w.data = inv * res
        w.data += a.harmonic_extension * w
        w.data += a.inner_solve * res

        err = Norm(w)
        # print ("|w| = ", err)
        if lasterr < 1e-8: break
        lasterr = err
        
        gfu.vec.data -= w
    else:
        newton_conv = False
        print ("Newton did not converge")

    # print ("step =", step, "tau =", tau.Get())
    if newton_conv:
        t = t + tau.Get()
        gfuold.vec.data = gfu.vec
    else:
        tau.Set(0.5*tau.Get())
        gfu.vec.data = gfuold.vec

    if step <= 8:
        tau.Set(1.1*tau.Get())        
    if step >= 10:
        tau.Set(0.8*tau.Get())
    if tau.Get()  > taumax:
        tau.Set(taumax)

    tau.Set(taumax)
    Redraw(blocking=True)
    scene_velhat.Redraw()

    epot = apot.Energy (gfu.vec)
    ekin = akin.Energy (gfu.vec)
    # print ("pot =", epot, ", ekin =", ekin, ", sum = ", epot+ekin)
    errE.append (epot+ekin)

    sleep(0.1)

from matplotlib.pyplot import plot

plot (errE)

[<matplotlib.lines.Line2D at 0x114259400>]