# Finite strain Elasticity bricks¶

This brick implements some classical hyperelastic constitutive law for large deformation elasticity.

## Some recalls on finite strain elasticity¶

Let \(\Omega\) be the reference configuration and \(\Omega_t\) the deformed configuration of an elastic media. Then for \(X \in \Omega\) we will denote by \(\Phi(x) = u(X) + X\) the deformation. the vector field \(u\) is the displacement with respect to the initial position.

The Cauchy-Green tensor is defined by

The deformation tensor (Green-Lagrange)

(In the case of linear elasticity, \({\nabla u^T}{\nabla u}\) is neglected).

One has

Both tensors \(E\) and \(C\) are used to describe finite strain elasticity constitutive laws.

### Main invariants and derivatives¶

The description of finite strain elasticity constitutive laws often requires the principal invariants of the deformation tensors:

\(i_1,i_2,i_3\) are the invariants of orders \(1,2\) and \(3\):

The derivatives of the invariants with respect to the tensor \(E\) in the direction \(H\) are:

We will write

Let us also recall that

The second derivatives of the invariants are fourth order tensors defined by

The notation \(A:B\) denotes the Frobenius product \(A:B = \displaystyle\sum_{ij}A_{ij}B_{ij}\). This product has the following properties:

Note also that

This property enables us to write the constitutive laws as a function of the Cauchy-Green tensor invariants, especially for the case of the generalized Blatz-Ko strain energy.

### Potential elastic energy and its derivative¶

The stress in the reference configuration can be describe by the second Piola-Kirchhoff stress tensor \({\hat{\hat{\sigma}}} = \nabla\Phi^{-1}\sigma\nabla\Phi^{-t}~\det \nabla\Phi\) where \(\sigma\) is the Cauchy stress tensor in the deformed configuration \(\Omega_t\). An hyper-elastic constitutive law is given by

where \({W}\) is the density of strain energy of the material. The total strain energy is given by

and the derivative of the energy in a direction \(v\) can be writen

because in particular

and \(A:B = A:(B+B^T)/2\) when A is symmetric which is the case for \({\hat{\hat{\sigma}}}\).

Another way is to consider the static equilibrium which can be written as follows in the reference configuration:

Integrating by parts, one obtains:

### Tangent matrix¶

The displacement \(u\) is fixed. In order to obtain the tangent matrix, one subsitutes \(u\) with \(u+h\)

and considers the linear part w.r.t. \(h\), which is

which is symmetric w.r.t. \(v\) and \(h\). It can be rewritten as

where \(\mathcal{A}\) is the symmetric \(3\times3\times3\times3\) tensor given by \(\mathcal{A}_{ijkl} = ((\frac{\partial^2 W}{\partial E^2})_{ijkl} + (\frac{\partial^2 W}{\partial E^2})_{ijlk})/2\).

### Some classical constitutive laws¶

`Linearized: Saint-Venant Kirchhoff law (small deformations)`

¶

`Three parameters Mooney-Rivlin law`

¶

Compressible material.

where \(c_1\), \(c_2\) and \(d_1\) are given coefficients and

and then

Incompressible material.

The incompressibility constraint \(i_3( C) = 1\) is handled with a Lagrange multiplier \(p\) (the pressure)

constraint: \(\sigma = -pI \Rightarrow {\hat{\hat{\sigma}}} = -p\nabla\Phi\nabla\Phi^{-T}\det\nabla\Phi\)

`Ciarlet-Geymonat law`

¶

with \(\lambda, \mu\) the Lame coefficients and \(\max(0,\frac{\mu}{2}-\frac{\lambda}{4})<a<\frac{\mu}{2}\) (see [ciarlet1988]).

`Generalized Blatz-Ko law`

¶

Since \(\frac{\partial}{\partial C} {W}(C) = \displaystyle\sum_{j}\frac{\partial W}{\partial i_j(C)} \frac{\partial i_j(C)}{\partial C}\), and \(\frac{\partial^2}{\partial C^2} {W}(C) = \displaystyle\sum_{j} \displaystyle\sum_{k} \frac{\partial^2 W}{\partial i_j(C) \partial i_k(C)} \frac{\partial i_k(C)}{\partial C} \otimes \frac{\partial i_j(C)}{\partial C} + \displaystyle\sum_{j} \frac{\partial W}{\partial i_j(C)} \frac{\partial^2 i_j(C)}{\partial C^2}\) we must compute the derivatives of the strain energy function with respect to the Cauchy-Green tensor invariants (we don’t need to compute the invariants derivatives with respect to \(E\) since \(\frac{\partial i_j}{\partial E}(C;H) = 2 \frac{\partial i_j}{\partial C}(C;H)\)) :

`Plane strain hyper-elasticity`

¶

All previous models are valid in volumic domains. Corresponding plane strain 2D models can be obtained by restricting the stress tensor and the fourth order tensor \(\mathcal{A}\) to their plane components.

## Add an nonlinear elasticity brick to a model¶

This brick represents a large strain elasticity problem. It is defined in the files `getfem/getfem_nonlinear_elasticity.h`

and `getfem/getfem_nonlinear_elasticity.cc`

