# Specifying 'F' for Mixed finite element for Nearly Incompressible Elasticity

Hi

I'm trying to modify source code that Corrado Maurini provided a while back. In it, the strain energy (W) was provided and 'F' was calculated automatically using F = derivative(

However, I would like to specify 'F' manually as I'd like to manipulate the 2PK stress tensor. It's a two-field mixed formulation with u (displacement) and p (pressure) as unknowns.

From literature, 'F' should be:

DPi(u,p).du = integral( [J(u)*p*C(u)^-1 + 2*(dW/dC) ] : dE(u) )dV + DPi_ext

DPi(u,p).dp = integral( [J(u) - 1]*dp )dV

My attempts result in either the nonlinear variational solver producing nans or PETSC ERRORs. Below, I've included the code. Could I please get some help as to what I'm doing wrong.

Many thanks

Pietro

=======

from dolfin import *

# Create mesh and define function space

mesh = UnitSquare(20,20)

#plot(mesh)

# Mark boundary subdomains

def left(x, on_boundary):

return x[0]==0 and on_boundary

b=Constant(

# Material parameters

Y, nu = 1.0, 0.49999

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

# Defintion of function spaces

porder=2

Pu = VectorFunctionS

Pp = FunctionSpace(mesh, "CG", porder-1) # space for pressure

V = MixedFunctionSp

ndim = Pu.cell().d # number of space dimensions

# Create test and trial functions

#

vq = TestFunction(V) # Testfunction in the mixed space

dup = TrialFunction(V) # Trial function in the mixed space to define the incremental solution

du, dp = split(dup)

up = Function(V) # Function in the mixed space to store the solution

dupsol = Function(V) # Function in the mixed space to store the solution of the linearized problem

p, u = split(up) # Function in each subspace to write the functional

q, v = split(vq)

# Boundary conditions

bcL= DirichletBC(

# Kinematics

I = Identity(ndim) # Identity tensor

F = I + grad(u) # Deformation gradient

C = transpose(F)*F # Right Cauchy-Green tensor

E = 0.5*( C - I ) # Green-Lagrange tensor

# Invariants of deformation tensors

Ic = tr(C)

J = det(F)

# Stored strain energy density (compressible neo-Hookean model, mixed formulation)

W = (mu/2)*(Ic - ndim) - mu*ln(J) + p*ln(J) - 1/(2*lmbda)*p**2

# Potential energy

#Pi_bulk = W*dx - inner(b,u)*dx

# First derivatives of the bulk and surface contributions to the potential energy

#F_bulk=

#======

#Specifying F Manually

#======

dE = 0.5 * ( grad(v).T*F + F.T*grad(v) )

CC = variable(C)

S = J*p*inv(C) + 2*diff(W,CC)

F_bulk = inner( S, dE )*dx + inner(b,v)*dx + (J-1)*q*dx

#======

## Second derivatives of the bulk and surface contributions to the potential energy

dF_bulk=

# Set-up the FEniCs build-in Newton

MyNonlinearVari

solver = NonlinearVariat

solver.

# Solve the problem

solver.solve()

# Plot and hold solution

plotu=plot(u, mode = "displacement"

#######

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Jan Domurath

- Solved:
- 2013-02-15

- Last query:
- 2013-02-15

- Last reply:
- 2013-02-15

Jan Domurath (jan-domurath) said : | #1 |

Hello Pietro,

I could be wrong here but I think you are confusing the automated

differentiation algorithm with

CC = variable(C)

S = J*p*inv(C) + 2*diff(W,CC)

diff(W, CC) tries to compute the derivative of W with respect to CC which is

zero, since W doesn't depend on CC.

You should define

C = variable(C)

before you use it for the first time. This way your energy density function is a

function of the "variable" C and diff(W, C) can compute the derivative with

respect to C.

So when you define the right Cauchy-Green tensor just add

C = variable(C)

directly after it or you can define it directly as a variable with

C = variable(

and then change the expression for the stress to

S = J*p*inv(C) + 2*diff(W,C)

this way it works for me.

Also have a look at the Cahn-Hilliard demo: http://

Cheers

Jan

Pietro Maximoff (segment-x) said : | #2 |

Hi Jan,

Thanks, this was indeed the problem and it works fine now.

Many thanks,

Pietro

Pietro Maximoff (segment-x) said : | #3 |

Thanks Jan Domurath, that solved my question.