Contents     Functions     Index     Previous Next     PDF

5.2  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.

5.2.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,x 0 0
0 N,y 0
0 0 N,z
0 N,z N,y
N,z 0 N,x
N,y N,x 0
] {
u
v
w
}     (5.1)
and stress by
{
x
y
z
yz
zx
xy
} =[
d1,1 N,x+d1,5 N,z+d1,6 N,y d1,2 N,y+d1,4 N,z+d1,6 N,x d1,3 N,z+d1,4 N,y+d1,5 N,x
d2,1 N,x+d2,5 N,z+d2,6 N,y d2,2 N,y+d2,4 N,z+d2,6 N,x d2,3 N,z+d2,4 N,y+d2,5 N,x
d3,1 N,x+d3,5 N,z+d3,6 N,y d3,2 N,y+d3,4 N,z+d3,6 N,x d3,3 N,z+d3,4 N,y+d3,5 N,x
d4,1 N,x+d4,5 N,z+d4,6 N,y d4,2 N,y+d4,4 N,z+d4,6 N,x d4,3 N,z+d4,4 N,y+d4,5 N,x
d5,1 N,x+d5,5 N,z+d5,6 N,y d5,2 N,y+d5,4 N,z+d5,6 N,x d5,3 N,z+d5,4 N,y+d5,5 N,x
d6,1 N,x+d6,5 N,z+d6,6 N,y d6,2 N,y+d6,4 N,z+d6,6 N,x d6,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. In particular volume elements inherited from MODULEF order shear stresses differently xy, yz, zx (these elements are obtained by setting p_solid integ value to zero. In fe_stres 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-
1
]
0
0
[
G 0 0
0 G 0
0 0 G
]
]     (5.2)

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

For constitutive law building, see p_solid.

5.2.2  2D elasticity

With m_elasticsubtype 4, p_solid deals with 2D mechanical volumes with strain defined by (see q4p constants)
{
x
y
xy
} =[
N,x 0
0 N,y
N,y N,x
] {
u
v
}     (5.3)
and stress by
{
x
y
xy
} =[
d1,1 N,x+d1,3 N,y d1,2 N,y+d1,3 N,x
d2,1 N,x+d2,3 N,y d2,2 N,y+d2,3 N,x
d3,1 N,x+d3,3 N,y d3,2 N,y+d3,3 N,x
] {
u
v
}     (5.4)

For isotropic plane stress (p_solid form=1), one has
D=
E
1-2
[
1 0
1 0
0 0
1-
2
]     (5.5)

For isotropic plane strain (p_solid form=0), one has
D=
E(1-
(1+)(1-2)
[
1
1-
0
1-
1 0
0 0
1-2
2(1-)
]     (5.6)

5.2.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
}     (5.7)

This replaces the earlier flui4 ... elements.

5.2.4  Classical lamination theory

5.2.5  Piezo-electric volumes

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,x 0 0 0
0 N,y 0 0
0 0 N,z 0
0 N,z N,y 0
N,z 0 N,x 0
N,y N,x 0 0
0 0 0 -N,x
0 0 0 -N,y
0 0 0 -N,z
] {
u
v
w
}     (5.8)
where is the electric potential (V).

The constitutive law associated with this strain state is given by
{
T
D
} = [
CE -eT
e S
]{
S
E
}     (5.9)
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 5.2.1 for formulas (there CE is noted D). Note that using -E rather than E makes the constitutive law symetric.

Alternatively, one can use the constitutive equations written in the following manner :
{
S
D
} = [
sE dT
d T
]{
T
E
}     (5.10)
In which sE is the matrix of mechanical compliances, [d] is the matrix of piezoeletric constants (m/V=Cb/N):
[d] = [
d11 d12 d13 d14 d15 d16
d21 d22 d23 d24 d25 d26
d31 d32 d33 d34 d35 d36
]     (5.11)

Matrices [e] and [d] are related through
[e] = [d] [ CE]     (5.12)

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
]     (5.13)

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

For this reason, the input value for the computation should be [T].
Also notice that usually relative permitivities are given in datasheets:
r =
0
    (5.15)
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] = [
11T 0 0
0 22T 0
0 0 33T
]     (5.16)

For PVDF, the matrix of piezoelectric constants is given by
[d] = [
0 0 0 0 0 0
0 0 0 0 0 0
d31 d32 d33 0 0 0
]     (5.17)

and for PZT materials :
[d] = [
0 0 0 0 d15 0
0 0 0 d24 0 0
d31 d32 d33 0 0 0
]     (5.18)

5.2.6  Geometric non-linearity for hyperelasticity

The principle of virtual work in non-linear total Lagrangian formulation for an hyperelastic medium is
ó
õ
 


W0
(0 u'', v) + ó
õ
 


W0
S : e = ó
õ
 


W0
f . v " v     (5.19)
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}     (5.20)
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 realy 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 use 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(FTF -I) and its variation
deij = ( FT dF ) Sym = ( Fki {N,j}T{ qk} ) Sym     (5.21)

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


