Hello every one hope you are fine.
I have the next question to ask to you about the formulation in multiphenics about non linear contact problem using unbiased Nitsche method.
The idea is to simulate the contact between two bodies considering Coulomb Friction Law and small strains. The bodies are assumed to be already in contact (gap function =0 ).
The formula used was taken from Finite Element Approximation of Contact and Friction in Elasticity (Chapter 10, page 245) from Franz Chouly, Patrick Hild, Yves Renard:
I am getting no convergence using this method:
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
Traceback (most recent call last):
File "/home/user/fenics_work/MODELO/multiphenics_nonlinear/example_two_rectangles_multiphenics.py", line 223, in <module>
solver.solve()
File "/home/user/.local/lib/python3.12/site-packages/multiphenics/nls/block_petsc_snes_solver.py", line 27, in solve
no_iterations, convergence_flag = PETScSNESSolver.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
***
*** https://fenicsproject.discourse.group/
***
*** 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.2.0.13.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------
The MWE uses legacy multiphenics inspired in the tutorials:
from dolfin import *
from multiphenics import *
# Solver parameters
snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver": {"linear_solver": "mumps",
"maximum_iterations": 20,
"report": True,
"error_on_nonconvergence": True}}
parameters["allow_extrapolation"] = True
# RectangleMesh(Point(-1, -1), Point(1.0, 1.0), 100, 50)
mesh = RectangleMesh(Point(-1.0, -2.0), Point(1.0, 2.0), 50, 100)
bdry = MeshFunction('size_t', mesh, 1, 0)
subdomains = MeshFunction("size_t", mesh, 2, 0)
dim = mesh.topology().dim()
# Surface subdomains
class Omega2(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.
class Omega1(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
# Boundary subdomains
class Top(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 2.0) and on_boundary)
class TRight(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 1.0) and between(x[1], (0.0, 2.0)) and on_boundary)
class TLeft(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], -1.0) and between(x[1], (0.0, 2.0)) and on_boundary)
class BRight(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 1.0) and between(x[1], (-2.0, 0.0)) and on_boundary)
class BLeft(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], -1.0) and between(x[1], (-2.0, 0.0)) and on_boundary)
class Bot(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], -2.0) and on_boundary)
# parametros físicos
E1 = 1.1E11 # ACERO GPa
E2 = 1.1E11
nu = 0.35
top = 1
bottom = 2
interface = 3
tright = 4
bright = 5
tleft = 6
bleft = 7
O1 = 1
O2 = 2
alfa = 20
Omega2().mark(subdomains, O2)
Omega1().mark(subdomains, O1)
Interface().mark(bdry, interface)
Top().mark(bdry, top)
TRight().mark(bdry, tright)
TLeft().mark(bdry, tleft)
BRight().mark(bdry, bright)
BLeft().mark(bdry, bleft)
Bot().mark(bdry, bottom)
def pos_f(u):
return 0.5*(abs(u) + u)
def neg_f(u):
return 0.5*(abs(u) - u)
def tnv(v, normal):
return v - dot(v, normal)*normal
def epsilon(u):
return sym(grad(u))
def D(E, nu, u):
mu = E/(2*(1+nu))
lamb = E*nu/((1-2*nu)*(1+nu))
return 2*mu*epsilon(u)+lamb*tr(epsilon(u))*Identity(dim)
def Pin(v, normal, E, nu):
jumpu = jump(v, normal) # scalar
sigmn = avg(dot(D(E, nu, v)*normal, normal)) # scalar
return gamN*jumpu - sigmn
def Pit(v, normal, E, nu):
jumpu = jump(tnv(v, normal)) # vector
sig_t = avg(tnv(D(E, nu, v)*normal, normal)) # vector
return gamN*jumpu - sig_t
def Pit_fric(v, normal, E, nu, fric):
# Coulomb friction law: - fric*|sig_n(u) | * u_t/|u_t|
jumpu = jump(tnv(v, normal)) # vector
sigmn = avg(abs(dot(D(E, nu, v)*normal, normal)))
v_t = avg(tnv(v, normal))
Coulomb = -fric*sigmn*v_t/sqrt(dot(v_t, v_t))
return gamN*jumpu - Coulomb
# print(domains.array())
omega1 = MeshRestriction(mesh, Omega1())
omega2 = MeshRestriction(mesh, Omega2())
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=bdry)
dS = Measure("dS", domain=mesh, subdomain_data=bdry)
Pv = VectorFunctionSpace(mesh, "CG", 1)
# u1 = u on Bottom domain
# u2 = u on Top domain
# u1, u2
Hh = BlockFunctionSpace([Pv, Pv],
restrict=[omega1, omega2])
sol = BlockFunction(Hh)
u1, u2, = block_split(sol)
test = BlockTestFunction(Hh)
v1, v2 = block_split(test)
dsol = BlockTrialFunction(Hh) # To use with Jacobian
bcUtleft = DirichletBC(Hh.sub(1), Constant((0, 0)), bdry, tleft)
bcUbleft = DirichletBC(Hh.sub(0), Constant((0, 0)), bdry, bleft)
bcs = BlockDirichletBC([bcUtleft, bcUbleft])
h = CellDiameter(mesh)
h_avg = (h("+") + h("-")) / 2.0
n = FacetNormal(mesh)
gamN = alfa/h_avg
# Define bilinear form tanto en el domino 1 como en el dominio 2
aO1 = inner(D(E1, nu, u1), epsilon(v1)) * dx(O1)
aO2 = inner(D(E2, nu, u2), epsilon(v2)) * dx(O2)
# Add Nitsche terms
############
aO1 += -0.5*inner(1/gamN*avg(D(E1, nu, u1)*n),
avg(D(E1, nu, v1)*n))*dS(interface)
aO2 += -0.5*inner(1/gamN*avg(D(E2, nu, u2)*n),
avg(D(E2, nu, v2)*n))*dS(interface)
############
aO1 += 0.5*inner(1/gamN*pos_f(Pin(u1, n, E1, nu)),
Pin(v1, n, E1, nu))*dS(interface)
aO2 += 0.5*inner(1/gamN*pos_f(Pin(u2, n, E2, nu)),
Pin(v2, n, E2, nu))*dS(interface)
############
# Coulomb friction law: - Fric*|sig_n(u)|* u_t/|u_t|
fric1 = 1E-10
fric2 = 1E-10
aO1 += 0.5 * inner(1/gamN*Pit_fric(u1, n, E1, nu, fric1),
Pit(v1, n, E1, nu))*dS(interface)
aO2 += 0.5 * inner(1/gamN*Pit_fric(u2, n, E2, nu, fric2),
Pit(v2, n, E2, nu))*dS(interface)
# Define linear form
loadf = Constant((10, 0))
LO1 = inner(loadf, v1) * ds(bright)
LO2 = inner(loadf, v2) * ds(tright)
F = [aO1 - LO1, aO2 - LO2]
J = block_derivative(F, sol, dsol)
# Solve
problem = BlockNonlinearProblem(F, sol, bcs, J)
solver = BlockPETScSNESSolver(problem)
solver.parameters.update(snes_solver_parameters["snes_solver"])
solver.solve()
u2h, u1h = block_split(sol)
Can anyone please could help me to solve this problem. I appreciate your answers.