. The function adding this brick to a model is

```
ind = getfem::add_nonlinear_elasticity_brick
(md, mim, varname, AHL, dataname, region = -1);
```

where `AHL`

is an object of type `getfem::abstract_hyperelastic_law`

which represents the considered hyperelastic law. It has to be chosen between:

```
getfem::SaintVenant_Kirchhoff_hyperelastic_law AHL;
getfem::Ciarlet_Geymonat_hyperelastic_law AHL;
getfem::Mooney_Rivlin_hyperelastic_law AHL(compressible, neohookean);
getfem::plane_strain_hyperelastic_law AHL(pAHL);
getfem::generalized_Blatz_Ko_hyperelastic_law AHL;
```

The Saint-Venant Kirchhoff law is a linearized law defined with the two Lame coefficients, Ciarlet Geymonat law is defined with the two Lame coefficients and an additional coefficient (\(\lambda, \mu, a\)).

The Mooney-Rivlin law accepts two optional flags, the first one determines if the material will be compressible (\(d_1 \neq 0\)) and the second one determines if the material is neo Hookean (\(c_2 = 0\)). Depending on these flags one to three coefficients may be necessary. By default it is defined as incompressible and non neo Hookean, thus it needs two material coefficients (\(c_1\), \(c_2\)). In this case, it is to be used with the large strain incompressibility condition.

The plane strain hyperelastic law takes a pointer on a hyperelastic law as a parameter and performs a 2D plane strain approximation.

`md`

is the model variable, `mim`

the integration method, `varname`

the string being the name of the variable on which the term is added, `dataname`

the string being the name of the data in the model representing the coefficients of the law (can be constant or describe on a finite element method) and `region`

is the region on which the term is considered (by default, all the mesh).

The program `nonlinear_elastostatic.cc`

in `tests`

directory and `demo_nonlinear_elasticity.m`

in `interface/tests/matlab`

directory are some examples of use of this brick with or without an incompressibility condition.

Note that the addition of a new hyperelastic constitutive law consists in furnishing the expression of the strain energy, the stress tensor and the derivative of the stress tensor. See the file `getfem/getfem_nonlinear_elasticity.cc`

for more details. In particular, expression of the invariants and their derivatives are available.

A function which computes the Von Mises or Tresca stresses is also available:

```
VM = compute_Von_Mises_or_Tresca
(md, varname, AHL, dataname, mf_vm, VM, tresca)
```

It returns a vector of the degrees of freedom of the Von Mises or Tresca stress on the finite element method mf_vm. `tresca`

is a boolean whose value should be `true`

for Tresca stress and `false`

for Von Mises stress.

## Add a large strain incompressibility brick to a model¶

This brick adds an incompressibility condition in a large strain problem of type

A Lagrange multiplier representing the pressure is introduced in a mixed formulation. The function adding this brick to a model is

```
ind = add_nonlinear_incompressibility_brick
(md, mim, varname, multname, region = -1)
```

where `md`

is the model, `mim`

the integration method, `varname`

the variable of the model on which the incompressibility condition is added, `multanme`

the multiplier variable corresponding to the pressure (be aware that at least a linear Ladyzhenskaja-Babuska-Brezzi inf-sup condition is satisfied between the f.e.m. of the variable and the one of the multiplier). `region`

is an optional parameter correponding to the mesh region on which the term is considered (by default, all the mesh).

## High-level generic assembly versions¶

The generic weak form language (GWFL) gives access to the hyperelastic potential and constitutive laws implemented in *GetFEM*. This allows to directly use them in the language, for instance using a generic assembly brick in a model or for interpolation of certain quantities (the stress for instance).

Here is the list of nonlinear operators in the language which can be useful for nonlinear elasticity:

```
Det(M) % determinant of the matrix M
Trace(M) % trace of the matrix M
Matrix_i2(M) % second invariant of M (in 3D): (sqr(Trace(m)) - Trace(m*m))/2
Matrix_j1(M) % modified first invariant of M: Trace(m)pow(Det(m),-1/3).
Matrix_j2(M) % modified second invariant of M: Matrix_I2(m)*pow(Det(m),-2/3).
Right_Cauchy_Green(F) % F' * F
Left_Cauchy_Green(F) % F * F'
Green_Lagrangian(F) % (F'F - Id(meshdim))/2
Cauchy_stress_from_PK2(sigma, Grad_u) % (Id+Grad_u)*sigma*(I+Grad_u')/det(I+Grad_u)
```

The potentials:

