Why FEniCS values are so different than theoretical ones?

Legacy DOLFIN
1

Hello everyone, I’m trying to adapt Hertzian contact 3D code given here into 2D. I am satisfied with the visualization results but the values are not matching with the theoretical values.

Maximum pressure FEniCS :  0.37048 Hertz:  1.39916
Applied force    FEniCS :  0.38579 Hertz:  0.02930

Here is MWE -

from dolfin import *
from mshr import *
import numpy as np

domain = Rectangle(Point(-1, -1),Point(1, 1))
mesh = generate_mesh(domain, 10)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1) and on_boundary
        
def bottom(x, on_boundary):
        return near(x[1], -1) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

R = 0.5
d = 0.02
obstacle = Expression("-d+(pow(x[0],2))/2/R", d=d, R=R, degree=2)

V = VectorFunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)

u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
gap = Function(V2, name="Gap")
p = Function(V0, name="Contact pressure")

bc =[DirichletBC(V, Constant((0., 0.)), bottom)]

E = Constant(10.)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
    return (x+abs(x))/2.

pen = Constant(1e4)
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[1], ppos(u[1]-obstacle))*ds(1)
J = derivative(form, u, du)

parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)

solver_parameters = {
    "newton_solver": {
        "linear_solver": "cg",
        "preconditioner" : "ilu",
        "relative_tolerance" : 1e-10,
        "absolute_tolerance" : 1e-15,
        "maximum_iterations" : 30
    }
}

form_compiler_parameters = {
    "cpp_optimize": True,
    "representation": "quadrature",
    "quadrature_degree": 2
}

solve(form == 0 , u , bc , J=J , solver_parameters= solver_parameters,form_compiler_parameters= form_compiler_parameters)

p.assign(-project(sigma(u)[1, 1], V0))
gap.assign(project(obstacle-u[1], V2))

a = sqrt(R*d)
F = 4/3.*float(E)/(1-float(nu)**2)*a*d
p0 = 3*F/(2*pi*a**2)
print("Maximum pressure FE: {0:8.5f} Hertz: {1:8.5f}".format(max(np.abs(p.vector().get_local())), p0))
print("Applied force    FE: {0:8.5f} Hertz: {1:8.5f}".format(4*assemble(p*ds(1)), F))

Any suggestions would be highly appreciable.
Thanks.

2

10*10 seems a bit low for the mesh maybe you should try to refine ( like mesh = generate_mesh(domain, 20) ) to see if you get closer.

if not there is definitely another problem

3

Thanks for the useful suggestion. It did work increasing mesh density, for 30, the FEM values increase a bit as follows -

Maximum pressure FEM:  1.04893 Hertz:  1.39916
Applied force    FEM:  0.28177 Hertz:  0.02930

But it didn’t get close to the calculated values on further increment.

I also changed the newton solver by adding the Custom Newton Solver by @nate and tried changing different linear solvers and preconditioners but it didn’t help much and the values were same as above.

Any further pointers are highly encouraged.

Here is the custom newton solver integrated MWE -

from dolfin import *
from mshr import *
import numpy as np

domain = Rectangle(Point(-1, -1),Point(1, 1))
mesh = generate_mesh(domain, 30)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1) and on_boundary
        
def bottom(x, on_boundary):
        return near(x[1], -1) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

R = 0.5
d = 0.02
obstacle = Expression("-d+(pow(x[0],2))/2/R", d=d, R=R, degree=2)

V = VectorFunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)

u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
gap = Function(V2, name="Gap")
p = Function(V0, name="Contact pressure")

bc =[DirichletBC(V, Constant((0., 0.)), bottom)]

class Problem(NonlinearProblem):
    def __init__(self, J, F, bcs):
        self.bilinear_form = J
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        assemble(self.linear_form, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self.bilinear_form, tensor=A)
        for bc in self.bcs:
            bc.apply(A)

class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("pc_type", "ilu")
        self.linear_solver().set_from_options()

E = Constant(10.)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
    return (x+abs(x))/2.

pen = Constant(1e4)
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[1], ppos(u[1]-obstacle))*ds(1)
J = derivative(form, u)

problem = Problem(J, form, bc)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())

p.assign(-project(sigma(u)[1, 1], V0))
gap.assign(project(obstacle-u[1], V2))

a = sqrt(R*d)
F = 4/3.*float(E)/(1-float(nu)**2)*a*d
p0 = 3*F/(2*pi*a**2)
print("Maximum pressure FEM: {0:8.5f} Hertz: {1:8.5f}".format(max(np.abs(p.vector().get_local())), p0))
print("Applied force    FEM: {0:8.5f} Hertz: {1:8.5f}".format(4*assemble(p*ds(1)), F))
4

This is a nonlinear problem, so in principle there may be more than one solution: is it possible that the nonlinear solver is converging to something which is still a solution, but is not the solution you were expecting?

