Contents     Functions         Previous Next     PDF Index

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 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)

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

{
єx 
єy 
єz 
γyz 
γzx 
γxy 
} [
 N,x0
 0N,y
 00N,z 
 0N,zN,y 
 N,z0N,x 
 N,yN,x
] {
u 
v 
w 
}     (6.1)

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

{
σx 
σy 
σz 
σyz 
σzx 
σxy 
} =[
 d1,1 N,x+d1,5 N,z+d1,6 N,yd1,2 N,y+d1,4 N,z+d1,6 N,xd1,3 N,z+d1,4 N,y+d1,5 N,x 
 d2,1 N,x+d2,5 N,z+d2,6 N,yd2,2 N,y+d2,4 N,z+d2,6 N,xd2,3 N,z+d2,4 N,y+d2,5 N,x 
 d3,1 N,x+d3,5 N,z+d3,6 N,yd3,2 N,y+d3,4 N,z+d3,6 N,xd3,3 N,z+d3,4 N,y+d3,5 N,x 
 d4,1 N,x+d4,5 N,z+d4,6 N,yd4,2 N,y+d4,4 N,z+d4,6 N,xd4,3 N,z+d4,4 N,y+d4,5 N,x 
 d5,1 N,x+d5,5 N,z+d5,6 N,yd5,2 N,y+d5,4 N,z+d5,6 N,xd5,3 N,z+d5,4 N,y+d5,5 N,x 
 d6,1 N,x+d6,5 N,z+d6,6 N,yd6,2 N,y+d6,4 N,z+d6,6 N,xd6,3 N,z+d6,4 N,y+d6,5 N,x 
] {
u 
v 
w 
}

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

D=[
 
E(1−ν)
(1+ν)(1−2ν)
[
1
ν
1−ν
ν
1−ν
 
 
ν
1−ν
1
ν
1−ν
 
 
ν
1−ν
ν
1−ν
]
 0
[
G00
0G0
00G
]
]     (6.2)

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

For orthotropic materials, the compliance is given by

{є} = [D]−1{σ}=[
1/E1
ν21
E2
ν31
E3
00
ν12
E1
1/E2
ν32
E3
00
ν13
E1
ν23
E2
1/E300
000
1
G23
00
0000
1
G31
 
00000
1
G12
 
]
{
σx 
σy 
σz 
σyz 
σzx 
σxy 
}
    (6.3)

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)

{
єx 
єy 
γxy 
} =[
 N,x
 0N,y 
 N,yN,x 
] {
u 
v 
}     (6.4)

and stress by

{
σ єx 
σ єy 
σ γxy 
} =[
 d1,1 N,x+d1,3 N,yd1,2 N,y+d1,3 N,x 
 d2,1 N,x+d2,3 N,yd2,2 N,y+d2,3 N,x 
 d3,1 N,x+d3,3 N,yd3,2 N,y+d3,3 N,x 
] {
u 
v 
}     (6.5)

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

D=
E
1−ν2
[
 1ν0
ν10
00
1−ν
2
 
]     (6.6)

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

D=
E(1−ν
(1+ν)(1−2ν)
[
 1
ν
1−ν
0
ν
1−ν
10
00
1−2ν
2(1−ν)
 
]     (6.7)

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

{
p,x 
p,y 
p,z 
} =[
 N,x 
 N,y 
 N,z 
] {
p 
}     (6.8)

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

     
 







σ11
σ22
σ33
τ23
τ13
τ12
 
 















C11C12C1300
 C22C2300
  C3300
   C440
 (s)  C55
     C66
  
 















є11
є22
є33
γ23
γ13
γ12
 
 







      (6.9)

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

     
  σ33=0 ⇒ є33=−
1
C33

C13є11+C23є22
,
      (6.10)

and for isotropic material,

     
  σ33=0 ⇒ є33=−
ν
1−ν

є1122
.
      (6.11)

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

     
 





σ11
σ22
σ12
σ23
σ13
 





=  






Q11Q1200
Q12Q2200
        00Q660
 000Q44
 0000Q55 
  
 












є11
є22
γ12
γ23
γ13
 





.
      (6.12)

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

     
Qij=





Cij− 
Ci3Cj3
C33
if i,j=1,2,
Cij if i,j=4,5,6.
 
      (6.13)

The reduced elastic law for an isotropic plate becomes,

     
 
 



σ11
σ22
τ12
 
 



E
(1−ν2)







1ν
ν1
 00
1−ν
2
 
 










є11
є22
γ12
 
 



,
      (6.14)

and

     
 
 



τ23
τ13
 
 



E
2(1+ν)



 1
0
  
 






 γ23
γ13
  
 



.
      (6.15)

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

     
є=









u1,1+x3β1,1
1
2
(u1,2+u2,1+x31,22,1))
1
2
1+w,1)
 u2,2+x3β2,2
