Contents     Functions              PDF Index

6.1  FEM problem formulations

• 3D elasticity
• 2D elasticity
• Acoustics
• Classical lamination theory
• Piezo-electric volumes
• Piezo-electric shells
• Geometric non-linearity
• Thermal pre-stress
• Hyperelasticity
• Gyroscopic effects
• Centrifugal follower forces
• Poroelastic materials
• Heat equation
  • Test case

This section gives a short theoretical reminder of supported FEM problems. The selection of the formulation for each element group is done through the material and element properties.

6.1.1  3D elasticity

Elements with a p_solid property entry with a non-zero integration rule are described under p_solid. They correspond exactly to the *b elements, which are now obsolete. These elements support 3D mechanics (DOFs .01 to .03 at each node) with full anisotropy, geometric non-linearity, integration rule selection, ... The elements have standard limitations. In particular they do not (yet)

• have any correction for shear locking found for high aspect ratios
• have any correction for dilatation locking found for nearly incompressible materials

With m_elastic subtypes 1 and 3, p_solid deals with 3D mechanics with strain defined by

where the engineering notation γ_yz=2є_yz, ... is used. Stress by

Note that the strain states are {є_x є_y є_z γ_yz γ_zx γ_xy} which may not be the convention of other software.

Note that NASTRAN, SAMCEF, ANSYS and MODULEF order shear stresses with σ_xy, σ_yz, σ_zx (MODULEF elements are obtained by setting p_solid integ value to zero). Abaqus uses σ_xy, σ_xz, σ_yz

In fe_stress the stress reordering can be accounted for by the definition of the proper TensorTopology matrix.

For isotropic materials

with at nominal G=E/(2(1+ν)).

For orthotropic materials, the compliance is given by

For constitutive law building, see p_solid.

6.1.2  2D elasticity

With m_elastic subtype 4, p_solid deals with 2D mechanical volumes with strain defined by (see q4p constants)

and stress by

For isotropic plane stress (p_solid form=1), one has

For isotropic plane strain (p_solid form=0), one has

6.1.3  Acoustics

With m_elastic subtype 2, p_solid deals with 2D and 3D acoustics (see flui4 constants) where 3D strain is given by

This replaces the earlier flui4 ... elements.

6.1.4  Classical lamination theory

Both isotropic and orthotropic materials are considered. In these cases, the general form of the 3D elastic material law is

Plate formulation consists in assuming one dimension, the thickness along x₃, negligible compared with the surface dimensions. Thus, vertical stress σ₃₃=0 on the bottom and upper faces, and assumed to be neglected throughout the thickness,

and for isotropic material,

By eliminating σ₃₃, the plate constitutive law is written, with engineering notations,

The reduced stiffness coefficients Q_ij (i,j = 1,2,4,5,6) are related to the 3D stiffness coefficients C_ij by

The reduced elastic law for an isotropic plate becomes,

and

Under Reissner-Mindlin's kinematic assumption the linearized strain tensor is

So, the strain vector is written,

with є^m the membrane, κ the curvature or bending, and γ the shear strains,

Note that the engineering notation with γ₁₂=u_1,2+u_2,1 is used here rather than the tensor notation with є₁₂=(u_1,2+u_2,1)/2 . Similarly κ₁₂=β_1,2+β_2,1, where a factor 1/2 would be needed for the tensor.

The plate formulation links the stress resultants, membrane forces N_αβ, bending moments M_αβ and shear forces Q_α3, to the strains, membrane є^m, bending κ and shearing γ,

The stress resultants are obtained by integrating the stresses through the thickness of the plate,

with α, β = 1, 2.

Therefore, the matrix extensional stiffness matrix [A], extension/bending coupling matrix [B], and the bending stiffness matrix [D] are calculated by integration over the thickness interval [hb ht]

An improvement of Mindlin's plate theory with tranverse shear consists in modifying the shear coefficients F_ij by

where k_ij are correction factors. Reddy's 3^rd order theory brings to k_ij=2/3. Very commonly, enriched 3^rd order theory are used, and k_ij are equal to 5/6 and give good results. For more details on the assessment of the correction factor, see [32].
For an isotropic symmetric plate (hb=−ht=h/2), the in-plane normal forces N₁₁, N₂₂ and shear force N₁₂ become