```
Saint_Venant_Kirchhoff_potential(Grad_u, [lambda; mu])
Plane_Strain_Saint_Venant_Kirchhoff_potential(Grad_u, [lambda; mu])
Generalized_Blatz_Ko_potential(Grad_u, [a;b;c;d;n])
Plane_Strain_Generalized_Blatz_Ko_potential(Grad_u, [a;b;c;d;n])
Ciarlet_Geymonat_potential(Grad_u, [lambda;mu;a])
Plane_Strain_Ciarlet_Geymonat_potential(Grad_u, [lambda;mu;a])
Incompressible_Mooney_Rivlin_potential(Grad_u, [c1;c2])
Plane_Strain_Incompressible_Mooney_Rivlin_potential(Grad_u, [c1;c2])
Compressible_Mooney_Rivlin_potential(Grad_u, [c1;c2;d1])
Plane_Strain_Compressible_Mooney_Rivlin_potential(Grad_u, [c1;c2;d1])
Incompressible_Neo_Hookean_potential(Grad_u, [c1])
Plane_Strain_Incompressible_Neo_Hookean_potential(Grad_u, [c1])
Compressible_Neo_Hookean_potential(Grad_u, [c1;d1])
Plane_Strain_Compressible_Neo_Hookean_potential(Grad_u, [c1;d1])
Compressible_Neo_Hookean_Bonet_potential(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Bonet_potential(Grad_u, [lambda;mu])
Compressible_Neo_Hookean_Ciarlet_potential(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Ciarlet_potential(Grad_u, [lambda;mu])
```

The second Piola-Kirchhoff stress tensors:

```
Saint_Venant_Kirchhoff_PK2(Grad_u, [lambda; mu])
Plane_Strain_Saint_Venant_Kirchhoff_PK2(Grad_u, [lambda; mu])
Generalized_Blatz_Ko_PK2(Grad_u, [a;b;c;d;n])
Plane_Strain_Generalized_Blatz_Ko_PK2(Grad_u, [a;b;c;d;n])
Ciarlet_Geymonat_PK2(Grad_u, [lambda;mu;a])
Plane_Strain_Ciarlet_Geymonat_PK2(Grad_u, [lambda;mu;a])
Incompressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2])
Plane_Strain_Incompressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2])
Compressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2;d1])
Plane_Strain_Compressible_Mooney_Rivlin_PK2(Grad_u, [c1;c2;d1])
Incompressible_Neo_Hookean_PK2(Grad_u, [c1])
Plane_Strain_Incompressible_Neo_Hookean_PK2(Grad_u, [c1])
Compressible_Neo_Hookean_PK2(Grad_u, [c1;d1])
Plane_Strain_Compressible_Neo_Hookean_PK2(Grad_u, [c1;d1])
Compressible_Neo_Hookean_Bonet_PK2(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Bonet_PK2(Grad_u, [lambda;mu])
Compressible_Neo_Hookean_Ciarlet_PK2(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Ciarlet_PK2(Grad_u, [lambda;mu])
```

Note that the derivatives with respect to the material parameters have not been implemented apart for the Saint Venant Kirchhoff hyperelastic law. Therefore, it is not possible to make the parameter depend on other variables of a model (derivatives are not necessary complicated to implement but for the moment, only a wrapper with old implementations has been written).

Note that the coupling of models is to be done at the weak formulation level. In a general way, it is recommended not to use the potential to define a problem. Main couplings cannot be obtained at the potential level. Thus the use of potential should be restricted to the actual computation of the potential.

An example of use to add a Saint Venant-Kirchhoff hyperelastic term to a variable `u`

in a model or a ga_workspace is given by the addition of the following assembly string:

```
"((Id(meshdim)+Grad_u)*(Saint_Venant_Kirchhoff_PK2(Grad_u,[lambda;mu]))):Grad_Test_u"
```

Note that in that case, `lambda`

and `mu`

have to be declared data of the model/ga_workspace. It is of course possible to replace them by explicit constants or expressions depending on several data.

Concerning the incompressible Mooney-Rivlin law, it has to be completed by an incompressibility term. For instance by adding the following incompressibility brick:

```
ind = add_finite_strain_incompressibility_brick(md, mim, varname, multname, region = -1);
```

This brick just adds the term `p*(1-Det(Id(meshdim)+Grad_u))`

if `p`

is the multiplier and `u`

the variable which represents the displacement.

The addition of an hyperelastic term to a model can also be done thanks to the following function:

```
ind = add_finite_strain_elasticity_brick(md, mim, lawname, varname, params,
region = size_type(-1));
```

where `md`

is the model, `mim`

the integration method, `varname`

the variable of the model representing the large strain displacement, `lawname`

is the constitutive law name which could be `Saint_Venant_Kirchhoff`

, `Generalized_Blatz_Ko`

, `Ciarlet_Geymonat`

, `Incompressible_Mooney_Rivlin`

, `Compressible_Mooney_Rivlin`

, `Incompressible_Neo_Hookean`

, `Compressible_Neo_Hookean`

, `Compressible_Neo_Hookean_Bonet`

or `Compressible_Neo_Hookean_Ciarlet`

. `params`

is a string representing the parameters of the law defined as a small vector or a vector field.

The Von Mises stress can be interpolated with the following function:

```
void compute_finite_strain_elasticity_Von_Mises(md, varname, lawname, params, mf_vm, VM,
rg=mesh_region::all_convexes());
```

where `md`

is the model, `varname`

the variable of the model representing the large strain displacement, `lawname`

is the constitutive law name (see previou brick), `params`

is a string representing the parameters of the law, `mf_vm`

a (preferably discontinuous) Lagrange finite element method on which the interpolation will be done and `VM`

a vector of type `model_real_plain_vector`

in which the interpolation will be stored.