Different results using different projections

1

Hi everyone,

I am solving a nonlinear elasticity for an incompressible solid where I have enforced incompressibility constraint (det(F) = 1) using a Lagrange multiplier p. My code seems to be working and my displacement values seem to be correct. However, there is a problem in the last four lines of the code:

In order to check satisfaction of incompressibility I am visualizing det(F) using paraview. If I project my obtained result in a function space W with CG(1) I can see that det(F) = 1 every where which seems correct. However, If I change CG(2), I get det(F) ~ 7 everywhere.
Is there any problem with projection? Also, when I look at solution for p, I see strong fluctuations. anyone has any idea about this issue?
Thanks

Here is my code:

from dolfin import *
import numpy as np
parameters[“form_compiler”][“quadrature_degree”] =3
parameters[“form_compiler”][“cpp_optimize”] =True
parameters[“form_compiler”][“optimize”] =True

ffc_options = {“optimize”: True,
“eliminate_zeros”: True,
“precompute_basis_const”: True,
“precompute_ip_const”: True,
“quadrature_degree”: 2}

Lx = 100
Ly = 10

mesh = RectangleMesh.create([Point(0,0),Point(Lx,Ly)],[40,4],CellType.Type.quadrilateral)

V1 = VectorElement(‘CG’, mesh.ufl_cell(), 2) # displacement finite element
V2 = FiniteElement(‘CG’, mesh.ufl_cell(), 2) # lagrange multiplier
V = FunctionSpace(mesh, MixedElement([V1, V2]))

vtest = TestFunction(V)
u_t = Function(V)
(u,p)=split(u_t)
(v_u,v_p)=split(vtest)

left = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 0.0)
bcl = DirichletBC(V.sub(0), Constant((0.,0.)), left)
bcs = [bcl]

mu = Constant(10**6)
B = Constant((0.0, -1.0))

d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr©
detF = det(F)

Sigma = mu * F - p * (inv(F).T)
Func = inner(grad(v_u),Sigma) * dx - dot(v_u,B)*dx + v_p * (detF - 1) *dx

solve(Func == 0, u_t, bcs,form_compiler_parameters=ffc_options,solver_parameters={“newton_solver”:{“relaxation_parameter”:0.5,“relative_tolerance”:1e-4 ,“maximum_iterations”:50,} })#form_compiler_parameters=ffc_options)
(u,p)=u_t.split()

file = File(“displacement.pvd”);
file << u;
file = File(“pressure.pvd”);
file << p;

#Here CG = 1, gives correct results bu CG =2 does not!!!
W=FunctionSpace(mesh,“CG”,2)# ???
projected_detF = project(det(I+grad(u)),W)
file = File(“detF.pvd”)
file << projected_detF

2

The first thing I’ll point out is that you shouldn’t use the same degree spaces for the Lagrange multiplier and displacement. If you want to use CG2 for displacement, CG1 would be a better choice for the Lagrange multiplier. (This corresponds to the classic Taylor–Hood element used for incompressible flow.) To this more deeply, look into the concept of “inf-sup stability” for discretizations of incompressible flows. (While this might not sound directly relevant to incompressible solids, note that Stokes flow is formally-identical to incompressible linear elasticity.)

However, even after correcting this issue, you’ve set the quadrature degree too low in the line

parameters["form_compiler"]["quadrature_degree"] = 3

If you increase the quadrature degree, you’ll start getting more reasonable results in projected_detF.

2 Likes
3

Thank for your response, it seems to be working now :). Do you know the reason for using CG1 for Lagrange multiplier and using CG2 for displacement? Is it because second order derivative of displacement appears in equilibrium equations and first order derivative of Lagrange multiplier?
If I want to couple equations with another physics, say electricity, where second order derivative of potential appears in equations do I have to use CG2?

Thanks for your help

4

The short, heuristic answer is that the number of pressure degrees of freedom is the number of constraints that you impose on the discrete problem for displacement, and you shouldn’t impose too many constraints on the displacement relative to the number of displacement degrees of freedom. Having an appropriate number of constraints in a discrete saddle point problem is a special case of something that numerical analysts call “inf-sup stability”, because the precise mathematical statement involves an infimum (inf) and a supremum (sup). Using CG k for displacement and CG (k-1) for pressure, with k \geq 2, is known as the “Taylor–Hood element”, and is provably inf-sup stable for the case of incompressible linear elasticity (which is formally-identical to Stokes flow).

I cover the basic theory of inf-sup stability in the notes for my graduate course on finite element analysis for coupled problems, and have a section specifically on incompressibility, although I do not delve into the proofs of inf-sup stability for specific element types (e.g., Taylor–Hood, Scott–Vogelius, etc.), as those can become quite technical.

3 Likes
5

Thank you very much. That was great help, I will look into your lecture notes.