I am trying to solve a small-strain linear elasticity problem, where I am attempting to obtain the stress expression from the Helmholtz free energy
$$\psi\,(\varepsilon) = \frac{1}{2} \kappa\, \mathrm{tr}\,(\varepsilon)^2 + \mu \, \mathrm{dev}\,(\varepsilon) : \mathrm{dev}\,(\varepsilon)$$
where the Green-Lagrange strain tensor is defined by
$$\varepsilon = \mathrm{sym}\,(\mathrm{grad}\,(u))$$
When I call solve, the nonlinear solver is activated, and cannot converge to a solution. Nevertheless, when I directly input the stress and not enclose eps in a variable(), I obtain the solution as normal (you can test by uncommenting the relevant lines below). What am I doing wrong?
Another question: When I look into the Fenics Book, I see some finite strain examples that build the problem by differentiating the free energy with respect to a corresponding strain, e.g. "Figure 17.1: UFL implementation of hyperelasticity equations with a Mooney-Rivlin material law.".
In the demos of the most recent version (2016.1), however, I did not come upon any problems that utilizes this technique. In the hyperelasticity demo, for example, the Gateaux derivative of the total potential energy form Pi = psi*dx ... is taken to obtain the bilinear form for the solution, rather than the bilinear form being constructed from the stress and the strain. Is this way more preferred now in Fenics programs?
from __future__ import print_function
from dolfin import *
# Sub domain for clamp at left end
def left(x, on_boundary):
return x[0] < 0.001 and on_boundary
# Sub domain for rotation at right end
def right(x, on_boundary):
return x[0] > 0.99 and on_boundary
# Load mesh and define function space
mesh = UnitCubeMesh(15, 15, 15)
# Define function space
V = VectorFunctionSpace(mesh, "CG", 1)
scalar_space = FunctionSpace(mesh,'CG',1)
tensor_space = TensorFunctionSpace(mesh, "CG", 1)
# Test and trial functions
v = TestFunction(V)
u = TrialFunction(V)
u_ = Function(V)
kappa = 10.
mu = 10.
# External forces (body and applied tractions
f = Constant((10.0, 0.0, 0.0))
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
force_boundary = AutoSubDomain(right)
force_boundary.mark(boundary_subdomains, 3)
# Define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)
# Set up boundary condition at left end
zero = Constant((0.0, 0.0, 0.0))
bcs = DirichletBC(V, zero, left)
eps = variable(sym(grad(u)))
psi = kappa / 2. * tr(eps) ** 2 + mu * inner(dev(eps), dev(eps))
sigma = diff(psi, eps)
# Uncommenting this yields the correct solution
# eps = sym(grad(u))
# sigma = kappa * tr(eps) * Identity(len(u)) + 2.0 * mu * dev(eps)
a = inner(sigma, sym(grad(v)))*dx
L = inner(f, v)*dss(3)
F = a - L
F = action(F, u_)
solve(F == 0, u_, bcs, solver_parameters={"newton_solver": {"linear_solver": "gmres"}},)
# Output displacement
xdmf_file = XDMFFile(mesh.mpi_comm(), "nonlinear_elasticity.xdmf")
u_.rename("u", "displacement")
xdmf_file.write(u_, 0.0)
