Hello people,
I am working on contact mechanics between hyperelastic and rigid bodies.
As the rigid bodies are always plane, and the hyperelastic bodies can have various shapes, the displacement caused by contact, is realized by a minimum displacement for the whole mesh. The Solver computes the mechanical field, when every meshpoint reached at least it’s minimum displacement.
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
# Check that DOLFIN is configured with PETSc and CGAL
if not has_petsc():
print ('DOLFIN must be compiled with PETSc to run this demo.')
exit(0)
#-----------------------------------------------------------------------------2D-domain and mesh
domain = Rectangle(Point(10, 0), Point(12, 2))
# Create mesh
mesh = generate_mesh(domain, 20)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
x = SpatialCoordinate(mesh) #Return symbolic physical coordinates for given mesh.
def eps(v):
return sym(as_tensor([[v[0].dx(0) , 0 , 0.5*(v[1].dx(0)+v[0].dx(1)) ],
[0 , v[0]/x[0] , 0 ],
[0.5*(v[1].dx(0)+v[0].dx(1)) , 0 , v[1].dx(1) ]]))
# Define functions
du = TrialFunction(V) # Incremetal displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, 0.0)) # Body force per unit volume
# Kinematics
I = Identity(3) # Identity tensor
F = I + eps(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 2) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx - dot(B, u)*dx
# Compute first variation of Pi (directional derivative about u in the
# direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of FCa
J = derivative(F, u, du)
# Symmetry condition (to block rigid body rotations)
tol = mesh.hmin()
def symmetry_line(x):
return abs(x[1]) < DOLFIN_EPS
bc = DirichletBC(V.sub(1), 0., symmetry_line,method="pointwise")
zmax=4
rmin=10.05
zmin=0
rmax=12.5+DOLFIN_EPS
# The displacement u must be such that the current configuration x+u
# remains in the box [xmin,xmax] x [umin,ymax]
constraint_max = Expression(("rmax - x[0]","zmax - x[1]"), rmax=rmax, zmax=zmax, degree=1)
umax = interpolate(constraint_max, V)
constraint_min = Expression(("rmin - x[0]","zmin - x[1]"), rmin=rmin, zmin=zmin, degree=1)
umin = interpolate(constraint_min, V)
##Set-up and calling of solver realized with common PETScSNESSolver.cpp.
# Define the solver parameters
snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver" : {"linear_solver" : "cg",
"maximum_iterations": 2000,
"report": True,
"error_on_nonconvergence": True,
}
}
# Set up the non-linear problem
problem = NonlinearVariationalProblem(F, u, bc, J=J)
problem.set_bounds(umin, umax)
# Set up the non-linear solver
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
#info(solver.parameters, True)
# Solve the problem
(iter, converged) = solver.solve()
# Check for convergence
if not converged:
warning("This demo is a complex nonlinear problem. Convergence is not guaranteed when modifying some parameters or using PETSC 3.2.")
# Save solution in VTK format
file = File("displacement.pvd")
file << u
# plot the current configuration
plot(mesh)
plot(u, mode="displacement", title="Displacement field")
plt.grid(b=None, which='major', axis='both')
plt.show()
For the given example, this works perfect up to a displacement of 0.05 m on the left side as you can see in the following picture.
With 0.06 m displacement, the mimimum displacement in absolute values gets bigger then the space between the outer meshpoint on the left side and first inner meshpoint. Once this point is reached, the simulation diverges. It is no matter weather I increase the displacement or decrease the mesh size.
The Error I get is
Solving nonlinear variational problem.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
Traceback (most recent call last):
File “…/exmple.py”, line 96, in
(iter, converged) = solver.solve()
RuntimeError:*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with PETScSNESSolver.
*** Reason: Solver did not converge.
*** Where: This error was encountered inside PETScSNESSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset: 640ef70247ae212e16ae9f24fc9bc603b506f78a
*** -------------------------------------------------------------------------
The reason for the divergation should be a division by zero, but I don’t really get the reason for this division. I also tried with an circular domain, where with nearly every amount of displacement, there are points on the border of the domain getting a minimum displacement near zero.
It works until a point, which is not part of the border changes the direction of its minimum displacement. But it is not important in which direction it changes…
On the left side, the error occurs with the change from a negative to the positive direction and on the right side the same error occurs with a change from the positive to negative direction of its minimum displacement.
Does anyone have an advice, how to fix this or where to look for a solution?