1
2
2+w,2)
(s) 
  
 









.  
      (6.16)

So, the strain vector is written,

     


є

=







є11m+x3κ11
є22m+x3κ22
γ12m+x3κ12
γ23
γ13
 








,  
      (6.17)

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

     
єm=



u1,1
u2,2
u1,2+u2,1
 
 



,  κ=



β1,1
β2,2
β1,22,1
 



,  γ=

β2+w,2
β1+w,1
 

,  
      (6.18)


Note that the engineering notation with γ12=u1,2+u2,1 is used here rather than the tensor notation with є12=(u1,2+u2,1)/2 . Similarly κ121,22,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 γ,

     
 



N
M
Q
 
 



=  


AB
BD
00F 
 






єm
κ
γ
 



.
      (6.19)

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

     
Nα β=
ht


hb
σα β dx3,   Mα β=
ht


hb
x3 σα β dx3,   Qα 3=
ht


hb
σα 3 dx3,  
      (6.20)


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]

Aij=
ht


hb
Qij dx3,
Bij=
ht


hb
x3 Qij dx3,
  
Dij=
ht


hb
x32 Qij dx3,
Fij=
ht


hb
Qij dx3
      (6.21)

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

     
   Hij=kijFij,      (6.22)

where kij are correction factors. Reddy's 3rd order theory brings to kij=2/3. Very commonly, enriched 3rd order theory are used, and kij 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 N11, N22 and shear force N12 become

     




N11
N22
N12
 



Eh
1−ν2






1ν0
 10
(s) 
1−ν
2
  









u1,1
u2,2
u1,2+u2,1
 



,  
      (6.23)

the 2 bending moments M11, M22 and twisting moment M12

     




M11
M22
M12
 



Eh3
12(1−ν2)






1ν0
 10
(s) 
1−ν
2
  









β1,1
β2,2
β1,22,1
 



,  
      (6.24)

and the out-of-plane shearing forces Q23 and Q13,

     


Q23
Q13
 

Eh
2(1+ν)


10
01
  



β2+w,2
β1+w,1
 

.  
      (6.25)

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

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

     
Aij=h[Qij]    Bij=x3(0)h [Qij]     Cij= (x3(0)2h+h3/12) [Qij]     Fij=h[Qij]        (6.26)

where is clearly appears that the constitutive matrix is a polynomial function of h, h3, x3(0)2h and x3(0)h. If the ply thickness is kept constant, the constitutive law is a polynomial function of 1,x3(0),x3(0)2.

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) :

{
S 
E 
} = {
єx 
єy 
єz 
γyz 
γzx 
γxy 
Ex 
Ey 
Ez 
} =[
 N,x00
 0N,y0
 00N,z
 0N,zN,y
 N,z0N,x
 N,yN,x0
 000N,x 
 000N,y 
 000N,z 
] {
u 
v 
w 
φ 
}     (6.27)

where φ is the electric potential (V).

The constitutive law associated with this strain state is given by

  {
T 
D 
} = [
CEeT 
e−єS 
]{
S 
E 
}     (6.28)

in which D is the electrical displacement vector (a density of charge in Cb/m2), T is the mechanical stress vector (N/m2). CE is the matrix of elastic constants at zero electric field (E=0, short-circuited condition, see section 6.1.1 for formulas (there CE 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 :

  {
S 
D 
} = [
sEdT 
dєT 
]{
T 
E 
}     (6.29)

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

[d] = [
d11d12d13d14d15d16 
 d21d22d23d24d25d26 
 d31d32d33d34d35d36 
 ]     (6.30)

Matrices [e] and [d] are related through

  [e] = [d] [ CE]     (6.31)

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

S] = [
є11Sє12Sє13S 
 є21Sє22Sє23S 
 є31Sє32Sє33S  
 ]     (6.32)

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 :

S]= [єT] − [d] [e]T     (6.33)

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

єr = 
є
є0
    (6.34)

є0 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

T] = [
є11T0
 0є22T
 00є33T 
 ]     (6.35)

For PVDF, the matrix of piezoelectric constants is given by

[d] = [
00000
 00000
 d31d32d3300
 ]     (6.36)

and for PZT materials :

[d] = [
0000d15
 000d240
 d31d32d3300
 ]     (6.37)

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    Ezi). Ezi is assumed to be constant over the thickness hi of the layer and is therefore given by Ezi=−Δ φi/hi 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