the 2 bending moments M₁₁, M₂₂ and twisting moment M₁₂

and the out-of-plane shearing forces Q₂₃ and Q₁₃,

One can notice that because the symmetry of plate, that means the reference plane is the mid-plane of the plate (x₃(0)=0) the extension/bending coupling matrix [B] is equal to zero.

Using expression (6.21) for a constant Q_ij, one sees that for a non-zero offset, one has

where is clearly appears that the constitutive matrix is a polynomial function of h, h³, x₃(0)²h and x₃(0)h. If the ply thickness is kept constant, the constitutive law is a polynomial function of 1,x₃(0),x₃(0)².

6.1.5  Piezo-electric volumes

A revised version of this information is available at http://www.sdtools.com/pdf/piezo.pdf. Missing PDF links will be found there.

The strain state associated with piezoelectric materials is described by the six classical mechanical strain components and the electrical field components. Following the IEEE standards on piezoelectricity and using matrix notations, S denotes the strain vector and E denotes the electric field vector (V/m) :

where φ is the electric potential (V).

The constitutive law associated with this strain state is given by

in which D is the electrical displacement vector (a density of charge in Cb/m²), T is the mechanical stress vector (N/m²). C^E is the matrix of elastic constants at zero electric field (E=0, short-circuited condition, see section 6.1.1 for formulas (there C^E is noted D). Note that using −E rather than E makes the constitutive law symmetric.

Alternatively, one can use the constitutive equations written in the following manner :

In which s^E is the matrix of mechanical compliances, [d] is the matrix of piezoelectric constants (m/V=Cb/N):

Matrices [e] and [d] are related through

Due to crystal symmetries, [d] may have only a few non-zero elements.

Matrix [є^S] is the matrix of dielectric constants (permittivities) under zero strain (constant volume) given by

It is more usual to find the value of є^T (Permittivity at zero stress) in the datasheet. These two values are related through the following relationship :

For this reason, the input value for the computation should be [є^T].
Also notice that usually relative permittivities are given in datasheets:

є₀ is the permittivity of vacuum (=8.854e-12 F/m)

The most widely used piezoelectric materials are PVDF and PZT. For both of these, matrix [є^T] takes the form

For PVDF, the matrix of piezoelectric constants is given by

and for PZT materials :

6.1.6  Piezo-electric shells

Shell strain is defined by the membrane, curvature and transverse shear as well as the electric field components. It is assumed that in each piezoelectric layer i=1...n, the electric field takes the form E^→= (0    0    E_zi). E_zi is assumed to be constant over the thickness h_i of the layer and is therefore given by E_zi=−Δ φ_i/h_i where Δ φ_i is the difference of potential between the electrodes at the top and bottom of the piezoelectric layer i. It is also assumed that the piezoelectric principal axes are parallel to the structural orthotropy axes.

The strain state of a piezoelectric shell takes the form

There are thus n additional degrees of freedom Δ φ_i, n being the number of piezoelectric layers in the laminate shell

The constitutive law associated to this strain state is given by :

where D_zi is the electric displacement in piezoelectric layer (assumed constant and in the z-direction), z_mi is the distance between the midplane of the shell and the midplane of piezoelectric layer i, and G_i, H_i are given by

where . denotes the direction of polarization. If the piezoelectric is used in extension mode, the polarization is in the z-direction, therefore H_i =0 and G_i ={

}_i . If the piezoelectric is used in shear mode, the polarization is in the x or y-direction, therefore G_i=0, and H_i = {0 e₁₅ }_i or H_i = {e₂₄ 0 }_i . It turns out however that the hypothesis of a uniform transverse shear strain distribution through the thickness is not satisfactory, a more elaborate shell element would be necessary. Shear actuation should therefore be used with caution.

[R_s]_i and [R]_i are rotation matrices associated to the angle θ of the piezoelectric layer.

6.1.7  Geometric non-linearity

The following gives the theory of large transformation problem implemented in OpenFEM function of_mk_pre.c Mecha3DInteg.
The principle of virtual work in non-linear total Lagrangian formulation for an hyperelastic medium is

with p the vector of initial position, x = p +u the current position, and u the displacement vector. The transformation is characterized by

where the N,j is the derivative of the shape functions with respect to Cartesian coordinates at the current integration point and q_i corresponds to field i (here translations) and element nodes. The notation is thus really valid within a single element and corresponds to the actual implementation of the element family in elem0 and of_mk. Note that in these functions, a reindexing vector is used to go from engineering ({e₁₁ e₂₂ e₃₃ 2e₂₃ 2e₃₁ 2e₁₂}) to tensor [e_ij] notations ind_ts_eg=[1 6 5;6 2 4;5 4 3];e_tensor=e_engineering(ind_ts_eg);. One can also simplify a number of computations using the fact that the contraction of a symmetric and non symmetric tensor is equal to the contraction of the symmetric tensor by the symmetric part of the non symmetric tensor.

One defines the Green-Lagrange strain tensor e=1/2(F^TF −I) and its variation

Thus the virtual work of internal loads (which corresponds to the residual in non-linear iterations) is given by

and the tangent stiffness matrix (its derivative with respect to the current position) can be written as

which using the notation u_i,j = {N_,j}^T{q_i} leads to

The term associated with stress at the current point is generally called geometric stiffness or pre-stress contribution.

In isotropic elasticity, the 2nd tensor of Piola-Kirchhoff stress is given by

the building of the constitutive law matrix D is performed in p_solid BuildConstit for isotropic, orthotropic and full anisotropic materials. of_mk_pre.c nonlin_elas then implements element level computations. For hyperelastic materials ∂² W/∂ e² is not constant and is computed at each integration point as implemented in hyper.c.

For a geometric non-linear static computation, a Newton solver will thus iterate with

where external forces f are assumed to be non following.

For an example see staticNewton.

6.1.8  Thermal pre-stress

The following gives the theory of the thermoelastic problem implemented in OpenFEM function of_mk_pre.c nonlin_elas.
In presence of a temperature difference, the thermal strain is given by [e_T] = [α] (T−T₀), where in general the thermal expansion matrix α is proportional to identity (isotropic expansion). The stress is found by computing the contribution of the mechanical deformation

This expression of the stress is then used in the equilibrium (6.44), the tangent matrix computation(6.48), or the Newton iteration (6.51). Note that the fixed contribution ∫_Ω0 (−C:e_T) : δ e can be considered as an internal load of thermal origin.

The modes of the heated structure can be computed with the tangent matrix.

An example of static thermal computation is given in ofdemos ThermalCube.

6.1.9  Hyperelasticity

The following gives the theory of the thermoelastic problem implemented in OpenFEM function hyper.c (called by of_mk.c MatrixIntegration).
For hyperelastic media S=∂ W/∂ e with W the hyperelastic energy. hyper.c currently supports Mooney-Rivlin materials for which the energy takes one of following forms

where (J₁,J₂,J₃) are the so-called reduced invariants of the Cauchy-Green tensor

linked to the classical invariants (I₁,I₂,I₃) by

where one recalls that

Note : this definition of energy based on reduced invariants is used to have the hydrostatic pressure given directly by p=−K(J₃−1) (K “bulk modulus”), and the third term of W is a penalty on incompressibility.

Hence, computing the corresponding tangent stiffness and residual operators will require the derivatives of the above invariants with respect to e (or C). In an orthonormal basis the first-order derivatives are given by:

where (C_ij⁻¹) denotes the coefficients of the inverse matrix of (C_ij). For second-order derivatives we have:

where the є_ijk coefficients are defined by

Note: when the strain components are seen as a column vector (“engineering strains”) in the form (e₁₁,e₂₂,e₃₃,2e₂₃,2e₃₁,2e₁₂)′, the last two terms of (6.59) thus correspond to the following 2 matrices

We finally use chain-rule differentiation to compute

Note that a factor 2 arise each time we differentiate the invariants with respect to e instead of C.

The specification of a material is given by specification of the derivatives of the energy with respect to invariants. The laws are implemented in the hyper.c EnPassiv function.

6.1.10  Gyroscopic effects

Written by Arnaud Sternchuss ECP/MSSMat.

In the fixed reference frame which is Galilean, the Eulerian speed of the particle in x whose initial position is p is

and its acceleration is

Ω is the rotation vector of the structure with

in a (x,y,z) orthonormal frame. The skew-symmetric matrix [Ω] is defined such that

The speed can be rewritten

and the acceleration becomes

In this expression appear

• the acceleration in the rotating frame ∂²u/∂ t²,
• the centrifugal acceleration a_g=[Ω]²(p+u),
• the Coriolis acceleration a_c=∂[Ω]/∂ t(p+u)+2[Ω]∂u/∂ t.

S₀^e is an element of the mesh of the initial configuration S₀ whose density is ρ₀. [N] is the matrix of shape functions on these elements, one defines the following elementary matrices

6.1.11  Centrifugal follower forces

This is the embryo of the theory for the future implementation of centrifugal follower forces.

where δ v_R designates the radial component (in deformed configuration) of δv. One assumes that the rotation axis is along e_z. Noting n_R = 1/R {x₁  x₂   0}^T, one then has

Thus the non-linear stiffness term is given by

One has dR=n_R· dx(= dx_R) and dδ v_R = dn_R·δ v, with

Thus, finally

Which gives

with α=1,2.

6.1.12  Poroelastic materials

The poroelastic formulation comes from [33], recalled and detailed in [34].

Domain and variables description:

Weak formulation, for harmonic time dependence at pulsation ω:

Matrix formulation, for harmonic time dependence at pulsation ω:

where the frequency-dependent matrices correspond to:

N.B. if the material of the solid phase is homogeneous, the frequency-dependent parameters can be eventually factorized from the matrices:

where the matrices marked with bars are frequency independent:

Material parameters:

Constants:

Poroelastic specific (frequency dependent) variables:

To add here:

• coupling conditions with poroelastic medium, elastic medium, acoustic medium
• dissipated power in medium

6.1.13  Heat equation

This section is based on an OpenFEM contribution by Bourquin Frédéric and Nassiopoulos Alexandre from Laboratoire Central des Ponts et Chaussées.

The variational form of the Heat equation is given by

with

• ρ the density, c the specific heat capacity.
• K the conductivity tensor of the material. The tensor K is symmetric, positive definite, and is often taken as diagonal. If conduction is isotropic, one can write K=k(x)Id where k(x) is called the (scalar) conductivity of the material.
• Acceptable loads and boundary conditions are
  • Internal heat source f
    
  • Prescribed temperature (Dirichlet condition, also called boundary condition of first kind)
    

    θ=θext    on    ∂Ω     (6.78)

    modeled using a DofSet case entry.
  • Prescribed heat flux g (Neumann condition, also called boundary condition of second kind)
    

    (Kgrad  θ)·

    → n

    =g    on    ∂Ω     (6.79)

    leading to a load applied on the surface modeled using a FVol case entry.
  • Exchange and heat flux (Fourier-Robin condition, also called boundary condition of third kind)
    

    (Kgrad  θ)·

    → n

    +α(θ−θext)=g    on    ∂Ω     (6.80)

    leading to a stiffness term (modeled using a group of surface elements with stiffness proportional to α) and a load on the associated surface proportional to g+αθ_ext (modeled using FVol case entries).

Test case

One considers a solid square prism of dimensions L_x,L_y, L_z in the three directions (Ox), (Oy) and (Oz) respectively. The solid is made of homogeneous isotropic material, and its conductivity tensor thus reduces to a constant k.

The faces, Γ_i (i=1..6, ∪_i=1⁶ Γ_i = ∂ Ω), are subject to the following boundary conditions and loads

• f=40 is a constant uniform internal heat source
• Γ₁  (x=0) : exchange & heat flux (Fourier-Robin) given by α=1,g₁=α θ_ext + α f L_x²/2k=25
  
• Γ₂  (x=L_x) : prescribed temperature : θ(L_x,y,z)=θ_ext=20
• Γ₃  (y=0), Γ₄  (y=L_y), Γ₅  (z=0), Γ₆  (z=L_z): exchange & heat flux g+α θ_ext =αθ_ext +α f/2 k (L_x²−x²)+g₁=25−x²/20
  

The problem can be solved by the method of separation of variables. It admits the solution

The resolution for this example can be found in demo/heat_equation.