Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

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

More
1 year 7 months ago #4478 by amit112amit
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.5*E/(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)
More
1 year 7 months 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
More
1 year 7 months ago #4483 by amit112amit
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.5*E/(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)
More
1 year 7 months 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
 
More
1 year 7 months ago #4488 by amit112amit
Thank you very much for your reply.
Time to create page: 0.174 seconds