Compression of elastic sphere

variational formulation
1

Hello,

I would like to simulate how an elastic sphere gets compressed when squished against two surfaces.
The goal is to get to know the force (and deformation work) needed to speeze the sphere into a certain shape
Here an illustration of what I mean:

BallDefo
BallDefo1920Ă—1440 227 KB

In practice I would probably only push from the top and have a 0 deformation dirichlet condition at the bottom. That should yield the same result using the doubled force at the top.

I know how to compress an elastic box but dont get how I can do the same to a sphere as the loadSurface should change with force as the contact area between the top surface at the sphere grows.

import fenics as fx
import mshr

domain = mshr.Box(fx.Point(0, 0, 0), fx.Point(10, 10, 100))
mesh = mshr.generate_mesh(domain, 64)

# Elastic variables
E = 300 #young modulus
nu = 0.5 #poisson ratio
G = E / (2*(1+nu))+3 #+3 to not be just on the border of physicality
mu = G # mu and lambda_ are the factors in the naviers equation
lambda_ = G*(E-2*G) / (3*G - E)

V = fx.VectorFunctionSpace(mesh, 'Lagrange', 1)

tol = 1e-6
def boundary(x, on_boundary):
    return on_boundary and (fx.near(x[2], 0, tol))

bc = fx.DirichletBC(V, fx.Constant((0, 0, 0)), boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(fx.nabla_grad(u) + fx.nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    #return lambda_*fx.nabla_div(u)*fx.Identity(d) + 2*mu*epsilon(u)
    return lambda_* fx.div(u) * fx.Identity(d) + 2*mu*epsilon(u)


# Define variational problem
u = fx.TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = fx.TestFunction(V)

# define where force applies
l = 10
class LoadSurface(fx.SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-6
        return on_boundary and fx.near(x[2], 100, tol) and \
            (x[0] >= 5 - l/2 - tol and x[0] <= 5 + l/2 + tol) and \
            (x[1] >= 5 - l/2 - tol and x[1] <= 5 + l/2 + tol)

bnd_markers = fx.MeshFunction('size_t', mesh, d-1, 0)
loadSurface = LoadSurface()
loadSurface.mark(bnd_markers, 1)
ds = fx.Measure('ds', domain=mesh, subdomain_data=bnd_markers)

force = 1e3/(l*l) 
T = fx.Constant((0, 0, -force))
a = fx.inner(sigma(u), epsilon(v))*fx.dx
L = fx.dot(T, v)*ds(1)

# Compute solution
u = fx.Function(V)
fx.parameters['ghost_mode'] = 'shared_facets' 
fx.solve(a == L, u, bc, solver_parameters={'linear_solver':'mumps'})

Is there a way to simulate that with a changing load surface or are the better methods?

2

Contact mechanics is a difficult topic, but the “entry level” method is to use a penalty regularization, where the boundary is subject to a traction that increases rapidly in proportion to how badly the contact constraint is violated. For contact of an elastic body with a rigid object that has a known, analytical geometry, this is relatively easy to implement using FEniCS. Here is a quick example of a 2D analogue to the problem you describe, albeit under displacement control (i.e., given displacement, find the force resultant) instead of load control:

from dolfin import *
from mshr import *
from matplotlib import pyplot as plt

# Create mesh of quarter circle (to take
# advantage of symmetry).
N = 32
R = 1.0
c = Circle(Point(0,0),R,4*N)
r1 = Rectangle(Point(-2*R,-2*R),Point(0,2*R))
r2 = Rectangle(Point(-2*R,-2*R),Point(2*R,0))
mesh = generate_mesh(c-r1-r2,N)
X = SpatialCoordinate(mesh)
h = CellDiameter(mesh)

# Linear elasticity, for simplicity.  Replace with a nonlinear
# model (e.g., neo-Hookean) if displacements are large.
V = VectorFunctionSpace(mesh,"CG",1)
u = Function(V)
K = Constant(1.0)
G = Constant(1.0)
I = Identity(2)
eps = sym(grad(u))
sigma = K*tr(eps)*I + 2*G*(eps-tr(eps)*I/3)

# Penalty force:
x = X + u
R = Constant(R)
E = 9*K*G/(3*K+G)
C_pen = Constant(1e3) # (Dimensionless constant)
k_pen = C_pen*E/h # (Typical choice of contact penalty)
chi_pen = conditional(gt(x[1],R),1,Constant(0))
t_pen = k_pen*chi_pen*as_vector([Constant(0),R-x[1]])

# Residual of variational problem:
v = TestFunction(V)
res = inner(sigma,grad(v))*dx - dot(t_pen,v)*ds

# Quarter-symmetry BCs, with a constant vertical displacement
# applied to the bottom.
bdry_disp = Constant(0)
bcs = [DirichletBC(V.sub(0),Constant(0),"x[0]<DOLFIN_EPS"),
       DirichletBC(V.sub(1),bdry_disp,"x[1]<DOLFIN_EPS")]

# Increase displacement and output force:
N_STEP = 7
MAX_DISP = 0.25
disp_inc = MAX_DISP/N_STEP
for i in range(0,N_STEP):
    # Increase the applied displacement and re-solve:
    bdry_disp.assign(Constant(float(bdry_disp)+disp_inc))
    solve(res==0,u,bcs=bcs)
    # Assemble net penalty force to get resultant:
    print("Displacement = "+str(float(bdry_disp))
          +" , Force = "+str(assemble(-t_pen[1]*ds)))
    # Plot results:
    plot(u,mode="displacement")
    plt.show()

Doing this correctly with load control (i.e., input force and get displacement) is a bit trickier; simply applying a uniform traction to the bottom would not be the correct symmetry condition. You would need to enforce a constraint that the vertical displacement on the bottom is a constant whose value becomes an unknown of the problem (and then apply any traction with right resultant on top of that).

I wrote up a clever answer to a question about this type of constraint some time ago, although I can’t seem to find the corresponding forum thread. I still have a local copy of the code, though, which I’ve reproduced here:

from dolfin import *
N = 64
mesh = UnitSquareMesh(N,N)
cell = mesh.ufl_cell()
uE = VectorElement("CG",cell,1)
u0E = FiniteElement("Real",cell,0)
W = FunctionSpace(mesh,MixedElement([uE,u0E]))
x = SpatialCoordinate(mesh)

# u is the displacement, while u0 is a global scalar unknown, which is the
# constant x-displacement over the right face of the domain.
u,u0 = TrialFunctions(W)
v,v0 = TestFunctions(W)

# Fix displacement at left face:
bc = DirichletBC(W.sub(0),Constant((0,0)),"x[0]<DOLFIN_EPS")

# Characteristic function for right face:
right_char = conditional(gt(x[0],Constant(1.0-DOLFIN_EPS)),1.0,Constant(0))

# Linear elasticity:
K = Constant(1.0)
G = Constant(1.0)
I = Identity(2)
def sigma(u):
    eps = sym(grad(u))
    return K*tr(eps)*I + 2*G*(eps-tr(eps)*I/3)
n = FacetNormal(mesh)
f = as_vector((Constant(3)*sin(pi*x[1]),Constant(0)))
C_pen = Constant(1e1)
h = CellDiameter(mesh)
# Nitsche-type enforcement of dot(u,n) = u0 on the right side:
res = inner(sigma(u),grad(v))*dx - inner(f,v)*dx \
    + right_char*((C_pen/h)*(dot(u,n)-u0)*(dot(v,n)-v0)
                  - dot(sigma(u)*n,n)*(dot(v,n)-v0)
                  - dot(sigma(v)*n,n)*(dot(u,n)-u0))*ds
w = Function(W)
solve(lhs(res)==rhs(res),w,[bc,])
u_out,u0_out = w.split()

# Check flatness error:
e = (dot(u_out,n)-u0_out)
print(sqrt(assemble(right_char*e*e*ds)))

u_out.rename("u","u")
File("u.pvd") << u_out
3 Likes
3

Thanks a lot for the detailed answer!
Your first code is exaclty what I’m looking for.

Why can you say that T=t_pen?
basically in 1D what you wrote is:
t_{pen}=c_{pen}\frac{E}{h}(R-x)
which kinda follows from the definition of Youngs Modulus:
E=\frac{\sigma}{\epsilon}=\frac{F/A}{\Delta l / l_0} \rightarrow f=F/A=\frac{E}{l_0}(\Delta l)
except for the forefactor.
And you insert this into the equation of
-\nabla \cdot \sigma=f=t_{pen}
where f is the force per unit volume but why is it proportional to some set constant like c_pen = 1e3?

Another question on top would be if you have to rewrite t_pen for 3D to something like t_{pen}=c_{pen}\frac{E}{h^2}(R-x) so the units work out as force per unit volume like wanted in 3.3.1 of the fenics book: https://link.springer.com/content/pdf/10.1007%2F978-3-319-52462-7.pdf

4

Notice that t_pen is integrated over a surface. Physically, it is a traction, not a body force. The choice of the penalty constant k_pen comes from several considerations:

  • The traction t_pen must have the correct units of force per unit area. This choice of units for t_pen holds in both 2D and 3D. (If you are confused by the fact that it is integrated over a 1D curve in 2D, notice that the force per unit volume in the interior problem is integrated over a 2D area; thus all terms in the problem are “missing” a factor of depth in 2D, and the the formulation as a whole remains dimensionally consistent. You could multiply the whole problem through by a constant length scale to compensate, but this would not change the result, so the depth can be taken as 1 without loss of generality.) The formula for k_pen should be the same in 2D and 3D.
  • For good conditioning of the linear system, the stiffness matrix entries from the penalty terms should have consistent scaling with respect to mesh size as the stiffness matrix entries with the interior equation.
  • The penalty factor k_pen should diverge to infinity as h\to 0. This will approach exact enforcement of the constraint as you refine the mesh.

The dimensionless constant c_pen represents a tradeoff between constraint enforcement and conditioning of the algebraic equation system for a fixed mesh size. Larger values will lead to less violation of the constraint, but will also make the system of equations more difficult to solve (which becomes especially significant when using iterative linear solvers, as one would typically use in large 3D problems). Attempting to enforce the constraint too strongly on a fixed mesh will typically also lead to large spatial oscillations in the penalty traction.

1 Like
5

Hi kamensky,
I’ve also tried to implement the Fenicsx Example of Hertzian contact as in: Hertzian contact with a rigid indenter using a penalty approach — Numerical tours of continuum mechanics using FEniCS master documentation using load control, but since I’m using dolfinx and in a 3D mesh, I have some difficulties in understanding your code. Did you maybe also make a example like this using dolfinx and willing to share it? That would be wonderful. Thank you!