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

7 months 2 weeks ago #4478
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()

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

mesh.Export('TestSphere.msh', 'Gmsh2 Format')

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

def C(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)

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)

7 months 2 weeks ago #4481 by joachim
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

7 months 2 weeks ago #4483
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()

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.5*E/(1 + nu)
lam = E * nu / ((1 + nu)*(1 - 2*nu))

def C(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)

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)

7 months 2 weeks ago #4484 by joachim
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.5*lam/nu) - 1)
return mu/2 * (Trace(C)-3) - mu/2*log(Det(C)) + lam/2*(sqrt(Det(C))-1)**2