W
S : e = ó
õ
 


W
{ qk}T{N,j} Fki Sij     (5.22)
and the tangent stiffness matrix (its derivative with respect to the current position) can be written as
KG= ó
õ
 


W
Sij uk,i vk,j + ó
õ
 


W
de :
2 W
e2
: e     (5.23)
which using the notation ui,j = {N,j}T{qi} leads to
KGe= ó
õ
 


W
{ qm} {N,l} æ
ç
ç
è		 Fmk
2 W
e2
ijkl Fni + Slj ö
÷
÷
ø		 {N,j} {dqn}     (5.24)
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(e) I + 2µ e     (5.25)
the building of the constitutive law matrix D is performed in p_solid BuildConstit for isotropic, orthotropic and full anisotropic materials. For hyperelastic materials 2 W/ e2 is not constant and is computed at each integration point.

For a geometric non-linear static computation one will thus solve
[K(qn)]{qn+1-qn} = ó
õ
 


W
f . dv - ó
õ
 


W0
S : e     (5.26)
where external forces f are assumed to be non following.

5.2.7  Hyperelasticity

For hyperelastic media S= W/ e with W the hyperelastic energy.
of_mk.c MatrixIntegration 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,     (5.27)
W = C1(J1-3) + C2(J2-3) + K(J3-1) - (C1 + 2C2 + K)ln(J3),     (5.28)

where (J1,J2,J3) are the so-called reduced invariants of the Cauchy-Green tensor
C=I+2e,     (5.29)
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
,     (5.30)
where one recalls that
I1=tr C, I2=
1
2
[(tr C)2-tr C2], I3=det C.     (5.31)

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
= I1ij-Cij,
I3
Cij
= I3 Cij-1,     (5.32)
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
= -ikjl+ijkl,
2 I3
Cij Ckl
= Cmn ikmjln,     (5.33)
where the ijk coefficients are defined by
ì
í
î
ijk =0 when 2 indices coincide
  =1 when (i,j,k) even permutation of (1,2,3)
  =-1 when (i,j,k) odd permutation of (1,2,3)
    (5.34)
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 (5.33) thus correspond to the following 2 matrices
æ
ç
ç
ç
ç
ç
ç
è
0 1 1 0 0 0
1 0 1 0 0 0
1 1 0 0 0 0
0 0 0 -1/2 0 0
0 0 0 0 -1/2 0
0 0 0 0 0 -1/2
ö
÷
÷
÷
÷
÷
÷
ø		 ,     (5.35)
æ
ç
ç
ç
ç
ç
ç
è
0 C33 C22 -C23 0 0
C33 0 C11 0 -C13 0
C22 C11 0 0 0 -C12
-C23 0 0 -C11/2 C12/2 C13/2
0 -C13 0 C12/2 -C22/2 C23/2
0 0 -C12 C13/2 C23/2 -C33/2
ö
÷
÷
÷
÷
÷
÷
ø		 .     (5.36)

We finally use chain-rule differentiation to compute
S =
W
e
=
 
k
W
Ik
Ik
e
,     (5.37)
2 W
e2
=
 
k
W
Ik
2 Ik
e2
+
 
k
 
l
2 W
Ik Il
Ik
e
Il
e
.     (5.38)

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.

5.2.8  Gyroscopic effects

Written by Arnaud Sternchüss 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
+WÙ(p+u)
and its acceleration is
2x
t2
=
2u
t2
+
W
t
Ù(p+u)+2WÙ
u
t
+WÙWÙ(p+u)
W is the rotation vector of the structure with
W= é
ê
ê
ë
x
y
z
ù
ú
ú
û
in a (x,y,z) orthonormal frame. The skew-symmetric matrix [W] is defined such that
[W] = [
0 -z y
z 0 -x
-y x 0
]
The speed can be rewritten
x
t
=
u
t
+[W](p+u)
and the acceleration becomes
2x
t2
=
2u
t2
+
[W]
t
(p+u)+2[W]
u
t
+[W]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
20 [N][W] [N] dS0e   gyroscopic coupling
[Kae] =
ó
õ
 


S0e
0 [N]
[W]
t
[N]  dS0e   centrifugal acceleration
[Kge] =
ó
õ
 


S0e
0 [N][W]2 [N]  dS0e   centrifugal softening/stiffening
    (5.39)

5.2.9  Centrifugal follower forces

This is the embrio of the theory for the future implementation of centrifugal follower forces.
W= ó
õ
 


W
2 R(x) vR,     (5.40)
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
vR= nR· v.     (5.41)

Thus the non-linear stiffness term is given by
-d W= - ó
õ
 


W
2 (dR vR + R d vR).     (5.42)
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= - ó
õ
 


W
2 (du1 v1 + du2 v2).     (5.43)

Which gives
du1 v1 + du2 v2= { q}T {N}{N}T {d q},     (5.44)
with =1,2.

©1991-2007 by SDTools
Previous Up Next