Hi,
I am trying to solve a finite elasticity problem using a Neo-hookean material in which a bar gets stiffer with time and has inhomogeneous stiffness in the domain (in the x-direction). I am applying constraints (1) ux = 0 on the left and right side of the bar (x = 0 and x = 10), (2) uy = 0 on the y = 1 plane and (3) uz = 0 on the z = 1 plane.
I am applying the inhomogenity by modifying the stiffness at the gauss points by defining a the stiffness as a "Quadrature" function and then modifying contents of the array (see the variable "Tact").
For some reason, the solution does not converged when I use quadratic element and when the inhomogeneity gets large (given by the parameter "stiffness_diff"). However, the solution converge without any issue when a linear element is used (even with very inhomogeneous stiffness)
Any idea of what is happening? Have attached the code - any help is greatly appreciated! Thanks!!
from dolfin import *
import numpy as np
import math as math
quad_degree = 2
# Optimization options for the form compiler
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# Create mesh and define function space
L = 10.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(L, 1.0, 1.0), 10, 2, 2)
dx = dolfin.dx(mesh, metadata = {"integration_order":quad_degree})
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = L)
top = CompiledSubDomain("near(x[2], side) && on_boundary", side = 1.0)
bot = CompiledSubDomain("near(x[2], side) && on_boundary", side = 0.0)
front = CompiledSubDomain("near(x[1], side) && on_boundary", side = 1.0)
back = CompiledSubDomain("near(x[1], side) && on_boundary", side = 0.0)
facetboundaries = FacetFunction("size_t", mesh)
facetboundaries.set_all(0)
right.mark(facetboundaries, 2)
left.mark(facetboundaries, 1)
top.mark(facetboundaries, 3)
front.mark(facetboundaries, 4)
W = VectorFunctionSpace(mesh, "Lagrange", 1)
Quads = FunctionSpace(mesh, "Quadrature", quad_degree)
bcright = DirichletBC(W.sub(0), Constant(("0")), facetboundaries, 2)
bctop = DirichletBC(W.sub(2), Constant(("0")), facetboundaries, 3)
bcfront = DirichletBC(W.sub(1), Constant(("0")), facetboundaries, 4)
bcleft = DirichletBC(W.sub(0), Constant(("0")), facetboundaries, 1)
bcs = [bcleft, bctop, bcfront, bcright]
# Define functions
du = TrialFunction(W) # Incremental displacement
v = TestFunction(W) # Test function
u = Function(W) # Displacement from previous iteration
Tact = Function(Quads)
d = u.geometric_dimension()
I = Identity(d)
F = I + grad(u)
J = det(F)
Cmat = F.T*F
Fdev = J**(-1.0/3.0)*F
Cdev = Fdev.T*Fdev
Kappa = Constant(1.0e5);
f0 = Constant(("1.0", "0.0", "0.0"))
W1 = Constant(1000.0)*(tr(Cdev) - Constant(3.0)) + Kappa/Constant(2.0)*(J - 1.0)**2.0
W2 = Tact/2.0 * (Cmat[0,0] - 1.0)
Wtotal = (W1 + W2)
# Total potential energy
Pi = Wtotal*dx
# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of F
Jac = derivative(F, u, du)
t_niter = 40
r_threshold = 1e-9
# Solve variational problem
dt = 2.0
t = 0.0
Tmax = 80e3
Tend = 100.0;
ngpt = len(Tact.vector().array()[:])
t_niter = 40
r_threshold = 1e-9
# Varying the inhomogenity of the stiffness
stiffness_diff = 50000
disp_vtk_file = File("nonlinearelasticity_block_disp_diag.pvd")
while(t < Tend):
# Modifying stiffness at the Gauss points
Tact.vector()[:] = np.ones(ngpt)*sin(t/Tend*math.pi)*80e3 + np.random.rand(1,ngpt)[0]*stiffness_diff
print "Max stiffness = ", max(Tact.vector().array()), " Min stiffness = ", min(Tact.vector().array())
solve(F == 0, u, bcs, J=Jac,form_compiler_parameters=ffc_options)
t = t + dt
disp_vtk_file << (u, t)
# Plot and hold solution
plot(u, interactive=True, mode="displacement")
