next up previous contents index
Next: II.2 The ``ClaDb'' class Up: II. Composite Reference Manual Previous: II.0 Introduction   Contents   Index

Subsections


II.1 Theoretical background

The purpose of this Chapter is to summarize the classical laminate theory, and to provide the information needed by the user of composite classes to understand a few conventions that have been assumed for the programming of Classical Laminate Analysis in FeResPost (axes, angles, numbering of layers,...).

The programmer will find a presentation of the classical laminate theory that follows closely what is programmed in C++ language in FeResPost. However those who are interested in studying the theory, or who are not familiar with it are referred to more extensive presentations of the classical laminate theory [Gay97,Pal99]. Only for the out-of-plane shear behavior of the laminate, is the presentation original, even though inspired by information found in [Sof04a].

The Chapter is organized as follows:


II.1.1 Conventions

Figure II.1.1 represents some of the conventions used for the definition of laminate in FeResPost. The laminate coordinate system is defined in such a way that the $ z$ axis is perpendicular to laminate surface and points from bottom to top surface. $ x$ and $ y$ vectors are parallel to the laminate plies. Plies are numbered from the bottom to the top surface. If $ k$ is the index of a ply, one considers that it is limited by $ z$ coordinates $ z_{k-1}$ and $ z_k$ . The origin of the coordinate system is located at mid laminate thickness so that if $ t$ is the laminate total thickness and the laminate has $ N$ plies, top surface of the laminate is located at $ z_N=+t/2$ and bottom surface at $ z_0=-t/2$ . (The plies are numbered from 1 to $ N$ .)

Figure II.1.1: Conventions for the numbering of plies in a laminate and for the definition of the laminate coordinate system. (A six plies laminate is represented.)
\begin{figure}\centerline{%
\input{claManual/theory/xfig/laminate_axes.pstex_t}}
\end{figure}

In the laminate, the plies are characterized by their material, thickness and by their angle in the laminate. Figure II.1.2 shows the convention for the orientation of a ply in a laminate. $ 1$ , $ 2$ and $ 3$ are the ply axes; $ x$ , $ y$ and $ z$ are the laminate axes. Of course, because only 2D laminates are considered here, one has always $ z=3$ . If $ \xi$ is the angle of the ply in the laminate, this angle is as represented in Figure II.1.2: a positive angle corresponds to a rotation from axis $ x$ to axis $ y$ . If the angle is $ 0^\circ$ , the first ply axis 1 is parallel to the first laminate axis $ x$ .

Figure II.1.2: Conventions for the orientation of plies in the laminate.)
\begin{figure}\centerline{%
\input{claManual/theory/xfig/ply_axes.pstex_t}}
\end{figure}


II.1.2 Rotation in XY plane and algebraic notations

One common operation in classical laminate analysis is to rotate vectors, tensors and matrices. One summarizes here the operations one uses in the rest of this Chapter and in FeResPost. This rotation is represented in Figure II.1.3.

Figure: Rotation of angle $ \theta$ around the origin in $ O_{12}$ plane.
\begin{figure}\centerline{%
\input{claManual/theory/xfig/rotation.pstex_t}}
\end{figure}

II.1.2.0.1 Rotation of base vectors

For such a rotation, the vectors $ e$ $ _x$ and $ e$ $ _y$ are expressed as a function of $ e$ $ _1$ and $ e$ $ _2$ as:

$\displaystyle \left(\begin{array}{c}
\mbox{\boldmath$e$}_x \\
\mbox{\boldmat...
...}{c}
\mbox{\boldmath$e$}_1 \\
\mbox{\boldmath$e$}_2
\end{array}\right)\ .
$

To simplify the notations, one introduces the symbols $ c=\cos\theta$ and $ s=\sin\theta$ . Also, one prefers to write the more general 3D version of the transformation:

$\displaystyle \left(\begin{array}{c} \mbox{\boldmath$e$}_x \\ \mbox{\boldmath$e...
...h$e$}_1 \\ \mbox{\boldmath$e$}_2 \\ \mbox{\boldmath$e$}_3 \end{array}\right)\ .$ (II.11)

The inverse relation corresponds to a rotation of angle $ -\theta$ and is obtained by changing the signs of the sinuses in the rotation matrix:

$\displaystyle \left(\begin{array}{c} \mbox{\boldmath$e$}_1 \\ \mbox{\boldmath$e...
...h$e$}_x \\ \mbox{\boldmath$e$}_y \\ \mbox{\boldmath$e$}_z \end{array}\right)\ .$ (II.12)

II.1.2.0.2 Transformation of vector and tensor components

The expressions (II.1.1) and (II.1.2) can be used to transform the components of vectors. For example:

$\displaystyle \left(\begin{array}{c} V_x \\ V_y \\ V_z \end{array}\right) = \le...
...ay}\right) \cdot \left(\begin{array}{c} V_1 \\ V_2 \\ V_3 \end{array}\right)\ .$ (II.13)

For the transformation of 2D tensors, the transformation matrix is used twice. For example, a Cauchy stress tensor is transformed as follows:

$\displaystyle \left(\begin{array}{ccc} \sigma_{xx} & \sigma_{xy} & \sigma_{zx} ...
...t(\begin{array}{ccc} c & -s & 0 \\ s & c & 0 \\ 0 & 0 & 1 \end{array}\right)\ .$ (II.14)

II.1.2.0.3 Matricial notations

As the Cauchy stress tensor is symmetric, expression (II.1.4) is more conveniently written in a matricial form as follows:

$\displaystyle \left(\begin{array}{c} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{zz} ...
...22} \\ \sigma_{33} \\ \tau_{23} \\ \tau_{13} \\ \tau_{12} \end{array}\right)\ .$ (II.15)

The same expression applies to the components of the strain tensor, which is also symmetric:

$\displaystyle \left(\begin{array}{c}
\epsilon_{xx} \\
\epsilon_{yy} \\
\ep...
...
\epsilon_{23} \\
\epsilon_{13} \\
\epsilon_{12}
\end{array}\right)\ .
$

However, unfortunately, the classical laminate analysis is universally written using angular shear components for the strain tensor:

$\displaystyle \gamma_{ij}=2\epsilon_{ij}\ (i\neq j)\ .
$

Using the angular components, the matricial expression to be used for the rotation becomes:

$\displaystyle \left(\begin{array}{c} \epsilon_{xx} \\ \epsilon_{yy} \\ \epsilon...
...epsilon_{33} \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \end{array}\right)\ .$ (II.16)

An interesting aspect of the transformations (II.1.5) and (II.1.6) is that one can apply the transformation separately on sub-groups of components: This contributes to justify some of the simplifications of the classical laminate analysis; among others, the decoupling of in-plane and flexural deformation of the laminate on one hand from the out-of-plane shear on the other hand. The third direction is systematically neglected: $ \sigma_{33}=\epsilon_{33}=0$ . The inverse of relation (II.1.7) is obviously;

$\displaystyle \left(\begin{array}{c} \sigma_{11} \\ \sigma_{22} \\ \tau_{12} \e...
...(\begin{array}{c} \sigma_{xx} \\ \sigma_{yy} \\ \tau_{xy} \end{array}\right)\ ,$ (II.110)

$\displaystyle \left(\begin{array}{c} \epsilon_{11} \\ \epsilon_{22} \\ \gamma_{...
...egin{array}{c} \epsilon_{xx} \\ \epsilon_{yy} \\ \gamma_{xy} \end{array}\right)$ (II.111)

II.1.2.0.4 Introduction of a short notation

In order to simplify the notations, one introduces the following notations:

$\displaystyle \left[T_+\right]_{(\theta)} =\left(\begin{array}{ccc} c^2 & s^2 & -2cs \\ s^2 & c^2 & 2cs \\ cs & -cs & c^2-s^2 \end{array}\right)\ ,$ (II.112)

$\displaystyle \left[T_-\right]_{(\theta)} =\left(\begin{array}{ccc} c^2 & s^2 & 2cs \\ s^2 & c^2 & -2cs \\ -cs & cs & c^2-s^2 \end{array}\right)\ ,$ (II.113)

$\displaystyle \left[T_+^\prime\right]_{(\theta)} =\left(\begin{array}{ccc} c^2 & s^2 & cs \\ s^2 & c^2 & -cs \\ -2cs & 2cs & c^2-s^2 \end{array}\right)\ ,$ (II.114)

$\displaystyle \left[T_-^\prime\right]_{(\theta)} =\left(\begin{array}{ccc} c^2 & s^2 & -cs \\ s^2 & c^2 & cs \\ 2cs & -2cs & c^2-s^2 \end{array}\right)\ ,$ (II.115)

$\displaystyle \left[S_-\right]_{(\theta)} =\left(\begin{array}{cc} c & -s \\ s & c \end{array}\right)\ ,$ (II.116)

$\displaystyle \left[S_+\right]_{(\theta)} =\left(\begin{array}{cc} c & s \\ -s & c \end{array}\right)\ .$ (II.117)

These matrices are not independent. For example:

$\displaystyle \left[T_-\right]_{(\theta)}
=\left(\left[T_+\right]_{(\theta)}\right)^{-1}
=\left[T_+\right]_{(-\theta)}\ ,
$

$\displaystyle \left[T_-\right]_{(\theta)}
=\left(\left[T_+\right]_{(\theta)}\right)^{-1}
=\left[T_+\right]_{(-\theta)}\ ,
$

$\displaystyle \left[T_+^\prime\right]_{(\theta)}
=\left(\left[T_+\right]_{(\theta)}\right)^{T}\ ,
$

$\displaystyle \left[S_-\right]_{(\theta)}
=\left(\left[S_+\right]_{(\theta)}\r...
...\right]_{(\theta)}\right)^{-1}
=\left(\left[S_+\right]_{(-\theta)}\right)\ .
$

The transformations of the components of strain tensor (II.1.8) and stress tensor (II.1.10) are then written:

$\displaystyle \left\{\epsilon\right\}_{xy} = \left[T_+^\prime\right]_{(\theta)} \cdot \left\{\epsilon\right\}_{12}\ ,$ (II.118)

$\displaystyle \left\{\sigma\right\}_{12} = \left[T_+\right]_{(\theta)} \cdot \left\{\sigma\right\}_{xy}\ ,$ (II.119)

$\displaystyle \left\{\epsilon\right\}_{12} = \left[T_-^\prime\right]_{(\theta)} \cdot \left\{\epsilon\right\}_{xy}\ ,$ (II.120)

$\displaystyle \left\{\sigma\right\}_{xy} = \left[T_-\right]_{(\theta)} \cdot \left\{\sigma\right\}_{12}\ .$ (II.121)

Similarly, for the out-of-plane shear stresses and strains one writes the following relations:

$\displaystyle \left\{\gamma_s\right\}_{xy} = \left[S_+\right]_{(\theta)} \cdot \left\{\gamma_s\right\}_{12}\ ,$ (II.122)

$\displaystyle \left\{\tau_s\right\}_{xy} = \left[S_+\right]_{(\theta)} \cdot \left\{\tau_s\right\}_{12}\ ,$ (II.123)

$\displaystyle \left\{\gamma_s\right\}_{12} = \left[S_-\right]_{(\theta)} \cdot \left\{\gamma_s\right\}_{xy}\ ,$ (II.124)

$\displaystyle \left\{\tau_s\right\}_{12} = \left[S_-\right]_{(\theta)} \cdot \left\{\tau_s\right\}_{xy}\ .$ (II.125)


II.1.3 Materials and plies

One summarizes in this section a few results that are commonly found in composite literature.

II.1.3.1 Plies

Each ply is defined by:

The orientation of the ply is given by an angle $ \theta$ .

II.1.3.2 Materials and constitutive equations

When a material is used in the definition of a laminate, assumptions are done about the axes defined in the laminate. Axes 1 and 2 are parallel to the laminate plane and axis 3 is orthogonal to the laminate.

The classical laminate analysis is based on the assumption that the relation between stress and strain tensors is linear. Then, as these two tensors are symmetric, a $ 6\times 6$ matrix contains all the elastic coefficients defining the material:

$\displaystyle \left(\begin{array}{c} \sigma_{11} \\  \sigma_{22} \\  \sigma_...
...33} \\  \gamma_{23} \\  \gamma_{31} \\  \gamma_{12} \end{array}\right) \ .$ (II.126)

One shows that, because the peculiar choice of angular strain tensor components, the matrix $ \left[C_{ijkl}\right]$ containing the elastic coefficients is symmetric. Therefore, the matrix has only 21 independent coefficients. $ \left[C_{ijkl}\right]$ is the stiffness matrix of the material.

Equation (II.1.26)can be reversed as follows:


In expression (II.1.27), one added the thermo-elastic and moisture expansion terms in previous expression. They are characterized by CTE and CME tensors noted $ \left\{\alpha_{kl}\right\}$ and $ \left\{\beta_{kl}\right\}$ respectively. Note that shear components of these two tensors are angular components. Practically, it does not matter much as most materials have zero shear components for CTE or CME tensors. $ \left[c_{ijkl}\right]$ is the compliance matrix of the material. Obviously $ \left[c_{ijkl}\right]=\left[C_{ijkl}\right]^{-1}$ . One often defines laminates with orthotropic materials: For an isotropic material, the definition of $ E$ and either $ G$ or $ \nu$ is sufficient to characterize the material. Then one has:

$\displaystyle E_1=E_2=E_3=E\ ,
$

$\displaystyle G_{12}=G_{23}=G_{13}=G\ ,
$

$\displaystyle \nu_{12}=\nu_{23}=\nu_{13}=\nu\ .
$

$ E$ , $ G$ and $ \nu$ satisfy the following relation:

$\displaystyle E-2G(1+\nu)=0\ .
$

Finally, one introduces shorter notations that allow to rewrite expressions (II.1.26) and (II.1.27) respectively as follows:

$\displaystyle \left\{\epsilon\right\} =  \left[c\right] \cdot \left\{\sigma\right\} +\Delta T \left\{\alpha\right\} +\Delta H \left\{\beta\right\} \ ,$ (II.127)

