3D nonlinear elasticity set one component of displacement as zero on a face

Main Forum Old-Forum - do not post here
1

Hi
I am trying to learn NGSPy. I am looking at one-eight of a thick spherical with internal pressure. I want to restrict the displacements on the y-z plane to have u_x = 0, in the z-x plane to have u_y = 0, and in the x-y plane to have u_z = 0.

Can I set the face to have Dirichlet Boundary condition? Does that mean all the three components of the displacement have to be specified on the boundaries marked as dirichlet? I was not able to find an example in the documentation that sets just one component of displacement to 0. Can you kindly give a sample how to do it?

I have attached the code for the geometry, the load and the FESpace:

from ngsolve import Draw, Redraw, Id, Grad, Trace, Det, Parameter, specialcf, H1
from netgen.csg import Plane, Sphere, Pnt, Vec, CSGeometry

Prepare mesh

yz_plane = Plane (Pnt(0, 0, 0), Vec(-1, 0, 0)).bc(“xface”)
zx_plane = Plane (Pnt(0, 0, 0), Vec(0, -1, 0)).bc(“yface”)
xy_plane = Plane (Pnt(0, 0, 0), Vec(0, 0, -1)).bc(“zface”)

inner_sphere = Sphere(Pnt(0, 0, 0), 5).maxh(0.5).bc(“inner”)
outer_sphere = Sphere(Pnt(0, 0, 0), 25).maxh(5).bc(“outer”)

one_eight_sphere = ((outer_sphere - inner_sphere) * yz_plane * zx_plane * xy_plane).mat(‘rubber’)

geo = CSGeometry()
geo.Add(one_eight_sphere)

mesh = geo.GenerateMesh()
mesh.Refine()
mesh.SecondOrder()

mesh.Export(‘TestSphere.msh’, ‘Gmsh2 Format’)

Prepare constitutive equation

E, nu = 1.794, 0.3
mu = 0.5E/(1 + nu)
lam = E * nu / ((1 + nu)(1 - 2*nu))

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

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

Prepare incremental load

factor = Parameter(0)
factor.Set(0.05)
pressure = factor * specialcf.normal(3)

Define the FESpace

fes = H1(mesh, order=2, dirichlet=‘xface|yface|zface’, dim=3)

2

Hi,

if you replacer H1(… dim=3) by VectorH1 (and you don’t give the dim-flag), you can use “dirichletx”, “dirichlety” and “dirichletz” flags to set Dirichlet-constraints on the individual Cartesian coordinates.

Joachim

3

Thank you very much for the reply @joachim. I was able to get the code to run. But I am confused about the output. In my understanding, I applied pressure on the inner surface of the sphere but the output shows that the highest magnitude of displacement is on the outer surface of the sphere?! Please help me understand whether I am applying pressure on the inner surface correctly.

a += Variation((-InnerProduct(pressure, u)).Compile()*ds(‘inner’))

The full working code is below:

from netgen.csg import CSGeometry, Plane, Pnt, Sphere, Vec
from ngsolve import (BilinearForm, Det, Grad, GridFunction, Id,
InnerProduct, Mesh, Norm, Parameter, SetVisualization, Trace, Variation, VectorH1,
ds, dx, specialcf)
from ngsolve.webgui import Draw

Prepare mesh

yz_plane = Plane (Pnt(0, 0, 0), Vec(-1, 0, 0)).bc(“xface”)
zx_plane = Plane (Pnt(0, 0, 0), Vec(0, -1, 0)).bc(“yface”)
xy_plane = Plane (Pnt(0, 0, 0), Vec(0, 0, -1)).bc(“zface”)

inner_sphere = Sphere(Pnt(0, 0, 0), 5).maxh(0.5).bc(“inner”)
outer_sphere = Sphere(Pnt(0, 0, 0), 25).maxh(5).bc(“outer”)

one_eight_sphere = ((outer_sphere - inner_sphere) * yz_plane * zx_plane * xy_plane).mat(‘rubber’)

geo = CSGeometry()
geo.Add(one_eight_sphere)

mesh = geo.GenerateMesh()
mesh.Refine()
mesh.SecondOrder()
#mesh.Export(‘TestSphere.msh’, ‘Gmsh2 Format’)
mesh = Mesh(mesh)

Prepare constitutive equation

E, nu = 1.794, 0.3
mu = 0.5E/(1 + nu)
lam = E * nu / ((1 + nu)(1 - 2*nu))

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

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

Prepare incremental load

factor = Parameter(0)
factor.Set(1.5)
pressure = factor * specialcf.normal(3)

Define the FESpace

fes = VectorH1(mesh, order=2, dirichletx=‘xface’, dirichlety=‘yface’, dirichletz=‘zface’)
u = fes.TrialFunction()

Set up the governing weak-form

a = BilinearForm(fes, symmetric=False)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(pressure, u)).Compile()*ds(‘inner’))

Set up the solver

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

def SolveNewton():
res = u.vec.CreateVector()

for it in range(5):
    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

Solve!

SolveNewton()
Draw(Norm(u), mesh, “Displacement”)
SetVisualization(deformation=True)

4

The setup is fine, but I think its a problem with this version of Neo-Hooke material law.
It looks good if I use the following one (but then we need loadsteps for convergence).

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

5

Thank you very much for your reply.