If the “Hertz” calculations are based on an analytical formula for the solution, write down that formula in ufl/dolfin.Expression, interpolate it on the FE grid, and use that function as the initial guess for the nonlinear solver:

  1. if the nonlinear iterations converge quickly to a stationary point, and that is the solution you were expecting, then you are in the case of my first sentence.
  2. if the nonlinear iterations diverge or converge to something else, even when starting close to a function that you expect to be a solution (your favorite one), then it’s likely that you have made en error in the implementation, for instance used wrong parameter values.
1 Like
5

I think that analytical solution for displacement will be a downward parabolic Function and therefore took the initial guess as -

u.interpolate(Expression(('0.', '-d+(pow(x[0],2))/2/R'), degree = 2, d = 0.02, R = 0.5))

But still the values remain same -

Maximum pressure FEM:  1.04893 Hertz:  1.39916
Applied force    FEM:  0.28177 Hertz:  0.02930

Although Newton solver converged in 13 iterations -

Newton iteration 0: r (abs) = 5.380e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 1.539e+00 (tol = 1.000e-10) r (rel) = 2.860e-01 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 1.723e+00 (tol = 1.000e-10) r (rel) = 3.203e-01 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.107e-01 (tol = 1.000e-10) r (rel) = 1.507e-01 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 7.671e-01 (tol = 1.000e-10) r (rel) = 1.426e-01 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 3.512e-01 (tol = 1.000e-10) r (rel) = 6.527e-02 (tol = 1.000e-09)
Newton iteration 6: r (abs) = 4.720e-01 (tol = 1.000e-10) r (rel) = 8.773e-02 (tol = 1.000e-09)
Newton iteration 7: r (abs) = 3.125e-01 (tol = 1.000e-10) r (rel) = 5.808e-02 (tol = 1.000e-09)
Newton iteration 8: r (abs) = 2.763e-01 (tol = 1.000e-10) r (rel) = 5.136e-02 (tol = 1.000e-09)
Newton iteration 9: r (abs) = 1.336e-01 (tol = 1.000e-10) r (rel) = 2.484e-02 (tol = 1.000e-09)
Newton iteration 10: r (abs) = 8.851e-02 (tol = 1.000e-10) r (rel) = 1.645e-02 (tol = 1.000e-09)
Newton iteration 11: r (abs) = 4.648e-02 (tol = 1.000e-10) r (rel) = 8.639e-03 (tol = 1.000e-09)
Newton iteration 12: r (abs) = 9.796e-09 (tol = 1.000e-10) r (rel) = 1.821e-09 (tol = 1.000e-09)
Newton iteration 13: r (abs) = 7.493e-14 (tol = 1.000e-10) r (rel) = 1.393e-14 (tol = 1.000e-09)
Newton solver finished in 13 iterations and 2387 linear solver iterations.
6

13 iterations and an initial residual of 5 can’t be boiled down to FE discretization errors. Either the u you are interpolating is not the actual solution the theory says, or you are not solving the problem you had in mind, or possibly both.

1 Like
7

Thanks for the insightful answer. It is true that I do not know the analytical solution to this problem and I am just trying to reproduce the results of Hertzian contact 3D to 2D problem in terms of matching the values to the Hertz calculation done by the author.

My last question regarding this problem is that are there any other method to check if the simulation is running right or not in the interest of displacement, strain and stress values ?

Once again thanks a lot for saving my time and effort.

8

Those are the only ones that come to my mind (i.e., a person who is not in the scientific field of contact problems), but you better ask your advisor or a senior colleague, since they may have more insight on the specific application.

1 Like
9

Hi,
why do you expect that 2D simulations will give the same results as the 3D case ?
In your setting, the indenter is a cylinder not a sphere, the analytical solution should be different than that of the sphere.

3 Likes
10

Thanks for pointing that out. I completely forgot about that.

PS: It is bit digressing from the question asked but I would also like to thank you for creating great tutorial series Numerical tours of continuum mechanics using FEniCS. Learnt a lot from them. Thanks for the effort.

Also, I have some questions -
Q.1

But indenter used here is obstacle = Expression("-d+(pow(x[0],2))/2/R", d=d, R=R, degree=2) and it seems parabolic to me which is contradictory from your remark.

Q.2 How did you manage to show the indenter in the tutorial ?

image
image1018×700 53.3 KB

Q.3 I’m after elasto-plastic contact and therefore thinking of integrating Hertzian contact with a rigid indenter using a penalty approach and Elasto-plastic analysis of a 2D von Mises material, if you have any suggestions on this approach, please correct me with what the right approach will be ?

11

Thanks. Regarding your questions
Q.1. Yes but in 3D it is a sphere, or more precisely, a paraboloid approximation of a part of a sphere owing to the large radius of curvature compared to the contact area. The analytical Hertz solution of the tutorial corresponds to the sphere case. What I was saying is that in 2D (non axisymmetric) you cannot hope to model the contact of a sphere on a half-space, it will necessarily be that of a cylinder (or parabolic cylinder with a similar approximation) and that you need to compare against the Hertz solution of a cylinder against a half-space.

Q.2 I think I used in Paraview something like Sources > Geometrical Shapes > Sphere and clipped some parts for better visualization

Q3. Yes you can try combining the two approaches

1 Like