No convergence for nonlinear elastic problem with spectral decomposition constitutive equation

dolfinx
1

Dear community,

I am solving a 2D nonlinear elastic problem using dolfinx. The spectral decomposition is employed to decompose the stress tensor into positive and negative parts. However, the issue is that when I setup the nonlinear Newton solver, the result cannot converge. If I switch the constitutive equation into the commonly used linear ones, I can solve the problem again. So I guess the problem is about the spectral decomposition part. Any ideas?

The minimal problem is attached in the following.

import dolfinx
from mpi4py import MPI 
from petsc4py import PETSc
import numpy as np 
import ufl
from ufl import sym, grad, tr, Identity, dev, inner, dx, ds, conditional, lt,gt,inv, derivative, as_tensor, sqrt
import time, os, sympy

n = 16
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n, n)

# define space
V = dolfinx.fem.FunctionSpace(mesh, ('CG', 1))
W = dolfinx.fem.VectorFunctionSpace(mesh, ('CG', 1))
WW = dolfinx.fem.FunctionSpace(mesh, ('CG', 1))
p, q = ufl.TrialFunction(V), ufl.TestFunction(V)
u, du, v = dolfinx.fem.Function(W), ufl.TrialFunction(W), ufl.TestFunction(W)

# the material parameters
lmbda = 121.1538e3
mu = 80.7692e3
K = lmbda + mu

# ==================== define the eigen decomposition methods ============================ #
# Compute the symbolic expression for eigenvalues by sympy
T = sympy.Matrix(2, 2, lambda i, j: sympy.Symbol('T[%d, %d]' % (i, j), symmetric =True, real=True))
eig_expr = T.eigenvects() 

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv = sympy.Matrix(eigv) # Column vector
eigw = sympy.Matrix(eigw) # Row has the elements
eigw = sympy.Matrix([[eigw[0],eigw[2]],[eigw[1],eigw[3]]])
eigv = sympy.Matrix([[eigv[0], 0.0],[0.0, eigv[1]]])

eigv_expr = list(map(str, eigv))
eigw_expr = list(map(str, eigw))

# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return list(map(eval, eigv_expr))
# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return list(map(eval, eigw_expr))

def epsilon(u):
    return sym(grad(u))

# Positive strain
def strn_p(u):
    T = sym(grad(u))
    v00, v01, v10, v11 = eigv(T)
    v00 = conditional(gt(v00,0.0),v00,0.0)
    v11 = conditional(gt(v11,0.0),v11,0.0)
    w00, w01, w10, w11 = eigw(T)
    wp = ([w00,w01],[w10,w11])
    wp = as_tensor(wp)
    vp = ([v00,v01],[v10,v11])
    vp = as_tensor(vp)
    return wp*vp*inv(wp)
    
# Negative strain
def strn_n(u):
    T = sym(grad(u))
    v00, v01, v10, v11 = eigv(T)
    v00 = conditional(lt(v00,0.0),v00,0.0)
    v11 = conditional(lt(v11,0.0),v11,0.0)
    w00, w01, w10, w11 = eigw(T)
    wn = ([w00,w01],[w10,w11])
    wn = as_tensor(wn)
    vn = ([v00,v01],[v10,v11])
    vn = as_tensor(vn)
    return wn*vn*inv(wn)

# Positive stress
def sigma_pos(u):
    brack_pos = (0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))
    sigma_pos = lmbda*brack_pos*Identity(2)+2*mu*strn_p(u)
    return sigma_pos

# Negative stress
def sigma_neg(u):
    brack_neg = (0.5*(-tr(epsilon(u))+abs(tr(epsilon(u)))))
    sigma_neg = -lmbda*brack_neg*Identity(2)+2*mu*strn_n(u)
    return sigma_neg

# select boundaries
def top(x):
    return np.isclose(x[1], 1.0)
def bottom(x):
    return np.isclose(x[1], 0)
class MyExpression:
    def __init__(self):
        self.t = 1e-2
    
    def eval(self):
        return np.array((self.t), dtype=PETSc.ScalarType)

# set bc
fdim = mesh.topology.dim - 1
top_boundary = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top)
top_load_dofs = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, top_boundary)
top_disp_dofs = dolfinx.fem.locate_dofs_topological(W.sub(1), fdim, top_boundary)
load = MyExpression()
top_bc_load = dolfinx.fem.dirichletbc(load.eval(), top_load_dofs, W.sub(0))
top_bc_disp = dolfinx.fem.dirichletbc(np.array((0.0), dtype=PETSc.ScalarType), top_disp_dofs, W.sub(1))

bottom_boundary = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom)
bottom_dofs = dolfinx.fem.locate_dofs_topological(W, fdim, bottom_boundary)
bottom_u = np.array((0,0), dtype=PETSc.ScalarType)
bottom_bc = dolfinx.fem.dirichletbc(bottom_u, bottom_dofs, W)

bc_u = [top_bc_load, top_bc_disp, bottom_bc]

# variational form
unew, uold = dolfinx.fem.Function(W), dolfinx.fem.Function(W)
F_u = inner(grad(v), sigma_pos(unew)+sigma_neg(unew))*dx
J_u = derivative(F_u, unew)
problem_u = dolfinx.fem.petsc.NonlinearProblem(F_u, unew, bc_u, J_u)
solver_u = dolfinx.nls.petsc.NewtonSolver(mesh.comm, problem_u)

# solver
solver_u.atol = 1e-10
solver_u.rtol = 1e-9
num_its, converged = solver_u.solve(unew)
print(np.linalg.norm(unew.x.array))
1 Like
3

Hello, I have the same problem, and I was wondering if you have found the solution. I think it is the sqrt or a power less than 1 in the functional, because if I change that it solves.

Best regards, Sara