Small contact area of force results in no displacement in an elastic body

Errors
1

Hello.

I am doing a simulation of elastic bodies and ran into some strange results.
This is my MRE:
I have an elastic bar (10x10x100) and am compressing it with a set force on the long axis.
Now I wrote code to apply my force not on the whole area (10x10) on the top but only a square in the middle of it with a side length of l.

Now when I set l < 2.82 I get a maximum displacement of 0.0.
I think that is an error. If I apply an absurd force of 1e90N to a 2x2 square on a 10x10 square, there should be some visible displacement but fenics says theres nothing.

Is that a mistake in my code and can I fix that somehow?

import fenics as fx
import mshr

#define the mesh in m
domain = mshr.Box(fx.Point(0, 0, 0), fx.Point(10, 10, 100))
mesh = mshr.generate_mesh(domain, 64)

# unscaled variables for an elastic body
E = 300
G = 103
mu = G
lambda_ = G*(E-2*G) / (3*G - E)

# this follows linearElasticity.ipynb of the fenics tutorial closely
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 = 2.0 #side length of the applied square
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 = 1e90/(l*l) #absurd 1e90 force
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.solve(a == L, u, bc, solver_parameters={'linear_solver':'mumps'})

V = fx.FunctionSpace(mesh, 'P', 1)
u_magnitude = fx.sqrt(fx.dot(u, u))
u_magnitude = fx.project(u_magnitude, V)
print('min/max u:',
      u_magnitude.vector().min(),
      u_magnitude.vector().max())

This prints me 0.0, 0.0.
Setting l=10 the maximum displacement prints the theoretically expected 0.003m (when appliying 1N Force) so the code works there.

MRE_pic
MRE_pic754×541 32.8 KB

Another problem is visible when setting l=5 (black square) and applying a force of 1e3. The bar is bend not only straight down but also to the side and the side seems to be random when retrying with different values for l.

Thanks in advance!

2

What is happening here is that the boundary facets are too large relative to the load surface, so none of them are being marked as belonging to the subdomain. What loadSurface.mark(bnd_markers, 1) is doing is looking at each exterior facet of the mesh and marking a facet with 1 if all of its vertices and its midpoint are in loadSurface. (One can optionally switch off the midpoint check with the keyword argument check_midpoint=False.) If loadSurface is too small relative to the facet size, no facets meet this criterion.

In general, you will see poor accuracy when attempting to integrate over subdomains that are not unions of elements or facets, since the integral over the marked elements/facets will differ from the integral over the desired subdomain. Even for larger values of l, which give nonzero displacements, you can see that there is a significant discrepancy between l*l and fx.assemble(fx.Constant(1)*ds(1)), unless l = 10.0 (in which case loadSurface is exactly a union of mesh facets).

The “random” bending observed for intermediate values of l is a consequence of the fact that the desired load surface is being approximated (poorly) by whatever collection of facets are fully contained in the subdomain loadSurface, and the force resultant is therefore off-center.

3 Likes
3

Thanks a lot for the explanation!

Is there anything I can do about that?
Increasing the mesh resolution to something like mesh = mshr.generate_mesh(domain, 512) yields better results for the discrepancy between l*l and fx.assemble(fx.Constant(1)*ds(1)). but there I quickly get problems with the solver like:

*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 76 (Error in external library).

And switchiing off the midpoint checking did not influence this discrepancy at all.

Should I use something else than loadSurface() or define a denser mesh where I want to apply the force?

4

The best practice would be to use a mesh that conforms to any discontinuities in your problem data. Even if you came up with a way to exactly integrate the discontinuous traction on the unfitted mesh, there would still be a jump in derivatives of the exact displacement solution cutting through element interiors, which the displacement function space would not approximate optimally. Building subdomains into meshes is outside the scope of what mshr can do in 3D, though, so you’d need to use some external mesh generator like Gmsh. @dokken wrote up some Gmsh tutorials on his personal webpage:

1 Like
5

Thanks a lot I will look into that!