3.2. Twisted beam

cover page picture of:

Nonlinear Solid Mechanics for Finite Element Analysis: Statics 1st Edition by Javier Bonet (Author), Antonio J. Gil (Author), Richard D. Wood (Author)

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
import ipywidgets as widgets

ea = { "euler_angles" : (-120,70,116) }

d = 0.01
box = Box ( (-d/2,-d/2,0), (d/2,d/2,0.1) ) + Box( (-d/2, -3*d/2,0.1), (d/2, 3*d/2, 0.1+d) )
box.faces.Min(Z).name = "bottom"
box.faces.Max(Z).name = "top"
Draw (box, **ea);

mesh = Mesh(OCCGeometry(box).GenerateMesh(maxh=0.005))
Draw (mesh, **ea);

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

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

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

A follower load is rotated with the deformation gradient \(F = I + \nabla u\).

loadfactor = Parameter(1)
force = loadfactor * CF ( (-y, x, 0) )

fes = H1(mesh, order=3, dirichlet="bottom", dim=mesh.dim)
u,v = fes.TnT()

a = BilinearForm(fes, symmetric=True)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += ((Id(3)+Grad(u.Trace()))*force)*v*ds("top")

gfu = GridFunction(fes)
gfu.vec[:] = 0

gfu_history = GridFunction(fes, multidim=0)
scene = Draw (gfu, deformation=True, **ea)
tw = widgets.Text(value='step = 0')
display(tw)

gfu.vec[:] = 0

numsteps = 30
for step in range(numsteps):
    loadfactor.Set(300*step/numsteps)
    solvers.Newton (a, gfu, printing=False, dampfactor=0.5)
    scene.Redraw()
    tw.value = 'step = '+str(step+1)+'/'+str(numsteps)
    gfu_history.AddMultiDimComponent(gfu.vec)

Draw (gfu_history, mesh, animate=True, min=0, max=0.04, autoscale=True, deformation=True, **ea);