Hyperelastic model with two different materials

General
1

hii all,

im trying to use the hyperplastic demo code and apply it on a model with 2 different materials and it seems to be failing on the boundary level.
as of now my code consist only the boundary and subdomain definition for the two materials with the material properties definition and it still failing.
any idea why?

import matplotlib.pyplot as plt
from dolfin import *
import fenics as fe
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

E, nu = 10.0, 0.3
# Create mesh and define function space
length, width, height = 1.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 10  # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 5, 5)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= length/2 + DOLFIN_EPS


class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > length/2 - DOLFIN_EPS


subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)


# Boundary conditions

class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, DOLFIN_EPS)


class right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1, DOLFIN_EPS)


facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
facets.set_all(0)
left().mark(facets, 1)
right().mark(facets, 2)

bc1 = DirichletBC(V.sub(0), Constant(0), facets, 1)
bc2 = DirichletBC(V.sub(1), Constant(0), facets, 1)
bc3 = DirichletBC(V.sub(2), Constant(0), facets, 1)
bc4 = DirichletBC(V.sub(0), Constant(5), facets, 2)
bc5 = DirichletBC(V.sub(1), Constant(5), facets, 2)
bc6 = DirichletBC(V.sub(2), Constant(5), facets, 2)
bcs = [bc1, bc2, bc3, bc4, bc5, bc6]


mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Define functions
du = TrialFunction(V)            # Incremental displacement
v = TestFunction(V)             # Test function
u = Function(V)                 # Displacement from previous iteration
B = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T = Constant((0.0, 0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = len(u)
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcs, J=J)

# Save solution in VTK format
file = File("displacement.pvd");
file << u;

# Plot solution
plot(u)
plt.show()

this is the error
Traceback (most recent call last):
File “/home/mirialex/PycharmProjects/fewmaterials/hyperwithfewmaterials.py”, line 91, in
solve(F == 0, u, bcs, J=J)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 233, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 314, in _solve_varproblem
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 NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------

Process finished with exit code 1

2

Keep in mind that hyperelasticity is a nonlinear problem which you solve with a Newton method which might not always converge. You seem to apply very large displacements on the boundary, first try to apply much lower displacements for improving the convergence behavior.

1 Like