{
єxx 
єyy 
2 єxy 
κxx 
κyy 
2 κxy 
γxz 
γyz
Ez1 
... 
Ezn  
}=[
 N,x00000...
 0N,y0000...0
 N,yN,x0000...
 0000N,x0...
 000N,y00...0
 000N,xN,y0...
 00N,x0N0...
 00N,yN00...
00000
1
h1
...
...............0...
1
hn
 
  
] {
u 
v 
w 
ru 
rw 
Δ φ1 
... 
Δ φn 
}     (6.38)

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 :

  {
N 
M 
Q 
Dz1 
... 
Dzn 
} = [
AB0G1T...GnT 
BD0zm1 G1T...zmn GnT 
00FH1T...HnT 
G1zm1 G1H1−є1...
.........0...
Gnzmn GnHn0...−єn
] {
є 
κ
γ 
Ez1 
... 
Ezn  
}     (6.39)

where Dzi is the electric displacement in piezoelectric layer (assumed constant and in the z-direction), zmi is the distance between the midplane of the shell and the midplane of piezoelectric layer i, and Gi, Hi are given by

Gi = {
e.1e.2
}i [Rs]i     (6.40)
Hi = {
e.4e.5 
}i [R]i     (6.41)

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

e31e320

}i . If the piezoelectric is used in shear mode, the polarization is in the x or y-direction, therefore Gi=0, and Hi = {0 e15 }i or Hi = {e24 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.

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

[Rs] = [
cos2 θsin2 θsinθ cosθ 
sin2 θcos2 θ− sinθ cos θ 
−2 sinθ cosθ2 sinθ cosθcos2 θ −  sin2 θ 
]     (6.42)
[R] = [
cosθ−sinθ 
sinθcosθ 
]     (6.43)

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

 
 


Ω0
 (ρ0 u″, δ v) + 
 


Ω0
 S : δ e = 
 


Ω0
 f . δ v   ∀ δ v     (6.44)

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

  Fi,j = I + ui,j = δij+{N,j}T{qi}     (6.45)

where the N,j is the derivative of the shape functions with respect to Cartesian coordinates at the current integration point and qi 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 ({e11 e22 e33 2e23 2e31 2e12}) to tensor [eij] 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(FTFI) and its variation

  deij = 
FT dF
Sym = 
Fki {N,j}T{qk}
Sym     (6.46)

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

 
 


Ω
 S : δ e = 
 


Ω
 {δ qk}T{N,jFki Sij     (6.47)

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

  KG=
 


Ω
 Sij δ uk,i ul,j + 
 


Ω
 de : 
2 W
∂ e2
 : δ e     (6.48)

which using the notation ui,j = {N,j}T{qi} leads to

  KGe=
 


Ω
 {δ qm} {N,l


Fmk
2 W
∂ e2
ijkl Fni + Slj


{N,j} {dqn}     (6.49)

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

  S = D:e(u) = 
2 W
∂ e2
:e(u) =  λ Tr(eI + 2µ e      (6.50)

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/∂ e2 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

  [K(qn)]{qn+1qn} =   R(qn) = 
 


Ω
 f . dv − 
 


Ω0
 S(qn) : δ e      (6.51)

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 [eT] = [α] (TT0), 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

  S = C:(e − eT) =  λ Tr(eI + 2µ e − (C:[α])(TT0)      (6.52)

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:eT) : δ 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

W = C1(J1−3) + C2(J2−3) + K(J3−1)2,     (6.53)
W = C1(J1−3) + C2(J2−3) + K(J3−1) − (C1 + 2C2 + K)ln(J3),     (6.54)

where (J1,J2,J3) are the so-called reduced invariants of the Cauchy-Green tensor

C=I+2e,     (6.55)

linked to the classical invariants (I1,I2,I3) by

J1=I1 I
1
3
 
3
,   J2=I2 I
2
3
 
3
,    J3=I
1
2
 
3
,     (6.56)

where one recalls that

I1=tr C,    I2=
1
2
[(tr C)2tr C2],    I3=det C.     (6.57)

Note : this definition of energy based on reduced invariants is used to have the hydrostatic pressure given directly by p=−K(J3−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:

∂ I1
∂ Cij
 = δij,   
∂ I2
∂ Cij
 = I1δijCij,   
∂ I3
∂ Cij
 = I3 Cij−1,     (6.58)

where (Cij−1) denotes the coefficients of the inverse matrix of (Cij). For second-order derivatives we have:

 
2 I1
∂ Cij∂ Ckl
 = 0,   
2 I2
∂ Cij∂ Ckl
 = −δikδjlijδkl,   
2 I3
∂ Cij∂ Ckl
 = Cmn єikmєjln,     (6.59)

where the єijk coefficients are defined by





    єijk=0when 2 indices coincide
 =1when (i,j,k) even permutation of (1,2,3)
 =−1when (i,j,k) odd permutation of (1,2,3)
    (6.60)

Note: when the strain components are seen as a column vector (“engineering strains”) in the form (e11,e22,e33,2e23,2e31,2e12)′, the last two terms of (6.59) thus correspond to the following 2 matrices








    011000
    101000
    110000
    000−1/200
    0000−1/20
    00000−1/2







,     (6.61)







    0C33C22C2300
    C330C110C130
    C22C11000C12
    −C2300C11/2C12/2C13/2
    0C130C12/2C22/2C23/2
    00C12C13/2C23/2C33/2







.     (6.62)

We finally use chain-rule differentiation to compute

S = 
∂ W
∂ e
 = 
 
k
 
∂ W
∂ Ik
 
∂ Ik
∂ e
,     (6.63)
2 W
∂ e2
 =
 
k
 
∂ W
∂ Ik
 
2 Ik
∂ e2
 
k
 
l
 
2 W
∂ Ik∂ Il
 
∂ Ik
∂ e
∂ Il
∂ e
.     (6.64)

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

x
∂ t
u
∂ t
+Ω∧(p+u)

and its acceleration is

2x
∂ t2
 =
2u
∂ t2
+
Ω
∂ t
∧(p+u)+2Ω
u
∂ t
+ΩΩ∧(p+u)

Ω is the rotation vector of the structure with

Ω=


ωx 
ωy 
ωz



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

[Ω] = [
0−ωzωy  
ωz0−ωx  
−ωyωx0
]

The speed can be rewritten

x
∂ t
u
∂ t
+[Ω](p+u

and the acceleration becomes

2x
∂ t2
 =
2u
∂ t2
+
∂[Ω]
∂ t
(p+u)+2[Ω]
u
∂ t
+[Ω]2(p+u)

In this expression appear

S0e is an element of the mesh of the initial configuration S0 whose density is ρ0. [N] is the matrix of shape functions on these elements, one defines the following elementary matrices

 
[Dge] =
 


S0e
 2ρ0 [N][Ω] [NdS0e   gyroscopic coupling
[Kae] =
 


S0e
 ρ0 [N]
∂[Ω]
∂ t
[N]  dS0e   centrifugal acceleration
[Kge] =
 


S0e
 ρ0 [N][Ω]2 [N]  dS0e   centrifugal softening/stiffening
    (6.65)

6.1.11  Centrifugal follower forces

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

δ Wω
 


Ω
ρ ω2 R(x) δ vR,     (6.66)

where δ vR designates the radial component (in deformed configuration) of δv. One assumes that the rotation axis is along ez. Noting nR = 1/R {x1  x2   0}T, one then has

δ vRnR·δ v.     (6.67)

Thus the non-linear stiffness term is given by

dδ Wω= − 
 


Ω
ρ ω2 (dR δ vR + R dδ vR).     (6.68)

One has dR=nR· dx(= dxR) and dδ vR = dnR·δ v, with

dnR=−
dR
R
nR + 
1
R
{dx1  dx2   0}T.

Thus, finally

dδ Wω= − 
 


Ω
ρ ω2 (du1 δ v1 + du2 δ v2).     (6.69)

Which gives

du1 δ v1 + du2 δ v2= {δ qα}T {N}{N}T {d qα},     (6.70)

with α=1,2.

6.1.12  Poroelastic materials

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

Domain and variables description:

ΩPoroelastic domain
∂ΩBounding surface of poroelastic domain
nUnit external normal of ∂Ω
uSolid phase displacement vector
uFFluid phase displacement vector uF = φ/ρ22ω2p − ρ1222u
pFluid phase pressure
σStress tensor of solid phase
σtTotal stress tensor of porous material σt=σ−φ(1+Q/R)pI

Weak formulation, for harmonic time dependence at pulsation ω:

 
 
 


Ω
σ(u) : є(δ u)  dΩ − ω2 
 


Ω
ρ  u.δ u  dΩ −
 


Ω
φ
α
∇ p.δ u  dΩ 
        −
 


Ω
φ


1+
Q
R



p∇.δ u  dΩ −
 


∂Ω
t(u).n).δ u  dS =0   ∀ δ u
    (6.71)
 
 
 


Ω
φ2
αρoω2
∇ p.∇δ p  dΩ −
 


Ω
φ2
R
p δ p  dΩ −
 


Ω
φ
α
 u.∇ δ p  dΩ 
        −
 


Ω
φ


1+
Q
R



δ p∇. u  dΩ −
 


∂Ω
φ(uFu).n δ p  dS = 0  ∀ δ p
    (6.72)

Matrix formulation, for harmonic time dependence at pulsation ω:

  [
K−ω2MC1C2
C1TC2T
1
ω2
FKp 
] {
u
p
} = {
Fst
Ff
}     (6.73)

where the frequency-dependent matrices correspond to:

 
 


Ω
σ(u):є(δ udΩ
⇒δ uT K u
 
 


Ω
ρ  u.δ u dΩ
⇒δ uT M u
 
 


Ω
φ2
αρo
∇ p.∇δ p
⇒δ pT Kp p 
 
 


Ω
φ2
R
p δ p
⇒δ pT F p 
 
 


Ω
φ
α
∇ p.δ u  dΩ
⇒δ uT C1 p 
 
 


Ω
φ


1+
Q
R



p∇.δ u  dΩ
⇒δ uT C2 p 
 
 


∂Ω
t(u).n).δ u  dS
⇒δ uT Fst
 
 


∂Ω
φ(uFu).n δ p  dS
⇒δ pT Ff
    (6.74)

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

  [
(1+iηs)K−ω2ρM
φ
α
C1− φ


1+
Q
R



C2 
φ
α
C1T−φ


1+ 
Q
R



C2T
1
ω2
φ2
R
F− 
φ2
αρo
Kp 
] {
u
p
} = {
Fst
Ff
}     (6.75)

where the matrices marked with bars are frequency independent:

 K=(1+iηs)KMM
C1=
φ
α
C1 
C2


1+
Q
R



C2
F=
φ2
R
F
Kp=
φ2
αρo
Kp 
    (6.76)

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
GShear modulus of the skeleton
νPoisson coefficient of the skeleton
ηsStructural loss factor of the skeleton
ρoFluid 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:

Po=1,01× 105 PaAmbient pressure
Pr=0.71Prandtl number

Poroelastic specific (frequency dependent) variables:

ρ11Apparent density of solid phase ρ11 = (1−φ)ρ−ρ12
ρ22Apparent density of fluid phase ρ22 = φρo−ρ12
ρ12Interaction apparent density ρ12=−φρo−1)
ρEffective density of solid phase ρ = ρ11 − (ρ12)222
ρ11Effective density of solid phase ρ11 = ρ11 + b/ iω
ρ22Effective density of fluid phase ρ22 = ρ22 + b/ iω
ρ12Interaction effective density ρ12 = ρ12 − b/ iω
bViscous damping coefficient b = φ2σ1 + i2ηρoω/ σ2Λ2φ2
γCoupling coefficient γ = φ(ρ1222 − Q/R)
QElastic coupling coefficient
   Biot formulation Q=1−φ−Kb/Ks/1−φ−Kb/KsKs/Kfφ Ks
   Approximation from Kb/Ks<<1 Q=(1−φ)Kf
RBulk modulus of air in fraction volume
   Biot formulation R=φ2Ks/1−φ−Kb/KsKs/Kf
   Approximation from Kb/Ks<<1 R=φKf
KbBulk modulus of porous material in vacuo Kb= 2G(1+ν)/ 3(1−2ν)
KsBulk modulus of elastic solid
   est. from Hashin-Shtrikman's upper bound Ks=1+2φ/1−φKb
KfEffective bulk modulus of air in pores Kf= Po/ 1 − γ −1/ γ α ′
α ′Function in Kf (Champoux-Allard model) α ′ = 1 + ωT/ 2iω(1+ iω/ ωT)1/2
ωTThermal characteristic frequencyωT= 16η/ PrΛ′2ρo

To add here:

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

  
 


Ω
 (ρ c θ)(vdx +
 


Ω
 (K grad  θ)(grad  vdx
 


∂Ω
 αθ v  dγ = 
 
 


Ω
 f v   dx +  
 


∂Ω
 (g+α θextv  dγ 
 
  ∀ v ∈ H1(Ω)
    (6.77)

with

Test case

One considers a solid square prism of dimensions Lx,Ly, Lz 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=16 Γi = ∂ Ω), are subject to the following boundary conditions and loads

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

 θ(x,y,z) =−
f
k
 x2 + θext + 
 f Lx2
2k
 = 
g(x)
α
= 25 − 
x2
20

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


Figure 6.1: Temperature distribution along the x-axis


©1991-2016 by SDTools
Previous Up Next