Convergence/Linearization Problem in Magnetoelasticity

Question

Hi Guys,

I'm currently working on nonlinear magnetoelasticity of magneto-sensitive elastomers.
The governing equations are the balance of linear momentum and the Gauss's law for magnetism, while the displacement field and the magnetic scalar potential are the unknowns.
I want to use a constitutive model which takes into account the effects of the magnetic saturation, but I'm facing a problem trying to solve the nonlinear variational problem.
It seems like Fenics cannot properly linearize a tanh expression in the definition of the magnetic induction vector. If i remove the tanh(...) term the code runs just fine.
But this tanh(...) term is essential to model the saturation behavior.

I have attached a minimum running example to reproduce the error.

from __future__ import print_function
from fenics import *
import numpy as np

# Turn on C++ optimization of generated code
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
parameters["form_compiler"]["quadrature_degree"] = 2
# Select the newer code generation backend 'uflacs'
parameters["form_compiler"]["representation"] = "uflacs"

# Material parameters
G = 0.1                     # N/mm^2        shear modulus
nu = 0.4                    # -             Poisson's ratio

mu0 = 4.0*pi * 1e-13                # N/mA^2        Magnetic Permeability Free Space
g1 = -1.0 * 1e-12                   # N/mA^2        Materialparameter
m0 = 0.4998 * 1e-6          # N/(mA*mm)
m1 = 309339.5               # mA/mm
c1 = 1250                   # -
c0 = 624.5*4.0*pi * 1e-13           # N/mA^2

# Create mesh
mesh = BoxMesh(Point(0, 0, 0), Point(10, 10, 10), 10, 10, 10)


# Build function spaces (Taylor-Hood)
U = VectorElement('P', tetrahedron, 2)
M = FiniteElement('P', tetrahedron, 1)
element = MixedElement([U, M])
W = FunctionSpace(mesh, element)



# Define boundary condition
tol = 1E-14

def left(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

def right(x, on_boundary):
    return on_boundary and abs(x[0]-10) < tol   

bc1 = DirichletBC(W.sub(0), Constant((0, 0, 0)), left)
bc2 = DirichletBC(W.sub(1), Constant(20*1e3), right)
bc3 = DirichletBC(W.sub(1), Constant(0), left)

# Collect boundary conditions
bcs = [bc1, bc2, bc3]


# Define trial and test functions
(v, q)  = TestFunctions(W)
dw = TrialFunction(W)
w = Function(W)

# Split system functions to access components
du, dm =split(dw)
u, m = split(w)

# Define stress vector
Tr = Constant((0.0, 0.0, 0.0))

# Define constants
G = Constant(G)
beta = Constant((2*nu)/(1-2*nu))
mu0 = Constant(mu0)
m0 = Constant(m0)
m1 = Constant(m1)
g1 = Constant(g1)
c1 = Constant(c1)
c0 = Constant(c0)

# Define strain measures

# Identity tensor
def I(u):
   return Identity(len(u))

# deformation gradient  
def F(u):
   return Identity(len(u)) + grad(u)

# right Cauchy-Green deformation tensor 
def C(u):
    return F(u).T*F(u)

# Determinate deformation gradient  
def J(u):
    return det(F(u))

# 1.Invariant   
def I1(u):
    return inner(C(u), I(u))    

# Define magnetic measures  

# Magnetic Field Strength
def H(m):
    return -grad(m)

# 4. Invariant  
def I4(m):
    return dot(H(m), H(m))

# Square root 4. Invariant  
def SI4(m):
    return (I4(m))**(0.5)

# Outer Product of Magnetic Field Strength
def HH(m):
    return outer(H(m), H(m))    

# Magnetic Induction Vector
def B(u, m):
    return - (I1(u) - 3.0)*g1*H(m) + 2.0*c0*H(m) - 2.0*c1*mu0*C(u)*H(m) +   tanh((SI4(m)/m1))*(m0/SI4(m))*H(m)


# Define stress measures    

# 1.PK stress tensor
def PK(u, m):
    return G*(F(u) - ((J(u))**(-beta))*inv(F(u).T)) + g1*I4(m)*F(u) + 2.0*c1*mu0*F(u)*HH(m) 

# Define nonlinear problem

R1 = inner(PK(u, m), grad(v))*dx - dot(Tr, v)*ds
R3 = dot(B(u, m), grad(q))*dx

R= R1 + R3
JJ = derivative(R, w, dw)

problem = NonlinearVariationalProblem(R, w, bcs, J=JJ)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-9
prm['newton_solver']['relative_tolerance'] = 1E-25
prm['newton_solver']['maximum_iterations'] = 15
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'ilu'

# Create VTK files for visualization output
file_u = File("Solution/SolutionU.pvd");
file_phi = File("Solution/SolutionPhi.pvd");

# Solve variational problem for timestep
solver.solve()

# Save solution in VTK format
u_1, m_1 = w.split(deepcopy=True)
file_u << (u_1)
file_phi << (m_1)

I will be grateful for any help.

Comments

Hi, what happens if you replace tanh with some polynomial approximation?

---

Hi,

when i replace the tanh with a taylor series approximation the problem still remains.
I'm not able to solve the nonlinear variational problem even with homogeneous Dirichlet BCs.
Fenics also seems to have problems to properly handle the square root expression:

# 4. Invariant  
def I4(m):
    return dot(H(m), H(m))

# Square root 4. Invariant  
def SI4(m):
    return (I4(m))**(0.5)