Finite strain Elasticity bricks
This brick implements some classical hyperelastic constitutive law for large deformation elasticity.
Some recalls on finite strain elasticity
Let be the reference configuration and the deformed configuration of an elastic media. Then for we will denote by the deformation. the vector field 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, is neglected).
One has
Both tensors and 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:
are the invariants of orders and :
The derivatives of the invariants with respect to the tensor in the direction are:
We will write
Let us also recall that
The second derivatives of the invariants are fourth order tensors defined by
The notation denotes the Frobenius product . 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 where is the Cauchy stress tensor in the deformed configuration . An hyper-elastic constitutive law is given by
where 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 can be writen
because in particular
and when A is symmetric which is the case for .
Another way is to consider the static equilibrium which can be written as follows in the reference configuration:
Integrating by parts, one obtains:
Some classical constitutive laws
Three parameters Mooney-Rivlin law
Compressible material.
where , and are given coefficients and
and then
Incompressible material.
The incompressibility constraint is handled with a Lagrange multiplier (the pressure)
constraint:
Ciarlet-Geymonat law
with the Lame coefficients and (see [ciarlet1988]).
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 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 ().
The Mooney-Rivlin law accepts two optional flags, the first one determines if the material will be compressible () and the second one determines if the material is neo Hookean (). 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 (, ). 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 decribe 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 weak form language 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_sigma(Grad_u, [lambda; mu])
Plane_Strain_Saint_Venant_Kirchhoff_sigma(Grad_u, [lambda; mu])
Generalized_Blatz_Ko_sigma(Grad_u, [a;b;c;d;n])
Plane_Strain_Generalized_Blatz_Ko_sigma(Grad_u, [a;b;c;d;n])
Ciarlet_Geymonat_sigma(Grad_u, [lambda;mu;a])
Plane_Strain_Ciarlet_Geymonat_sigma(Grad_u, [lambda;mu;a])
Incompressible_Mooney_Rivlin_sigma(Grad_u, [c1;c2])
Plane_Strain_Incompressible_Mooney_Rivlin_sigma(Grad_u, [c1;c2])
Compressible_Mooney_Rivlin_sigma(Grad_u, [c1;c2;d1])
Plane_Strain_Compressible_Mooney_Rivlin_sigma(Grad_u, [c1;c2;d1])
Incompressible_Neo_Hookean_sigma(Grad_u, [c1])
Plane_Strain_Incompressible_Neo_Hookean_sigma(Grad_u, [c1])
Compressible_Neo_Hookean_sigma(Grad_u, [c1;d1])
Plane_Strain_Compressible_Neo_Hookean_sigma(Grad_u, [c1;d1])
Compressible_Neo_Hookean_Bonet_sigma(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Bonet_sigma(Grad_u, [lambda;mu])
Compressible_Neo_Hookean_Ciarlet_sigma(Grad_u, [lambda;mu])
Plane_Strain_Compressible_Neo_Hookean_Ciarlet_sigma(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. Two reasons are first that the second order derivative of the potential (necessary to obtain the tangent system) can be very complicated and non-optimized and 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_sigma(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.