Hey!
I try to solve the hyperelastic uniaxial stress example as a nonlinearproblem in a mathematically correct form. It shows the same strain as in the analytical proove in an order of 1e-4 for all pressures. Here, the code is given for incompressible solid, means nu=.5
For stability, a MixedElement function space was used with the hydrostatic pressure, that should result in 0.
While setting the external pressure up to 5e8 the problem converges. Even while plotting strain over different stresses, the curve looks as it should. With higher pressure, resulting in a strain > 3.6 the problem not converges anymore.
I hope, anyone can give me a hint, why the problem results in the following and a solution cannot be found.
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 7e99df3564aea9dc05729ddb64498ae47b6bc15d
the code:
from dolfin import *
boundary_vals = {
"hold": 1,
"free": 2,
"pressure": 3,
"mirrorY": 4, # hold y
"mirrorX": 5, # hold x
"mirrorZ": 6, # hold Z
"cells": 0 # internal cells
}
def solve_elasticity(
p: float = 4.0e6,
mu: float = 362.2e3,
K: float = 0.,
nu: float = .5,
resolution: int = 10,
fem_rel_tol: float = 1e-6,
) -> int:
# === input parameters after nondimensionalisation
p_inner = -p/mu
# assuming l, rho = 1., 1.
dt = Constant(mu**(1/2))
mu = 1.
# mesh and faccet functions
mesh = UnitCubeMesh(resolution, resolution, resolution)
facet_function = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
x0 = AutoSubDomain(lambda x: near(x[0], 0))
x1 = AutoSubDomain(lambda x: near(x[0], 1))
y0 = AutoSubDomain(lambda x: near(x[1], 0))
x0.mark(facet_function, boundary_vals['mirrorX'])
y0.mark(facet_function, boundary_vals['mirrorY'])
z0 = AutoSubDomain(lambda x: near(x[2], 0))
z0.mark(facet_function, boundary_vals['mirrorZ'])
x1.mark(facet_function, boundary_vals['pressure'])
mesh = facet_function.mesh()
# ==========================================================================
vtk_resuls_file = File('test_u.pvd')
# Build function space
element_v = VectorElement("Lagrange", mesh.ufl_cell(), 1)
element_s = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
mixed_element = MixedElement([element_v, element_v, element_s])
W = FunctionSpace(mesh, mixed_element)
# ## == Prepare Dirichlet BCs =============================================
bcmX = DirichletBC(
W.sub(0).sub(0), Constant(0.0), facet_function,
boundary_vals["mirrorX"])
bcmY = DirichletBC(
W.sub(0).sub(1), Constant(0.0), facet_function,
boundary_vals["mirrorY"])
bcmZ = DirichletBC(
W.sub(0).sub(2), Constant(0.0), facet_function,
boundary_vals["mirrorZ"])
bcs = [bcmY, bcmX, bcmZ]
# dimensionalized output Functions
w_write = Function(W)
write_u, v_write, p_write = w_write.split()
write_u.rename("u", "displacement")
dx = Measure("dx")
ds = Measure("ds", subdomain_data=facet_function, subdomain_id=2)
# Limit quadrature degree
dx = dx(degree=4)
ds = ds(degree=4)
I_ = Identity(W.mesh().geometry().dim())
w = Function(W)
u, v, p = split(w)
w0 = Function(W)
u0, v0, p0 = split(w0)
_u, _v, _p = TestFunctions(W)
dw = TrialFunction(W)
P, pp, T, S = stress(u, p, I_, K, mu)
P0, pp0, T0, S0 = stress(u0, p0, I_, K, mu)
F1a = - inner(v, _u) * dx
F2a = inner(P, grad(_v)) * dx
Fp = pp * _p * dx
Fv = (1.0 / dt) * inner(u - u0, _u) * dx + F1a
Fa = (1.0 / dt) * inner(v - v0, _v) * dx + F2a
# Traction at Neumann boundary
dss = ds(subdomain_data=facet_function)
F = I_ + grad(u)
facet_normal = FacetNormal(mesh)
bF_magnitude = Constant(0.0)
bF = det(F) * dot(bF_magnitude * facet_normal, inv(F)) # 1PK on boundary
FF = inner(bF, _v) * dss(boundary_vals['pressure'])
F_ = Fv + Fa + Fp + FF
Jac = derivative(F_, w, dw)
# ## == Initialize solver =================================================
problem = NonlinearVariationalProblem(F_, w, bcs=bcs, J=Jac)
solver = NonlinearVariationalSolver(problem)
solver.parameters['symmetric'] = False # default value
solver.parameters['newton_solver']['relative_tolerance'] = fem_rel_tol
solver.parameters['newton_solver']['maximum_iterations'] = 50
solver.parameters['newton_solver']['krylov_solver']['relative_tolerance'] = fem_rel_tol
solver.parameters['newton_solver']['krylov_solver']['absolute_tolerance'] = fem_rel_tol
solver.parameters['newton_solver']['linear_solver'] = 'superlu_dist'
# Extract solution components
u, v, p = w.split()
bF_magnitude.assign(Constant(p_inner))
# w0.assign(w)
solver.solve()
write_u.assign(project(
u, solver_type="cg", preconditioner_type="amg"))
vtk_resuls_file << (write_u)
def stress(u, p, I_, K, mu):
"""Returns Piola-Kirchhoff stresses and hydrostatic pressure
for given u, p."""
F = I_ + grad(u)
J = det(F)
B = F * F.T # left Cauchy–Green deformation tensor
sigma = -p*I_ + mu * (B) # cauchy stress
T = J**(-1) * sigma
P = J * T * inv(F).T # 1st Piola-Kirchhoff stress
S = J * dot(dot(inv(F), T), inv(F).T) # 2nd Piola-Kirchhoff stress
# hydrostatic pressure from \partial \psi / \partial J
hydrostat_p = - Constant(1.0) + J
return P, hydrostat_p, T, S
if __name__ == '__main__':
solve_elasticity(p=5e8) # works
solve_elasticity(p=6e8) # not converges
using SNES solver, the same happens.
solver.parameters['nonlinear_solver'] = 'snes'
solver.parameters['snes_solver']['linear_solver'] = 'superlu_dist'
