6.1 FEM problem formulations#
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 nonzero 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 nonlinearity, 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
(6.1) 
where the engineering notation γ_{yz}=2є_{yz}, ... is used. Stress by
(6.2) 
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
(6.3) 
with at nominal G=E/(2(1+ν)). For isotropic materials, interpolation of ρ,η,E,ν,G,α with temperature is supported.
For orthotropic materials, the compliance is given by
(6.4) 
For constitutive law building, see p_solid. Material orientation can be interpolated by defining v1 and v2 fields in the InfoAtNode. Interpolation of non isotropic material properties was only implemented for of_mk >= 1.236.
6.1.2 2D elasticity#
With m_elastic subtype 4, p_solid deals with 2D mechanical volumes with strain defined by (see q4p constants)
(6.5) 
and stress by
(6.6) 
For isotropic plane stress (p_solid form=1), one has
(6.7) 
For isotropic plane strain (p_solid form=0), one has
(6.8) 
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
(6.9) 
This replaces the earlier flui4 ... elements.
The mass and stiffness matrices are given by
(6.10) 
(6.11) 
The source associated with a enforced velocity on a surface
(6.12) 
When an impedance Z=ρ C R(1+iη) is considered on a surface, the associated viscous damping matrix is given by
(6.13) 
6.1.4 Piezoelectric 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) :
(6.14) 
where φ is the electric potential (V).
The constitutive law associated with this strain state is given by
(6.15) 
in which D is the electrical displacement vector (a density of charge in Cb/m^{2}), T is the mechanical stress vector (N/m^{2}). C^{E} is the matrix of elastic constants at zero electric field (E=0, shortcircuited 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 :
(6.16) 
In which s^{E} is the matrix of mechanical compliances, [d] is the matrix of piezoelectric constants (m/V=Cb/N):
(6.17) 
Matrices [e] and [d] are related through
(6.18) 
Due to crystal symmetries, [d] may have only a few nonzero elements.
Matrix [є^{S}] is the matrix of dielectric constants (permittivities) under zero strain (constant volume) given by
(6.19) 
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 :
(6.20) 
For this reason, the input value for the computation should be
[є^{T}].
Also notice that usually relative permittivities are given in datasheets:
(6.21) 
є_{0} is the permittivity of vacuum (=8.854e12 F/m)
The most widely used piezoelectric materials are PVDF and PZT. For both of these, matrix [є^{T}] takes the form
(6.22) 
For PVDF, the matrix of piezoelectric constants is given by
(6.23) 
and for PZT materials :
(6.24) 
6.1.5 Piezoelectric shells#
A revised version of this information is available at http://www.sdtools.com/pdf/piezo.pdf.
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
(6.25) 
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 :
(6.26) 
where D_{zi} is the electric displacement in piezoelectric layer (assumed constant and in the zdirection), 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
(6.27) 
(6.28) 
where . denotes the direction of polarization. If the piezoelectric is used in extension mode, the polarization is in the zdirection, therefore H_{i} =0 and G_{i} ={
e_{31}  e_{32}  0 
}_{i} . If the piezoelectric is used in shear mode, the polarization is in the x or ydirection, therefore G_{i}=0, and H_{i} = {0 e_{15} }_{i} or H_{i} = {e_{24} 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.29) 
(6.30) 
6.1.6 Geometric nonlinearity#
The following gives the theory of large transformation problem implemented in OpenFEM function of_mk_pre.c Mecha3DInteg.
The principle of virtual work in nonlinear total Lagrangian formulation for an hyperelastic medium is
(6.31) 
with p the vector of initial position, x = p +u the current position, and u the displacement vector. The transformation is characterized by
(6.32) 
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_{11} e_{22} e_{33} 2e_{23} 2e_{31} 2e_{12}}) 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 GreenLagrange strain tensor e=1/2(F^{T}F −I) and its variation
(6.33) 
Thus the virtual work of internal loads (which corresponds to the residual in nonlinear iterations) is given by
(6.34) 
and the tangent stiffness matrix (its derivative with respect to the current position) can be written as
(6.35) 
which using the notation u_{i,j} = {N_{,j}}^{T}{q_{i}} leads to
(6.36) 
The term associated with stress at the current point is generally called geometric stiffness or prestress contribution. For implementation, the variable names are d2wde2, Sigma and the large displacement computation (F_{mk} ∂^{2} W/∂ e^{2}_{ijkl} F_{ni} + S_{lj}) has a reference implementation in elem0('LdDD'). The result is called dd in the code.
In isotropic elasticity, the 2nd tensor of PiolaKirchhoff stress is given by
(6.37) 
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 ∂^{2} W/∂ e^{2} is not constant and is computed at each integration point as implemented in hyper.c.
For a geometric nonlinear static computation, a Newton solver will thus iterate with
(6.38) 
where external forces f are assumed to be non following.
For an example see staticNewton.
6.1.7 Thermal prestress#
Note that more recent developments are found in SDTnlsim, see sdtweb('hyper3D'). 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_{0}), 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
(6.39) 
This expression of the stress is then used in the equilibrium (6.31), the tangent matrix computation(6.35), or the Newton iteration (6.38). 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.8 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 MooneyRivlin materials for which the energy takes one of following forms
(6.40) 
(6.41) 
where (J_{1},J_{2},J_{3}) are the socalled reduced invariants of the CauchyGreen tensor
(6.42) 
linked to the classical invariants (I_{1},I_{2},I_{3}) by
(6.43) 
where one recalls that
(6.44) 
Note : this definition of energy based on reduced invariants is used to have the hydrostatic pressure given directly by p=−K(J_{3}−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 firstorder derivatives are given by:
(6.45) 
where (C_{ij}^{−1}) denotes the coefficients of the inverse matrix of (C_{ij}). For secondorder derivatives we have:
(6.46) 
where the є_{ijk} coefficients are defined by
(6.47) 
Note: when the strain components are seen as a column vector ("engineering strains") in the form (e_{11},e_{22},e_{33},2e_{23},2e_{31},2e_{12})′, the last two terms of (6.46) thus correspond to the following 2 matrices
(6.48) 
(6.49) 
We finally use chainrule differentiation to compute
(6.50) 
(6.51) 
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.9 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
(6.52) 
and its acceleration is
(6.53) 
Ω is the rotation vector of the structure with
(6.54) 
in a (x,y,z) orthonormal frame. The skewsymmetric matrix [Ω] is defined such that
(6.55) 
The speed can be rewritten
(6.56) 
and the acceleration becomes
(6.57) 
In this expression appear
 the acceleration in the rotating frame ,
 the centrifugal acceleration ,
 the Coriolis acceleration .
S_{0}^{e} is an element of the mesh of the initial configuration S_{0} whose density is ρ_{0}. [N] is the matrix of shape functions on these elements, one defines the following elementary matrices
(6.58) 
The traditional fe_mknl MatType in SDT are 7 for gyroscopic coupling and 8 for centrifugal softening.
6.1.10 Centrifugal follower forces#
This is the embryo of the theory for the future implementation of centrifugal follower forces.
(6.59) 
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_{1} x_{2} 0}^{T}, one then has
(6.60) 
Thus the nonlinear stiffness term is given by
(6.61) 
One has dR=n_{R}· dx(= dx_{R}) and dδ v_{R} = dn_{R}·δ v, with
dn_{R}=− 
 n_{R} + 
 {dx_{1} dx_{2} 0}^{T}. 
