5. Elastic Pendulum

\[ \textit{Here is presented a simple example of an elastic pendulum} \]

from ngsolve import *
from ngsolve.webgui import Draw

The equation of Elastodynamics is:

\[\begin{align*} \rho \ddot u &= \nabla \cdot \sigma + f^{\text{ext}} \\ \sigma &= C : E \\ E &= \frac{1}{2} (\nabla u + \nabla u^T + \nabla u ^T\nabla u) \end{align*}\]

where \(\rho\) is the density, \(u\) is the displacement, \(\sigma\) is the stress, \(C\) is the elasticity tensor, \(\varepsilon\) is the strain, and \(f^{\text{ext}}\) is the external force.

Suppose that at \(t=0\) the displacement and the velocity are zero, moreover the external force is the gravity force.

from netgen.occ import *
bar = Rectangle(0.1,1).Face()
bar.edges.Max(Y).name="top"
bar.faces.name="bar"
bar.faces.maxh=0.05

hole = Circle((0.05, 0), 0.1).Face()
hole.faces.name="hole"

circ = Circle((0.05, 0), 0.2).Face()
circ.faces.name="circ"
circ.faces.maxh=0.1

shape = Glue([bar - (circ), circ - hole])

shape = shape.Rotate(Axis((0.05, 1, 0), (0, 0, 1)), 180)
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1))
mesh.Curve(4)
Draw (mesh);

Given the disolacement, we define the Right Cauchy-Green tensor as:

\[\begin{align*} C = (I + \nabla u)^T (I + \nabla u) \end{align*}\]

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

The Lamé parameters \(\lambda\) and \(\mu\) are derived from the Young’s modulus \(E\) and the Poisson’s ratio \(\nu\) as follows:

# E, nu = 210, 0.2
E, nu = 80, 0.4  # rubber

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

We can define the hyperelastic material strain energy density function as:

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

The upper boundary is free to swing, to do so we control the \(x\) and \(y\) mean value of the upper boundary using a Number space.

V = VectorH1(mesh, order=2)
Q = NumberSpace(mesh, definedon=mesh.Boundaries('top'))
fes = V * Q**2
(u,q), (v, p) = fes.TnT()

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

gfut = GridFunction(V, multidim=0)

# define the needed GridFunctiond required in the scheme
gfu = GridFunction(fes)
gfv = GridFunction(fes)
gfa = GridFunction(fes)

gfuold = GridFunction(fes)
gfvold = GridFunction(fes)
gfaold = GridFunction(fes)

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

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

Here the unknowns are the displacement \(u^{n+1}\) and the acceleration \(a^{n+1}\), the velocity are to be determined via the time-stepping scheme.

bfa = BilinearForm(fes)
bfa += Variation(NeoHooke(C(u))*dx).Compile() 
bfa += (InnerProduct(u, p) + InnerProduct(v, q)) * ds('top') # we add the constraints

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*}\]

That becomes:

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

tau = 0.025
tend = 10
force = CF((0, -1))


vel_new = 2/tau * (u-gfuold.components[0]) - gfvold.components[0]
acc_new = 2/tau * (vel_new-gfvold.components[0]) - gfaold.components[0]

# need to add to the bilinear form since it depends on the current valurs of the GridFunctions

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

from ngsolve.solvers import Newton
gfu.vec[:] = 0
gfut.AddMultiDimComponent(gfu.components[0].vec)
#scene = Draw(gfu.components[0], mesh, "deformation", deformation=True)  
t = 0
i =1
with TaskManager():
    while t < tend:
        i += 1
        t += tau
        Newton(a=bfa, u=gfu, printing=False, inverse="sparsecholesky")

        if i % 5 == 0:
            #scene.Redraw()
            gfut.AddMultiDimComponent(gfu.components[0].vec)
        gfv.vec[:] = 2/tau * (gfu.vec-gfuold.vec) - gfvold.vec
        gfa.vec[:] = 2/tau * (gfv.vec-gfvold.vec) - gfaold.vec

        gfuold.vec[:] = gfu.vec
        gfvold.vec[:] = gfv.vec
        gfaold.vec[:] = gfa.vec

settings = {"Multidim": {
    "speed" : 10
}}

#Draw(gfut, mesh, animate=True, min=0, 
#     max=0.5, autoscale=True, deformation=True)#, settings = settings);
Draw(gfut, mesh, interpolate_multidim=True, deformation = True, animate=True, autoscale = False, min = 0, max = 2,  settings = settings);