$\displaystyle \left\{\sigma\right\} =  \left[C\right] \cdot \left\{\epsilon...
...ft\{\alpha\right\} -\Delta H \left[C\right] \cdot \left\{\beta\right\} \ .$ (II.128)

One introduces also the ``Mechanical Strain Tensor'' estimated as follows:

$\displaystyle \left\{\epsilon^{\text{Mech}}\right\} =  \left[c\right] \cdot \left\{\sigma\right\} \ .$ (II.129)

This new strain tensor differs from the one defined by (II.1.29) by the fact that no thermo-elastic or hygro-elastic contribution is taken into account to estimate its components. It is the strain that corresponds to the actual material stress, when no thermo-elastic or hygro-elastic expansion is considered. This ``Mechanical Strain Tensor'' is also sometimes called ``Equivalent Strain Tensor''.

II.1.3.3 In-plane properties

One considers the properties of the ply in a plane parallel to the laminate. Then the constitutive equation (II.1.28) reduces to:

$\displaystyle \left(\begin{array}{c} \epsilon_{11} \\  \epsilon_{22} \\  \ga...
...{array}{c} \beta_{11} \\  \beta_{22} \\  \beta_{12} \end{array}\right) \ .$ (II.130)

The indices in this notation are integers and indicate that the corresponding properties are given in ply coordinate system. The equation (II.1.32) is written more shortly as follows:

$\displaystyle \left\{\epsilon\right\}_{\text{ply}} =  \left[c\right]_{\text{p...
...{\alpha\right\}_{\text{ply}} +\Delta H \left\{\beta\right\}_{\text{ply}} \ .$ (II.131)

One introduces in (II.1.33) the material in-plane compliance matrix $ \left[c\right]_{\text{ply}}$ . In order to avoid too complicated notations, one uses the same notations as for the $ 6\times 6$ full material compliance matrix introduced in (II.1.29). This will be done systematically for the in-plane matricial and vectorial quantities in the rest of the document ( $ \left[c\right]$ , $ \left[C\right]$ , $ \left\{\alpha\right\}$ , $ \left\{\beta\right\}$ , $ \left\{\epsilon\right\}$ ,...

The inverse of expression (II.1.33) is noted:

$\displaystyle \left\{\sigma\right\}_{\text{ply}} =  \left[C\right]_{\text{ply...
...a H \left[C\right]_{\text{ply}} \cdot \left\{\beta\right\}_{\text{ply}} \ .$ (II.132)

In (II.1.34) one introduces the in-plane stiffness matrix $ \left[C\right]_{\text{ply}}=\left[c\right]_{\text{ply}}^{-1}$ .

Plies are characterized by their orientation in the laminate. Let $ \xi$ be the angle of the ply in the laminate axes. Then, the laminate axes are obtained by rotating the ply axes by an angle $ -\xi$ . Equations (II.1.33) and (II.1.34) are expressed in the laminate coordinate system as follows:

$\displaystyle \left\{\epsilon\right\}_{\text{lam}}
=
\left[T_+^\prime\righ...
...\left[T_+^\prime\right]_{(-\xi)}\cdot\left\{\beta\right\}_{\text{ply}}
\ ,
$

\begin{multline}
\left\{\sigma\right\}_{\text{lam}}
=
\left[T_-\right]_{(...
...t{ply}}
\cdot\left\{\beta\right\}_{\text{ply}}\nonumber
\ .
\end{multline}

This leads to the new expression in laminate axes:

$\displaystyle \left\{\epsilon\right\}_{\text{lam}}
=
\left[c\right]_{\text...
...a\right\}_{\text{lam}}
+\Delta H
\left\{\beta\right\}_{\text{lam}}
\ ,
$

$\displaystyle \left\{\sigma\right\}_{\text{lam}}
=
\left[C\right]_{\text{l...
...left[C\right]_{\text{lam}}
\cdot
\left\{\beta\right\}_{\text{lam}}
\ ,
$

where one introduces new notations for in-plane ply properties rotated by an angle $ -\xi$ (in laminate axes):

$\displaystyle \left[c\right]_{\text{lam}}= \left[T_+^\prime\right]_{(-\xi)}\cdot\left[c\right]_{\text{ply}} \cdot\left[T_+\right]_{(-\xi)}\ ,$ (II.133)

$\displaystyle \left\{\alpha\right\}_{\text{lam}}= \left[T_+^\prime\right]_{(-\xi)}\cdot\left\{\alpha\right\}_{\text{ply}}\ ,$ (II.134)

$\displaystyle \left\{\beta\right\}_{\text{lam}}= \left[T_+^\prime\right]_{(-\xi)}\cdot\left\{\beta\right\}_{\text{ply}}\ ,$ (II.135)

$\displaystyle \left[C\right]_{\text{lam}}= \left[T_-\right]_{(-\xi)}\cdot\left[C\right]_{\text{ply}} \cdot\left[T_-^\prime\right]_{(-\xi)}\ ,$ (II.136)

$\displaystyle \left[C\right]_{\text{lam}}\cdot\left\{\alpha\right\}_{\text{lam}...
...i)}\cdot\left[C\right]_{\text{ply}} \cdot\left\{\alpha\right\}_{\text{ply}}\ ,$ (II.137)

$\displaystyle \left[C\right]_{\text{lam}}\cdot\left\{\beta\right\}_{\text{lam}}...
...xi)}\cdot\left[C\right]_{\text{ply}} \cdot\left\{\beta\right\}_{\text{ply}}\ .$ (II.138)

When a matrix is transformed as in (II.1.35) or a vector as in (II.1.36), one says that they are rotated with $ \left[T_+^\prime\right]_{(-\xi)}$ rotation matrix.

II.1.3.4 Out-of-plane shear properties

One makes developments similar to those in the previous section. The out-of-plane shear constitutive equations are written as follows:

$\displaystyle \left\{\tau_s\right\}_{\text{ply}} =  \left[g\right]_{\text{ply...
...H \left[g\right]_{\text{ply}} \cdot \left\{\beta_s\right\}_{\text{ply}} \ ,$ (II.139)

$\displaystyle \left\{\gamma_s\right\}_{\text{ply}} =  \left[g^{-1}\right]_{\t...
...pha_s\right\}_{\text{ply}} +\Delta H \left\{\beta_s\right\}_{\text{ply}} \ .$ (II.140)

If $ \xi$ is the angle of the ply in the laminate, the previous relations can be written in laminate axes by rotating them by an angle $ -\xi$ . For example:

$\displaystyle \left\{\gamma_s\right\}_{\text{lam}}
=
\left[S_+\right]_{(-\xi)}
\cdot
\left\{\gamma_s\right\}_{\text{ply}}\ ,
$

$\displaystyle \left\{\tau_s\right\}_{\text{lam}}
=
\left[S_+\right]_{(-\xi)}
\cdot
\left\{\tau_s\right\}_{\text{ply}}\ ,
$

$\displaystyle \left\{\alpha_s\right\}_{\text{lam}}
=
\left[S_+\right]_{(-\xi)}
\cdot
\left\{\alpha_s\right\}_{\text{ply}}\ ,
$

$\displaystyle \left\{\beta_s\right\}_{\text{lam}}
=
\left[S_+\right]_{(-\xi)}
\cdot
\left\{\beta_s\right\}_{\text{ply}}\ ,
$

$\displaystyle \left\{\gamma_s\right\}_{\text{ply}}
=
\left[S_-\right]_{(-\...
...
\left[S_+\right]_{(\xi)}
\cdot
\left\{\gamma_s\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\tau_s\right\}_{\text{ply}}
=
\left[S_-\right]_{(-\xi...
...
\left[S_+\right]_{(\xi)}
\cdot
\left\{\tau_s\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\alpha_s\right\}_{\text{ply}}
=
\left[S_-\right]_{(-\...
...
\left[S_+\right]_{(\xi)}
\cdot
\left\{\alpha_s\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\beta_s\right\}_{\text{ply}}
=
\left[S_-\right]_{(-\x...
...
\left[S_+\right]_{(\xi)}
\cdot
\left\{\beta_s\right\}_{\text{lam}}\ .
$

Then, one makes consecutive transformations of relations (II.1.41) as follows:

\begin{multline}
\left[S_+\right]_{(\xi)}
\cdot
\left\{\tau_s\right\}_{\te...
...)}
\cdot
\left\{\beta_s\right\}_{\text{lam}}\nonumber
\ ,
\end{multline}

\begin{multline}
\left\{\tau_s\right\}_{\text{lam}}
=
\left[S_+\right]_{(...
...)}
\cdot
\left\{\beta_s\right\}_{\text{lam}}\nonumber
\ ,
\end{multline}

$\displaystyle \left\{\tau_s\right\}_{\text{lam}}
=
\left[g\right]_{\text{l...
...ft[g\right]_{\text{lam}}
\cdot
\left\{\beta_s\right\}_{\text{lam}}
\ ,
$

where one introduced:

$\displaystyle \left[g\right]_{\text{lam}}
=
\left[S_+\right]_{(-\xi)}
\cdot
\left[g\right]_{\text{ply}}
\cdot
\left[S_+\right]_{(\xi)}
\ .
$

One says that tensor $ \left[g\right]$ is rotated by matrix $ \left[S_+\right]_{(-\xi)}$ which corresponds to the expression of the shear stiffness tensor in a new coordinate system obtained by rotating the previous one by an angle $ -\xi$ .

The transformation of the out-of-plane shear compliance tensor by the same angle $ -\xi$ is made with the same expression as for the stiffness tensor:

$\displaystyle \left[g^{-1}\right]_{\text{lam}}
=
\left[S_+\right]_{(-\xi)}...
...
\left[g^{-1}\right]_{\text{ply}}
\cdot
\left[S_+\right]_{(\xi)}
\ .
$


II.1.4 Thickness and mass of laminate

The total laminate thickness is the sum of the thickness of each of its plies:

$\displaystyle t_{\text{lam}}=\sum_{k=1}^N e_k\ .
$

Correspondingly the surfacic mass is given by:

$\displaystyle sm_{\text{lam}}=\sum_{k=1}^N \rho e_k\ .
$

And the laminate average density is:

$\displaystyle \rho_{\text{lam}}=sm_{\text{lam}}/t_{\text{lam}}\ .
$


II.1.5 In-plane and flexural laminate behavior

The classical laminate analysis is based on the assumption that in-plane and flexural behavior of the laminate is not related to out-of-plane shear loading. The corresponding laminate properties can be studied separately. The same remark is true for the load response calculation. In this section, the in-plane and flexural behavior of laminates are studied.

In this section the thermal and moisture expansions are not taken into account. The out-of-plane shear properties and loading of laminates is also discussed in a separate section. One summarizes the results of classical laminate analysis. The reader shall refer to the literature if more information on the developments that lead to these results are needed. In this section, the different equations are written in laminate axes and the corresponding indices are noted $ x$ and $ y$ .

Laminate compliance and stiffness matrices relate the in-plane forces and bending moments on one hand to the average strain and curvatures on the other hand. Those different quantities are defined as follows:

Note that average strain tensor, as well as the true tensor are not ``real'' tensors because their shear components (i.e. non-diagonal components are angular components.)

The relations between the four tensors are then given by two equations:

$\displaystyle \left\{N\right\}_{\text{lam}}=
\left[A\right]_{\text{lam}}\cdot\...
...lam}}
+\left[B\right]_{\text{lam}}\cdot\left\{\kappa\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{M\right\}_{\text{lam}}=
\left[B\right]_{\text{lam}}\cdot\...
...lam}}
+\left[D\right]_{\text{lam}}\cdot\left\{\kappa\right\}_{\text{lam}}\ .
$

One defines below the different matrices and vectors introduced in these equations: All the new matrices and vectors are obtained by summation of the ply contributions. In order to obtain the ply stiffness matrix in laminate axes $ \left[C\right]_{\text{lam}}^k$ and the ply thermo-elastic CTE coefficients in thermo-elastic axes $ \left\{\alpha\right\}_{\text{lam}}^k$ , one uses the transformations (II.1.36) and (II.1.38) respectively. Note however that if a ply is characterized by an orientation $ \theta$ wrt to laminate axes, the rotation of ply properties must be of an angle $ -\theta$ .

The laminate compliance matrices $ \left[a\right]_{\text{lam}}$ , $ \left[b\right]_{\text{lam}}$ and $ \left[d\right]_{\text{lam}}$ are obtained by inversion of the $ 6\times 6$ $ \left[ABBD\right]_{\text{lam}}$ matrix:

$\displaystyle \left[\begin{array}{cc} \left[a\right]_{\text{lam}} & \left[b\rig...
...[B\right]_{\text{lam}} & \left[D\right]_{\text{lam}} \end{array}\right]^{-1}\ ,$ (II.142)

Then the average laminate strain and its curvature tensor can be calculated as follows:

$\displaystyle \left\{\epsilon^0\right\}_{\text{lam}}=
\left[a\right]_{\text{la...
...text{lam}}
+\left[b\right]_{\text{lam}}\cdot\left\{M\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\kappa\right\}_{\text{lam}}=
{\left[b\right]}^{\rm T}_{\t...
...text{lam}}
+\left[d\right]_{\text{lam}}\cdot\left\{M\right\}_{\text{lam}}\ .
$

One often calculates equivalent moduli corresponding to the calculated stiffness matrices $ A$ and $ D$ . This is done as follows (we follow the expressions presented in [Pal99]):


II.1.6 Out-of-plane shear of laminate

One presents one version of the out-of-plane shear theory for laminates based on information found in Chapter 13 of [Sof04a]. Only, one presents here a more general version of the calculation that takes $ xy$ and $ yz$ components of the out-of-plane shear stress into account at the same time.

II.1.6.1 Out-of-plane shear equilibrium equations

In Chapter 13 of [Sof04a] one considers the equilibrium in direction $ x$ of a small portion of the material (Figure II.1.4) of lengths $ dx$ and $ dz$ respectively:

$\displaystyle \cfrac{\partial\tau_{xz}}{\partial z}+ \cfrac{\partial\sigma_{xx}}{\partial x}=0\ .$ (II.143)

Similarly, the equilibrium of a portion $ dx$ of the full laminate is given globally by the expression:

$\displaystyle Q_{xz}-\cfrac{\partial M_{xx}}{\partial x}=0\ .$ (II.144)

Figure II.1.4: Out-of-plane XZ shear equilibrium in a laminate (from [Sof04a]).
\begin{figure}\centerline{%
\input{claManual/theory/xfig/transverseShear.pstex_t}}
\end{figure}

Then, in Chapter 13 of [Sof04a], developments are done to calculate the relations between $ Q_{xz}$ and $ \tau_{xz}$ . All the developments are based on the local $ x$ equilibrium relation.

In this document, a more general presentation of the out-of-plane shear behavior of laminates is done. The $ x$ and $ y$ components of in-plane local equilibrium are written as follows:

$\displaystyle \sigma_{xx,x}+\tau_{xy,y}+\tau_{xz,z}=0\ ,$ (II.145)

$\displaystyle \tau_{xy,x}+\sigma_{yy,y}+\tau_{yz,z}=0\ .$ (II.146)

Correspondingly, a global equilibrium is expressed by the two equations:

$\displaystyle M_{xx,x}+M_{xy,y}=Q_{xz}\ ,$ (II.147)

$\displaystyle M_{xy,x}+M_{yy,y}=Q_{yz}\ .$ (II.148)

Those equations shall be developed and will ultimately allow the calculation of $ \tau_{xz}$ and $ \tau_{yz}$ from the global shear $ Q_{xz}$ and $ Q_{yz}$ .

II.1.6.2 Triangular distribution of in-plane stresses

In most expressions below, the components of tensors are expressed in laminate axes. Therefore, the ``lam'' underscore is often added to the different quantities used in the equations.

First, one calculates the components of Cauchy stress tensor. However, a few simplifying assumptions shall be done. The strain tensor components are calculated from the laminate average strain tensor and curvature as follows:

$\displaystyle \left\{\epsilon\right\}_{\text{lam}}=
\left\{\epsilon^0\right\}_{\text{lam}}
+z\left\{\kappa\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\epsilon^0\right\}_{\text{lam}}=
\left[a\right]_{\text{la...
...text{lam}}
+\left[b\right]_{\text{lam}}\cdot\left\{M\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\kappa\right\}_{\text{lam}}=
{\left[b\right]}^{\rm T}_{\t...
...text{lam}}
+\left[d\right]_{\text{lam}}\cdot\left\{M\right\}_{\text{lam}}\ .
$

(The thermo-elastic contributions have been neglected.) In most out-of-plane shear theories presented in the literature, one assumes a decoupling between in-plane load response and out-of-plane shear response. This allows us to neglect a few terms in the equations:

$\displaystyle \left\{\epsilon^0\right\}_{\text{lam}}\approx
\left[b\right]_{\text{lam}}\cdot\left\{M\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\kappa\right\}_{\text{lam}}\approx
\left[d\right]_{\text{lam}}\cdot\left\{M\right\}_{\text{lam}}\ .
$

One then writes a simple expression of the in-plane laminate deformation tensor:

$\displaystyle \left\{\epsilon\right\}_{\text{lam}}=
\left\{\epsilon^0\right\}_{\text{lam}}
+z\left\{\kappa\right\}_{\text{lam}}\ .
$

Then, the components of Cauchy stress tensor are given by:

$\displaystyle \left\{\sigma\right\}_{\text{lam}}(x,y,z)= \left[C\right]_{\text{...
...+z\left[d\right]_{\text{lam}}\right) \cdot\left\{M\right\}_{\text{lam}}(x,y)\ .$ (II.149)

In this last expression, the $ 3\times 3$ matrix $ \left[C\right]_{\text{lam}}(z)$ corresponds to the plies in-plane moduli expressed in laminate axes. It depends on $ z$ because the components generally change from one ply to another. However, one shall assume that the components of the moduli matrix are constant in each ply.

Note that, in the local and global equilibrium relations (II.1.47) to (II.1.50), only partial derivatives of bending moments and Cauchy stress tensor components appear. One assumes the decoupling between the out-of-plane shear behavior and the absolute bending in laminate. However, as shown by expressions (II.1.49) and (II.1.50), the out-of-plane shear is related to the gradient of bending moment. One derives equation (II.1.51) wrt to $ x$ and $ y$ :

$\displaystyle \left\{\begin{array}{c}
\sigma_{xx,x}(x,y,z) \\
\sigma_{yy,x}(...
...y) \\
M_{yy,x}(x,y) \\
M_{xy,x}(x,y)
\end{array}\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\begin{array}{c}
\sigma_{xx,y}(x,y,z) \\
\sigma_{yy,y}(...
...y) \\
M_{yy,y}(x,y) \\
M_{xy,y}(x,y)
\end{array}\right\}_{\text{lam}}\ .
$

At this point, one no longer needs to assume a dependence of the gradient of bending moments wrt $ x$ and $ y$ . The same is true for the gradient of Cauchy stress tensor. One also introduces a new notation:

$\displaystyle \left[F\right]_{\text{lam}}(z)=
\left[C\right]_{\text{lam}}(z)\cdot
\left(\left[b\right]_{\text{lam}}
+z\left[d\right]_{\text{lam}}\right)\ .
$

Then, the components of Cauchy stress tensor gradient are obtained from the components of bending moments gradient with the following expression:

$\displaystyle \left\{\begin{array}{c} \sigma_{xx,x}(z) \\ \sigma_{yy,x}(z) \\ \...
... \\ \hline M_{xx,y} \\ M_{yy,y} \\ M_{xy,y} \end{array}\right\}_{\text{lam}}\ .$ (II.150)


II.1.6.3 Out-of-plane shear stress partial derivative equations

Note that the global equilibrium equation (II.1.49) and (II.1.50) do not contain the components $ M_{xx,y}$ and $ M_{yy,x}$ of the bending moments tensor. Similarly, the local equilibrium equations do not contain the components $ \sigma_{xx,y}$ and $ \sigma_{yy,x}$ of the Cauchy stress tensor. Then, these components can be considered as nil without modifying the result of the developments. The corresponding lines and columns could be removed from the equations (II.1.52).

Actually, one can do better than that. The local equilibrium equations (II.1.47) and (II.1.48) are rewritten as follows:

$\displaystyle \left\{\begin{array}{c} \tau_{xz,z}(z) \\ \tau_{yz,z}(z) \end{arr...
...y}(z) \\ \sigma_{yy,y}(z) \\ \tau_{xy,y}(z) \end{array}\right\}_{\text{lam}}\ .$ (II.151)

Similarly, one writes:

$\displaystyle \left\{\begin{array}{c} M_{xx,x} \\ M_{yy,x} \\ M_{xy,x} \\ \hlin...
..._{xx,x} \\ M_{xy,x} \\ M_{yy,y} \\ M_{xy,y} \end{array}\right\}_{\text{lam}}\ .$ (II.152)

The substitution of (II.1.52) and (II.1.54) in (II.1.53) leads to the following expression:


One would like to eliminate the four $ x$ and $ y$ partial derivative of bending moment tensor components in the previous expression. For this, one uses the global equilibrium equations (II.1.49) and (II.1.50). This leaves some arbitrary choice in the determination of dependence wrt out-of-plane shear. For example:

$\displaystyle M_{xx,x}$ $\displaystyle =\mu_xQ_{xz}\ ,$    
$\displaystyle M_{xy,y}$ $\displaystyle =(1-\mu_x)Q_{xz}\ ,$    
$\displaystyle M_{yy,y}$ $\displaystyle =\mu_yQ_{yz}\ ,$    
$\displaystyle M_{xy,x}$ $\displaystyle =(1-\mu_y)Q_{yz}\ .$    

$\displaystyle \left\{\begin{array}{c}
M_{xx,x} \\
M_{xy,x} \\
M_{yy,y} \\ ...
...ft\{\begin{array}{c}
Q_{xz} \\
Q_{yz}
\end{array}\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\begin{array}{c} M_{xx,x} \\ M_{yy,x} \\ M_{xy,x} \\ \hline M_{xx,y} \\ M_{yy,y} \\ M_{xy,y} \end{array}\right\}_{\text{lam}}$ $\displaystyle = \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 &...
...ot \left\{\begin{array}{c} Q_{xz} \\ Q_{yz} \end{array}\right\}_{\text{lam}}\ ,$    
  $\displaystyle = \left[\begin{array}{cc} \mu_x & 0 \\ 0 & 0 \\ 0 & (1-\mu_x) \\ ...
...ot \left\{\begin{array}{c} Q_{xz} \\ Q_{yz} \end{array}\right\}_{\text{lam}}\ .$    

This allows to find a new expression of the relation between bending moment gradients and out-of-plane shear stress. One first calculates a new matrix as follows:

$\displaystyle \left[X\right]_{\text{lam}}(z)
=-
\left[\begin{array}{rrr\vert ...
... \\
\hline
0 & 0 \\
0 & \mu_y \\
(1-\mu_y) & 0
\end{array}\right]\ .
$

The choice $ \mu_x=\mu_y=1/2$ gives more symmetry to the relation between $ \left\{\tau\right\}_{\text{lam}}(z)$ and $ \left\{Q\right\}_{\text{lam}}(z)$ . This choice leads to:

$\displaystyle \left[X\right]_{\text{lam}}(z)
=-\cfrac{1}{2}
\left[\begin{arra...
...0 \\
0 & 1 \\
\hline
0 & 0 \\
0 & 1 \\
1 & 0
\end{array}\right]\ .
$

The choice $ \mu_x=\mu_y=1$ is the default choice in FeResPost. The same choice seems to have been done in other software, like ESAComp. In the rest of the document, the following notations are used:

$\displaystyle \left[P_1\right]= \left[\begin{array}{rrr\vert rrr} 1 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 & 0 \\ \end{array}\right]\ ,$ (II.153)

$\displaystyle \left[P_2(\mu_x,\mu_y)\right]= \left[\begin{array}{cc} \mu_x & 0 ...
...& (1-\mu_x) \\ \hline 0 & 0 \\ 0 & \mu_y \\ (1-\mu_y) & 0 \end{array}\right]\ .$ (II.154)

$ \left[X\right]_{\text{lam}}(z)$ is a $ 2\times 2$ matrix that relates the out-of-plane shear stress components partial derivatives wrt $ z$ to the out-of-plane shear force components:

$\displaystyle \left\{\begin{array}{c} \tau_{xz,z}(z) \\ \tau_{yz,z}(z) \end{arr...
...ot \left\{\begin{array}{c} Q_{xz} \\ Q_{yz} \end{array}\right\}_{\text{lam}}\ .$ (II.155)

Note that by replacing the partial derivative of bending moments by out-of-plane shear force components, one has done an approximation that involves much arbitrariness. Actually, Expression (II.1.55) shows that out-of-plane shear stresses in the laminate are related to partial derivative of bending moments, but not directly to out-of-plane shear force components. The approximations that have been done allow us to define a theory that is more in line with what is usually presented in the literature.

Its main advantage is that the decoupling between in-plane and flexural load response on one side and out-of-plane shear response on the other side can be completed. Moreover, estimates of the partial derivatives of bending moments are not always available. For example, no finite element results corresponding to these unknown are commonly available.

Its main disadvantage is that the laminate out-of-plane shear equations lose their objectivity wrt rotations around axis $ z$ as illustrated by the example described in section IV.5.5. This example also allows to estimate the effects of the approximation on the precision of results given by the theory.

The matrix $ \left[X\right]_{\text{lam}}(z)$ depends on $ z$ for two reasons: because of the triangular distribution of strains through the thickness, and because material moduli depend on plies material and orientation. In a given ply of index $ k$ , one has:

$\displaystyle \left[X^k\right]_{\text{lam}}(z)
=\left[X_0^k\right]_{\text{lam}}
+z\left[X_1^k\right]_{\text{lam}}\ ,
$

in which the components of the two matrices $ \left[X_0^k\right]_{\text{lam}}$ and $ \left[X_1^k\right]_{\text{lam}}$ are constant. Similarly one can write a polynomial expression for $ \left[F\right]_{\text{lam}}(z)$ if one splits the definition by plies:

$\displaystyle \left[F^k\right]_{\text{lam}}(z)$ $\displaystyle =\left[C^k\right]_{\text{lam}}\cdot \left(\left[b\right]_{\text{lam}} +z\left[d\right]_{\text{lam}}\right)\ ,$    
  $\displaystyle =\left(\left[C^k\right]_{\text{lam}}\cdot\left[b\right]_{\text{la...
...+\left(\left[C^k\right]_{\text{lam}}\cdot\left[d\right]_{\text{lam}}\right)z\ ,$    
  $\displaystyle =\left[F_0^k\right]_{\text{lam}} +\left[F_1^k\right]_{\text{lam}}z\ .$    

Of course, one has the two relations:

$\displaystyle \left[X_0^k\right]_{\text{lam}}=
-
\left[P_1\right]
\cdot\left...
...ht]_{\text{lam}}
\end{array}\right]
\cdot
\left[P_2(\mu_x,\mu_y)\right]\ ,
$

$\displaystyle \left[X_1^k\right]_{\text{lam}}=
-
\left[P_1\right]
\cdot\left...
...ht]_{\text{lam}}
\end{array}\right]
\cdot
\left[P_2(\mu_x,\mu_y)\right]\ .
$

II.1.6.4 Integration of out-of-plane shear stress equation

The out-of-plane shear stress components are obtained by integration of expression (II.1.58) along the thickness. This leads to the following expression:

$\displaystyle \left\{\begin{array}{c}
\tau_{xz}(z) \\
\tau_{yz}(z)
\end{arr...
...c}
\tau_{xz}(z_0) \\
\tau_{yz}(z_0)
\end{array}\right\}_{\text{lam}}
\ .
$

One assumes zero shear stress along the bottom surface of the laminate. This corresponds to a free surface, or at least to a surface that receives no contact forces in direction $ x$ and $ y$ . This assumption leads to the following expression:

$\displaystyle \left\{\begin{array}{c} \tau_{xz}(z) \\ \tau_{yz}(z) \end{array}\right\}_{\text{lam}}$ $\displaystyle =\left(\int_{z_0}^z\left[X\right]_{\text{lam}}(z)\hspace{0.4em}\t...
...t \left\{\begin{array}{c} Q_{xz} \\ Q_{yz} \end{array}\right\}_{\text{lam}} \ ,$    
  $\displaystyle =\left[Y\right]_{\text{lam}}(z) \cdot \left\{\begin{array}{c} Q_{xz} \\ Q_{yz} \end{array}\right\}_{\text{lam}} \ ,$ (II.156)

where a new matrix has been introduced:

$\displaystyle \left[Y\right]_{\text{lam}}(z)
=\int_{z_0}^z\left[X\right]_{\text{lam}}(z)\hspace{0.4em}\text{d}z\ .
$

An explicit expression of the integrated matrix is calculated ply-by-ply, from bottom layer to top layer. If $ z_{k-1}\leq z \leq z_k$ :

$\displaystyle \left[Y\right]_{\text{lam}}(z)$ $\displaystyle =\left[Y^k\right]_{\text{lam}}(z)\ ,$    
  $\displaystyle =\int_{z_0}^z\left[X^k\right]_{\text{lam}}(z)\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\left[Y^k\right]_{\text{lam}}(z_{k-1}) +\int_{z_{k-1}}^z\left[X^k\right]_{\text{lam}}(z)\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\left[Y^k\right]_{\text{lam}}(z_{k-1}) +\int_{z_{k-1}}^z\left[X_...
...right]_{\text{lam}} +z\left[X_1^k\right]_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\left[Y^k\right]_{\text{lam}}(z_{k-1}) +\left[X_0^k\right]_{\tex...
...ft(z-z_{k-1}\right) +\left[X_1^k\right]_{\text{lam}}\cfrac{z^2-z_{k-1}^2}{2}\ ,$    
  $\displaystyle =\left[Y_0^k\right]_{\text{lam}} +\left[Y_1^k\right]_{\text{lam}}z +\left[Y_2^k\right]_{\text{lam}}z^2\ .$ (II.157)

In expression (II.1.60), one introduced new symbols that are calculated as follows:

$\displaystyle \left[Y_0^k\right]_{\text{lam}}
=\left[Y^k\right]_{\text{lam}}(z...
..._{\text{lam}}z_{k-1}
-\left[X_1^k\right]_{\text{lam}}\cfrac{z_{k-1}^2}{2}\ ,
$

$\displaystyle \left[Y_1^k\right]_{\text{lam}}
=\left[X_0^k\right]_{\text{lam}}\ ,
$

$\displaystyle \left[Y_2^k\right]_{\text{lam}}
=\cfrac{\left[X_1^k\right]_{\text{lam}}}{2}\ .
$

Note that the expression above involve the a priori unknown quantity $ \left[Y^k\right]_{\text{lam}}(z_{k-1})$ . To calculate this expression, one uses the continuity of $ \left[Y^k\right]_{\text{lam}}(z)$ across ply interfaces:

$\displaystyle \left[Y^k\right]_{\text{lam}}(z_{k-1})
=\left[Y^{k-1}\right]_{\text{lam}}(z_{k-1})\ ,
$

This relation corresponds to the continuity of out-of-plane shear stress at ply interfaces. One develops the relation as follows.

$\displaystyle \left[Y_0^k\right]_{\text{lam}}$ $\displaystyle =\left[Y^{k-1}\right]_{\text{lam}}(z_{k-1}) -\left[X_0^k\right]_{\text{lam}}z_{k-1} -\left[X_1^k\right]_{\text{lam}}\cfrac{z_{k-1}^2}{2}\ ,$    
  $\displaystyle =\left[Y_0^{k-1}\right]_{\text{lam}} +\left[Y_1^{k-1}\right]_{\te...
...ht]_{\text{lam}}z_{k-1} -\left[X_1^k\right]_{\text{lam}}\cfrac{z_{k-1}^2}{2}\ ,$    
  $\displaystyle =\left[Y_0^{k-1}\right]_{\text{lam}} +\left(\left[Y_1^{k-1}\right...
...ht]_{\text{lam}} -\cfrac{\left[X_1^k\right]_{\text{lam}}}{2}\right)z_{k-1}^2\ .$ (II.158)

The last line of this development allows to calculate recursively the components of $ \left[Y_0^k\right]_{\text{lam}}$ from bottom ply to top ply. For bottom ply, the condition $ \tau_{xz}(z_0)=\tau_{yz}(z_0)=0$ leads to the following expressions:

$\displaystyle \left[Y^1\right]_{\text{lam}}(z_0)
=\left[Y_0^1\right]_{\text{la...
...left[Y_1^1\right]_{\text{lam}}z_0
+\left[Y_2^1\right]_{\text{lam}}z_0^2=0\ .
$

$\displaystyle \left[Y_0^1\right]_{\text{lam}}
=-\left[Y_1^1\right]_{\text{lam}}z_0
-\left[Y_2^1\right]_{\text{lam}}z_0^2\ .
$

Then, it becomes possible to calculate recursively the $ \left[Y_0^i\right]_{\text{lam}}$ matrices.

One checks easily that the condition $ \tau_{xz}(z_0)=\tau_{yz}(z_0)=0$ ensures also that $ \tau_{xz}(z_N)=\tau_{yz}(z_N)=0$ . Indeed, one has:

$\displaystyle \left[Y\right]_{\text{lam}}(z_N)$ $\displaystyle =\int_{z_0}^{z_N}\left[X\right]_{\text{lam}}(z)\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =- \left[P_1\right] \cdot \int_{z_0}^{z_N} \left[\begin{array}{c\...
...nd{array}\right] \hspace{0.4em}\text{d}z \cdot \left[P_2(\mu_x,\mu_y)\right]\ ,$    
  $\displaystyle =- \left[P_1\right] \cdot \left[\begin{array}{c\vert c} \int_{z_0...
...space{0.4em}\text{d}z \end{array}\right] \cdot \left[P_2(\mu_x,\mu_y)\right]\ .$    

The last line of previous equation contains twice the integral of $ \left[F\right]_{\text{lam}}(z)$ along the laminate thickness. One develops this integral as follows:

$\displaystyle \int_{z_0}^{z_N}\left[F\right]_{\text{lam}}(z)\hspace{0.4em}\text{d}z$ $\displaystyle =\int_{z_0}^{z_N}\left[C\right]_{\text{lam}}(z) \cdot\left(\left[...
...ht]_{\text{lam}} +z\left[d\right]_{\text{lam}}\right)\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\left(\int_{z_0}^{z_N}\left[C\right]_{\text{lam}}(z)\hspace{0.4e...
...\text{lam}}(z)\hspace{0.4em}\text{d}z\right)\cdot\left[d\right]_{\text{lam}}\ ,$    
  $\displaystyle =\left[A\right]_{\text{lam}}\cdot\left[b\right]_{\text{lam}} +\left[B\right]_{\text{lam}}\cdot\left[d\right]_{\text{lam}}\ .$    

On the other hand, equation (II.1.44) allows to write:

$\displaystyle \left[\begin{array}{cc} \left[A\right] & \left[B\right] \\ \left[...
... \left[b\right] \\ {\left[b\right]}^{\rm T} & \left[d\right] \end{array}\right]$ $\displaystyle = \left[\begin{array}{cc} \left[A\right]\cdot\left[a\right] +\lef...
...ht]\cdot\left[b\right] +\left[D\right]\cdot\left[d\right] \end{array}\right]\ ,$    
  $\displaystyle = \left[\begin{array}{cc} \left[I\right] & \left[0\right] \\ \left[0\right] & \left[I\right] \end{array}\right]\ .$    

(The ``lam'' subscript has been omitted for shortness sake.) The identification of the right upper corner of the last expression with the integration of $ \left[F\right]_{\text{lam}}(z)$ along the laminate thickness shows that this integral must be zero. Consequently, one also has:

$\displaystyle \left[Y\right]_{\text{lam}}(z_N)$ $\displaystyle =0\ ,$    
$\displaystyle \tau_{xz}(z_N)$ $\displaystyle =0\ ,$    
$\displaystyle \tau_{yz}(z_N)$ $\displaystyle =0\ .$    

It is interesting to remark that the ply out-of-plane shear moduli have not been used in the calculations to obtain (II.1.59). The out-of-plane shear stresses depend only on out-of-plane shear forces and ply in-plane material properties. One shows in section II.1.6.5 that on the other hand, the calculation of out-of-plane shear strains caused by out-of-plane shear forces requires the knowledge of ply out-of-plane material constants.


II.1.6.5 Out-of-plane laminate shear stiffness

One assumes a linear relation between out-of-plane shear components of strain tensor and the corresponding components of Cauchy stress tensor:

$\displaystyle \left\{\begin{array}{c}
\gamma_{xz}(z) \\
\gamma_{yz}(z)
\end...
...ay}{c}
\tau_{xz}(z) \\
\tau_{yz}(z)
\end{array}\right\}_{\text{lam}}
\ .
$

To this relation corresponds a relation between the average out-of-plane shear strains and the out-of-plane shear force:

$\displaystyle \left\{\begin{array}{c}
\Gamma_{xz} \\
\Gamma_{yz}
\end{array...
...\{\begin{array}{c}
Q_{xz} \\
Q_{yz}
\end{array}\right\}_{\text{lam}}
\ .
$

II.1.6.5.1 Introduction of new notations

In the definition of loadings, the out-of-plane components of shear force $ Q$ can be replaced by average out-of-plane shear stress $ \overline{\tau}$ . Then the conversion between these two types of components is done simply by multiplication or division by laminate total thickness $ t$ :

$\displaystyle \overline{\tau}=\cfrac{Q}{t}\ .
$

One introduces notations that simplifies the writing of equations:

$\displaystyle \left\{\Gamma_s\right\}_{\text{*}}
=\left\{\begin{array}{c}
\Gamma_{xz} \\
\Gamma_{yz}
\end{array}\right\}_{\text{*}}
\ ,
$

$\displaystyle \left\{Q_s\right\}_{\text{*}}
=\left\{\begin{array}{c}
Q_{xz} \\
Q_{yz}
\end{array}\right\}_{\text{*}}
\ ,
$

$\displaystyle \left\{\overline{\tau_s}\right\}_{\text{*}}
=\left\{\begin{array...
...ine{\tau}_{xz} \\
\overline{\tau}_{yz}
\end{array}\right\}_{\text{*}}
\ .
$

In these expressions the $ *$ subscripts can be replaced by a symbol specific to the coordinate system in which the components of the vector are expressed (for example "load", "ply", "lam"...).

II.1.6.5.2 Energetic approach

The components of matrix $ \left[g\right]_{\text{lam}}(z)$ are easily obtained from the orientation and material of plies. The components of $ \left[G\right]_{\text{lam}}$ are obtained by a calculation of out-of-plane shear strain surface energy. One first calculates an estimate of this surfacic energy using the local expression of shear strains:

$\displaystyle W_{\text{shear}}^{\text{local}}$ $\displaystyle =\frac{1}{2}\int_{z_0}^{z_N} {\left\{\gamma_s(z)\right\}}^{\rm T}...
...xt{lam}} \cdot \left\{\tau_s(z)\right\}_{\text{lam}} \hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\frac{1}{2}\int_{z_0}^{z_N} {\left\{\tau_s(z)\right\}}^{\rm T}_{...
...lam}}(z) \cdot \left\{\tau_s(z)\right\}_{\text{lam}} \hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\frac{1}{2}\int_{z_0}^{z_N} {\left\{Q_s\right\}}^{\rm T}_{\text{...
...{\text{lam}}(z) \cdot \left\{Q_s\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\frac{1}{2} {\left\{Q_s\right\}}^{\rm T}_{\text{lam}} \cdot \lef...
...ce{0.4em}\text{d}z\right]_{\text{lam}} \cdot \left\{Q_s\right\}_{\text{lam}}\ .$ (II.159)

similarly, the surfacic energy can be estimated from the out-of-plane shear global equation:

$\displaystyle W_{\text{shear}}^{\text{global}}$ $\displaystyle =\frac{1}{2} {\left\{\Gamma_s(z)\right\}}^{\rm T}_{\text{lam}} \cdot \left\{Q_s\right\}_{\text{lam}} \ ,$    
  $\displaystyle =\frac{1}{2} {\left\{Q_s\right\}}^{\rm T}_{\text{lam}} \cdot \left[G\right]^{-1}_{\text{lam}} \cdot \left\{Q_s\right\}_{\text{lam}} \ .$ (II.160)

As there is only one surfacic energy, $ W_{\text{shear}}^{\text{global}}
=W_{\text{shear}}^{\text{local}}$ and

$\displaystyle \left[G\right]^{-1}_{\text{lam}}$ $\displaystyle =\int_{z_0}^{z_N} {\left[Y\right]}^{\rm T}_{\text{lam}}(z) \cdot ...
...{\text{lam}}(z) \cdot \left[Y\right]_{\text{lam}}(z) \hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\int_{z_0}^{z_N} \left[U\right]_{\text{lam}}(z) \hspace{0.4em}\text{d}z\ .$    

Here again, the integration can be calculated ply-by-ply. More precisely, one calculates on ply $ k$ :

$\displaystyle \left[U^k\right]_{\text{lam}}(z)$ $\displaystyle ={\left[Y^k\right]}^{\rm T}_{\text{lam}}(z) \cdot \left[g^k\right]^{-1}_{\text{lam}} \cdot \left[Y^k\right]_{\text{lam}}(z)\ ,$    
  $\displaystyle =\left[U_0^k\right]_{\text{lam}} +\left[U_1^k\right]_{\text{lam}}...
...}z^2 +\left[U_3^k\right]_{\text{lam}}z^3 +\left[U_4^k\right]_{\text{lam}}z^4\ .$    

where

$\displaystyle \left[U_0^k\right]_{\text{lam}}$ $\displaystyle ={\left[Y_0^k\right]}^{\rm T}_{\text{lam}} \cdot \left[g^k\right]^{-1}_{\text{lam}} \cdot \left[Y_0^k\right]_{\text{lam}}\ ,$    
$\displaystyle \left[U_1^k\right]_{\text{lam}}$ $\displaystyle ={\left[Y_0^k\right]}^{\rm T}_{\text{lam}} \cdot \left[g^k\right]...
...dot \left[g^k\right]^{-1}_{\text{lam}} \cdot \left[Y_0^k\right]_{\text{lam}}\ ,$    
$\displaystyle \left[U_2^k\right]_{\text{lam}}$ $\displaystyle ={\left[Y_0^k\right]}^{\rm T}_{\text{lam}} \cdot \left[g^k\right]...
...dot \left[g^k\right]^{-1}_{\text{lam}} \cdot \left[Y_0^k\right]_{\text{lam}}\ ,$    
$\displaystyle \left[U_3^k\right]_{\text{lam}}$ $\displaystyle ={\left[Y_1^k\right]}^{\rm T}_{\text{lam}} \cdot \left[g^k\right]...
...dot \left[g^k\right]^{-1}_{\text{lam}} \cdot \left[Y_1^k\right]_{\text{lam}}\ ,$    
$\displaystyle \left[U_4^k\right]_{\text{lam}}$ $\displaystyle ={\left[Y_2^k\right]}^{\rm T}_{\text{lam}} \cdot \left[g^k\right]^{-1}_{\text{lam}} \cdot \left[Y_2^k\right]_{\text{lam}}\ .$    

Then the integral above develops as follows:

$\displaystyle \left[G\right]^{-1}_{\text{lam}}$ $\displaystyle =\int_{z_0}^{z_N} \left[U\right]_{\text{lam}}(z) \hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\int_{z_{k-1}}^{z_k} \left[U^k\right]_{\text{lam}}(z) \hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\int_{z_{k-1}}^{z_k}\left( \sum_{i=0}^4\left[U_i^k\right]_{\text{lam}}z^i \right)\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left(\sum_{i=0}^4\left[U_i^k\right]_{\text{lam}} \int_{z_{k-1}}^{z_k}z^i\hspace{0.4em}\text{d}z\right)\ ,$    
  $\displaystyle =\sum_{k=1}^N\sum_{i=0}^4\left(\left[U_i^k\right]_{\text{lam}} \frac{z_k^{i+1}-z_{k-1}^{i+1}}{i+1}\right)\ .$ (II.161)

One notes the stiffness matrix $ \left[G\right]_{\text{lam}}$ and the compliance matrix $ \left[G^{-1}\right]_{\text{lam}}$ . Note that once the laminate out-of-plane shear stiffness and compliance matrices are known, the laminate out-of-plane shear equivalent moduli are calculated from the components of the compliance matrix with the following expressions:

$\displaystyle G_{xz}=\cfrac{1}{t\cdot\left[G^{-1}\right]_{11}}\ ,
$

$\displaystyle G_{yz}=\cfrac{1}{t\cdot\left[G^{-1}\right]_{22}}\ ,
$

in which the $ \left[G^{-1}\right]_{\text{lam}}$ matrix has first been rotated into the appropriate axes.

II.1.6.6 Calculation algorithm for shear stiffness

One describes below the calculation sequences that is used to calculate the laminate out-of-plane shear stiffness properties, and the out-of-plane shear stresses related to a given loading of the laminate.

The calculation sequence is described below. It involves two loops on the laminate layers.

  1. Calculate laminate in-plane and flexural properties. This is necessary because one needs the matrices $ \left[b\right]_{\text{lam}}$ and $ \left[d\right]_{\text{lam}}$ to calculate out-of-plane shear properties.
  2. One initializes the $ 2\times 2$ matrix $ \left[G\right]^{-1}_{\text{lam}}$ to zero.
  3. Then for each layer $ k$ with $ 1\leq k\leq N$ , one performs the following sequence of operations:
    1. One estimates the $ 3\times 3$ matrix of in-plane stiffness coefficients in laminate axes $ \left[C^k\right]_{\text{lam}}$ . For other calculations, one also need properties like the laminate thickness $ t$ and the positions $ z_k$ of different layer interfaces.
    2. This matrix is used to calculate the two matrices $ \left[F_0^k\right]_{\text{lam}}$ and $ \left[F_1^k\right]_{\text{lam}}$ . (See section II.1.6.3 for more details.) One has:

      $\displaystyle \left[F^k\right]_{\text{lam}}(z)
=\left[F_0^k\right]_{\text{lam}}
+z\left[F_1^k\right]_{\text{lam}}\ .
$

    3. Then, one calculates two other $ 2\times 2$ matrices $ \left[X_0^k\right]_{\text{lam}}$ and $ \left[X_1^k\right]_{\text{lam}}$ . (See section II.1.6.3.) One has:

      $\displaystyle \left[X^k\right]_{\text{lam}}(z)
=\left[X_0^k\right]_{\text{lam}}
+z\left[X_1^k\right]_{\text{lam}}\ .
$

    4. Then one calculates the $ \left[Y_i^k\right]_{\text{lam}}$ $ 2\times 2$ matrices:

      $\displaystyle \left[Y_0^k\right]_{\text{lam}}
=\left[Y_0^{k-1}\right]_{\text{l...
..._{\text{lam}}
-\cfrac{\left[X_1^k\right]_{\text{lam}}}{2}\right)z_{k-1}^2\ ,
$

      $\displaystyle \left[Y_1^k\right]_{\text{lam}}
=\left[X_0^k\right]_{\text{lam}}\ ,
$

      $\displaystyle \left[Y_2^k\right]_{\text{lam}}
=\cfrac{\left[X_1^k\right]_{\text{lam}}}{2}\ .
$

      As the expression of $ \left[Y_0^k\right]_{\text{lam}}$ is recursive, one needs another expression for the first value. The expression is:

      $\displaystyle \left[Y_0^1\right]_{\text{lam}}$ $\displaystyle =-\left[Y_1^1\right]_{\text{lam}}z_0 -\left[Y_2^1\right]_{\text{lam}}z_0^2\ .$    
        $\displaystyle =-\left[X_0^1\right]_{\text{lam}}z_0 -\cfrac{1}{2}\left[X_1^1\right]_{\text{lam}}z_0^2\ .$    

      These matrices allow to calculate the out-of-plane shear stress from the global out-of-plane shear force:

      $\displaystyle \left[Y^k\right]_{\text{lam}}(z)
=\left[Y_0^k\right]_{\text{lam}}
+\left[Y_1^k\right]_{\text{lam}}z
+\left[Y_2^k\right]_{\text{lam}}z^2\ ,
$

      $\displaystyle \left\{\tau_s(z)\right\}_{\text{lam}}
=\left[Y\right]_{\text{lam}}(z)
\cdot
\left\{Q_s\right\}_{\text{lam}}
\ .
$

      Actually, one is interested in stresses in ply axes, rather than in laminate axes. If $ \xi$ is the orientation of the ply in laminate axes, then:

      $\displaystyle \left\{\tau_s(z)\right\}_{\text{ply}}
=\left[S_+\right]_{(\xi)}
\cdot
\left\{\tau_s(z)\right\}_{\text{lam}}
\ .
$

      Then, one has simply:

      $\displaystyle \left\{\tau_s(z)\right\}_{\text{ply}}$ $\displaystyle =\left[S_+\right]_{(\xi)} \cdot \left[Y\right]_{\text{lam}}(z) \cdot \left\{Q_s\right\}_{\text{lam}} \ ,$    
        $\displaystyle = \left[Y\right]_{\text{ply}}^{\text{lam}}(z) \cdot \left\{Q_s\right\}_{\text{lam}} \ .$    

      One stores a $ 2\times 2$ matrix $ \left[Y\right]_{\text{ply}}^{\text{lam}}$ for each station through laminate thickness where out-of-plane shear stress might be requested. Actually it is done at top, mid and bottom surfaces in each ply. This means that $ 3N$ $ 2\times 2$ matrices $ \left[Y^j\right]_{\text{ply}}^{\text{lam}}$ are stored in the ClaLam object.

    5. One calculates the $ \left[U_i^k\right]_{\text{lam}}$ $ 2\times 2$ matrix. (See the end of section II.1.6.5 for the expressions to be used.) Then to $ \left[G\right]^{-1}_{\text{lam}}$ , one adds one term:

      $\displaystyle \left[G\right]^{-1}_{\text{lam}}=\left[G\right]^{-1}_{\text{lam}}...
...left[U_i^k\right]_{\text{lam}}
\frac{z_k^{i+1}-z_{k-1}^{i+1}}{i+1}\right)\ .
$

  4. At the end of the loop on layers, the shear stiffness matrix $ \left[G\right]_{\text{lam}}$ is calculated by inversion of $ \left[G\right]^{-1}_{\text{lam}}$ .

One also defines an out-of-plane shear compliance matrix calculated as follows:

$\displaystyle \left[g\right]_{\text{lam}}=\left[G\right]^{-1}_{\text{lam}}.
$

This matrix allows to calculate the laminate out-of-plane shear moduli:

$\displaystyle G_{xz}=\cfrac{1}{t\cdot g_{\text{11}}},
$

$\displaystyle G_{yz}=\cfrac{1}{t\cdot g_{\text{22}}},
$

Note that the values calculated above do not correspond to an out-of-plane shear stiffness of a material equivalent to the defined laminate. To convince yourself of this you can define a laminate with a single ply of orthotropic material. Then, you will observe that $ G_{xz}=G_{yz}=5G/6$ in which $ G$ is material shear modulus. (The usual $ 5/6$ factor in shell theory is recovered.)


II.1.7 CTE and CME calculations

One assumes a linear dependence of the temperature on the location through the laminate thickness:

$\displaystyle T(z)=T_0+zT_{,z}\ .$ (II.162)

Similarly, the water uptake depends linearly on $ z$ :

$\displaystyle H(z)=H_0+zH_{,z}\ .$ (II.163)

The calculation of laminate response to hygrometric loading is very similar to its response to thermo-elastic loading. Therefore, the following developments are done for thermo-elastic loading only. Later, they are transposed to hygrothermal solicitations.


II.1.7.1 In-plane and flexural thermo-elastic behavior

One calculates the stresses induced in plies for a thermo-elastic loading assuming that the material strain components are all constrained to zero. Equation (II.1.34) becomes:

$\displaystyle \left\{\sigma\right\}_{\text{ply}}
=
-\Delta T
\left[C\right]_{\text{ply}}
\cdot
\left\{\alpha\right\}_{\text{ply}}
\ .
$

In laminate axes, the equation is rewritten:

$\displaystyle \left\{\sigma\right\}_{\text{lam}}
=
-\Delta T
\left[C\right]_{\text{lam}}
\cdot
\left\{\alpha\right\}_{\text{lam}}
\ .
$

One substitutes in the equation the assumed temperature profile:

$\displaystyle \left\{\sigma\right\}_{\text{lam}}
=
-[(T_0-T_{\text{ref}})+zT...
...
\left[C\right]_{\text{lam}}
\cdot
\left\{\alpha\right\}_{\text{lam}}
\ .
$

The corresponding laminate in-plane force tensor is obtained by integrating the Cauchy stress tensor along the thickness:

$\displaystyle \left\{N\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}\left\{\sigma\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =-\int_{-h/2}^{h/2}[(T_0-T_{\text{ref}})+zT_{,z}] \left[C\right]_{\text{lam}}\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =-\left(\int_{-h/2}^{h/2}\left[C\right]_{\text{lam}}\cdot \left\{...
...\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\right)T_{,z}\ ,$    
  $\displaystyle =-\left\{\alpha Eh\right\}_{\text{lam}}(T_0-T_{\text{ref}}) -\left\{\alpha Eh^2\right\}_{\text{lam}}T_{,z}\ .$    

In the previous expression, two new symbols have been introduced that are calculated as follows:

$\displaystyle \left\{\alpha Eh\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}\left[C\right]_{\text{lam}}\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[C\right]_{\text{lam}}^k \cdot\left\{\alpha\right\}_{\text{lam}}^k \left(z_k-z_{k-1}\right)\ .$ (II.164)

$\displaystyle \left\{\alpha Eh^2\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}z\left[C\right]_{\text{lam}}\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[C\right]_{\text{lam}}^k \cdot\left\{\alpha\right\}_{\text{lam}}^k \left(\frac{z_k^2-z_{k-1}^2}{2}\right)\ .$ (II.165)

Similarly the bending moment tensor is obtained by integrating the Cauchy stress tensor multiplied by $ z$ along the thickness:

$\displaystyle \left\{M\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}z\left\{\sigma\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =-\int_{-h/2}^{h/2}z[(T_0-T_{\text{ref}})+zT_{,z}] \left[C\right]_{\text{lam}}\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =-\left(\int_{-h/2}^{h/2}z\left[C\right]_{\text{lam}}\cdot \left\...
...\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\right)T_{,z}\ ,$    
  $\displaystyle =-\left\{\alpha Eh^2\right\}_{\text{lam}}(T_0-T_{\text{ref}}) -\left\{\alpha Eh^3\right\}_{\text{lam}}T_{,z}\ .$    

In the previous expression, one new symbol has been introduced:

$\displaystyle \left\{\alpha Eh^3\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}z^2\left[C\right]_{\text{lam}}\cdot \left\{\alpha\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[C\right]_{\text{lam}}^k \cdot\left\{\alpha\right\}_{\text{lam}}^k \left(\frac{z_k^3-z_{k-1}^3}{3}\right)\ .$ (II.166)

Because of the linearity of all the equations, the thermo-elastic loading may be considered as an additional loading applied to the laminate, and if one considers an additional imposition of average in-plane strain and of a curvature, the laminate in-plane forces and bending moments are given by:


Using relation (II.1.44), the previous expression is reversed as follows:

\begin{multline}
\left\{\begin{array}{c}
\left\{\epsilon^0\right\}_{\text{lam}...
...^3\right\}_{\text{lam}}
\end{array}\right\}T_{,z}
\nonumber\ .
\end{multline}

In the last expression, four new quantities can be identified:

$\displaystyle \left\{\alpha_0^\epsilon\right\}_{\text{lam}} =\left[a\right]_{\t...
...xt{lam}} +\left[b\right]_{\text{lam}}\left\{\alpha Eh^2\right\}_{\text{lam}}\ ,$ (II.167)

$\displaystyle \left\{\alpha_1^\epsilon\right\}_{\text{lam}} =\left[a\right]_{\t...
...xt{lam}} +\left[b\right]_{\text{lam}}\left\{\alpha Eh^3\right\}_{\text{lam}}\ ,$ (II.168)

$\displaystyle \left\{\alpha_0^\kappa\right\}_{\text{lam}} ={\left[b\right]}^{\r...
...xt{lam}} +\left[d\right]_{\text{lam}}\left\{\alpha Eh^2\right\}_{\text{lam}}\ ,$ (II.169)

$\displaystyle \left\{\alpha_1^\kappa\right\}_{\text{lam}} ={\left[b\right]}^{\r...
...xt{lam}} +\left[d\right]_{\text{lam}}\left\{\alpha Eh^3\right\}_{\text{lam}}\ .$ (II.170)

So that finally, the ``compliance'' equation is:

\begin{multline}
\left\{\begin{array}{c}
\left\{\epsilon^0\right\}_{\text{lam}...
..._1^\kappa\right\}_{\text{lam}}
\end{array}\right\}T_{,z}
\ .
\end{multline}


II.1.7.2 Out-of-plane shear thermo-elastic behavior

Starting with the out-of-plane shear constitutive equation (II.1.42) and of the expression defining the out-of-plane shear force vectors $ \left\{Q_s\right\}_{\text{lam}}$ one makes developments similar to those of section II.1.7.1 and defines the following quantities:

$\displaystyle \left\{\alpha_sGh\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}\left[g\right]_{\text{lam}}\cdot \left\{\alpha_s\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[g\right]_{\text{lam}}^k \cdot\left\{\alpha_s\right\}_{\text{lam}}^k \left(z_k-z_{k-1}\right)\ .$ (II.171)

$\displaystyle \left\{\alpha_sGh^2\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}z\left[g\right]_{\text{lam}}\cdot \left\{\alpha_s\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[g\right]_{\text{lam}}^k \cdot\left\{\alpha_s\right\}_{\text{lam}}^k \left(\frac{z_k^2-z_{k-1}^2}{2}\right)\ .$ (II.172)

They are used in the expression:

$\displaystyle \left\{Q\right\}_{\text{lam}} = \left[G\right]_{\text{lam}} \cdot...
...t{lam}}(T_0-T_{\text{ref}}) -\left\{\alpha_sGh^2\right\}_{\text{lam}}T_{,z} \ .$ (II.173)

Correspondingly, one estimates the laminate out-of-plane shear CTE vectors:

$\displaystyle \left\{\alpha_0^s\right\}_{\text{lam}} =\left[G\right]_{\text{lam}}^{-1}\cdot\left\{\alpha Gh\right\}_{\text{lam}}$ (II.174)

$\displaystyle \left\{\alpha_1^s\right\}_{\text{lam}} =\left[G\right]_{\text{lam}}^{-1}\cdot\left\{\alpha Gh^2\right\}_{\text{lam}}$ (II.175)

These two expressions allow to write the expression of the inverse of (II.1.78):

$\displaystyle \left[\Gamma\right]_{\text{lam}}=\left[G\right]_{\text{lam}}^{-1}...
...text{lam}}(T_0-T_{\text{ref}}) +\left\{\alpha_1^s\right\}_{\text{lam}}T_{,z}\ .$ (II.176)


II.1.7.3 Hygrometric behavior of laminates

One transposes below the results for thermo-elastic behavior. It is done simply by replacing the CTE $ \alpha$ by the CME $ \beta$ in the definitions and equations.

$\displaystyle \left\{\beta Eh\right\}_{\text{lam}}$ $\displaystyle =\sum_{k=1}^N\left[C\right]_{\text{lam}}^k \cdot\left\{\beta\right\}_{\text{lam}}^k \left(z_k-z_{k-1}\right)\ ,$ (II.177)

$\displaystyle \left\{\beta Eh^2\right\}_{\text{lam}}$ $\displaystyle =\sum_{k=1}^N\left[C\right]_{\text{lam}}^k \cdot\left\{\beta\right\}_{\text{lam}}^k \left(\frac{z_k^2-z_{k-1}^2}{2}\right)\ ,$ (II.178)

$\displaystyle \left\{\beta Eh^3\right\}_{\text{lam}}$ $\displaystyle =\sum_{k=1}^N\left[C\right]_{\text{lam}}^k \cdot\left\{\beta\right\}_{\text{lam}}^k \left(\frac{z_k^3-z_{k-1}^3}{3}\right)\ ,$ (II.179)


$\displaystyle \left\{\beta_0^\epsilon\right\}_{\text{lam}} =\left[a\right]_{\te...
...ext{lam}} +\left[b\right]_{\text{lam}}\left\{\beta Eh^2\right\}_{\text{lam}}\ ,$ (II.180)

$\displaystyle \left\{\beta_1^\epsilon\right\}_{\text{lam}} =\left[a\right]_{\te...
...ext{lam}} +\left[b\right]_{\text{lam}}\left\{\beta Eh^3\right\}_{\text{lam}}\ ,$ (II.181)

$\displaystyle \left\{\beta_0^\kappa\right\}_{\text{lam}} ={\left[b\right]}^{\rm...
...ext{lam}} +\left[d\right]_{\text{lam}}\left\{\beta Eh^2\right\}_{\text{lam}}\ ,$ (II.182)

$\displaystyle \left\{\beta_1^\kappa\right\}_{\text{lam}} ={\left[b\right]}^{\rm...
...ext{lam}} +\left[d\right]_{\text{lam}}\left\{\beta Eh^3\right\}_{\text{lam}}\ ,$ (II.183)

\begin{multline}
\left\{\begin{array}{c}
\left\{\epsilon^0\right\}_{\text{lam}...
..._1^\kappa\right\}_{\text{lam}}
\end{array}\right\}H_{,z}
\ ,
\end{multline}

$\displaystyle \left\{\beta_sGh\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}\left[g\right]_{\text{lam}}\cdot \left\{\beta_s\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[g\right]_{\text{lam}}^k \cdot\left\{\beta_s\right\}_{\text{lam}}^k \left(z_k-z_{k-1}\right)\ ,$ (II.184)

$\displaystyle \left\{\beta_sGh^2\right\}_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{h/2}z\left[g\right]_{\text{lam}}\cdot \left\{\beta_s\right\}_{\text{lam}}\hspace{0.4em}\text{d}z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left[g\right]_{\text{lam}}^k \cdot\left\{\beta_s\right\}_{\text{lam}}^k \left(\frac{z_k^2-z_{k-1}^2}{2}\right)\ ,$ (II.185)

$\displaystyle \left\{Q\right\}_{\text{lam}} = \left[G\right]_{\text{lam}} \cdot...
...xt{lam}}(H_0-H_{\text{ref}}) -\left\{\beta_sGh^2\right\}_{\text{lam}}H_{,z} \ ,$ (II.186)

$\displaystyle \left\{\beta_0^s\right\}_{\text{lam}} =\left[G\right]_{\text{lam}}^{-1}\cdot\left\{\beta Gh\right\}_{\text{lam}}$ (II.187)

$\displaystyle \left\{\beta_1^s\right\}_{\text{lam}} =\left[G\right]_{\text{lam}}^{-1}\cdot\left\{\beta Gh^2\right\}_{\text{lam}}$ (II.188)

$\displaystyle \left[\Gamma\right]_{\text{lam}}=\left[G\right]_{\text{lam}}^{-1}...
...\text{lam}}(H_0-H_{\text{ref}}) +\left\{\beta_1^s\right\}_{\text{lam}}H_{,z}\ .$ (II.189)


II.1.7.4 Full sets of equations

Finally, the full set of constitutive equation written with stiffness matrices looks like:



The two previous expressions are inversed as follows:




II.1.8 Calculation of load response

One always considered a decoupling of in-plane and flexural of laminates on one side, and the out-of-plane shear of laminates on the other side. These two aspects are discussed in sections II.1.8.1 and II.1.8.2 respectively.


II.1.8.1 In-plane and flexural response

Beside thermo-elastic or hygro-elastic loading, the composite classes of FeResPost allows the definition different types of mechanical loads:

The type of loading is specified component-by-component. This means that a single loading may have some components imposed as normal forces and bending moments, with other components imposed as average strains, and other components as average stresses or flexural stresses. The mechanical part of loading is also characterized by a direction $ \lambda$ wrt laminate axes. The subscript ``load'' indicates that the components are given in loading axes. One explains below how the laminate response is calculated.

  1. The solver first checks if average or flexural stresses are imposed. If such components of the loading are found, they are converted to in-plane forces and bending moments with the following equations:

    $\displaystyle \left\{N\right\}_{\text{load}}
=h\left\{\sigma^0\right\}_{\text{load}}\ ,
$

    $\displaystyle \left\{M\right\}_{\text{load}}
=\frac{h^2}{6}\left\{\sigma^f\right\}_{\text{load}}\ ,
$

    in which $ h$ is the laminate thickness.

  2. The mechanical part of loading is characterized by a direction wrt laminate axes. This direction is given by an angle $ \lambda$ . In order to have laminate properties and loading given in the same coordinate system, the laminate stiffness matrices and CTE vectors are calculated in this new coordinate system. (It is more convenient for the elimination of components imposed as average strains or curvatures.) More precisely, the stiffness matrices and CTE vector are rotated with the following expressions:

    $\displaystyle \left[A\right]_{\text{load}}=
\left[T_-\right]_{(\lambda)}\cdot\left[A\right]_{\text{lam}}
\cdot\left[T_-^\prime\right]_{(\lambda)}\ ,
$

    $\displaystyle \left[B\right]_{\text{load}}=
\left[T_-\right]_{(\lambda)}\cdot\left[B\right]_{\text{lam}}
\cdot\left[T_-^\prime\right]_{(\lambda)}\ ,
$

    $\displaystyle \left[D\right]_{\text{load}}=
\left[T_-\right]_{(\lambda)}\cdot\left[D\right]_{\text{lam}}
\cdot\left[T_-^\prime\right]_{(\lambda)}\ ,
$

    $\displaystyle \left\{\alpha Eh\right\}_{\text{load}}
=
\left[T_-\right]_{(\lambda)}
\cdot
\left\{\alpha Eh\right\}_{\text{lam}}\ ,
$

    $\displaystyle \ldots
$

    $\displaystyle \left\{\beta Eh^3\right\}_{\text{load}}
=
\left[T_-\right]_{(\lambda)}
\cdot
\left\{\beta Eh^3\right\}_{\text{lam}}\ .
$

    The calculation of CTE and CME related quantities is done only if the corresponding temperature or moisture contributions have been defined in the loading. The system of equations looks like:

    \begin{multline}
\left[\begin{array}{cc}
\left[A\right]_{\text{load}} & \left[...
...h^3\right\}_{\text{load}}
\end{array}\right\}H_{,z}\ .\nonumber
\end{multline}

    (Here again the CTE and CME related terms are optional.) Actually, one can write a single set of 6 equations with 6 unknowns. The general form of this system is

    $\displaystyle \sum_{j=1}^n K_{ij}u_j=b_i,\hspace{5mm}i=1...n.
$

  3. Now, one considers a case in which one component of vector $ u_j$ is constrained to be a certain value. For example $ u_k=a$ . This equation replaces the $ k^{\text{th}}$ equation of the system:

    $\displaystyle \sum_{j=1}^n K_{ij}u_j$ $\displaystyle = b_i, \hspace{5mm}i=1...n,\ i\neq k$    
    $\displaystyle u_k$ $\displaystyle =a\ .$    

    The unknown $ u_k$ can be easily eliminated from the linear

    $\displaystyle \sum_{j=1,j\neq k}^n K_{ij}u_j$ $\displaystyle = b_i-K_{ik}a, \hspace{5mm}i=1...n,\ i\neq k$    
    $\displaystyle u_k$ $\displaystyle =a\ .$    

    The first line above corresponds to a new linear system of $ n-1$ equations with $ n-1$ unknowns. The set of two lines define the algebraic operations that are performed in FeResPost when one imposes an average strain or curvature component.

    Actually, the operation can be simplified. It is sufficient to replace line $ k$ in the linear system of equations by the constraint equation $ u_k=a$ and perform the ``usual'' Gaussian elimination to solve the linear system of equations.

  4. When all the components of loading imposed as average strains or curvature have been eliminated from the linear system, a classical Gaussian elimination algorithm calculates the other unknowns of the system.

    Then the components of tensors $ \left\{\epsilon^0\right\}$ and $ \left\{\kappa\right\}$ are known in loading axes.

  5. The normal forces and bending moments are then calculated in loading axes with the following equations:

    \begin{multline}
\left\{N\right\}_{\text{load}}=
\left[A\right]_{\text{load}}\...
...})
-\left\{\beta Eh^2\right\}_{\text{load}}H_{,z}
\nonumber\ ,
\end{multline}

    \begin{multline}
\left\{M\right\}_{\text{load}}=
\left[B\right]_{\text{load}}\...
...})
-\left\{\beta Eh^3\right\}_{\text{load}}H_{,z}
\nonumber\ .
\end{multline}

    (The CTE and CME related terms are optional.)

  6. If $ \lambda$ is the angle characterizing the loading orientation wrt laminate axes a rotation of $ -\lambda$ of the two vectors gives the average strain and curvature tensors in laminate axes: $ \left\{\epsilon^0\right\}_{\text{lam}}$ and $ \left\{\kappa\right\}_{\text{lam}}$ .

    $\displaystyle \left\{\epsilon^0\right\}_{\text{lam}}=
\left[T_+^\prime\right]_{(-\lambda)}
\cdot\left\{\epsilon^0\right\}_{\text{load}}\ ,
$

    $\displaystyle \left\{\kappa\right\}_{\text{lam}}=
\left[T_+^\prime\right]_{(-\lambda)}
\cdot\left\{\kappa\right\}_{\text{load}}\ .
$

    Similarly, the normal forces and bending moments components are re-expressed in laminate axes:

    $\displaystyle \left\{N\right\}_{\text{lam}}=
\left[T_-\right]_{(-\lambda)}
\cdot\left\{N\right\}_{\text{load}}\ ,
$

    $\displaystyle \left\{M\right\}_{\text{lam}}=
\left[T_-\right]_{(-\lambda)}
\cdot\left\{M\right\}_{\text{load}}\ .
$

  7. For each ply, one calculates (if required) the stresses and strains as follows:
    1. One rotates the laminate average strain and curvature tensors to obtain them in ply axes. If the ply is characterized by an angle $ \xi$ wrt laminate axes, the two tensors are rotated by the same angle $ \xi$ :

      $\displaystyle \left\{\epsilon^0\right\}_{\text{ply}}=
\left[T_+^\prime\right]_{(\xi)}
\cdot\left\{\epsilon^0\right\}_{\text{lam}}\ ,
$

      $\displaystyle \left\{\kappa\right\}_{\text{ply}}=
\left[T_+^\prime\right]_{(\xi)}
\cdot\left\{\kappa\right\}_{\text{lam}}\ .
$

      Note that, even though the components of these two tensors are now given in one of the plies coordinate system, they correspond to strain or curvature of the laminate at mid-thickness.
    2. At the different stations through the thickness at which strains and stresses are required, the strain components are calculated with:

      $\displaystyle \left\{\epsilon\right\}_{\text{ply}}(z)=
\left\{\epsilon^0\right\}_{\text{ply}}
+z\left\{\kappa\right\}_{\text{ply}}\ .
$

      In FeResPost $ z$ may have the values $ z_{\text{inf}}$ , $ z_{\text{mid}}$ , $ z_{\text{sup}}$ . Then the stress components are given by:


      (Here again the CTE and CME related terms are optional.) A peculiar version of the ply strain tensor that corresponds to ply stresses, but without thermo-elastic or moisture contribution is calculated as follows:

      $\displaystyle \left\{\epsilon^{\text{Mech}}\right\}_{\text{ply}}(z)
=\left[c\right]_{\text{ply}}
\cdot\left\{\sigma\right\}_{\text{ply}}(z)\ .
$

      This version of the strain tensor is called the ``Mechanical Strain Tensor'' or ``Equivalent Strain Tensor''. This is the version of the strain tensors that is used for the strain failure criteria. Note however that a ``Total Strain'' version of the criteria is proposed as well.
At the end of the calculations, the laminate object which has been used to perform those calculations stores a few results: The ply results may be used later to calculate failure indices or reserve factors.


II.1.8.2 Out-of-plane shear response

Some of the quantities calculated above, and stored in the ClaLam object are used to estimate laminate shear loading response.

The different steps of the calculation are described below:

  1. The first step of the calculation is to resolve the loading in out-of-plane shear forces in loading axes $ \left\{Q\right\}_{\text{load}}$ (or $ \left\{\Gamma\right\}_{\text{load}}$ ). For this, one proceeds as in section II.1.5, but with the following differences: the conversion of average out-of-plane shear strain to out-of-plane shear force components requires the knowledge of out-of-plane shear stiffness matrix in loading axes. This one is readily obtained by transforming the corresponding matrix in laminate axes:

    $\displaystyle \left[G\right]_{\text{load}}
=
\left[S_+\right]_{(\lambda)}
\cdot
\left[G\right]_{\text{lam}}
\cdot
\left[S_-\right]_{(\lambda)}\ .
$

  2. The out-of-plane shear loading can be expressed by specifying the out-of-plane shear forces, or the out-of-plane average shear strain, or a combination of the two. In all cases, the components are specified in loading axes.

    If out-of-plane average shear forces are specified, the resolution of the following linear system of equations allows to calculate the corresponding out-of-plane shear strains:

    \begin{multline}
\left[G\right]_{\text{load}}
\cdot
\left\{\Gamma\right\}_{\t...
...)
+\left\{\beta_sGh^2\right\}_{\text{load}}H_{,z}
\ .\nonumber
\end{multline}

    (The CTE and CME related terms are optional.) The resolution of this equation is done following the same approach as for the in-plane and bending loading. One performs a Gaussian elimination in a $ 2\times 2$ matrix. Constraints can be imposed if out-of-plane shear strains are specified for some components of the loading instead of out-of-plane lineic shear force.

  3. At this stage, whatever the type of loading applied to the laminate, $ \left\{\Gamma\right\}_{\text{load}}$ is known. One can obtain the lineic out-of-plane shear forces with

    \begin{multline}
\left[Q\right]_{\text{load}}
=\left[G\right]_{\text{load}}\cd...
...}})
-\left\{\beta_1^s\right\}_{\text{load}}H_{,z}
\ .\nonumber
\end{multline}

    (The CTE and CME related terms are optional.) Once $ \left\{Q\right\}_{\text{load}}$ and $ \left\{\Gamma\right\}_{\text{load}}$ are known, the corresponding loading in laminate axes is obtained with:

    $\displaystyle \left\{Q\right\}_{\text{lam}}
=
\left[S_-\right]_{(\lambda)}
\cdot
\left\{Q\right\}_{\text{load}}\ ,
$

    $\displaystyle \left\{\Gamma\right\}_{\text{lam}}
=
\left[S_-\right]_{(\lambda)}
\cdot
\left\{\Gamma\right\}_{\text{load}}\ .
$

  4. The ply out-of-plane shear stress components are calculated at the different requested locations by:

    $\displaystyle \left\{\tau^j_s(z)\right\}_{\text{ply}}
=\left[Y^j(z)\right]_{\text{ply}}^{\text{lam}}
\cdot
\left\{Q_s\right\}_{\text{lam}}
\ ,
$

    in which the matrix $ \left[Y^j\right]_{\text{ply}}^{\text{lam}}$ is relative to the station at which the stress is requested.

  5. Finally, for the stations where out-of-plane shear stresses have been calculated the out-of-plain shear strain is also calculated using the corresponding ply material coefficients:

    \begin{multline}
\left\{\gamma_s(z)\right\}_{\text{ply}}
=
\left[g\right]^{-...
...\}_{\text{ply}}(H_0-H_{\text{ref}}+zH_{,z})\right)\ .
\nonumber
\end{multline}

    (The CTE and CME related terms are optional. One takes benefits of the ``decoupling of out-of-plane shear'' assumption.)


II.1.8.3 Out-of-plane T/C deformation

The Classical Lamination Theory is based on the assumption that $ \sigma_{33}=0$ . Consequently, $ \epsilon_{33}$ and $ \epsilon^{\text{Mech}}_{33}$ are generally not zero. These strain tensor components can be estimated from (II.1.27):

$\displaystyle \epsilon^{\text{Mech}}_{33}=
c_{3311}\sigma_{11}+c_{3322}\sigma_...
...3}\sigma_{33}+
c_{3323}\sigma_{23}+c_{3331}\sigma_{31}+c_{3312}\sigma_{12}\ ,
$

$\displaystyle \epsilon_{33}=\epsilon^{\text{Mech}}_{33}+\alpha_{33}\Delta T
+\beta_{33}\Delta H\ .
$


II.1.9 Accelerating the calculation of load response

In section II.1.8, one explained the different steps to solve the laminate load response equation, and estimate ply stresses and strains. The number of operations involved in these calculations is very important, and when repetitive calculation of laminate load response is done, the computation time can increase unacceptably. This is the case, for example, when laminate load response analysis is performed on loads extracted from finite element model results.

However, many operations described in section II.1.8 will be the same for each different loading. This means that these operations could be done only once for the different laminate loads considered in the analysis. We investigate in this section the possible accelerations of laminate load response analysis.


II.1.9.1 Calculation of laminate loads and strains

We explain in this section how the laminate load response calculation can be accelerated. In particular, when the calculation of laminate load response is done repititively with similar loading, the benefit of simplifying the sequence of operations to estimate laminate loading becomes obvious.

One explains the calculation of laminate in-plane and flexural load response in section II.1.8.1. The sequence of operations results in the building of $ 6\times 6$ matrix $ K$ that depends on laminate definition and loading angle. For all the calculations done with a common laminate and loading angle, the operations can be simplified as follows:

  1. The matrix $ K$ is assembled.
  2. The components of loading that are specified as in-plane strain, or curvature lead to the imposition of constraints on matrix $ K$ . For example, if one imposes $ u_k=a$ , it is sufficient to replace the elements of line $ k$ in matrix $ K$ by 0, except $ K_{kk}$ =1. When all the constraints have been imposed, one obtains a new $ 6\times 6$ matrix that we call $ K_{\text{constr}}$ .
  3. Finally, this matrix is inversed, and on obtained the $ 6\times 6$ matrix $ K^{-1}_{\text{constr}}$ .
This $ K^{-1}_{\text{constr}}$ matrix is the same for all the loadings that we apply with on the same laminate, with the same loading angle $ \theta$ , and with the same constraints. When calculating load response, this matrix is used as follows:
  1. For each laminate load, one assembles the 6 components vector $ b_{\text{constr}}$ as explained in section II.1.8.1. The components usually correspond to shell forces or moments, but the ``constrained components'' (components specified as strains or curvature) are replaced by components of shell in-plane strains or curvatures.
  2. Shell in-plane strains and curvatures in loading axes are obtained by calculating the following matricial product:

    $\displaystyle \left\{\begin{array}{c}
\left\{\epsilon^0\right\}_{\text{load}}...
...
\end{array}\right\}
=
K^{-1}_{\text{constr}}\cdot b_{\text{constr}}\ .
$

  3. Then, the shell strains and curvature components can be expressed in laminate axes by performing the following matricial operation:

    \begin{multline}
\left\{\begin{array}{c}
\left\{\epsilon^0\right\}_{\text{la...
...
\cdot
K^{-1}_{\text{constr}}
\cdot
b_{\text{constr}}\ .
\end{multline}

Finally, the last expression leads to the definition of a new $ 6\times 6$ matrix

$\displaystyle K^{-1}_{\text{constr, lam}}
=
\left[\begin{array}{cc}
\left...
...\right]_{(-\lambda)}
\end{array}\right]
\cdot
K^{-1}_{\text{constr}}\ ,
$

that allows to write

$\displaystyle \left\{\begin{array}{c} \left\{\epsilon^0\right\}_{\text{lam}} \...
...}} \end{array}\right\} =K^{-1}_{\text{constr, lam}}\cdot b_{\text{constr}}\ .$ (II.190)

This matrix can be constructed once and for all for a given laminate, loading angle and loading characteristics.

A similar approach can be used for the simplification of the laminate out-of-plane shear response calculation:

  1. The construction of a $ 2\times 2$ $ K^{\prime-1}_{\text{constr}}$ matrix allows to write:

    $\displaystyle \left[\Gamma\right]_{\text{load}}=K^{\prime-1}_{\text{constr}}\cdot b^\prime_{\text{constr}}\ ,
$

    in which $ b_{\text{constr}}$ corresponds to the out-of-plane shear loads (forces or strains).
  2. Then

    $\displaystyle \left[\Gamma\right]_{\text{lam}}
=
\left[S_-\right]_{(\lambda)}
\cdot
\left[\Gamma\right]_{\text{load}}\ .
$

  3. And this leads to the definition of a new $ 2\times 2$ matrix that allows to write:

    $\displaystyle K^{\prime-1}_{\text{constr, lam}}
=
\left[S_-\right]_{(\lambda)}
\cdot
K^{\prime-1}_{\text{constr}}\ ,
$

    $\displaystyle \left[\Gamma\right]_{\text{lam}}=K^{\prime-1}_{\text{constr, lam}}\cdot b^\prime_{\text{constr}}\ .$ (II.191)

Finally, now that the laminate strains and curvatures have been estimated in laminate axes, the corresponding laminate forces and moments are estimated as follows:

\begin{multline}
\left\{\begin{array}{c}
\left\{N\right\}_{\text{lam}} \\
...
...eta Eh^3\right\}_{\text{lam}}
\end{array}\right\}H_{,z}\ ,
\end{multline}

\begin{multline}
\left[Q\right]_{\text{lam}}
=
\left[G\right]_{\text{lam}}...
...{ref}})
-\left\{\beta_1^s\right\}_{\text{lam}}H_{,z}
\ .
\end{multline}

One can substitute (II.1.103) and (II.1.104) in expressions (II.1.105) and (II.1.106). This leads to the following expressions:

\begin{multline}
\left\{\begin{array}{c}
\left\{N\right\}_{\text{lam}} \\
...
...y}\right]
\cdot
\left\{b^{FM}_{\text{constr}}\right\}\ ,
\end{multline}

\begin{multline}
\left[Q\right]_{\text{lam}}
=
\left[\begin{array}{c\vert ...
...rray}\right]
\cdot
\left\{b^Q_{\text{constr}}\right\}\ .
\end{multline}

From equations (II.1.103), (II.1.104), (II.1.107) and (II.1.108), one identifies four matrices and four vectors that allow the calculation of laminate stress/strain state in laminate axes: The four matrices depend on laminate definition, angle $ \theta$ of the loading wrt laminate axes, and on the set of components that are constrained to strain or curvature values. The four vectors also depend on the values of the particular loading that is examined. This means that the four matrices can be calculated once and for all the particular loadings. On the other-hand the four vectors must be re-estimated for each element defining the load.


II.1.9.2 Calculation of plies stresses and strains

The calculation of ply stresses and strains from the laminate loads is easily accelerated. Indeed, let us consider a laminate loading corresponding to:

These quantities correspond to 12 real values that characterize entirely the laminate loading. For a given ply, the stresses can be calculated at a given height ($ z$ value) from these 12 real values. All the calculations are linear.

If one defines a vector with 12 components as follows:

$\displaystyle \left\{LLd\right\}_{\text{lam}}=
\left\{\begin{array}{c}
\lef...
...f}} \\
T_{,z} \\
H_0-H_{\text{ref}} \\
H_{,z}
\end{array}\right\},
$

That contains all the laminate loading, there must be a matrix $ \left[K_\sigma\right]_{\text{ply}}(z_i)$ that allows to calculate ply stresses as follows:

$\displaystyle \left\{\sigma\right\}_{\text{ply}}(z_i)=\left[K_\sigma\right]_{\text{ply}}(z_i)
\cdot\left\{LLd\right\}_{\text{lam}}\ .
$

Matrix $ \left[K_\sigma\right]_{\text{ply}}(z_i)$ is a $ \times 12$ matrix that depends only on laminate definition. This means that this matrix can be calculated once and for all when laminate is created in the database.

Similarly one can also define matrices $ \left[K_\epsilon\right]_{\text{ply}}(z_i)$ and $ \left[K^{Mech}_\epsilon\right]_{\text{ply}}(z_i)$ for the calculations of $ \left\{\epsilon\right\}_{\text{ply}}(z_i)$ and $ \left\{\epsilon^{\text{Mech}}\right\}_{\text{ply}}(z_i)$ respectively. We explain here, how these three matrices can be constructed.

II.1.9.2.1 Loading in ply axes

The first step of ply stresses or strain calculations consists in expressing the laminate loading in ply axes. The following operations are performed:

$\displaystyle \left\{\epsilon^0\right\}_{\text{ply}}=
\left[T_+^\prime\right]_{(\xi)}
\cdot\left\{\epsilon^0\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\kappa\right\}_{\text{ply}}=
\left[T_+^\prime\right]_{(\xi)}
\cdot\left\{\kappa\right\}_{\text{lam}}\ ,
$

$\displaystyle \left\{\tau^j_s(z)\right\}_{\text{ply}}
=\left[Y^j(z)\right]_{\text{ply}}^{\text{lam}}
\cdot
\left\{Q_s\right\}_{\text{lam}}\ .
$

The four real values corresponding to laminate temperature and moisture loading are not affected by the modification of coordinate system. The three relations above allow us to write the following equation:

$\displaystyle \left\{\begin{array}{c} \left\{\epsilon^0\right\}_{\text{ply}} \...
...t{ref}} \\  T_{,z} \\  H_0-H_{\text{ref}} \\  H_{,z} \end{array}\right\}\ .$ (II.192)

Note that the out-of-plane shear response is now expressed as out-of-plane shear stresses at the specified height in selected ply and no longer as laminate out-of-plane shear forces.

The strains, temperatures and moisture at height $ z$ in selected ply is easily obtained with the following expressions:

$\displaystyle \left\{\epsilon\right\}_{\text{ply}}(z)=
\left\{\epsilon^0\right\}_{\text{ply}}
+z\left\{\kappa\right\}_{\text{ply}}\ ,
$

$\displaystyle T(z)-T_{\text{ref}}=T_0-T_{\text{ref}}+zT_{,z}\ ,
$

$\displaystyle H(z)-H_{\text{ref}}=H_0-H_{\text{ref}}+zH_{,z}\ .
$

The combination of these three expressions in a single matricial expression gives:

$\displaystyle \left\{\begin{array}{c} \left\{\epsilon\right\}_{\text{ply}}(z) ...
...t{ref}} \\  T_{,z} \\  H_0-H_{\text{ref}} \\  H_{,z} \end{array}\right\}\ .$ (II.193)

( $ \left[\delta_3\right]$ and $ \left[\delta_2\right]$ are $ 3\times 3$ and $ 2\times 2$ unit matrices respectively.) Then, considering equation (II.1.101), one writes:

$\displaystyle \left\{\begin{array}{c} \left\{\sigma\right\}_{\text{ply}}(z) \\...
...ply}} \\  T(z)-T_{\text{ref}} \\  H(z)-H_{\text{ref}} \end{array}\right\}\ .$ (II.194)

The vector at left hand side of expression (II.1.111) contains all the ply stress components. One can remove the two lower lines of the equation as follows:

$\displaystyle \left\{\begin{array}{c} \left\{\sigma\right\}_{\text{ply}}(z) \\...
...ply}} \\  T(z)-T_{\text{ref}} \\  H(z)-H_{\text{ref}} \end{array}\right\}\ .$ (II.195)

Then, the components can be reordered as follows:

$\displaystyle \left\{\sigma\right\}_{\text{mat}}(z) = \left(\begin{array}{c}...
...y}}(z) \\  \left\{\tau^j_s(z)\right\}_{\text{ply}} \\  \end{array}\right\}\ .$ (II.196)

One will also use:

$\displaystyle \left(\begin{array}{c} \left\{\sigma\right\}_{\text{mat}}(z) \\ ...
...ply}} \\  T(z)-T_{\text{ref}} \\  H(z)-H_{\text{ref}} \end{array}\right\}\ .$ (II.197)

One uses (II.1.27) to estimate the strain tensor:

$\displaystyle \left\{\epsilon\right\}_{\text{mat}}(z) = \left(\begin{array}{c...
... T(z)-T_{\text{ref}} \\  \hline H(z)-H_{\text{ref}}  \end{array}\right) \ .$ (II.198)

The vector that appears in right-hand-side of the previous expression The so-called ``mechanical strain tensor'' is given by:

$\displaystyle \left\{\epsilon\right\}^{\text{Mech}}_{\text{mat}}(z) = \left(\...
...  \left[c_{ijkl}\right] \cdot \left\{\sigma\right\}_{\text{mat}}(z) \\  \ .$ (II.199)

All the operations (II.1.109) to (II.1.116) reduce to matricial products. The characteristics of the matrices used in these operations are summarized in Table II.1.1. Two of these matrices are ``re-ordering'' metrices, and do not depend on the ply material or position accross laminate thickness. The other matrices depend on the ply material. Only two of the matrices depends on height $ z$ .


Table II.1.1: Summary of the matrices that have been defined for the acceleration of ply stresses and strains calculation.
Matrix Reference Size Dependence on
$ \left[K_1\right]$ (II.1.109) $ 12\times 12$ Ply angle, material and height $ z$
$ \left[K_2\right]$ (II.1.110) $ 7\times 12$ Height $ z$
$ \left[K_3\right]$ (II.1.111) $ 7\times 7$ Ply material
$ \left[K_4\right]$ (II.1.112) $ 5\times 7$ Ply material
$ \left[K_5\right]$ (II.1.113) $ 6\times 5$ --
$ \left[K_6\right]$ (II.1.114) $ 8\times 7$ --
$ \left[K_7\right]$ (II.1.115) $ 6\times 8$ Ply material
$ \left[K_8\right]$ (II.1.116) $ 6\times 6$ Ply material

In the end, one writes:

$\displaystyle \left\{\sigma\right\}_{\text{mat}}(z)$ $\displaystyle = \left[K_5\right] \cdot \left[K_4\right] \cdot \left[K_2\right](z) \cdot \left[K_1\right](z) \cdot \left\{LLd\right\}_{\text{lam}}\ ,$    
  $\displaystyle = \left[K_{\sigma}\right](z) \cdot \left\{LLd\right\}_{\text{lam}}\ .$ (II.1100)

$\displaystyle \left\{\epsilon\right\}_{\text{mat}}(z)$ $\displaystyle = \left[K_7\right] \cdot \left[K_6\right] \cdot \left[K_3\ri...
...ight](z) \cdot \left[K_1\right](z) \cdot \left\{LLd\right\}_{\text{lam}}\ ,$    
  $\displaystyle = \left[K_{\epsilon}\right](z) \cdot \left\{LLd\right\}_{\text{lam}}\ .$ (II.1101)

$\displaystyle \left\{\epsilon\right\}^{\text{Mech}}_{\text{mat}}(z)$ $\displaystyle = \left[K_8\right] \cdot \left\{\sigma\right\}_{\text{mat}}(z)\ ,$    
  $\displaystyle = \left[K_8\right] \cdot \left[K_{\sigma}\right](z) \cdot \left\{LLd\right\}_{\text{lam}}\ ,$    
  $\displaystyle = \left[K^{\text{Mech}}_{\epsilon}\right](z) \cdot \left\{LLd\right\}_{\text{lam}}\ .$ (II.1102)

We have shown that if one wishes to calculate stresses and strains in plies from the laminate loading $ \left\{LLd\right\}_{\text{lam}}$ , at the three ``bot'' ``mid'' and ``sup'' heights of laminate plies, one needs to calculate $ 3\times 3=9$ matrices per ply. Each of the matrices has $ 6\times 12$ size and depends only on the laminate definition. This means that these matrices can be calculated once and for all when laminate is defined.

Note that the acceleration matrices for the ply stresses and strains calculation must be re-estimated each time the laminate or one of its materials is modified. This is the reason why method ``reInitAllPliesAccelMatrices'' has been added to the ``ClaLam'' and ``ClaDb'' classes.


II.1.10 Failure theories

When stresses and strains have been calculated in plies (or some of the plies), the failure indices can be estimated too. One presents below the different failure theories that are proposed in FeResPost, and how these failure theories can be used to estimate laminate reserve factors.

In this section, one conventionally uses integer subscripts to denote that tensor components are given in ply axes and Roman subscripts to indicate principal components. Often, only in-plane components of ply stress or strain tensors are used to estimate criteria. Then, the principal components are estimated as follows:

$\displaystyle \sigma_{\text{I,II}} =\frac{1}{2}\left(\sigma_{11}+\sigma_{22}\r...
...)\pm \sqrt{\frac{1}{4}\left(\sigma_{11}+\sigma_{22}\right)^2 +\tau_{12}^2}\ .$ (II.1103)

In the rest of this section all the allowables used in failure criteria have positive values. Even the compressive allowables are $ \ge 0$ .

Table II.1.2 summarizes the criteria available in FeResPost. For each criterion, the Table provides:

  1. A String corresponding to the argument that identifies the selected criterion when derivation is asked.
  2. A description of the type of material (metallic or isotropic, unidirectional tape, fabric,...).
  3. A reference to the section in which the criterion is presented and discussed.
  4. Specification whether an equivalent stress for this criterion can be derived in FeResPost or not.
  5. One specifies whether a failure index can be derived with FeResPost. Generally, the failure index is calculated according to the ``usual'' definition in litterature. when no such standard failure index definition is available, one provides a default definition which corresponds to the inverse of the reserve factor calculated with $ FoS=1$ .
  6. One specifies whether a reserve factor and/or strength ratio can be calculated with FeResPost:


Table II.1.2: Summary of failure criteria available in FeResPost and information on what can be calculated.
Criterion Material Section Derived    
Name Type Number Stress F.I. R.F. and S. R.
``Tresca2D'' metallic II.1.10.1 yes yes yes
``VonMises2D'' metallic II.1.10.2 yes yes yes
``VonMises3D'' metallic II.1.10.3 yes yes yes
``MaxStress'' tape or fabric II.1.10.4 no yes yes
``MaxStress3D'' tape or fabric II.1.10.5 no yes yes
``MaxStrain'' tape or fabric II.1.10.6 no yes yes
``MaxStrain3D'' tape or fabric II.1.10.7 no yes yes
``CombStrain2D'' tape or fabric II.1.10.8 no yes yes
``MaxTotalStrain'' tape or fabric II.1.10.6 no yes yes
``MaxTotalStrain3D'' tape or fabric II.1.10.7 no yes yes
``CombTotalStrain2D'' tape or fabric II.1.10.8 no yes yes
``TsaiHill'' fabric II.1.10.9 no yes yes
``TsaiHill_b'' fabric II.1.10.10 no yes yes
``TsaiHill_c'' fabric II.1.10.11 no yes yes
``TsaiHill3D'' fabric II.1.10.12 no yes yes
``TsaiHill3D_b'' fabric II.1.10.13 no yes yes
``TsaiWu'' fabric II.1.10.14 no yes yes
``TsaiWu3D'' fabric II.1.10.15 no yes yes
``Hoffman'' fabric II.1.10.16 no yes yes
``Puck'' tape II.1.10.17 no yes yes
``Puck_b'' tape II.1.10.17 no yes yes
``Puck_c'' tape II.1.10.17 no yes yes
``Hashin'' tape II.1.10.18 no yes yes
``Hashin_b'' tape II.1.10.18 no yes yes
``Hashin_c'' tape II.1.10.18 no yes yes
``Hashin3D'' tape II.1.10.19 no yes yes
``Hashin3D_b'' tape II.1.10.19 no yes yes
``Hashin3D_c'' tape II.1.10.19 no yes yes
``YamadaSun'' tape II.1.10.20 no yes yes
``YamadaSun_b'' fabric II.1.10.21 no yes yes
``Honey3D'' honeycomb II.1.10.22 no yes yes
``HoneyShear'' honeycomb II.1.10.23 no yes yes
``HoneyShear_b'' honeycomb II.1.10.24 yes yes yes
``Ilss'' all II.1.10.25 yes yes yes
``Ilss_b'' all II.1.10.25 yes yes yes

Note that many of the criteria presented here are particular cases of a general quadratic criterion that requires first the calculation of a failure index:

\begin{multline}
f=F_{11}\sigma_{11}^2+F_{22}\sigma_{22}^2+F_{33}\sigma_{33}^2...
...a_{33}
+F_{4}\sigma_{23}+F_{5}\sigma_{31}+F_{6}\sigma_{12}\ ,
\end{multline}

then a test is done on the calculated value:

$\displaystyle f\leq 1\ .
$

Several failure theories discussed below are obtained by expressing the coefficients in the expressions above by expressions depending on the material allowables. Also for the Tsai-Wu criteria discussed in section II.1.10.14 and II.1.10.15, the parameters $ F_{ij}$ are directly characterized for the material.

Note that the 2D criteria defined in this section often correspond to the failure criteria defined in ESAComp.


Table II.1.3: Correspondence between the ESAComp failure criteria and the failure criteria defined in CLA module.
ESAComp failure criterion CLA criterion ID section
``Maximum Shear Stress (Tresca)'' ``Tresca2D'' II.1.10.1
``Von Mises'' ``VonMises2D'' II.1.10.2
``Maximum Strain'' (in ply axes) ``MaxStrain'' II.1.10.6
``Maximum Stress'' (in ply axes) ``MaxStress'' II.1.10.4
``Tsai-Wu'' ``TsaiWu'' II.1.10.14
``Tsai-Hill'' ``TsaiHill'' II.1.10.9
``Hoffman'' ``Hoffman'' II.1.10.16
``Simple Puck'' ``Puck'' II.1.10.17
``Modified Puck'' ``Puck_b'' II.1.10.17
``Hashin'' ``Hashin'' II.1.10.18


II.1.10.1 Tresca criterion (2D)

Using the stress tensor components a scalar equivalent shear stress is given by:

$\displaystyle \tau_{\text{eq}}^{\text{Tresca2D}}=\max\left( \frac{1}{2}\left\v...
...I}}\right\vert, \frac{1}{2}\left\vert\sigma_{\text{II}}\right\vert \right)\ .$ (II.1104)

This equivalent shear stress allows to define a Tresca failure index as follows:

$\displaystyle FI^{\text{Tresca2D}}=\frac{2\cdot\tau_{\text{eq}}^{\text{Tresca2D}}} {\sigma_t^{\text{all}}}\ ,$ (II.1105)

The reserve factor is given by:

$\displaystyle RF^{\text{Tresca2D}}=\frac{1}{FoS\cdot FI^{\text{Tresca2D}}} =\frac{\sigma_t^{\text{all}}} {FoS\cdot 2\cdot\tau_{\text{eq}}^{\text{Tresca2D}}}\ ,$ (II.1106)

and the strength ratio by:

$\displaystyle SR^{\text{Tresca2D}}=FoS\cdot FI^{\text{Tresca2D}} =\frac{FoS\cdot 2\cdot\tau_{\text{eq}}^{\text{Tresca2D}}} {\sigma_t^{\text{all}}}\ .$ (II.1107)

This criterion is referred to as ``Tresca2D'' criterion in FeResPost.


II.1.10.2 Von Mises criterion (2D)

Using the stress tensor components a scalar equivalent shear stress is given by:

$\displaystyle \sigma_{\text{eq}}^{\text{VonMises2D}}=\sqrt{\sigma_{11}^2 +\sigma_{22}^2-\sigma_{11}\sigma_{22}+3\tau_{12}^2}\ .$ (II.1108)

The corresponding failure index is:

$\displaystyle FI^{\text{VonMises2D}}=\frac{\sigma_{\text{eq}}^{\text{VonMises2D}}} {\sigma_t^{\text{all}}}\ ,$ (II.1109)

The reserve factor is given by:

$\displaystyle RF^{\text{VonMises2D}}=\frac{1}{FoS\cdot FI^{\text{VonMises2D}}}...
...rac{\sigma_t^{\text{all}}} {FoS\cdot\sigma_{\text{eq}}^{\text{VonMises2D}}}\ ,$ (II.1110)

and the strength ratio by:

$\displaystyle SR^{\text{VonMises2D}}=FoS\cdot FI^{\text{VonMises2D}} =\frac{FoS\cdot\sigma_{\text{eq}}^{\text{VonMises2D}}} {\sigma_t^{\text{all}}}\ .$ (II.1111)

This criterion is referred to as ``VonMises2D'' criterion in FeResPost.


II.1.10.3 Von Mises criterion (3D)

Using the stress tensor components a scalar equivalent shear stress is given by:

$\displaystyle \sigma_{\text{eq}}^{\text{VonMises3D}}=\sqrt{\sigma_{11}^2 +\sig...
...\sigma_{33}-\sigma_{33}\sigma_{11} +3\tau_{12}^2+3\tau_{23}^2+3\tau_{31}^2}\ .$ (II.1112)

As for the 2D version, the corresponding failure index is given by II.1.127, the reserve factor by II.1.128 and the strength ratio by II.1.129. This criterion is referred to as ``VonMises3D'' criterion in FeResPost.


II.1.10.4 Maximum stress criterion

The failure index is calculated as follows:

$\displaystyle FI^{\text{MaxStress}}=\max\left( \left\vert\frac{\sigma_{11}}{\s...
...ert, \left\vert\frac{\tau_{12}}{\tau_{12}^{\text{all}}}\right\vert \right)\ ,$ (II.1113)

in which the $ \sigma_{11}^{\text{all}}$ and $ \sigma_{22}^{\text{all}}$ allowables depend on the sign of $ \sigma_{11}$ and $ \sigma_{22}$ respectively. (Generally, tensile and compressive allowables are different for orthotropic materials.) The reserve factor is calculated as follows:

$\displaystyle RF^{\text{MaxStress}}=\frac{1}{FoS\cdot FI^{\text{MaxStress}}}\ ,$ (II.1114)

and the strength ratio is given by:

$\displaystyle SR^{\text{MaxStress}}=FoS\cdot FI^{\text{MaxStress}}\ .$ (II.1115)

This criterion is referred to as ``MaxStress'' criterion in FeResPost.


II.1.10.5 Maximum stress criterion (3D)

The failure index is calculated as follows:

$\displaystyle FI^{\text{MaxStress3D}}=\max\left( \left\vert\frac{\sigma_{11}}{...
...ert, \left\vert\frac{\tau_{31}}{\tau_{31}^{\text{all}}}\right\vert \right)\ ,$ (II.1116)

in which the $ \sigma_{11}^{\text{all}}$ , $ \sigma_{22}^{\text{all}}$ and $ \sigma_{33}^{\text{all}}$ allowables depend on the sign of $ \sigma_{11}$ , $ \sigma_{22}$ and $ \sigma_{33}$ respectively. (Generally, tensile and compressive allowables are different for orthotropic materials.) The reserve factor is calculated as follows:

$\displaystyle RF^{\text{MaxStress3D}}=\frac{1}{FoS\cdot FI^{\text{MaxStress3D}}}\ ,$ (II.1117)

and the strength ratio is given by:

$\displaystyle SR^{\text{MaxStress3D}}=FoS\cdot FI^{\text{MaxStress3D}}\ .$ (II.1118)

This criterion is referred to as ``MaxStress3D'' criterion in FeResPost.


II.1.10.6 Maximum strain criteria (2D)

The criterion is very similar to maximum stress criterion, except that it is calculated from the mechanical (or equivalent) strain tensor components. The failure index is calculated as follows:

$\displaystyle FI^{\text{MaxStrain}}=\max\left( \left\vert\frac{\epsilon_{11}^{...
...rac{\gamma_{12}^{\text{Mech}}}{\gamma_{12}^{\text{all}}}\right\vert \right)\ ,$ (II.1119)

in which the $ \epsilon_{11}^{\text{all}}$ and $ \epsilon_{22}^{\text{all}}$ allowables depend on the sign of $ \epsilon_{11}^{\text{Mech}}$ and $ \epsilon_{22}^{\text{Mech}}$ respectively. The reserve factor is calculated as follows:

$\displaystyle RF^{\text{MaxStrain}}=\frac{1}{FoS\cdot FI^{\text{MaxStrain}}}\ ,$ (II.1120)

and the strength ratio is given by:

$\displaystyle SR^{\text{MaxStrain}}=FoS\cdot FI^{\text{MaxStrain}}\ .$ (II.1121)

This criterion is referred to as ``MaxStrain'' criterion in FeResPost.

FeResPost proposes a second version of the criterion where the ``total'' strain tensor $ \left\{\epsilon\right\} $ is used instead of the mechanical strain tensor $ \left\{\epsilon^{\text{Mech}}\right\}$ . This version of the criterion is referred to as ``MaxTotalStrain'' criterion in FeResPost.


II.1.10.7 Maximum strain criterion (3D)

The criterion is very similar to maximum stress criterion, except that it is calculated from the mechanical (or equivalent) strain tensor components. The failure index is calculated as follows:

$\displaystyle FI^{\text{MaxStrain3D}}=\max\left( \left\vert\frac{\epsilon^{\te...
...rac{\gamma^{\text{Mech}}_{31}}{\gamma_{31}^{\text{all}}}\right\vert \right)\ .$ (II.1122)

in which the $ \epsilon_{11}^{\text{all}}$ , $ \epsilon_{22}^{\text{all}}$ and $ \epsilon_{33}^{\text{all}}$ allowables depend on the sign of the corresponding mechanical strain tensor component.

The reserve factor is calculated as follows:

$\displaystyle RF^{\text{MaxStrain3D}}=\frac{1}{FoS\cdot FI^{\text{MaxStrain3D}}}\ ,$ (II.1123)

and the strength ratio is given by

$\displaystyle SR^{\text{MaxStrain3D}}=FoS\cdot FI^{\text{MaxStrain3D}}\ .$ (II.1124)

This criterion is referred to as ``MaxStrain3D'' criterion in FeResPost.

FeResPost proposes a second version of the criterion where the ``total'' strain tensor $ \left\{\epsilon\right\} $ is used instead of the mechanical strain tensor $ \left\{\epsilon^{\text{Mech}}\right\}$ . This version of the criterion is referred to as ``MaxTotalStrain3D'' criterion in FeResPost.


II.1.10.8 Combined strain criterion (2D)

The criterion is a strain criterion that uses a combination of several components of the strain tensor. This criterion can be considered as a Tsai-type criterion adapted to strain tensor. The combined strain failure index is calculated as follows:

$\displaystyle FI^{\text{CombStrain2D}}=\max\left[ \sqrt{\left(\frac{\epsilon^{...
...frac{\gamma^{\text{Mech}}_{12}}{\gamma_{12}^{\text{all}}}\right)^2} \right]\ ,$ (II.1125)

in which the $ \epsilon_{11}^{\text{all}}$ and $ \epsilon_{22}^{\text{all}}$ allowables depend on the sign of $ \epsilon^{\text{Mech}}_{11}$ and $ \epsilon^{\text{Mech}}_{22}$ respectively. (The criterion is calculated from the mechanical or equivalent strain tensor components.) The reserve factor is calculated as follows:

$\displaystyle RF^{\text{CombStrain2D}}=\frac{1}{FoS\cdot FI^{\text{CombStrain2D}}}\ ,$ (II.1126)

and the strength ratio is given by:

$\displaystyle SR^{\text{CombStrain2D}}=FoS\cdot FI^{\text{CombStrain2D}}\ .$ (II.1127)

This criterion is referred to as ``CombStrain2D'' criterion in FeResPost.

FeResPost proposes a second version of the criterion where the ``total'' strain tensor $ \left\{\epsilon\right\} $ is used instead of the mechanical strain tensor $ \left\{\epsilon^{\text{Mech}}\right\}$ . This version of the criterion is referred to as ``CombTotalStrain2D'' criterion in FeResPost.


II.1.10.9 Tsai-Hill criterion

The Tsai-Hill criterion is a quadratic criterion. The 2D version of this criterion has a failure index calculated as follows:

$\displaystyle FI^{\text{TsaiHill}}= \left(\frac{\sigma_{11}}{\sigma_{11}^{\tex...
...2 -\frac{\sigma_{11}\sigma_{22}}{\left(\sigma_{11}^{\text{all}}\right)^2} \ .$ (II.1128)

Here again, the allowables $ \sigma_{11}^{\text{all}}$ and $ \sigma_{22}^{\text{all}}$ depend on the signs of $ \sigma_{11}$ and $ \sigma_{22}$ respectively. The Tsai-Hill failure index depends quadratically on the different components of the stress tensor. Therefore, the reserve factor is calculated as follows:

$\displaystyle RF^{\text{TsaiHill}}=\frac{1}{FoS\sqrt{FI^{\text{TsaiHill}}}}\ ,$ (II.1129)

and the strength ratio is given by:

$\displaystyle SR^{\text{TsaiHill}}=FoS\sqrt{FI^{\text{TsaiHill}}}\ .$ (II.1130)

This criterion is referred to as ``TsaiHill'' criterion in FeResPost.


II.1.10.10 Tsai-Hill criterion (version b)

This criterion is very similar to the one described in section II.1.10.10. It differs by the fact that only the tensile allowables $ \sigma_{11}^{\text{t,all}}$ and $ \sigma_{22}^{\text{t,all}}$ are considered for the calculations.

The calculation of reserve factor is as for the more classical Tsai-Hill criterion (II.1.147) and (II.1.148). This criterion is referred to as ``TsaiHill_b'' criterion in FeResPost.


II.1.10.11 Tsai-Hill criterion (version c)

This criterion is very similar to the one described in section II.1.10.9. The criterion has a failure index calculated as follows:

$\displaystyle FI^{\text{TsaiHill\_c}}= \left(\frac{\sigma_{11}}{\sigma_{11}^{\...
...rac{\sigma_{11}\sigma_{22}}{\left(\sigma_{11}^{\prime\text{all}}\right)^2} \ .$ (II.1131)

The criterion differs from the previous one by the fact that the allowable $ \sigma_{11}^{\prime\text{all}}$ is used in the calculation. This allowable is set to $ \sigma_{11}^{\text{t,all}}$ or $ \sigma_{11}^{\text{c,all}}$ depending on the sign of $ \sigma_{11}\sigma_{22}$ .

The calculation of reserve factor is as for the more classical Tsai-Hill criterion (II.1.147) and (II.1.148). This criterion is referred to as ``TsaiHill_c'' criterion in FeResPost.


II.1.10.12 Tsai-Hill criterion (3D)

The Tsai-Hill criterion is a quadratic criterion. The 3D version of this criterion has a failure index calculated as follows:

\begin{multline}
FI^{\text{TsaiHill3D}}=
\left(\frac{\sigma_{11}}{\sigma_{11...
...\text{all}}\right)^2}
\right)\sigma_{33}\sigma_{11}
\ .
\end{multline}

The allowables $ \sigma_{11}^{\text{all}}$ , $ \sigma_{22}^{\text{all}}$ and $ \sigma_{33}^{\text{all}}$ depend on the signs of $ \sigma_{11}$ , $ \sigma_{22}$ and $ \sigma_{33}$ respectively. The Tsai-Hill failure index depends quadratically on the different components of the stress tensor. The reserve factor is calculated as follows:

$\displaystyle RF^{\text{TsaiHill3D}}=\frac{1}{FoS\sqrt{FI^{\text{TsaiHill3D}}}}\ ,$ (II.1132)

and the strength ratio is given by:

$\displaystyle SR^{\text{TsaiHill3D}}=FoS\sqrt{FI^{\text{TsaiHill3D}}}\ .$ (II.1133)

This criterion is referred to as ``TsaiHill3D'' criterion in FeResPost.


II.1.10.13 Tsai-Hill criterion (3D version b)

This criterion is very similar to the one described in section II.1.10.13. It differs by the fact that only the tensile allowables $ \sigma_{11}^{\text{t,all}}$ , $ \sigma_{22}^{\text{t,all}}$ and $ \sigma_{33}^{\text{t,all}}$ are considered for the calculations.

The calculation of reserve factor is as for the more classical Tsai-Hill criterion (II.1.151) and (II.1.152). This criterion is referred to as ``TsaiHill3D_b'' criterion in FeResPost.


II.1.10.14 Tsai-Wu criterion

The Tsai-Wu criterion is a quadratic criterion. The 2D version of this criterion has a failure index calculated as follows:

\begin{multline}
FI^{\text{TsaiWu}}=
\frac{\sigma_{11}^2}{\sigma_{11}^{\text...
...gma_{22}
+2F_{12}^{\text{TW}}\sigma_{11}\sigma_{22}
\ .
\end{multline}

In this expression, $ F_{12}^{\text{TW}}$ is a material parameter to be obtained by characterization tests. This parameter must satisfy the following relation:

$\displaystyle \left(F_{12}^{\text{TW}}\right)^2<
\cfrac{1}{\sigma_{11}^{\text...
...1}^{\text{c,all}}
\sigma_{22}^{\text{t,all}}\sigma_{22}^{\text{c,all}}}\ .
$

Its units are, for example, $ \left[{\text{Pa}}^{-2}\right]$ . Sometimes the corresponding dimensionless parameter is used instead:

$\displaystyle F_{12}^{*}=F_{12}^{\text{TW}}
\sqrt{\sigma_{11}^{\text{t,all}}\...
...1}^{\text{c,all}}
\sigma_{22}^{\text{t,all}}\sigma_{22}^{\text{c,all}}}\ .
$

This dimensionless parameter must satisfy the relation:

$\displaystyle -1\leq F_{12}^{*}\leq 1\ .
$

The value $ F_{12}^{*}=-1/2$ corresponds to a generalized Von Mises criterion. The value

$\displaystyle F_{12}^{\text{TW}}
=\cfrac{-1}{2\sigma_{11}^{\text{t,all}}\sigma_{11}^{\text{c,all}}}
$

leads to the Hoffman criterion discussed in section II.1.10.16.

In some of the terms of expression (II.1.153) the components of Cauchy stress tensor appear linearly, and other terms, they appear quadratically. Therefore, the reserve factor expression is a little more complicated. FeResPost calculates it as follows:

$\displaystyle a=\frac{\sigma_{11}^2}{\sigma_{11}^{\text{t,all}}
\sigma_{11}^{...
...}}{\tau_{12}^{\text{all}}}\right)^2
+2F_{12}^{TW}\sigma_{11}\sigma_{22}\ ,
$

$\displaystyle b=\left(
\frac{1}{\sigma_{11}^{\text{t,all}}}
-\frac{1}{\sigm...
...{\text{t,all}}}
-\frac{1}{\sigma_{22}^{\text{c,all}}}\right)\sigma_{22}\ ,
$

$\displaystyle c=-1\ ,
$

$\displaystyle RF^{\text{TsaiWu}}=\frac{-b+\sqrt{b^2-4ac}}{2\cdot FoS\cdot a}\ ,$ (II.1134)

and the strength ratio is given by:

$\displaystyle SR^{\text{TsaiWu}}=\frac{1}{RF^{\text{TsaiWu}}}\ .$ (II.1135)

This criterion is referred to as ``TsaiWu'' criterion in FeResPost.


II.1.10.15 Tsai-Wu criterion (3D)

The 3D version of Tsai-Wu failure criterion leads to the following expression of the failure index:

\begin{multline}
FI^{\text{TsaiWu3D}}=
\frac{\sigma_{11}^2}{\sigma_{11}^{\te...
...gma_{33}
+2F_{31}^{\text{TW}}\sigma_{33}\sigma_{11}
\ .
\end{multline}

The values of $ F_{23}^{\text{TW}}$ and $ F_{31}^{\text{TW}}$ are submitted to the same limitations as $ F_{12}^{\text{TW}}$ in section II.1.10.14. The RF calculation is done as follows:

\begin{multline}
a=\frac{\sigma_{11}^2}{\sigma_{11}^{\text{t,all}}
\sigma_{1...
...{33}
+2F_{31}^{\text{TW}}\sigma_{33}\sigma_{11} \nonumber \ ,
\end{multline}

$\displaystyle b=\left(
\frac{1}{\sigma_{11}^{\text{t,all}}}
-\frac{1}{\sigm...
...{\text{t,all}}}
-\frac{1}{\sigma_{33}^{\text{c,all}}}\right)\sigma_{33}\ ,
$

$\displaystyle c=-1\ ,
$

$\displaystyle RF^{\text{TsaiWu3D}}=\frac{-b+\sqrt{b^2-4ac}}{2\cdot FoS\cdot a}\ ,$ (II.1136)

and the strength ratio is given by:

$\displaystyle SR^{\text{TsaiWu3D}}=\frac{1}{RF^{\text{TsaiWu3D}}}\ .$ (II.1137)

This criterion is referred to as ``TsaiWu3D'' criterion in FeResPost.


II.1.10.16 Hoffman criterion

The Hoffman criterion is very similar to the Tsai-Wu criterion. Only the last term of failure index is different:


In some of the terms of previous expression the components of Cauchy stress tensor appear linearly, and other terms, they appear quadratically. Therefore, the reserve factor expression is a little more complicated. FeResPost calculates it as follows:

$\displaystyle a=\frac{\sigma_{11}^2}{\sigma_{11}^{\text{t,all}}
\sigma_{11}^{...
...{11}\sigma_{22}}
{\sigma_{11}^{\text{t,all}}\sigma_{11}^{\text{c,all}}}\ ,
$

$\displaystyle b=\left(
\frac{1}{\sigma_{11}^{\text{t,all}}}
-\frac{1}{\sigm...
...{\text{t,all}}}
-\frac{1}{\sigma_{22}^{\text{c,all}}}\right)\sigma_{22}\ ,
$

$\displaystyle c=-1\ ,
$

$\displaystyle RF^{\text{Hoffman}}=\frac{-b+\sqrt{b^2-4ac}}{2\cdot FoS\cdot a}\ ,$ (II.1138)

and the strength ratio is given by:

$\displaystyle SR^{\text{Hoffman}}=\frac{1}{RF^{\text{Hoffman}}}\ .$ (II.1139)

This criterion is referred to as ``Hoffman'' criterion in FeResPost.


II.1.10.17 Puck criteria

This criterion is adapted to the justification of laminates with unidirectional plies. This criterion distinguishes two failure mode: one fiber failure mode in direction 1 and one matrix failure mode. One distinguishes three versions of the Puck criterion.

The first version of Puck failure index is calculated as follows:

$\displaystyle FI^{\text{Puck}}=\max\left( \left\vert\cfrac{\sigma_{11}}{\sigma...
...ight)^2 +\left(\cfrac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2} \right)\ ,$ (II.1140)

in which the allowables $ \sigma_{11}^{\text{all}}$ and $ \sigma_{22}^{\text{all}}$ depend on the signs of $ \sigma_{11}$ and $ \sigma_{22}$ respectively. The reserve factor is simply given by:

$\displaystyle RF^{\text{Puck}}=\cfrac{1}{FoS\cdot\max\left( \left\vert\cfrac{\...
...ght)^2 +\left(\cfrac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2} \right)}\ ,$ (II.1141)

and the strength ratio is given by:

$\displaystyle SR^{\text{Puck}}=FoS\cdot\max\left( \left\vert\cfrac{\sigma_{11}...
...ight)^2 +\left(\cfrac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2} \right)\ .$ (II.1142)

This criterion is referred to as ``Puck'' criterion in FeResPost.

A modified version of Puck criterion is defined as follows:

$\displaystyle FI^{\text{Puck\_b}}=\max\left( \left\vert\cfrac{\sigma_{11}}{\si...
...{22}^{\text{t,all}}} -\cfrac{1}{\sigma_{22}^{\text{c,all}}}\right) \right)\ .$ (II.1143)

The reserve factor is calculated as follows:

$\displaystyle RF_a=\cfrac{\sigma_{11}^{\text{all}}}{FoS\cdot\sigma_{11}}\ ,
$

$\displaystyle a=\cfrac{\sigma_{22}^2}{\sigma_{22}^{\text{c,all}}\sigma_{22}^{\text{t,all}}}
+\left(\cfrac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2\ ,
$

$\displaystyle b=\sigma_{22}\left(\frac{1}{\sigma_{22}^{\text{t,all}}}
-\frac{1}{\sigma_{22}^{\text{c,all}}}\right)\ ,
$

$\displaystyle c=-1\ ,
$

$\displaystyle RF_b=\frac{-b+\sqrt{b^2-4ac}}{2\cdot FoS\cdot a}\ ,
$

$\displaystyle RF^{\text{Puck\_b}}=\min\left(RF_a,RF_b\right)\ .$ (II.1144)

The calculation of failure index is based on the calculation of the reserve factor:

$\displaystyle FI^{\text{Puck\_b}}=\cfrac{1}{RF^{\text{Puck\_2}}}\ ,$ (II.1145)

in which the safety factor used in the calculation of the RF is 1. The advantage of this new expression is that the failure index is proportional to the components of stress tensor. This criterion is referred to as ``Puck_b'' criterion in FeResPost.

For the strength ratio, one uses the same expression, but the safety factor is not set to 1:

$\displaystyle SR^{\text{Puck\_b}}=\cfrac{1}{RF^{\text{Puck\_2}}}\ .$ (II.1146)

Sometimes, an additional term is given in the expression corresponding to the fiber failure and the modified version of Puck criterion is defined as follows:

$\displaystyle FI^{\text{Puck\_c}}=\max\left( \left\vert\cfrac{\sigma_{11}}{\si...
...{22}^{\text{t,all}}} -\cfrac{1}{\sigma_{22}^{\text{c,all}}}\right) \right)\ .$ (II.1147)

The calculation of reserve factor is done as for version ``b'' of Puck criterion, but one uses a modified expression for $ a$ parameter:

$\displaystyle a=\left(\cfrac{\sigma_{11}}{2\sigma_{11}^{\text{t,all}}}\right)^2...
...text{t,all}}}
+\left(\cfrac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2\ .
$

This criterion is referred to as ``Puck_c'' criterion in FeResPost.


II.1.10.18 Hashin criteria

This criterion is meant to be used for uni-directional materials. Direction 1 is assumed to be direction of fibers. One first presents the way the reserve factor is calculated:

Finally, the reserve factor is given by:

$\displaystyle RF^{\text{Hashin}}=\min\left(RF_b,RF_c\right)$ (II.1148)

Here again, the calculation of failure index is based on the calculation of the reserve factor:

$\displaystyle FI^{\text{Hashin}}=\cfrac{1}{RF^{\text{Hashin}}}\ ,$ (II.1149)

in which the safety factor used in the calculation of the RF is 1, and the strength ratio is calculated as

$\displaystyle SR^{\text{Hashin}}=\cfrac{1}{RF^{\text{Hashin}}}\ ,$ (II.1150)

in which one keeps the value of safety factor.

The criterion presented above is referred to by ``Hashin'' criterion in FeResPost. Correspondingly, one defines version ``Hashin_b'' in which only the fiber failure is checked and ``Hashin_c'' in which only the matrix failure is checked. (These correspond to the values $ RF_b$ and $ RF_c$ calculated above.)


II.1.10.19 Hashin criteria (3D)

A 3D version of the criterion defined in section II.1.10.18 is defined as follows

Finally, the reserve factor is given by:

$\displaystyle RF^{\text{Hashin3D}}=\min\left(RF_b,RF_c\right)$ (II.1151)

Here again, the calculation of failure index is based on the calculation of the reserve factor:

$\displaystyle FI^{\text{Hashin3D}}=\cfrac{1}{RF^{\text{Hashin3D}}}\ ,$ (II.1152)

in which the safety factor used in the calculation of the RF is set to 1, and the strength ratio is calculated as

$\displaystyle SR^{\text{Hashin3D}}=\cfrac{1}{RF^{\text{Hashin3D}}}\ ,$ (II.1153)

in which one keeps the value of safety factor.

The criterion presented above is referred to by ``Hashin3D'' criterion in FeResPost. Correspondingly, one defines version ``Hashin3D_b'' in which only the fiber failure is checked and ``Hashin3D_c'' in which only the matrix failure is checked. (These correspond to the values $ RF_b$ and $ RF_c$ calculated above.)


II.1.10.20 Yamada-Sun criterion

The Yamada-Sun criterion is a kind of Tsai criterion adapted to tape (unidirectional) materials. Its failure index is calculated as follows:

$\displaystyle FI^{\text{YamadaSun}}= \left(\frac{\sigma_{11}}{\sigma_{11}^{\text{all}}}\right)^2 +\left(\frac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2 \ .$ (II.1154)

The allowable $ \sigma_{11}^{\text{all}}$ depends on the sign of $ \sigma_{11}$ . The reserve factor is calculated as follows:

$\displaystyle RF^{\text{YamadaSun}} =\frac{1}{FoS\sqrt{FI^{\text{YamadaSun}}}}\ ,$ (II.1155)

and the strength ratio is given by:

$\displaystyle SR^{\text{YamadaSun}}=FoS\sqrt{FI^{\text{YamadaSun}}}\ .$ (II.1156)

This criterion is referred to as ``YamadaSun'' in FeResPost.


II.1.10.21 Yamada-Sun criterion (version b)

A second version of Yamada-Sun criterion more adapted to fabrics is proposed. The ``tape'' version of Yamada-Sun criterion is calculated in two directions, and the worst direction is considered for failure. The failure index is calculated as follows:

$\displaystyle FI^{\text{YamadaSun\_b}}= {\text{max}}\left[ \left(\frac{\sigma...
...\right)^2 +\left(\frac{\tau_{12}}{\tau_{12}^{\text{all}}}\right)^2 \right]\ .$ (II.1157)

The allowables $ \sigma_{11}^{\text{all}}$ and $ \sigma_{22}^{\text{all}}$ depends on the signs of $ \sigma_{11}$ and $ \sigma_{22}$ respectively. The reserve factor is calculated as follows:

$\displaystyle RF^{\text{YamadaSun\_b}} =\frac{1}{FoS\sqrt{FI^{\text{YamadaSun\_b}}}}\ ,$ (II.1158)

and the strength ratio is given by:

$\displaystyle SR^{\text{YamadaSun\_b}}=FoS\sqrt{FI^{\text{YamadaSun\_b}}}\ .$ (II.1159)

This criterion is referred to as ``YamadaSun_b'' in FeResPost.


II.1.10.22 3D honeycomb criterion

A general honeycomb criterion uses the out-of-plane tension/compression and the two out-of-plane shear components of the Cauchy stress tensor. The criterion read as follows:

$\displaystyle FI^{\text{Honey3D}}=\sqrt{ \left(\frac{\sigma_{33}}{\sigma_{33}^...
...all}}}\right)^2 +\left(\frac{\tau_{31}}{\tau_{31}^{\text{all}}}\right)^2} \ .$ (II.1160)

The allowable $ \sigma_{33}^{\text{all}}$ depends on the sign of $ \sigma_{33}$ (Generally, the compressive allowable is significantly smaller than the tensile one). The honeycomb material is generally defined in such a way that the allowable $ \tau_{31}^{\text{all}}=\tau_L$ is in ribbon direction (longitudinal allowable) and $ {\tau_{23}^{\text{all}}}$ is the transverse allowable $ \tau_W$ .

The reserve factor is calculated as follows:

$\displaystyle RF^{\text{Honey3D}}=\frac{1}{FoS\cdot FI^{\text{Honey3D}}}\ ,$ (II.1161)

and the strength ratio is given by:

$\displaystyle SR^{\text{Honey3D}}=FoS\cdot FI^{\text{Honey3D}}\ .$ (II.1162)

This criterion is referred to as ``Honey3D'' in FeResPost.


II.1.10.23 Honeycomb shear criterion

Depending of the modeling, the component $ \sigma_{33}$ of Cauchy stress tensor is sometimes zero. Then a simplified ``shear'' criterion is often used:

$\displaystyle FI^{\text{HoneyShear}}=\sqrt{ \left(\frac{\tau_{23}}{\tau_{23}^{...
...all}}}\right)^2 +\left(\frac{\tau_{31}}{\tau_{31}^{\text{all}}}\right)^2} \ .$ (II.1163)

The reserve factor is calculated as follows:

$\displaystyle RF^{\text{HoneyShear}} =\frac{1}{FoS\cdot FI^{\text{HoneyShear}}}\ ,$ (II.1164)

and the strength ratio is given by:

$\displaystyle SR^{\text{HoneyShear}} =FoS\cdot FI^{\text{HoneyShear}}\ .$ (II.1165)

This criterion is referred to as ``HoneyShear'' in FeResPost.


II.1.10.24 Honeycomb simplified shear criterion

Sometimes, one simplifies the criterion described in section II.1.10.23 by using the smallest shear allowable $ \tau_W$ . This new criterion is referred to as ``HoneyShear_b''.

As a single allowable is used, an equivalent shear stress can be defined:

$\displaystyle \tau_{\text{shear}}=\sqrt{\tau_{13}^2+\tau_{23}^2}\ .$ (II.1166)

The failure index is then calculated as follows:

$\displaystyle FI^{\text{HoneyShear\_b}}= \frac{\tau_{\text{shear}}}{\tau_{23}^{\text{all}}} \ .$ (II.1167)

The reserve factor is calculated as follows:

$\displaystyle RF^{\text{HoneyShear\_b}} =\frac{1}{FoS\cdot FI^{\text{HoneyShear\_b}}}\ ,$ (II.1168)

and the strength ratio is given by:

$\displaystyle SR^{\text{HoneyShear\_b}}=FoS\cdot FI^{\text{HoneyShear\_b}}\ .$ (II.1169)

This criterion is referred to as ``HoneyShear_b'' in FeResPost.


II.1.10.25 Inter-laminar shear criterion

The inter-laminar shear criterion is based on the comparison of inter-laminar shear stress with resin shear allowable:

$\displaystyle FI^{\text{Ilss}}=\cfrac{\tau_{\text{inter-laminar}}} {\tau_{\text{inter-laminar}}^{\text{all}}}\ ,$ (II.1170)

In which the inter-laminar shear stress is a scalar stress calculated as follows:

$\displaystyle \tau_{\text{inter-laminar}}=\sqrt{\tau_{xz}^2+\tau_{yz}^2}\ .
$

The reserve factor is of course:

$\displaystyle RF^{\text{Ilss}} =\frac{1}{FoS\cdot FI^{\text{Ilss}}}\ ,$ (II.1171)

and the strength ratio is given by:

$\displaystyle SR^{\text{Ilss}}=FoS\cdot FI^{\text{Ilss}}\ .$ (II.1172)

This criterion is different from the other criteria in this that it does not use ply material allowables. Instead, the ply ilss allowable or the laminate ilss allowable is used.

Note also that FeResPost calculates the inter-laminar shear criterion on the lower face of plies only. This criterion is referred to as ``Ilss'' in FeResPost.

A second version of the inter-laminar shear stress criterion, referred to as ``Ilss_b'' criterion in FeResPost, is defined in FeResPost. This criterion differs from the more usual ``Ilss'' criterion by the fact that its calculation is not limited to the lower face of laminate plies. It can also be calculated at mid-ply or ply upper face. This can be handy when one wishes to evaluate the ILSS criterion on finite element results that are extracted at mid ply thickness only. (Actually, the calculation of ILSS failure criterion using stresses extracted from FE results, or calculated from FE shell forces and moments, is the reason that has justified the introduction of this second version of ILSS criterion.)

II.1.11 Temperature diffusion in laminates

The thermal conservation equation in a solid material is written as follows:

$\displaystyle \rho C_p \cfrac{\partial T}{\partial t}=-\nabla\cdot$$\displaystyle \mbox{\boldmath$q$}$$\displaystyle +r\ ,$ (II.1173)

where $ T$ is the temperature, $ C_p$ is the heat capacity of the material and vector $ q$ is the heat flux. Note that $ C_p$ and $ \rho$ are used for transient thermal calculations only. Generally, the heat flux is related to the gradient of temperature by Fourrier's law:

$\displaystyle \mbox{\boldmath$q$}$$\displaystyle =-$$\displaystyle \mbox{\boldmath$\lambda^T$}$$\displaystyle \cdot\nabla T\ ,$ (II.1174)

where $ \lambda$ $ ^T$ is the tensor of thermal conductivity coefficients.

When thermal conductivity calculations are performed with laminates, two homogenized quantities must first be calculated:

These two quantities are obtained by integrating material properties along laminate thickness.

II.1.11.1 Material thermal parameters

Two scalar parameters influence the transient thermal behavior of laminates: the density $ \rho$ and the heat capacity (per unit of mass) $ C_p$ . As these two parameters are scalar, their characteristics do not depend on the possible anisotropy of the material.

On the other hand, the thermal conductivity $ \lambda$ $ ^T$ is a $ 2^{\text{nd}}$ order tensor. Generally, for an anisotropic material, the tensor may be written as follows:

$\displaystyle \left(\begin{array}{c}
q_1 \\
q_2 \\
q_3
\end{array}\right)...
...31} & \lambda^T_{32} & \lambda^T_{33}
\end{array}\right)
\cdot
\nabla T\ .
$

The tensor $ \lambda^T$ is symmetric, so that only 6 components must be defined:

$\displaystyle \left(\begin{array}{c} q_1 \\ q_2 \\ q_3 \end{array}\right) =-\le...
...a^T_{31} & \lambda^T_{23} & \lambda^T_{33} \end{array}\right) \cdot \nabla T\ .$ (II.1175)

For an orthotropic material, the previous equations reduces to:

$\displaystyle \left(\begin{array}{c} q_1 \\ q_2 \\ q_3 \end{array}\right) =-\le...
...ambda^T_{22} & 0 \\ 0 & 0 & \lambda^T_{33} \end{array}\right) \cdot \nabla T\ .$ (II.1176)

Then, only three material parameters define the thermal conductivity. Finally, for an isotropic material the thermal conductivity is defined by a single parameter $ \lambda^T$ :

$\displaystyle \left(\begin{array}{c} q_1 \\ q_2 \\ q_3 \end{array}\right) =-\le...
... \\ 0 & \lambda^T & 0 \\ 0 & 0 & \lambda^T \end{array}\right) \cdot \nabla T\ .$ (II.1177)

II.1.11.2 In-plane and out-of-plane components

As has been done for the motion equations, one assumes a decoupling of in-plane laminate thermal conductivity and out-of-plane conductivity. Therefore, the thermal flux is separated into an in-plane flux:

$\displaystyle \left\{q_{12}\right\}
=\left(\begin{array}{c}
q_1 \\
q_2
\end{array}\right)\ ,
$

and the out-of-plane component $ q_3$ .

Correspondingly, the tensor of thermal conductivity coefficients is separated into an in-plane conductivity tensor:

$\displaystyle \left[\lambda^T\right]_{12}
=\left(\begin{array}{cc}
\lambda^T_...
...& \lambda^T_{12} \\
\lambda^T_{12} & \lambda^T_{22}
\end{array}\right)\ ,
$

and the out-of-plane conductivity $ \lambda^T_{33}$ . The ``shear'' components $ \lambda^T_{31}$ and $ \lambda^T_{23}$ are neglected in the homogenization theory. (This means that out-of-plane and in-plane conductivities are decoupled.)

II.1.11.3 In-plane rotations of vectorial and tensorial properties

The material scalar properties are left unmodified by in-plane rotations. The same is true for the out-of-plane quantities $ q_3$ and $ \lambda^T_{33}$ . The transformation of components for $ \left\{q\right\}$ and $ \left[\lambda^T\right]$ is noted with the transformation matrices defined in expressions (II.1.12) to (II.1.17):

$\displaystyle \left[\lambda^T\right]_{12} = \left[S_-\right]_{(\theta)} \cdot \left[\lambda^T\right]_{xy} \cdot \left[S_+\right]_{(\theta)}\ ,$ (II.1178)

$\displaystyle \left[\lambda^T\right]_{xy} = \left[S_+\right]_{(\theta)} \cdot \left[\lambda^T\right]_{12} \cdot \left[S_-\right]_{(\theta)}\ ,$ (II.1179)

$\displaystyle \left\{q\right\}_{xy} = \left[S_+\right]_{(\theta)} \cdot \left\{q\right\}_{12}\ ,$ (II.1180)

$\displaystyle \left\{q\right\}_{12} = \left[S_-\right]_{(\theta)} \cdot \left\{q\right\}_{xy}\ .$ (II.1181)

II.1.11.4 Integration along the laminate thickness

One considers separately the laminate thermal in-plane conductivity, thermal out-of-plane conductivity and thermal capacity.

II.1.11.4.1 In-plane conductivity

To calculate the laminate in-plane thermal conductivity properties, one assumes that temperature is constant along the laminate thickness. Consequently, the temperature gradient does not depend on $ z$ . The thermal flux $ \left\{q\right\}_{xy}$ however depends on $ z$ because the thermal conductivity does:

$\displaystyle \left\{q\right\}_{xy}(z) =\left[\lambda^T\right]_{xy}(z) \cdot\left\{\nabla T\right\}_{xy}\ .$ (II.1182)

The laminate in-plane thermal conductivity is calculated as follows:

$\displaystyle \left\{Q\right\}_{xy}$ $\displaystyle =\int_{-h/2}^{-h/2}\left\{q\right\}_{xy}(z)$   d$\displaystyle z\ ,$    
  $\displaystyle =\int_{-h/2}^{-h/2}\left[\lambda^T\right]_{xy}(z) \cdot\left\{\nabla T\right\}_{xy}$   d$\displaystyle z\ ,$    
  $\displaystyle =\left(\int_{-h/2}^{-h/2}\left[\lambda^T\right]_{xy}(z)\hspace{0.4em}\text{d}z\right) \cdot\left\{\nabla T\right\}_{xy} \ ,$    
  $\displaystyle =\left(\int_{-h/2}^{-h/2} \left[S_+\right]_{(\theta)}(z) \cdot \l...
...\theta)}(z)\hspace{0.4em}\text{d}z\right) \cdot\left\{\nabla T\right\}_{xy} \ ,$    
  $\displaystyle =\left[\Lambda^T\right]_{xy}\cdot\left\{\nabla T\right\}_{xy} \ .$    

In the previous equation, one introduced the laminate in-plane thermal conductivity:

$\displaystyle \left[\Lambda^T\right]_{xy}$ $\displaystyle =\int_{-h/2}^{-h/2} \left[S_+\right]_{(\theta)}(z) \cdot \left[\lambda^T\right]_{12}(z) \cdot \left[S_-\right]_{(\theta)}(z)$   d$\displaystyle z\ ,$    
  $\displaystyle =\sum_{k=1}^N \left[S_+\right]_{(\theta_k)} \cdot \left[\lambda_k^T\right]_{12} \cdot \left[S_-\right]_{(\theta_k)}\ e_k\ ,$ (II.1183)

where $ e_k$ is the thickness of ply $ k$ .

II.1.11.4.2 Out-of-plane conductivity

One assumes that out-of-plane thermal flux is constant across laminate thickness. Then, as thermal conductivity depends on $ z$ , so will the out-of-plane gradient of temperature:

$\displaystyle \nabla T_z=\cfrac{q_z}{\lambda^T_{33}}\ .$ (II.1184)

The integration across the thickness gives the difference of temperature between upper and lower laminate surfaces:

$\displaystyle T_{\text{sup}}-T_{\text{inf}}$ $\displaystyle =\int_{-h/2}^{-h/2}\nabla T_z$   d$\displaystyle z\ ,$    
  $\displaystyle =\int_{-h/2}^{-h/2}\cfrac{q_z}{\lambda^T_{33}}$   d$\displaystyle z\ ,$    
  $\displaystyle =\left(\int_{-h/2}^{-h/2}\cfrac{1}{\lambda^T_{33}}\hspace{0.4em}\text{d}z\right) q_z\ ,$    
  $\displaystyle =R^T_{33}\ q_z\ .$    

In previous expression, one introduced the out-of-plane thermal resistance:

$\displaystyle R^T_{33}$ $\displaystyle =\int_{-h/2}^{-h/2}\cfrac{1}{\lambda^T_{33}}$   d$\displaystyle z\ ,$    
  $\displaystyle =\sum_{k=1}^N\cfrac{1}{(\lambda^T_{33})_k}\ e_k\ .$ (II.1185)

II.1.11.4.3 Thermal capacity

To estimate the laminate thermal capacity, one again assumes a temperature constant across the laminate thickness. Then, the heat energy stored per unit of surface is:

$\displaystyle U$ $\displaystyle =\int_{-h/2}^{-h/2}\rho C_p T$   d$\displaystyle z\ ,$    
  $\displaystyle =\left(\int_{-h/2}^{-h/2}\rho C_p\hspace{0.4em}\text{d}z\right) T\ ,$    
  $\displaystyle =\left(\rho C_p h\right)_{\text{lam}} T\ .$    

In previous expression, one introduced the surfacic thermal capacity

$\displaystyle \left(\rho C_p h\right)_{\text{lam}}$ $\displaystyle =\int_{-h/2}^{-h/2}\rho C_p$   d$\displaystyle z\ ,$    
  $\displaystyle =\sum_{k=1}^N\left(\rho C_p\right)_k\ e_k\ .$ (II.1186)

Note that in Nastran, when a thermal material MAT4 or MAT5 is defined, the density $ \rho$ and the heat capacity per unit of mass $ C_p$ are defined separately. So it is the responsibility of the user to select appropriate values for these two quantities. Also, in Nastran, the thickness is defined separately in the PSHELL property card. (PCOMP or PCOMPG cards do not accept thermal materials.)

II.1.12 Moisture diffusion in laminates

The moisture conservation equation in a solid material is written as follows:

$\displaystyle \cfrac{\partial M}{\partial t}=-\nabla\cdot$$\displaystyle \mbox{\boldmath$q$}$$\displaystyle _M\ ,$ (II.1187)

where $ M$ is the moisture content (for example in [kg/m$ ^3$ ]) and vector $ q$ $ _M$ is the massic flux of moisture (in [kg/m$ ^2$ s]). When several materials are present in a structure, as often in laminates, it is more practical to work with moisture percentage $ H$ (in [%w]). The moisture percentage is related to moisture content by the following expression:

$\displaystyle M=c\cdot H\ ,$ (II.1188)

where $ c$ is the moisture coefficient in [kg/(m$ ^3$ %w)].

Generally, the moisture flux is related to the gradient of moisture by Fick's law, and the moisture diffusion equation can be written:

$\displaystyle \cfrac{\partial H}{\partial t} =\nabla\cdot\left(\mbox{\boldmath$\lambda$}^H\cdot\nabla H\right)\ .$ (II.1189)

here $ \lambda$ $ ^H$ is the tensor of moisture conductivity coefficients. This expression is very similar to the equation of thermal diffusion: As one works with moisture percentages, nothing equivalent to the laminate thermal capacity $ \left(\rho C_p h\right)_{\text{lam}}$ is defined for moisture.


II.1.13 Units

All the quantities introduced in this Chapter have been given without dimensions. Since version 3.0.1, units can be attributed to all the CLA quantities defined in FeResPost, except of course the dimensionless quantities. This allows the user to express all the CLA quantities in a units system compatible, for example, with the unit system used for finite element modeling.

An engineer should be able to figure out the units of the different quantities introduced in this Chapter from the definition given for the quantities or from the expressions used for their calculations. However, we think it might be useful to remind the units of the different quantities to avoid ambiguities.

In the rest of the section, one assumes a consistent set of units compatible with MKS system is used. This is what we recommend for FeResPost, as well as finite element models. The default ``base'' units are:

One notes that units have been defined for force and energy, despite the fact that they can be derived from mass, length and time units. Practically, the definition of additional base units for the force and energy are of interest for mechanical engineers.

All the other units of FeResPost are obtained by combining these base units. The most important ones are summarized in Tables II.1.4, II.1.5 and II.1.6.

If one of the base units above is modified, the user is responsible for modifying the derived units coherently. For example, we expressed the moisture content in [%w]. This influences the units of moisture content $ H$ as well as the units of coefficients of moisture expansion $ \beta$ .


Table II.1.4: Units of the different quantities related to loading, stresses or load responses.
Quantities Symbols Units
Strains $ \epsilon$ , $ \gamma$ or $ \Gamma$ [L/L] or [-]
Stresses $ \sigma$ or $ \tau$ [F/L$ ^2$ ]
Curvatures $ \kappa$ [1/L]
Forces $ \{N\}$ or $ \{Q\}$ [F/L]
Moments $ \{M\}$ [FL/L] or [F]
Temperatures $ T$ [T]
Moistures $ H$ [W] (always [%w])
Failure indices $ FI$ []
Reserve factors $ RF$ []


Table II.1.5: Units of the different quantities related to material properties.
Quantities Symbols Units
Materials stiffnesses or moduli $ C_{ijkl}$ , $ [C]$ , $ E$ or $ G$ [F/L$ ^2$ ]
Materials compliance matrices $ c_{ijkl}$ , $ [c]$ , $ [g]$ [L$ ^2$ /F] ]
Poisson coefficients $ \nu_{ij}$ [-]
Thermal conductivity $ \lambda^T_{ij}$ [E/(LT)]
Coefficients of thermal expansion $ \alpha$ [L/(LT)] or [1/T]
Moisture conductivity $ \lambda^H_{ij}$ [1/(Lt)]
Coefficients of moisture expansion $ \beta$ [L/(LW)] or [1/W]
Density $ \rho$ [M/L$ ^3$ ]
Heat specific capacity $ C_p$ [E/(MT)]
Coefficients of quadratic failure criteria $ F_i$ [L$ ^2$ /F]
Coefficients of quadratic failure criteria $ F_{ij}$ [L$ ^4$ /F$ ^2$ ]


Table II.1.6: Units of the different quantities related to laminate properties.
Quantities Symbols Units
Thicknesses $ t$ [L]
Membrane stiffness matrix $ [A]$ [F/L]
Membrane-bending coupling stiffness matrix $ [B]$ [F]
Bending stiffness matrix $ [D]$ [FL]
Out-of-plane shear stiffness matrix $ [G]$ [F/L]
Membrane compliance matrix $ [a]$ [L/F]
Membrane-bending coupling compliance matrix $ [b]$ [1/F]
Bending compliance matrix $ [d]$ [1/(FL)]
Out-of-plane shear compliance matrix $ [g]$ [L/F]
In-plane thermal conductivity matrix $ \left[\Lambda^T\right]$ [E/(tT)]
Out-of-plane thermal resistance $ R^T_{33}$ [L$ ^2$ Tt/E]
Surfacic heat capacity $ \left(\rho C_p h\right)_{\text{lam}}$ [E/L$ ^2$ /T]
In-plane moisture conductivity matrix $ \left[\Lambda^H\right]$ [1/t]
Out-of-plane moisture resistance $ R^H_{33}$ [L$ ^2$ t]

Other units systems can be used with the CLA classes. Each CLA object has an attribute corresponding to the units system in which all its characteristics are defined.

The units of CLA object can be obtained with ``getUnits'' method that returns a Hash containing pairs of Strings. The first String describes the kind of dimension: ``L'' for length, ``M'' for mass, ``t'' for time, ``T'' for temperature, ``E'' for energy, ``W'' for moisture and ``F'' for force. The second String corresponds to one of the unit listed above for the selected quantity.

The units of a CLA object can also be modified by calling ``setUnits'' or ``changeUnits'' methods. The ``setUnits'' method change the units attributed to an object without modifying the values of the different quantities defining the object. (No unit conversion is done.) The ``changeUnits'' method performs the units conversion. Both ``setUnits'' and ``changeUnits'' methods has an Hash argument corresponding to the value returned by ``getUnits'' method.


next up previous contents index
Next: II.2 The ``ClaDb'' class Up: II. Composite Reference Manual Previous: II.0 Introduction   Contents   Index
FeResPost 2017-05-28