Thus, finally
(6.62) 
Which gives
(6.63) 
with α=1,2.
6.1.11 Poroelastic materials#
The poroelastic formulation comes from [40], recalled and detailed in [41].
Domain and variables description:
Poroelastic domain  
Bounding surface of poroelastic domain  
Unit external normal of  
Solid phase displacement vector  
Fluid phase displacement vector  
Fluid phase pressure  
Stress tensor of solid phase  
Total stress tensor of porous material 
Weak formulation, for harmonic time dependence at pulsation ω:
(6.64) 
(6.65) 
Matrix formulation, for harmonic time dependence at pulsation ω:
(6.66) 
where the frequencydependent matrices correspond to:
(6.67) 
N.B. if the material of the solid phase is homogeneous, the frequencydependent parameters can be eventually factorized from the matrices:
(6.68) 
where the matrices marked with bars are frequency independent:
(6.69) 
Material parameters:
φ  Porosity of the porous material 
σ  Resistivity of the porous material 
α_{∞}  Tortuosity of the porous material 
Λ  Viscous characteristic length of the porous material 
Λ′  Thermal characteristic length of the skeleton 
ρ  Density of the skeleton 
G  Shear modulus of the skeleton 
ν  Poisson coefficient of the skeleton 
η_{s}  Structural loss factor of the skeleton 
ρ_{o}  Fluid density 
γ  Heat capacity ratio of fluid (=1.4 for air) 
η  Shear viscosity of fluid (=1.84×10^{−5} kg m^{−1} s^{−1} for air) 
Constants:
P_{o}=1,01× 10^{5} Pa  Ambient pressure 
Pr=0.71  Prandtl number 
Poroelastic specific (frequency dependent) variables:
Apparent density of solid phase  
Apparent density of fluid phase  
Interaction apparent density  
Effective density of solid phase  
Effective density of solid phase  
Effective density of fluid phase  
Interaction effective density  
Viscous damping coefficient  
Coupling coefficient  
Elastic coupling coefficient  
Biot formulation  
Approximation from  
Bulk modulus of air in fraction volume  
Biot formulation  
Approximation from  
Bulk modulus of porous material in vacuo  
Bulk modulus of elastic solid  
est. from HashinShtrikman's upper bound  
Effective bulk modulus of air in pores  
Function in (ChampouxAllard model)  
Thermal characteristic frequency 
To add here:
 coupling conditions with poroelastic medium, elastic medium, acoustic medium
 dissipated power in medium
6.1.12 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
(6.70) 
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)
modeled using a DofSet case entry.(6.71)  Prescribed heat flux g (Neumann condition, also called boundary condition of second kind)
leading to a load applied on the surface modeled using a FVol case entry.(6.72)  Exchange and heat flux (FourierRobin condition, also called boundary condition of third kind)
(6.73) 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).

Internal heat source f
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, , are subject to the following boundary conditions and loads
 f=40 is a constant uniform internal heat source
 Γ_{1} (x=0) : exchange & heat flux (FourierRobin) given by α=1,g_{1}=α θ_{ext} + α f L_{x}^{2}/2k=25
 Γ_{2} (x=L_{x}) : prescribed temperature : θ(L_{x},y,z)=θ_{ext}=20
 Γ_{3} (y=0), Γ_{4} (y=L_{y}), Γ_{5} (z=0), Γ_{6} (z=L_{z}): exchange & heat flux g+α θ_{ext} =αθ_{ext} +α f/2 k (L_{x}^{2}−x^{2})+g_{1}=25−x^{2}/20
The problem can be solved by the method of separation of variables. It admits the solution
θ(x,y,z) =− 
 x^{2} + θ_{ext} + 
 = 
 = 25 − 

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