Hello,
I'm quite new to FEniCS and so far I'm very impressed, especially with the documentation and ease of use of the software.
As part of my PhD project I'm trying to develop a model of soft tissue swelling. The material is hyperelastic and the swelling is modelled using a multiplicative split of the deformation gradient into swelling and elastic components. At this stage the swelling is defined by a simple function of time using an "Expression" (later it may depend on other variables or an external subroutine).
The model converges quickly when only standard boundary conditions such as forces and displacements are applied. However, when I apply the swelling deformation the model fails to converge unless pseudo-timestepping with relatively small time steps is used. Therefore, I suspect that there is either a problem with the linearisation of the swelling part of the model (the Jacobian) or that the swelling introduces significant nonlinearities that can only be overcome by taking small time steps.
My question is: do I need to modify the Jacobian to account for the Expression that is changing with time or is this done automatically by FEniCS? In other words, is there any way that I can improve the convergence of the solution.
Here is an example of the code I have so far:
from dolfin import *
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True,
"quadrature_degree": 1}
# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume
T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary
# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
# Define Dirichlet boundary
l = Expression(("0.0", "0.0", "0.0"))
r = Expression(("0.0", "0.0", "0.0"))
bcl = DirichletBC(V, l, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
# Swelling
Fs = Expression(("1 + scale*t"), t=0.0, scale = 0.5)
Fsinv = 1 / Fs
Js = Fs ** 3
# Kinematics
F = I + grad(u) # Deformation gradient
Fe = F * Fsinv # Elastic deformation gradient
Ce = Fe.T * Fe # Elastic right Cauchy-Green tensor
# Invariants of deformation tensors
Ice = variable(tr(Ce))
Je = variable(det(Fe))
# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ice - 3) - mu*ln(Je) + (lmbda/2)*(ln(Je))**2
# Elastic second Piola-Kirchhoff stress tensor
Se = 2*diff(psi, Ice)*I + Je*Je*diff(psi, Je)*inv(Ce)
# Total second Piola-Kirchhoff stress
S = Js * Fsinv * Se * Fsinv
# First Piola-Kirchhoff stress tensor
P = F*S
# Variational problem
F = inner(P, grad(v))*dx - inner(B, v)*dx - inner(T, v)*ds
# Compute Jacobian of F
J = derivative(F, u, du)
# Timestepping
t = 0.0; dt = 0.1; T = 1.0
inc = 0
file = File("displacement.pvd");
file << u
while t < T:
t += dt
inc += 1
# Update swelling deformation
Fs.t = t
# Solve variational problem
solve(F == 0, u, bcs, J=J,
form_compiler_parameters=ffc_options)
# Save solution in VTK format
file << u;
# Plot and hold solution
plot(u, mode = "displacement", interactive = True)
