Samer Adeeb

Constitutive Laws: Linear Elastic Materials

Matrix of Material Properties of Linear Elastic Materials:

A linear elastic material is a material that exhibits a linear relationship between the components of the stress tensor and the components of the strain tensor. A linear elastic material constitutive law, under the assumption of small strains, is fully represented by a linear map C_s between the stress matrix \sigma\in\mathbb{M}^3 and the infinitesimal strain matrix \varepsilon\in\mathbb{M}^3 such that C_s:\sigma\mapsto\varepsilon. In general, the stress matrix has 3\times 3=9 components and the strain matrix has 3\times 3=9 components. A linear relationship between these components would be represented by a fourth-order tensor with 3\times 3\times 3\times 3=81 compnents. These components are unique for a particular linear elastic material, and thus, they are considered to be material constants. By adopting the Einstein summation convention the linear map between the stress and the strain tensors can be written in component form as

    \[ \sigma_{ij}=\sum_{k,l=1}^3 \left(C_{ijkl}\varepsilon_{kl}\right)=C_{ijkl}\varepsilon_{kl} \]

The components of the fourth-order tensor C_s can be visualized in the following form:

    \[ C_s=\left( \begin{array}{ccc} \left(\begin{array}{ccc} C_{1111} & C_{1112} & C_{1113}\\ C_{1121} & C_{1122} & C_{1123}\\ C_{1131} & C_{1132} & C_{1133} \end{array}\right) & \left(\begin{array}{ccc} C_{1211} & C_{1212} & C_{1213}\\ C_{1221} & C_{1222} & C_{1223}\\ C_{1231} & C_{1232} & C_{1233} \end{array}\right) & \left(\begin{array}{ccc} C_{1311} & C_{1312} & C_{1313}\\ C_{1321} & C_{1322} & C_{1323}\\ C_{1331} & C_{1332} & C_{1333} \end{array}\right) \\ \left(\begin{array}{ccc} C_{2111} & C_{2112} & C_{2113}\\ C_{2121} & C_{2122} & C_{2123}\\ C_{2131} & C_{2132} & C_{2133} \end{array}\right) & \left(\begin{array}{ccc} C_{2211} & C_{2212} & C_{2213}\\ C_{2221} & C_{2222} & C_{2223}\\ C_{2231} & C_{2232} & C_{2233} \end{array}\right) & \left(\begin{array}{ccc} C_{2311} & C_{2312} & C_{2313}\\ C_{2321} & C_{2322} & C_{2323}\\ C_{2331} & C_{2332} & C_{2333} \end{array}\right) \\ \left(\begin{array}{ccc} C_{3111} & C_{3112} & C_{3113}\\ C_{3121} & C_{3122} & C_{3123}\\ C_{3131} & C_{3132} & C_{3133} \end{array}\right) & \left(\begin{array}{ccc} C_{3211} & C_{3212} & C_{3213}\\ C_{3221} & C_{3222} & C_{3223}\\ C_{3231} & C_{3232} & C_{3233} \end{array}\right) & \left(\begin{array}{ccc} C_{3311} & C_{3312} & C_{3313}\\ C_{3321} & C_{3322} & C_{3323}\\ C_{3331} & C_{3332} & C_{3333} \end{array}\right) \end{array} \right) \]

For example, the elements of the matrix on the top left of the above equation can be viewed as the matrix of coefficients that are used to calculate the stress component \sigma_{11} as follows:

    \[ \sigma_{11}=C_{1111}\varepsilon_{11}+ C_{1112}\varepsilon_{12} + C_{1113}\varepsilon_{13}+C_{1121}\varepsilon_{21} +C_{1122}\varepsilon_{22} + C_{1123}\varepsilon_{23}+C_{1131}\varepsilon_{31} + C_{1132}\varepsilon_{32} + C_{1133}\varepsilon_{33} \]

Reductions due to Symmetries of the Stress and the Strain:

Because the stress matrix is symmetric, i.e., \sigma_{ij}=\sigma_{ji}, then this leads to the following relationship among the components of C_s

    \[ C_{ijkl}=C_{jikl} \]

This can be further explained using an example. Assume a strain state defined by a strain matrix with zero values except for the strain \varepsilon_{22}. Then, \sigma_{12}=C_{1222}\varepsilon_{22} and \sigma_{21}=C_{2122}\varepsilon_{22}. Since \sigma_{21}=\sigma_{12}, then, C_{1222}=C_{2122}. This applies to all the coefficients in the form of C_{ijkl} and C_{jikl}.

Similarly, because the strain matrix is symmetric, i.e., \varepsilon_{ij}=\varepsilon_{ji}, then this leads to the following relationship among the components of C_s

    \[ C_{ijkl}=C_{ijlk} \]

This can be further explained using an example. Assume a strain state defined by a strain matrix with zero values except for the strain components \varepsilon_{12}=\varepsilon_{21}. Then, \sigma_{23}=C_{2312}\varepsilon_{12}+C_{2321}\varepsilon_{21}. Since \varepsilon_{12}=\varepsilon_{21} we have: \sigma_{23}=\frac{C_{2312}+C_{2321}}{2}\varepsilon_{12}+\frac{C_{2312}+C_{2321}}{2}\varepsilon_{21}, i.e., the coefficients C_{2312} and C_{2321} can be chosen to be equal. This applies to all the coefficients in the form of C_{ijkl} and C_{ijlk}.

Therefore, because of the symmetries of the stress and the strain matrices, the number of independent coefficients in the matrix C_s can be reduced to 6\times 6=36 independent coefficients. This simplifies the relationship which can now be represented in matrix form by assuming that the stress and the strain matrices are vectors in \mathbb{R}^6 with the following forms:

    \[ \sigma=\left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)\qquad \varepsilon=\left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)= \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\\gamma_{12}\\\gamma_{13}\\\gamma_{23} \end{array}\right) \]

Notice that traditionally, the engineering shear strains \gamma_{12}, \gamma_{13}, and \gamma_{23} are used in the constitutive relationships and these are equal to double the corresponding shear strains. The linear map between the stress and the strain can be represented by a matrix C\in\mathbb{M}^6. The relationship between the vector representation of the stress and the vector representation of the strain has the following form:

    \[ \sigma=C\varepsilon\qquad \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)= \left(\begin{array}{cccccc} C_{1111}&C_{1122}&C_{1133}&C_{1112}&C_{1113}&C_{1123}\\ C_{2211}&C_{2222}&C_{2233}&C_{2212}&C_{2213}&C_{2223}\\ C_{3311}&C_{3322}&C_{3333}&C_{3312}&C_{3313}&C_{3323}\\ C_{1211}&C_{1222}&C_{1233}&C_{1212}&C_{1213}&C_{1123}\\ C_{1311}&C_{1322}&C_{1333}&C_{1312}&C_{1313}&C_{1323}\\ C_{2311}&C_{2322}&C_{2333}&C_{2312}&C_{2313}&C_{2323} \end{array} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \]

The above relationship is described by 36 independent coefficients. So far, there are no restrictions imposed on C for it to be symmetric. However, in the next section we will show that the restriction imposed by the preservation of energy leads to the symmetry of C.

Reductions due to Preservation of Energy:

For a linear elastic material, the energy stored during loading and/or the energy released during unloading have to be independent of the path of loading and/or unloading. Thus, for such a material, a potential energy function \overline{U} has to exist such that the stress is derived from this potential energy function as follows:

    \[ \sigma_{ij}=\frac{\partial \overline{U}}{\partial \varepsilon_{ij}} \]

See chapter 6 in the book “Elasticity: Theory, Applications, and Numerics” by Martin H. Sadd for a detailed description of this concept.
The stress component \sigma_{ij} can be obtained using the following relationship:

    \[ \sigma_{ij}=\sum_{k,l=1}^3 C_{ijkl}\varepsilon_{kl}\Rightarrow C_{ijkl}=\frac{\partial \sigma_{ij}}{\partial \varepsilon_{kl}}=\frac{\partial^2\overline{U}}{\partial \varepsilon_{ij}\partial \varepsilon_{kl}} \]

The equality of the second partial derivatives of \overline{U} leads to the symmetry of C as follows:

    \[ \frac{\partial^2\overline{U}}{\partial \varepsilon_{ij}\partial \varepsilon_{kl}}=\frac{\partial^2\overline{U}}{\partial \varepsilon_{kl}\partial \varepsilon_{ij}}\Rightarrow C_{ijkl}=C_{klij} \]

In matrix form, the relationship between the vector representation of the stress and the vector representation of the strain has the following form with 21 independent coefficients:

    \[ \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)= \left( \begin{matrix} C_{1111} & C_{1122} & C_{1133} & C_{1112} &C_{1113} & C_{1123}\\ & C_{2222} & C_{2233} & C_{2212}&C_{2213}& C_{2223} \\ & & C_{3333} & C_{3312}&C_{3313}& C_{3323} \\ & & & C_{1212}&C_{1213}& C_{1223} \\ & \multicolumn{2}{c}{\text{sym.}} & & C_{1313}& C_{1323} \\ & & & & & C_{2323} \end{matrix} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \]

The symmetry of the coefficient matrix will is crucial so that the material model would not predict energy creation or dissipation during loading cycles. This can be illustrated with an example in which the relationship between the in-plane strain and in-plane stress of a material is given by the following non-symmetric relationship:

    \[ \left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\end{array}\right)=\frac{1}{E}\left(\begin{matrix}1 & -0.3\\-0.4 & 1\end{matrix}\right) \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\end{array}\right) \]

Assume that this material is loaded such that \sigma_{11} is increased gradually from 0 to 100MPa. Then, \sigma_{11} is kept constant and \sigma_{22} is increased gradually from 0 to 100MPa. Then, assume that both stresses are decreases simultaneously from 100MPa to 0. Since this is an elastic material, no energy should be created or dissipated through the entire loading process since the final state of stress is identical to the initial state of stress even though the paths of loading and unloading are different. The energy in a path from state a to state b can be calculated using the equation:

    \[ \overline{U}=\int_a^b \! \sigma_{11} \, \mathrm{d}\varepsilon_{11}+\int_a^b \! \sigma_{22} \, \mathrm{d}\varepsilon_{22} \]

To simplify the integrals, we will use the following relationships:

    \[ \mathrm{d}\varepsilon_{11}=\frac{\mathrm{d}\sigma_{11}}{E}-0.3\frac{\mathrm{d}\sigma_{22}}{E} \]

    \[ \mathrm{d}\varepsilon_{22}=-0.4\frac{\mathrm{d}\sigma_{11}}{E}+1\frac{\mathrm{d}\sigma_{22}}{E} \]

Therefore, the integrals have the following form:

    \[ \overline{U}=\int_a^b \! \frac{\sigma_{11}}{E} \, \mathrm{d}\sigma_{11}+\int_a^b \! -0.3\frac{\sigma_{11}}{E} \, \mathrm{d}\sigma_{22}+\int_a^b \! -0.4\frac{\sigma_{22}}{E} \, \mathrm{d}\sigma_{11}+\int_a^b \! \frac{\sigma_{22}}{E} \, \mathrm{d}\sigma_{22} \]

The loading can be divided into three paths. In the first path, \sigma_{11} increases from 0 to 100MPa while \sigma_{22} stays constant at zero. Therefore, the energy can be calculated as follows:

    \[ \overline{U}_{path 1}=\int_0^{100} \! \frac{\sigma_{11}}{E} \, \mathrm{d}\sigma_{11}+0+0+0=\frac{100^2}{2E} \]

In the second path, \sigma_{11} stays constant at 100MPa while \sigma_{22} increases from 0 to 100MPa. Therefore, the energy in the second path can be calculated as follows:

    \[ \overline{U}_{path 2}=0+\int_0^{100} \! -0.3\frac{\sigma_{11}}{E} \, \mathrm{d}\sigma_{22}+0+\int_0^{100} \! \frac{\sigma_{22}}{E} \, \mathrm{d}\sigma_{22}=(-0.3)\frac{100}{E}{100}+\frac{100^2}{2E} \]

In the third path, \sigma_{11}=\sigma_{22}, \mathrm{d}\sigma_{11}=\mathrm{d}\sigma_{22} and both change from 100MPa to 0. Therefore, the energy in the third path can be calculated as follows:

    \[ \overline{U}_{path 3}=\int_{100}^{0} \! \frac{\sigma_{11}}{E} \, \mathrm{d}\sigma_{11}+\int_{100}^{0} \! -0.3\frac{\sigma_{11}}{E} \, \mathrm{d}\sigma_{22}+\int_{100}^{0} \! -0.4\frac{\sigma_{22}}{E} \, \mathrm{d}\sigma_{11}+\int_{100}^{0} \! \frac{\sigma_{22}}{E} \, \mathrm{d}\sigma_{22}=-\frac{100^2}{2E}+0.3\frac{100^2}{2E}+0.4\frac{100^2}{2E}-\frac{100^2}{2E} \]

Therefore, the total energy after the three paths is equal to:

    \[ \overline{U}=\overline{U}_{path 1}+\overline{U}_{path 2}+\overline{U}_{path 3}=(0.4-0.3)\frac{100^2}{2E} \]

which is contrary to the expectation that the material model preserves energy. The above loading cycle shows a dissipation of energy. The amount of energy dissipated is related to the difference between the off-diagonal components in the matrix of material coefficients. Had the matrix of coefficients been symmetric, the total energy stored or generated during the loading cycle would have been zero. Can you identify a different loading cycle in which the material would generate energy?

Reductions due to Material Symmetry:

The number of constants can be further reduced, depending on the symmetries of the material. If a material possesses a certain symmetry, this implies that the components of C_s do not change upon coordinate transformations within the specified symmetry. Since C_s is a fourth-order tensor, upon changing the coordinate system using an orthogonal matrix Q\in\mathbb{M}^3, the new material constants have the following form:

    \[ C'_{ijkl}=Q_{im}Q_{jn}Q_{ko}Q_{lp}C_{mnop} \]

If a material’s behaviour is independent of a certain change of the coordinate system described by a particular orthogonal matrix Q, then the components of C'_s should be equal to those of C_s. Such orthogonal matrix Q is termed a material symmetry. For example, if a material response is symmetric around a plane, then a reflection matrix Q_1\in\mathbb{M}^3 across this plane would produce the same components for the coefficients matrix. If a material response is the same within any rotation inside a plane, then a rotation Q_2\in\mathbb{M}^3 within that plane would produce the same components for C'_s as those for C_s. The number of the independent constants of C_s can be further reduced from 21, depending on the level of anisotropy of the material. In the following section, the different levels of anisotropy are described, and the physical interpretation of some of the material constants is given.

Orthotropic Linear Elastic Materials:

In an orthotropic linear elastic material, the 21 material constants are reduced to 9 material constants fully describing the relationship between the stress and the small strain. A proof of this argument is illustrated using the mathematica code found here. Assuming that the basis vectors e_1, e_2, and e_3 are the axis of orthotropy, then the relationship between the stress and the strain is given by:

(1)   \begin{equation*} \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)= \left( \begin{matrix} C_{1111} & C_{1122} & C_{1133} & 0 &0 & 0\\ & C_{2222} & C_{2233} & 0&0& 0 \\ & & C_{3333} & 0&0& 0\\ & & & C_{1212}&0& 0 \\ & \multicolumn{2}{c}{\text{sym.}} & & C_{1313}& 0 \\ & & & & & C_{2323} \end{matrix} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \end{equation*}

Notice that symmetry around each plane leads to the normal stresses being independent of the shear strains and vice versa. In other words, if an orthotropic material is pulled along a direction perpendicular to a plane of symmetry, then the associated shear strains are zero. For orthotropic materials, the relationship between the small strains and the stresses is traditionally given as function of nine physical material constants E_{11}, E_{22}, E_{33}, \nu_{12}, \nu_{13}, \nu_{23}, G_{12}, G_{13, and G_{23} where:
E_{ii}: The ratio of the uniaxial stress to uniaxial strain in the direction of the basis vector e_i.
\nu_{ij}: Negative the ratio between the axial strain in the direction of the basis vector e_j to the axial strain in the direction of the basis vector e_i when the material is uniaxially stressed in the direction of the basis vector e_i.
G_{ij}: The ratio between the shear stress \sigma_{ij} to the engineering shear strain \gamma_{ij}=2\varepsilon_{ij} in the plane defined by the directions of e_i and e_j.

The relationship can be written in terms of the above physical material constants as follows:

    \[ \begin{split} \varepsilon_{11} & =\frac{\sigma_{11}}{E_{11}}-\nu_{21}\frac{\sigma_{22}}{E_{22}}-\nu_{31}\frac{\sigma_{33}}{E_{33}}\\ \varepsilon_{22} & =-\nu_{12}\frac{\sigma_{11}}{E_{11}}+\frac{\sigma_{22}}{E_{22}}-\nu_{32}\frac{\sigma_{33}}{E_{33}}\\ \varepsilon_{33} & =-\nu_{13}\frac{\sigma_{11}}{E_{11}}-\nu_{23}\frac{\sigma_{22}}{E_{22}}+\frac{\sigma_{33}}{E_{33}}\\ 2\varepsilon_{12} & =\gamma_{12}=\frac{\sigma_{12}}{G_{12}}\\ 2\varepsilon_{13} & =\gamma_{13}=\frac{\sigma_{13}}{G_{13}}\\ 2\varepsilon_{23} & =\gamma_{23}=\frac{\sigma_{23}}{G_{23}} \end{split} \]

The relationship has the following matrix form:

(2)   \begin{equation*} \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)= \left( \begin{matrix} \frac{1}{E_{11}} & -\frac{\nu_{21}}{E_{22}} & -\frac{\nu_{31}}{E_{33}} & 0 &0 & 0\\ -\frac{\nu_{12}}{E_{11}}& \frac{1}{E_{22}} & -\frac{\nu_{32}}{E_{33}} & 0&0& 0 \\ -\frac{\nu_{13}}{E_{11}}& -\frac{\nu_{23}}{E_{22}}& \frac{1}{E_{33}} & 0&0& 0 \\ 0 &0 &0 & \frac{1}{G_{12}}&0& 0 \\ 0 & 0& 0 &0 & \frac{1}{G_{13}}& 0 \\ 0&0 &0 & 0& 0& \frac{1}{G_{23}} \end{matrix} \right) \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right) \end{equation*}

Notice that symmetry of the matrix C implies that 

    \[ \frac{\nu_{ij}}{E_{ii}}=\frac{\nu_{ji}}{E_{jj}} \]

Transversely Isotropic Linear Elastic Materials:

In a transversely isotropic linear elastic material, five material constants describe the relationship between the stress and the small strain. A proof of this argument is illustrated using the mathematica code found here. Assuming full isotropy in the plane formed by the basis vectors e_1 and e_2. i.e., the material response is identical given any rotation around e_3. Also, assuming that the plane formed by the vectors e_1 and e_2 is a plane of symmetry, then, the relationship between the stress and the strain is given by five material constants:

(3)   \begin{equation*} \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)= \left( \begin{matrix} C_{1111} & C_{1122} & C_{1133} & 0 &0 & 0\\ & C_{1111} & C_{1133} & 0&0& 0 \\ & & C_{3333} & 0&0& 0\\ & & & \frac{C_{1111}-C_{1122}}{2}&0& 0 \\ & \multicolumn{2}{c}{\text{sym.}} & & C_{1313}& 0 \\ & & & & & C_{1313} \end{matrix} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \end{equation*}

For transversely isotropic materials, the relationship between the small strains and the stresses are traditionally given as function of five physical material constants E, \nu, E_{33}, \nu_{13}, and G_{13} where:
E: The ratio of the uniaxial stress to uniaxial strain in any direction in the isotropy plane.
\nu: Negative the ratio between the transverse strain to the axial strain when the material is uniaxially stressed in any direction in the isotropy plane.
E_{33}: The ratio of the uniaxial stress to uniaxial strain in the direction perpendicular to the isotropy plane.
\nu_{13}: Negative the ratio between the axial strain \varepsilon_{33} to the axial strain \varepsilon_{11} when the material is stressed uniaxially in the direction of e_1
G_{13}: The ratio between the shear stress \sigma_{13} to the engineering shear strain \gamma_{13}=2\varepsilon_{13} in the plane defined by the directions of e_1 and e_3.
In this case, the relationship between the strains and the stresses can be written as:

    \[ \begin{split} \varepsilon_{11} & =\frac{\sigma_{11}}{E}-\nu\frac{\sigma_{22}}{E}-\nu_{31}\frac{\sigma_{33}}{E_{33}}\\ \varepsilon_{22} & =-\nu\frac{\sigma_{11}}{E}+\frac{\sigma_{22}}{E}-\nu_{31}\frac{\sigma_{33}}{E_{33}}\\ \varepsilon_{33} & =-\nu_{13}\frac{\sigma_{11}}{E}-\nu_{13}\frac{\sigma_{22}}{E}+\frac{\sigma_{33}}{E_{33}}\\ 2\varepsilon_{12} & =\gamma_{12}=\frac{2(1+\nu)}{E}\sigma_{12}\\ 2\varepsilon_{13} & =\gamma_{13}=\frac{\sigma_{13}}{G_{13}}\\ 2\varepsilon_{23} & =\gamma_{23}=\frac{\sigma_{23}}{G_{13}} \end{split} \]

The relationship has the following matrix form:

(4)   \begin{equation*} \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)= \left( \begin{matrix} \frac{1}{E} & -\frac{\nu}{E} & -\frac{\nu_{31}}{E_{33}} & 0 &0 & 0\\ -\frac{\nu}{E}}& \frac{1}{E} & -\frac{\nu_{31}}{E_{33}} & 0&0& 0 \\ -\frac{\nu_{13}}{E}& -\frac{\nu_{13}}{E}& \frac{1}{E_{33}} & 0&0& 0 \\ 0 & 0 & 0 & \frac{2(1+\nu)}{E}&0& 0 \\ 0 & 0 & 0 &0 & \frac{1}{G_{13}}& 0 \\ 0 & 0 &0 &0 & 0 & \frac{1}{G_{13}} \end{matrix} \right) \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right) \end{equation*}

The symmetry of the matrix C implies that

    \[ \frac{\nu_{13}}{E}=\frac{\nu_{31}}{E_{33}} \]

Isotropic Linear Elastic Materials:

In an isotropic linear elastic material, the response of the material to stress is independent of the coordinate system chosen. In other words, the relationship between the stress and the strain is independent of the coordinate system. This assumption reduces the number of material constants that describe the behaviour of the material in any direction to only two. A proof of this argument is illustrated using the mathematica code found here. The relationship between the stress and the strain components has the following form independent of the chosen coordinate system:

(5)   \begin{equation*} \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)= \left( \begin{matrix} C_{1111} & C_{1122} & C_{1122} & 0 &0 & 0\\ & C_{1111} & C_{1122} & 0&0& 0 \\ & & C_{1111} & 0&0& 0\\ & & & \frac{C_{1111}-C_{1122}}{2}&0& 0 \\ & \multicolumn{2}{c}{\text{sym.}} & & \frac{C_{1111}-C_{1122}}{2}& 0 \\ & & & & & \frac{C_{1111}-C_{1122}}{2} \end{matrix} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \end{equation*}

The two constants that are generally used are the Young’s modulus E and Poisson’s ratio \nu, with the following physical definitions:
Young’s modulus E: The slope of the stress-strain curve during uniaxial tensile tests of a sample of a linear elastic isotropic material.
Poisson’s ratio \nu: Negative the ratio of the transverse strain (perpendicular to the applied load) to the axial strain (in the direction of the applied load) during uniaxial tensile tests of a sample of a linear elastic isotropic material.

For such materials, the strains are given by the following equations

    \[ \begin{split} \varepsilon_{11} & =\frac{\sigma_{11}}{E}-\nu\frac{\sigma_{22}}{E}-\nu\frac{\sigma_{33}}{E}\\ \varepsilon_{22} & =-\nu\frac{\sigma_{11}}{E}+\frac{\sigma_{22}}{E}-\nu\frac{\sigma_{33}}{E}\\ \varepsilon_{33} & =-\nu\frac{\sigma_{11}}{E}-\nu\frac{\sigma_{22}}{E}+\frac{\sigma_{33}}{E}\\ 2\varepsilon_{12} & =\gamma_{12}=\frac{2(1+\nu)}{E}\sigma_{12}\\ 2\varepsilon_{13} & =\gamma_{13}=\frac{2(1+\nu)}{E}\sigma_{13}\\ 2\varepsilon_{23} & =\gamma_{23}=\frac{2(1+\nu)}{E}\sigma_{23} \end{split} \]

The relationship has the following matrix form:

(6)   \begin{equation*} \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)= \left( \begin{matrix} \frac{1}{E} & -\frac{\nu}{E} & -\frac{\nu}{E} & 0 &0 & 0\\ -\frac{\nu}{E}}& \frac{1}{E} & -\frac{\nu}{E} & 0&0& 0 \\ -\frac{\nu}{E}& -\frac{\nu}{E}& \frac{1}{E} & 0&0& 0 \\ 0 &0 &0 & \frac{1}{G}&0& 0 \\ 0 & 0 & 0& 0 & \frac{1}{G}& 0 \\ 0 & 0 &0 & 0 & 0 & \frac{1}{G} \end{matrix} \right) \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right) \end{equation*}

Where:

    \[ G=\frac{E}{2(1+\nu)} \]

The inverse of the above relationship has the following form:

(7)   \begin{equation*} \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right) = \frac{E}{(1-2\nu)(1+\nu)}\left( \begin{matrix} (1-\nu) & \nu & \nu & 0 &0 & 0\\ \nu & (1-\nu) & \nu & 0&0& 0 \\ \nu & \nu & (1-\nu) & 0&0& 0 \\ 0 &0 & 0& \frac{1-2\nu}{2}&0& 0 \\ 0 & 0 & 0& 0& \frac{1-2\nu}{2}& 0 \\ 0 & 0 & 0 &0 &0 & \frac{1-2\nu}{2} \end{matrix} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \end{equation*}

Note that equations 1, 2, 3, and 4 are only valid in particular coordinate systems while equations 5, 6, and 7 can be used in any coordinate system.

Alternative Material Constants:

The stress-strain relationships (Equations 6 and 7) rely on the definitions of E and \nu. However, other material constants can be used. The wikipedia page “Elastic Modulus” offers an extensive list of such constants. In particular, the shear modulus G, the bulk modulus K, and the Lamé’s constants \lambda and \mu are sometimes used. The following are their respective definitions:

Shear Modulus:

The shear modulus is the slope of the shear stress versus the engineering shear strain of a linear elastic isotropic material. It is related to E and \nu as follows:

    \[ G=\frac{E}{2(1+\nu)} \]

The above relationship can be deduced for an isotropic material as follows. We first use Equation 6 but assume that G, E, and \nu are independent of each other. Assume a state of pure shear stress with \sigma_{12}=\sigma_{21} being the only non-zero stress component. Using the isotropic linear elastic constitutive law equation (Equation 6), then the only non-zero strain components are 2\varepsilon_{12}=2\varepsilon_{21}=\frac{\sigma_{12}}{G}. Therefore, the stress and strain matrices are:

    \[ \sigma=\left(\begin{matrix} 0 & \sigma_{12} & 0\\ \sigma_{12} & 0 & 0\\ 0 & 0 & 0 \end{matrix}\right)\qquad \varepsilon=\left(\begin{matrix} 0 & \varepsilon_{12} & 0\\ \varepsilon_{12} & 0 & 0\\ 0 & 0 & 0 \end{matrix}\right) \]

We now apply a 45^{\circ} coordinate transformation in the plane described by e_1 and e_2 using the rotation matrix Q:

    \[ Q=\left(\begin{matrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0\\ -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0\\ 0 & 0 & 1 \end{matrix}\right)\]

In the new coordinate system, \sigma' and \varepsilon' have the following forms:

    \[\begin{split} \sigma'=Q\sigma Q^T=\left(\begin{matrix} \sigma_{12} & 0 & 0\\ 0 & -\sigma_{12} & 0\\ 0 & 0 & 0 \end{matrix}\right)=\left(\begin{matrix} \sigma'_{11} & 0 & 0\\ 0 & \sigma'_{22} & 0\\ 0 & 0 & 0 \end{matrix}\right) \\ \varepsilon'=Q\varepsilon Q^T =\left(\begin{matrix} \varepsilon_{12} & 0 & 0\\ 0 & -\varepsilon_{12} & 0\\ 0 & 0 & 0 \end{matrix}\right)=\left(\begin{matrix} \varepsilon'_{11} & 0 & 0\\ 0 & \varepsilon'_{22} & 0\\ 0 & 0 & 0 \end{matrix}\right) \end{split} \]

Since the material is isotropic, the relationship shown in Equation 6 applies in the new coordinate system as well. Therefore:

    \[ \varepsilon'_{11}=\frac{\sigma'_{11}}{E}-\nu\frac{\sigma'_{22}}{E}\Rightarrow \varepsilon_{12}=\frac{\sigma_{12}}{E}(1+\nu) \]

However, we already know from applying the constitutive equation in the original coordinate system that:

    \[ \varepsilon_{12}=\frac{\sigma_{12}}{2G} \]

However, since \varepsilon'_{11}=\varepsilon_{12} we get the required relationship: G=\frac{E}{2(1+\nu)}.

Bulk Modulus:

The bulk modulus is the slope of the applied hydrostatic stress versus the volumetric strain (\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}) of a linear elastic isotropic material. It is related to E and \nu as follows:

    \[ K=\frac{E}{3(1-2\nu)} \]

The relationship for the bulk modulus can be obtained as follows. First we assume a hydrostatic state of stress \sigma=pI where I is the identity matrix, and p\in\mathbb{R} is a real number representing the applied pressure. In this case, when adopting Equation 6, the developed strains are \varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}=p\frac{1-2\nu}{E}. Therefore,

    \[ \varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}=p \frac{3(1-2\nu)}{E} \]

Therefore, the slope of the relationship between p and the volumetric strain \varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33} is equal to \frac{E}{3(1-2\nu)}.

Lamé’s constants:

These constants are defined as follows:

    \[ \lambda=\frac{E\nu}{(1+\nu)(1-2\nu)}\qquad \mu=G=\frac{E}{2(1+\nu)} \]

These two constants provide a simplified way of writing 7 as follows.

(8)   \begin{equation*} \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right) = \left( \begin{matrix} \lambda+2\mu & \lambda & \lambda & 0 &0 & 0\\ \lambda & \lambda+2\mu & \lambda & 0&0& 0 \\ \lambda & \lambda & \lambda+2\mu & 0&0& 0 \\ 0 &0 & 0& \mu &0& 0 \\ 0 & 0 & 0& 0& \mu & 0 \\ 0 & 0 & 0 &0 &0 & \mu \end{matrix} \right) \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right) \end{equation*}

or in simple component form

    \[ \sigma_{ij}=2\mu\varepsilon_{ij}+\sum_{k=1}^3\lambda \varepsilon_{kk}\delta_{ij} \]

 

Restrictions on the Elastic Moduli:

The strain energy density per unit volume \overline{U} for a linear elastic material under a general state of stress is equal to:

    \[ \overline{U}=\frac{1}{2}\sum_{i,j=1}^3\sigma_{ij}\varepsilon_{ij}=\frac{1}{2}\left(\sigma_{11}\varepsilon_{11}+\sigma_{22}\varepsilon_{22}+\sigma_{33}\varepsilon_{33}+2\sigma_{12}\varepsilon_{12}+2\sigma_{13}\varepsilon_{13}+2\sigma_{23}\varepsilon_{23}\right) \]

As shown above, since each of the stress and strain matrices are symmetric, they can be represented as vectors in \mathbb{R}^6. The dot product can then be used to find a simplified expression for the strain energy density function:

    \[ \overline{U}=\frac{1}{2}(\sigma\cdot\varepsilon)=\frac{1}{2}(\varepsilon\cdot\sigma)=\frac{1}{2}(\varepsilon\cdot C\varepsilon) \]

We can use two physical restrictions on the materials to impose mathematical restrictions on the entries of the matrix C. These two restrictions are that the material is stable and deformable. Stability means that whenever a state of strain is applied, energy is stored in the material. Stability is also equivalent to the statement that the material deformation increases in the direction of the increasing stress. The equivalent mathematical restriction on C is:

    \[ \forall \varepsilon:\overline{U}= \frac{1}{2}(\varepsilon\cdot C\varepsilon)>0 \]

Since C is symmetric, then C has to be positive definite. This implies that all its eigenvalues have to be positive. This argument will be used to find restrictions on the values of the elastic moduli described above.

Restrictions on the Elastic Moduli of Linear Elastic Isotropic Materials:

The distinct eigenvalues of the matrix C in Equation 6 are:

    \[ \frac{1-\nu}{E},\qquad\frac{1+\nu}{E},\qquad \frac{2(1+\nu)}{E} \]

For these to be positive, the following inequalities have to be satisfied:

    \[ E>0 \qquad -1<\nu<\frac{1}{2}  \]

The following mathematica code presents the calculations above:
View Mathematica Code

G = Ee/2/(1 + Nu);
Cc = {{1/Ee, -Nu/Ee, -Nu/Ee, 0, 0, 0}, {-Nu/Ee, 1/Ee, -Nu/Ee, 0, 0,0}, {-Nu/Ee, -Nu/Ee, 1/Ee, 0, 0, 0}, {0, 0, 0, 1/G, 0, 0}, {0, 0,0, 0, 1/G, 0}, {0, 0, 0, 0, 0, 1/G}};
Cc // MatrixForm
a = Eigenvalues[Cc] 
Reduce[{a[[1]]>0, a[[2]]>0}, {Ee, Nu}]

Materials possessing a negative Poisson’s ratio are common, see for example the article here or the Wikipedia article on Auxetics (materials possessing negative Poisson’s ratio).
 
A more intuitive method to find the same restrictions is by assuming different states of stresses and ensuring that the energy stored per unit volume under these states of stresses is positive. The first state is a uniaxial state of stress where \sigma_{11} is the only non-zero stress component. Using Equation 6, the stresses and the strains have the following forms:

    \[ \sigma=\left(\begin{matrix} \sigma_{11} & 0 & 0\\0&0&0\\0&0&0 \end{matrix}\right)\qquad \varepsilon=\left(\begin{matrix} \frac{\sigma_{11}}{E} & 0 & 0\\0&-\frac{\nu\sigma_{11}}{E}&0\\0&0&-\frac{\nu\sigma_{11}}{E} \end{matrix}\right) \]

The restriction that the energy stored has to be greater than zero leads to E>0 as follows:

    \[ \overline{U}=\frac{1}{2}\sigma_{11}\frac{\sigma_{11}}{E}=\frac{\sigma_{11}^2}{2E}>0\Rightarrow E>0 \]

The second state is a shear state of stress where \sigma_{12}=\sigma_{21} are the only non-zero stress components. Using Equation 6, the stresses and the strains have the following forms:

    \[ \sigma=\left(\begin{matrix} 0& \sigma_{12} & 0\\\sigma_{12}&0&0\\0&0&0 \end{matrix}\right)\qquad \varepsilon=\left(\begin{matrix} 0 & \frac{\sigma_{12}}{2G} & 0\\\frac{\sigma_{12}}{2G} &0&0\\0&0&0 \end{matrix}\right) \]

The restriction that the energy stored has to be greater than zero leads to \nu>-1 as follows:

    \[ \overline{U}=\frac{1}{2}\left(2\sigma_{12}\frac{\sigma_{12}}{2G}\right)=\frac{\sigma_{12}^2}{2G}>0\Rightarrow G>0\Rightarrow \frac{E}{2(1+\nu)}>0\Rightarrow\nu>-1 \]

The third state is a state of hydrostatic stress where \sigma_{11}=\sigma_{22}=\sigma_{33}=p\in\mathbb{R} are the only non-zero stress components. Using Equation 6, the stresses and the strains have the following forms:

    \[ \sigma=\left(\begin{matrix} p& 0 & 0\\0&p&0\\0&0&p \end{matrix}\right)\qquad \varepsilon=\frac{(1-2\nu)p}{E}\left(\begin{matrix} 1 & 0 & 0\\0&1&0\\0&0&1 \end{matrix}\right) \]

The restriction that the energy stored has to be greater than zero leads to \nu<0.5 as follows:

    \[ \overline{U}=\frac{1}{2}\left(3p\frac{(1-2\nu)p}{E}\right)=\frac{3p^2(1-2\nu)}{2E}>0\Rightarrow \nu<\frac{1}{2} \]

Restrictions on the Elastic Moduli of Linear Elastic Orthotropic Materials:

For orthotropic linear elastic materials, the positive definiteness of C imply the following inequalities:

    \[ \begin{split} E_{11}, E_{22}, E_{33}, G_{12}, G_{13}, G_{23} >0 \\ |\nu_{ij}|<\left(\frac{E_{ii}}{E_{jj}}\right)^{\frac{1}{2}}\qquad i\neq j\\ 1-\nu_{12}\nu_{21}-\nu_{13}\nu_{31}-\nu_{23}\nu_{32}-2\nu_{21}\nu_{13}\nu_{32}>0 \end{split} \]

The above equations show that values for the different orthotropic Poisson’s ratios are not bounded by –1 and 0.5. In fact, orthotropic materials (typically, fibrous biological materials like ligaments) can exhibit values beyond those bounds. For example, ligaments can exhibit a Poisson’s ratio that is higher than 0.5.

Important Note:

The above restrictions ensure that, physically, a material is stable (possessing a rising stress strain curve) and that any state of stress will cause nonzero strain values inside the material. It is important to realize that it cannot be argued that E < 0 and/or that \nu=0.5 are physically impossible. They are simply unacceptable constants for linear elastic materials with a positive definite constitutive model. Materials that become unstable with loading (e.g., decrease in stress accompanying an increase in strain) or incompressible materials (materials exhibiting no strain for a stress state \sigma=pI\in\mathbb{M}^3, p\in\mathbb{R} may exist but cannot be modelled using the linear elastic constitutive model

 

Plane Isotropic Linear Elastic Materials Constitutive Laws:

In many solid mechanics problems, it is possible, due to loading and geometry symmetries, to simplify particular problems to be solved in \mathbb{R}^2 rather than \mathbb{R}^3. In such cases, the matrix of elastic moduli shown in Equations 6 and 7 can be reduced from a 6\times6 matrix to a 3\times3 or 4\times4 matrix that accounts for such symmetries. In the following, the material is assumed to be isotropic linear elastic and the measure of strain is the infinitesimal strain tensor.

Plane Stress:

The plane stress assumption is used when the non-zero stress components are \sigma_{11}, \sigma_{22}, and \sigma_{12}=\sigma_{21}. Usually, this assumption is used for situations when a material is allowed to freely contract and/or expand in one direction (called the thickness direction) while the loads are applied only in the plane perpendicular to the thickness direction. In this situation, the stress matrix has the following form:

    \[ \sigma=\left( \begin{matrix}\sigma_{11}&\sigma_{12} & 0 \\ \sigma_{12}& \sigma_{22} & 0\\0 &0&0 \end{matrix} \right) \]

Using Equation 6, \varepsilon_{13}=\varepsilon_{23}=\varepsilon_{31}=\varepsilon_{32}=0 and the strain matrix has the following form:

    \[ \varepsilon=\left( \begin{matrix}\varepsilon_{11}&\varepsilon_{12} & 0 \\ \varepsilon_{12}& \varepsilon_{22} & 0\\0 &0&\varepsilon_{33} \end{matrix} \right) \]

The relationship between the stress and strain components shown in Equation 6 can now be reduced to only account for the stresses that are not necessarily equal to zero:

(9)   \begin{equation*} \left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right)=\frac{1}{E}\left(\begin{matrix}1&-\nu&0\\-\nu & 1 & 0\\0 & 0 & 2(1+\nu)\end{matrix}\right)\left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right) \end{equation*}

The inverse of this relationship has the following form:

(10)   \begin{equation*} \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right)=\frac{E}{(1+\nu)}\left(\begin{matrix}\frac{1}{1-\nu}&\frac{\nu}{1-\nu}&0\\\frac{\nu}{1-\nu} & \frac{1}{1-\nu} & 0\\0 & 0 & \frac{1}{2}\end{matrix}\right)\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right) \end{equation*}

In the plane stress assumption, the material is allowed to expand and contract freely in the third direction. The corresponding value of the axial strain in the thickness direction can be obtained using the equation:

    \[ \varepsilon_{33}=-\frac{\nu}{E}(\sigma_{11}+\sigma_{22}) \]

Plane Strain:

The plane strain assumption is used when the non-zero strain components are \varepsilon_{11}, \varepsilon_{22}, and \varepsilon_{12}=\varepsilon_{21}. Usually, this assumption is used for situations when a material is restrained from contracting or expanding in one direction (called the thickness direction). The loads are applied in the plane perpendicular to the thickness direction and a restraining stress is applied in the thickness direction. In this situation, the strain matrix has the following form:

    \[ \varepsilon=\left( \begin{matrix}\varepsilon_{11}&\varepsilon_{12} & 0 \\ \varepsilon_{12}& \varepsilon_{22} & 0\\0 &0&0 \end{matrix} \right) \]

Using Equation 7, \sigma_{13}=\sigma_{23}=\sigma_{31}=\sigma_{32}=0 and the stress matrix has the following form:

    \[ \sigma=\left( \begin{matrix}\sigma_{11}&\sigma_{12} & 0 \\ \sigma_{12}& \sigma_{22} & 0\\0 &0&\sigma_{33} \end{matrix} \right) \]

The relationship between the stress and strain components shown in Equation 7 can now be reduced to only account for the strains that are not necessarily equal to zero:

(11)   \begin{equation*} \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right)=\frac{E}{(1+\nu)(1-2\nu)}\left(\begin{matrix}1-\nu &\nu&0\\\nu & 1-\nu & 0\\0 & 0 & \frac{1-2\nu}{2}\end{matrix}\right)\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right) \end{equation*}

The inverse of this relationship has the following form:

(12)   \begin{equation*} \left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\2\varepsilon_{12}\end{array}\right)=\frac{1+\nu}{E}\left(\begin{matrix}1-\nu&-\nu &0\\-\nu & 1-\nu& 0\\0 & 0 & 2\end{matrix}\right)\left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right) \end{equation*}

In the plane strain assumption, the material is restrained from expanding or contracting in the third direction. The restraining axial stress in the thickness direction can be obtained using the equation:

    \[ \sigma_{33}=\nu(\sigma_{11}+\sigma_{22}) \]

Axisymmetry:

Axisymmetry is used for modelling cylindrical shapes when the loads and boundary conditions are also cylindrical. By considering an cross section where e_2 is aligned with the axis of symmetry, the shear stresses \sigma_{13}=\sigma_{31} and \sigma_{23}=\sigma_{32} can be assumed to be zero. In this situation, the stress matrix has the following form:

    \[ \sigma=\left( \begin{matrix}\sigma_{11}&\sigma_{12} & 0 \\ \sigma_{12}& \sigma_{22} & 0\\0 &0&\sigma_{33} \end{matrix} \right) \]

Using Equation 6, \varepsilon_{13}=\varepsilon_{23}=\varepsilon_{31}=\varepsilon_{32}=0 and the strain matrix has the following form:

    \[ \varepsilon=\left( \begin{matrix}\varepsilon_{11}&\varepsilon_{12} & 0 \\ \varepsilon_{12}& \varepsilon_{22} & 0\\0 &0&\varepsilon_{33} \end{matrix} \right) \]

The axisymmetry assumption leads to the following relationship between the displacement u_1, the radius r and the strain component \varepsilon_{33} as follows:

    \[ \varepsilon_{33}=\frac{u_1}{r} \]

The relationship between the stress and strain components shown in Equation 6 can now be reduced to only account for the stresses that are not necessarily equal to zero:

    \[ \left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\end{array}\right)=\frac{1}{E}\left(\begin{matrix}1&-\nu&-\nu&0\\-\nu & 1 & -\nu&0\\-\nu & -\nu & 1 &0\\0 & 0 &0& 2(1+\nu)\end{matrix}\right)\left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\end{array}\right) \]

The inverse of this relationship has the following form:

    \[ \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\end{array}\right)=\frac{E}{(1+\nu)(1-2\nu)}\left(\begin{matrix}1-\nu &\nu&\nu&0\\\nu & 1-\nu &\nu& 0\\\nu & \nu& 1-\nu &0\\0 &0& 0 & \frac{1-2\nu}{2}\end{matrix}\right)\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\end{array}\right) \]

 

Examples and Problems:

Example 1:

The state of strain inside a continuum is represented by the small strain matrix:

    \[ \varepsilon=\left(\begin{matrix}0.002 & 0.001 & 0.002\\ 0.001 &-0.003& 0.001\\0.002 & 0.001 & 0.0015\end{matrix}\right) \]

Assume that the material is steel with Young’s modulus of 210 GPa and Poisson’s ratio of 0.3. Assume that the material follows the von Mises yield criterion with \sigma_{max} = 750 MPa. Find the principal strains and the principal strains directions, the principal stresses and their directions, and comment on whether this state of strain is possible before yielding of the material.

Solution:

The principal strains (eigenvalues) are

    \[ \varepsilon_1=0.004\qquad \varepsilon_2= -0.0033\qquad \varepsilon_3= -0.00026  \]

and the corresponding principal directions (eigenvectors) of the strains are:

    \[ ev_{\varepsilon_1}=\left(\begin{array}{c}0.7324\\ 0.1964\\ 0.6519\end{array}\right)\qquad  ev_{\varepsilon_2}=\left(\begin{array}{c}0.1282\\-0.9801\\0.1512\end{array}\right)\qquad  ev_{\varepsilon_3}=\left(\begin{array}{c}-0.6687\\ 0.02718\\ 0.7431\end{array}\right) \]

To apply the linear elastic constitutive relations, the strain matrix will be converted into a vector of strains, as follows:

    \[ \varepsilon=\left(\begin{array}{c}\varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23}\end{array}\right)=\left(\begin{array}{c}0.002\\-0.003\\0.0015\\0.002\\0.004\\0.002\end{array}\right) \]

Using Equation 7, the corresponding stress vector is:

    \[ \sigma=\left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right)=\left(\begin{array}{c}383.66\\-424.04\\302.89\\161.54\\323.08\\161.54\end{array}\right)MPa \]

The corresponding stress matrix is:

    \[ \sigma=\left(\begin{matrix}383.66&161.54&323.08\\161.54&-424.04&161.54\\323.08&161.54&302.89\end{matrix}\right)MPa \]

The principal stresses (eigenvalues) of the stress matrix are:

    \[\sigma_1=714.53MPa \qqquad \sigma_2=-470.1MPa \qquad \sigma_3=18.064MPa \]

the principal directions (the corresponding eigenvalues of the stress matrix) are:

    \[ ev_{\sigma_1}=\left(\begin{array}{c}0.7324\\ 0.1964\\ 0.6519\end{array}\right)\qquad  ev_{\sigma_2}=\left(\begin{array}{c}0.1282\\-0.9801\\0.1512\end{array}\right)\qquad  ev_{\sigma_3}=\left(\begin{array}{c}-0.6687\\ 0.02718\\ 0.7431\end{array}\right) \]

Since the material is an isotropic linear elastic material, the eigenvectors of the stress and the strain matrix coincide. The maximum principal stress is 714.53 MPa; however, the von Mises stress is equal to 1031.19 MPa, which is higher than \sigma_{max}=750 MPa. This indicates that the material cannot sustain the stress matrix above, and in fact, the strain state above is impossible to occur before yielding or failure of the steel metal.

View Mathematica Code
str={{0.002,0.001,0.002},{0.001,-0.003,0.001},{0.002,0.001,0.0015}};
a=Eigensystem[str];
a//MatrixForm  
str//MatrixForm  
Ee=210000;
Nu=0.3;
G=Ee/2/(1+Nu);
Strainvector={str[[1,1]],str[[2,2]],str[[3,3]],2*str[[1,2]],2*str[[1,3]],2*str[[2,3]]};
Cc={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
Dd=FullSimplify[Inverse[Cc]];
stressvector=Dd.Strainvector
Stressmatrix={{stressvector[[1]],stressvector[[4]],stressvector[[5]]},{stressvector[[4]],stressvector[[2]],stressvector[[6]]},{stressvector[[5]],stressvector[[6]],stressvector[[3]]}};
Stressmatrix//MatrixForm  
b=Eigensystem[Stressmatrix];
b//MatrixForm  
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-
sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
VonMises[Stressmatrix]

 

Example 2:

In a confined compression test of a linear elastic isotropic material, the confining stress was -500MPa while the vertical stress was -100 MPa. The axial strains measured along the directions of the vectors e_1, e_2, and e_3 were -0.003, 0.001, and -0.003, respectively. Find Young’s modulus and Poisson’s ratio.
Example2

Solution:

Using the constitutive relationships for linear elastic isotropic materials:

    \[ \varepsilon_{11}=-\frac{500}{E}+\nu\frac{100}{E}+\nu\frac{500}{E}=-0.003\Rightarrow -500+600\nu=-0.003E \]

    \[ \varepsilon_{22}=\nu\frac{500}{E}-\frac{100}{E}+\nu\frac{500}{E}=0.001\Rightarrow -100+1000\nu=0.001E \]

Solving the two equations:

    \[ \nu=0.222\qquad E=122,222.00MPa \]

 

Example 3:

Consider a transversely isotropic linear elastic piece of material with the isotropy plane being the plane of e_1 and e_2. The values of Young’s modulus and Poisson’s ratio in the plane of isotropy are equal to 100 MPa and 0.2, respectively. In a uniaxial tension test, when the material was stretched in the direction of e_1 by a stress of a value 1 MPa the strain measured in the direction of the vector e_3 was recorded to be –0.0005. In another uniaxial tension test, when the material was stretched in the direction of e_3 by a stress of a value 1 MPa, the strain measured in the direction of e_3 was 0.002. Find an estimate for the material constants E_{33}, \nu_{13}, and \nu_{31}.

Solution:

The Young’s modulus in the direction of e_3 can be obtained from the second experiment:

    \[ \varepsilon_{33}=\frac{\sigma_{33}}{E_{33}}\Rightarrow E_{33}=\frac{1}{0.002}=500MPa \]

From the first experiment, the value of Poisson’s ratio \nu_{13} can be obtained as follows:

    \[\varepsilon_{33}=-\nu_{13}\frac{\sigma_{11}}{E_{11}}\Rightarrow \nu_{13}=\frac{0.0005\times 100}{1}=0.05 \]

The value of Poisson’s ratio \nu_{31} can be obtained from the symmetry of the coefficients matrix:

    \[ \frac{\nu_{13}}{E_{11}}=\frac{\nu_{31}}{E_{33}}\Rightarrow \nu_{31}=\frac{0.05\times 500}{100}=0.25 \]

 

Example 4:

Consider an orthotropic material whose major axes of orthotropy lie in the directions of the vectors:

    \[ e'_1=\left(\begin{array}{c}\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}\\0\end{array}\right)\qquad  e'_2=\left(\begin{array}{c}-\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}\\0\end{array}\right) \qquad e'_3=\left(\begin{array}{c}0\\0\\1\end{array}\right) \]

The material properties in the coordinate systems defined by the orthotropy axes are: E_{11}=100MPa, E_{22}=200MPa, E_{33}=350MPa, \nu_{31}=0.2, \nu_{32}=0.1, \nu_{21}=0.05, G_{12}=50MPa, G_{13}=70MPa, and G_{23}=90MPa. Consider the coordinate system defined by the basis set B=\{e_1,e_2,e_3\} where

    \[ e_1=\left(\begin{array}{c}1\\0\\0\end{array}\right)\qquad  e_2=\left(\begin{array}{c}0\\1\\0\end{array}\right) \qquad e_3=\left(\begin{array}{c}0\\0\\1\end{array}\right) \]

If the material is in a state of uniaxial stress such that the only nonzero stress component is \sigma_{11}=1MPa, find the engineering strain matrix in the coordinate systems defined by the basis sets B and B'=\{e'_1,e'_2, e'_3\}.

Solution:

Since the material is orthotropic, the coefficients matrix C depends on the coordinate system involved. The values of the material coefficients are given in the coordinate system defined by B'. Therefore, it is best to apply the constitutive relationship in the coordinate system of B'.
The stress matrix in the coordinate system of B is given by:

    \[ \sigma=\left(\begin{matrix}1&0&0\\0&0&0\\0&0&0\end{matrix}\right) \]

The rotation matrix required for the coordinate transformation is:

    \[ Q=\left(\begin{matrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\0&0&1\end{matrix}\right) \]

The stress matrix in the coordinate system defined by B' is given by:

    \[ \sigma'=Q\sigma Q^T=\left(\begin{matrix}\frac{1}{2}&-\frac{1}{2}&0\\-\frac{1}{2}&\frac{1}{2}&0\\0&0&0\end{matrix}\right) \]

The relationship between the vector representations of the stress and the strain in the coordinate system of B' can be used to find the strains in the coordinate system of B' (See Equation 2) as follows:

    \[ \left(\begin{array}{c} \varepsilon'_{11}\\\varepsilon'_{22}\\\varepsilon'_{33}\\2\varepsilon'_{12}\\2\varepsilon'_{13}\\2\varepsilon'_{23} \end{array}\right)= \left( \begin{matrix} \frac{1}{100} & -\frac{0.05}{200}} & -\frac{0.2}{350} & 0 &0 & 0\\ -\frac{0.05}{200}& \frac{1}{200} & -\frac{0.1}{350} & 0&0& 0 \\ -\frac{0.2}{350}& -\frac{0.1}{350}& \frac{1}{350} & 0&0& 0 \\ 0 &0 &0 & \frac{1}{50}&0& 0 \\ 0 & 0& 0 &0 & \frac{1}{70}& 0 \\ 0&0 &0 & 0& 0& \frac{1}{90} \end{matrix} \right) \left(\begin{array}{c} \sigma'_{11}\\\sigma'_{22}\\\sigma'_{33}\\\sigma'_{12}\\\sigma'_{13}\\\sigma'_{23}\end{array}\right)=\left(\begin{array}{c}0.00488\\0.00238\\-0.00043\\-0.01\\0\\0\end{array}\right) \]

Therefore, the strain matrix in the coordinate system of B' is given by:

    \[ \varepsilon'=\left(\begin{matrix}0.00488&-0.005&0\\-0.005&0.00238&0\\0&0&-0.00043\end{matrix}\right) \]

On the other hand, the strain matrix \varepsilon in the coordinate system of B is given by:

    \[ \varepsilon'=Q\varepsilon Q^T\Rightarrow \varepsilon=Q^T\varepsilon' Q=\left(\begin{matrix}0.00863&0.00125&0\\0.00125&-0.001375&0\\0&0&-0.00043\end{matrix}\right) \]

It is very important to notice that, while a uniaxial state of stress was applied in the coordinate system defined by B, shear strain was developed in the same coordinate system. This is because the uniaxial state of stress was applied in a coordinate system different from the orthotropy axes coordinate system!
View Mathematica Code

E11=100;E22=200;E33=350;G12=50;G13=70;G23=90;
nu31=0.2;nu32=0.1;nu21=0.05;
nu13=nu31/E33*E11;nu12=nu21/E22*E11;nu23=nu32/E33*E22;
Ss={{1/E11,-nu21/E22,-nu31/E33,0,0,0} ,{-nu12/E11,1/E22,-nu32/E33,0,0,0} ,{-nu13/E11,-nu23/E22,1/E22,0,0,0},{0,0,0,1/G12,0,0},{0,0,0,0,1/G13,0},
{0,0,0,0,0,1/G23}};
s={{1,0,0},{0,0,0},{0,0,0}};
Q=RotationMatrix[-45Degree,{0,0,1}]  
sdash=Q.s.Transpose[Q];
sdash//MatrixForm  
mat=sdash  
sigmavector={mat[[1,1]],mat[[2,2]],mat[[3,3]],mat[[1,2]],mat[[1,3]],mat[[2,3]]}
sigmavector//MatrixForm
strvector=Ss.sigmavector;
strvector//MatrixForm
strdash={{strvector[[1]],strvector[[4]]/2,strvector[[5]]/2},{strvector[[4]]/2,strvector[[2]],strvector[[6]]/2},{strvector[[5]]/2,strvector[[6]]/2,strvector[[3]]}} ;
strdash//MatrixForm
str=Transpose[Q].strdash.Q;
str//MatrixForm
 

Example 5:

The in-plane strain components \varepsilon_{11}, \varepsilon_{22}, and \gamma_{12} of a specimen are 0.001, -0.002, and 0.003, respectively. If Young’s modulus is 10,000MPa and Poisson’s ratio is 0.2, find the stress and strain matrices if the material is in the state of plane strain and if the material is in the state of plane stress.

Solution:

For the state of plane strain, the strain matrix has the following form:

    \[ \varepsilon=\left(\begin{matrix}0.001&0.0015&0\\0.0015&-0.002&0\\0&0&0\end{matrix}\right) \]

Equation 11 can be used to find the stress components \sigma_{11}, \sigma_{22}, and \sigma_{12} as follows:

    \[ \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right)=\left(\begin{matrix}11111.1&2777.78 & 0\\2777.78&11111.1&0\\0&0&4166.67\end{matrix}\right)\left(\begin{array}{c}0.001\\-0.002\\0.003\end{array}\right)=\left(\begin{array}{c}5.56\\-19.44\\12.5\end{array}\right)MPa \]

In a plane strain condition, the movement in the direction of e_3 is restrained, and thus, a restraining stress component \sigma_{33} develops according to the relationship:

    \[ \sigma_{33}=\nu(\sigma_{11}+\sigma_{22})=0.2(5.56-19.4)=-2.78MPa \]

Thus, the stress matrix has the following form:

    \[ \sigma=\left(\begin{matrix}5.56&12.5&0\\12.5&-19.4&0\\0&0&-2.78\end{matrix}\right) \]

View Mathematica Code for the State of Plane Strain:
Ee=10000;Nu=0.2;
Cc=Ee/(1+Nu)*{{(1-Nu)/(1-2Nu),Nu/(1-2Nu),0} ,{Nu/(1-2Nu),(1-Nu)/(1-2Nu),0},{0,0,1/2}};
Cc//MatrixForm  
strvector={0.001,-0.002,0.003};
stressvector=Cc.strvector;
stressvector//MatrixForm  
s=Table[0,{i,1,3},{j,1,3}];
s[[1,1]]=stressvector[[1]];
s[[2,2]]=stressvector[[2]];
s[[1,2]]=stressvector[[3]];
s[[2,1]]=s[[1,2]];
s[[3,3]]=Nu (s[[1,1]]+s[[2,2]]);
s//MatrixForm

In a plane stress condition, a nonzero strain component \varepsilon_{33} can develop, and thus, the strain matrix has the form:

    \[ \varepsilon=\left(\begin{matrix}0.001&0.0015&0\\0.0015&-0.002&0\\0&0&\varepsilon_{33}\end{matrix}\right) \]

The relationship between the stresses and the strains in a plane stress condition (Equation 10) can be utilized to find the values of the non-zero stress components, as follows:

    \[ \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right)=\left(\begin{matrix}10416.7&2083.33 & 0\\2083.33&10416.7&0\\0&0&4166.67\end{matrix}\right)\left(\begin{array}{c}0.001\\-0.002\\0.003\end{array}\right)=\left(\begin{array}{c}6.25\\-18.75\\12.5\end{array}\right)MPa \]

Thus, the stress matrix has the following form:

    \[ \sigma=\left(\begin{matrix}6.25&12.5&0\\12.5&-18.75&0\\0&0&0\end{matrix}\right) \]

In a plane stress condition, \varepsilon_{33} develops according to the relationship:

    \[ \varepsilon_{33}=-\frac{\nu}{E}(\sigma_{11}+\sigma_{22})=-\frac{0.2}{10000}(6.25 - 18.75)=0.00025 \]

View Mathematica Code for the State of Plane Stress:
Ee=10000;Nu=0.2;
Cc=Ee/(1+Nu)*{{1/(1-Nu),Nu/(1-Nu),0},{Nu/(1-Nu),1/(1-Nu),0},{0,0,1/2}};
Cc//MatrixForm  
strvector={0.001,-0.002,0.003};
stressvector=Cc.strvector;
stressvector//MatrixForm  
s=Table[0,{i,1,3},{j,1,3}];
s[[1,1]]=stressvector[[1]];
s[[2,2]]=stressvector[[2]];
s[[1,2]]=stressvector[[3]];
s[[2,1]]=s[[1,2]];
s[[3,3]]=0;
s//MatrixForm  
str=Table[0,{i,1,3} ,{j,1,3}];
str[[1,1]]=strvector[[1]];
str[[2,2]]=strvector[[2]];
str[[1,2]]=str[[2,1]]=strvector[[3]]/2;
str[[3,3]]=-Nu/Ee*(s[[1,1]]+s[[2,2]]);
str//MatrixForm
 

Example 6:

The position function in dimensions of m of a plate adheres to the following form:

    \[ \begin{split} x_1 & =X_1+0.001X_2+0.0002X_1^2\\ x_2 & =1.001X_2-0.0002X_2^2\\ x_3 & =X_3 \end{split} \]

If the plate has dimensions 2m in the direction of e_1, and 1m in the direction of e_2, with the bottom left corner coinciding with the origin of the coordinate system, then:

  • Show that this is a state of plane strain.
  • Find the engineering small strain matrix, and the Green strain matrix and show that the engineering small strain is adequate to describe the motion.
  • Draw the vector plot of the displacement function.
  • If Young’s modulus and Poisson’s ratio are equal to 210 GPa and 0.3, respectively, then find the equilibrium body forces to which this plate is subjected.
  • Draw the contour plot of the von Mises stress and the normal out-of-plane stress component.
Solution:

Unlike the previous examples, this example shows a nonlinear displacement function resulting in a strain field that varies according to position. I.e., every material point will have a different local strain matrix. The displacement function for this plate has the following form:

    \[ u=x-X=\left(\begin{array}{c}0.0002X_1^2+0.001X_2\\0.001X_2-0.0002X_2^2\\0\end{array}\right) \]

The following figures shows the vector plot of the displacement function as drawn by Mathematica:
Example5
The gradient of the displacement tensor is:

    \[ \nabla u=\left(\begin{matrix}0.0004X_1 & 0.001 & 0\\ 0 & 0.001-0.0004X_2 & 0\\0&0&0\end{matrix}\right) \]

The engineering and the Green strain matrices are:

    \[ \begin{split} \varepsilon_{eng} & = \frac{1}{2}\left(\nabla u + \nabla u^T\right)=\left(\begin{matrix}0.0004X_1&0.0005&0\\0.0005&0.001-0.0004X_2&0\\0&0&0\end{matrix}\right)\\ \varepsilon_{Green} & =\frac{1}{2}\left(\nabla u + \nabla u^T+\nabla u^T\nabla u\right)=\left(\begin{matrix}0.0004X_1+8\times 10^{-8}X_1^2&0.0005+2\times 10^{-7}X_1&0\\0.0005+2\times 10^{-7}X_1&0.001-0.0004X_2+8\times 10^{-8}X_2^2&0\\0&0&0\end{matrix}\right) \end{split} \]

The difference between the engineering and the Green strain matrices is very small and can be ignored. Thus, the engineering strain can be used to adequately describe the motion (i.e., the displacement is not associated with large rotations). In addition, the displacement component u_3=0 and the other displacement components are independent of the coordinate X_3. Therefore, the strain components \varepsilon_{33}=\varepsilon_{13}=\varepsilon_{23}=0, i.e., a state of plane strain.

In order to calculate the stresses, the vector representation of the strain will be used:

    \[ \varepsilon=\left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)=\left(\begin{array}{c} 0.0004X_1\\0.001-0.0004X_2\\0\\0.001\\0\\0 \end{array}\right) \]

Using the constitutive equation for the linear elastic isotropic material (Equation 7), the vector representation of the stress is:

    \[ \sigma=\left(\begin{array}{c} 121.154+113.077X_1-48.4615X_2\\282.692+48.4615X_1-113.077X_2\\121.154+48.4615X_1+48.4615X_2\\80.7692\\0\\0 \end{array}\right) MPa \]

Since the displacement is small, the coordinates x_i can be used instead of the coordinates X_i to find the equilibrium body forces:

    \[\begin{split} \rho b_1 & =-\sum_{j=1}^3\frac{\partial \sigma_{j1}}{\partial x_j}=-\frac{\partial \sigma_{11}}{\partial x_1}-\frac{\partial \sigma_{21}}{\partial x_2}-\frac{\partial \sigma_{31}}{\partial x_3}=-113.077 MN/m^3\\ \rho b_2 & =-\sum_{j=1}^3\frac{\partial \sigma_{j2}}{\partial x_j}=-\frac{\partial \sigma_{12}}{\partial x_1}-\frac{\partial \sigma_{22}}{\partial x_2}-\frac{\partial \sigma_{32}}{\partial x_3}=113.077 MN/m^3\\ \rho b_3 & =-\sum_{j=1}^3\frac{\partial \sigma_{j3}}{\partial x_j}=-\frac{\partial \sigma_{13}}{\partial x_1}-\frac{\partial \sigma_{23}}{\partial x_2}-\frac{\partial \sigma_{33}}{\partial x_3}=0 MN/m^3 \end{split} \]

The von Mises stress:

    \[ \sigma_{vM}=\sqrt{45665.7+4175.15(-2.5+X_1)X_1 - 20875.7X_2+4175.15X_1X_2+4175.15X_2^2} \]

The contour plot of the von Mises stress shows that the areas with the maximum von Mises stresses are the right and left bottom corners of the plate. Compare with the vector plot of the displacement (shown above) and note that the maximum displacements are not necessarily correlated with the maximum strain or stress. Since the plate is under a state of plane strain, the restraining stress component \sigma_{33} is non-zero. The contour plot of \sigma_{33} is shown below as well (Units are in MPa).

Example5bcolour

View the Mathematica code
Clear[Ee,Nu]  
x1=X1+0.001*X2+0.0002*X1^2;
x2=1.001*X2-0.0002*X2^2;
x3=X3;
X={X1,X2,X3};
x={x1,x2,x3};
u=x-X;
u//MatrixForm  
uplane=Table[u[[i]],{i,1,2}]  
VectorPlot[uplane,{X1,0,2},{X2,0,1},AspectRatio->Automatic,PlotLabel->”u”,VectorStyle->Black]  
gradu=Table[D[u[[i]],X[[j]]],{i,1,3},{j,1,3}];
gradu//MatrixForm  
strGreen=FullSimplify[1/2*(gradu+Transpose[gradu]+Transpose[gradu].gradu)];
Expand[strGreen]//MatrixForm  
str=FullSimplify[1/2*(gradu+Transpose[gradu])];
str//MatrixForm  
Strainvector={str[[1,1]],str[[2,2]],str[[3,3]],2*str[[1,2]],2*str[[1,3]],2*str[[2,3]]};
Strainvector//MatrixForm  
Ee=210000;
Nu=0.3;
G=Ee/2/(1+Nu);
Cc={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
Dd=FullSimplify[Inverse[Cc]];
stressvector=FullSimplify[Chop[Dd.Strainvector]];
stressvector//MatrixForm  
Stressmatrix={{stressvector[[1]],stressvector[[4]],stressvector[[5]]},{stressvector[[4]],stressvector[[2]],stressvector[[6]]},{stressvector[[5]],stressvector[[6]],stressvector[[3]]}};
Stressmatrix=FullSimplify[Chop[Stressmatrix]];
Stressmatrix//MatrixForm  
X={X1,X2,X3};
s=Stressmatrix;
-Sum[D[s[[i,1]],X[[i]]],{i,1,3}]  
-Sum[D[s[[i,2]],X[[i]]],{i,1,3}]  
-Sum[D[s[[i,3]],X[[i]]],{i,1,3}]  
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-
sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
a=FullSimplify[Chop[VonMises[Stressmatrix]]]  
ContourPlot[a,{X1,0,2},{X2,0,1},AspectRatio->Automatic,ContourLabels->All,PlotLabel->”vonMises”,ColorFunction->”GrayTones”] ContourPlot[Stressmatrix[[3,3]],{X1,0,2},{X2,0,1},AspectRatio->Automatic,ContourLabels->All,PlotLabel->”s33″,ColorFunction->”GrayTones”]
 

Example 7:

The in plane displacements of a plate of dimensions of 5m and 1m in the directions of e_1 and e_2, respectively, with the bottom left corner situated at the origin are according to the relationship:

    \[\begin{split} u_1 & = 0.02X1 + 0.02X_1X_2\\ u_2 & = 0.0015X_2X_1^2+0.02X_2 \end{split} \]

If the plate is in the state of plane stress, draw the contour plot of the change of thickness of the plate, assuming that Young’s modulus is 1000 MPa, Poisson’s ratio is 0.4, and the thickness of the plate is 1mm.

Solution:

Since the plate is in a state of plane stress, \sigma_{33}=  0 while a strain compoenent \varepsilon_{33} develops according to the equation:

    \[ \varepsilon_{33}=-\frac{\nu}{E}(\sigma_{11}+\sigma_{22}) \]

The in plane strain components can be obtained from the in plane gradient of the displacement function, as follows:

    \[ \nabla u=\left(\begin{matrix} 0.02+0.02X_2 & 0.02X_1\\ 0.003X_1X_2 & 0.02+0.0015X_1^2\end{matrix}\right) \]

    \[ \left(\begin{matrix} \varepsilon_{11} & \varepsilon_{12}\\ \varepsilon_{21} & \varepsilon_{22}\end{matrix}\right)=\frac{1}{2}(\nabla u + \nabla u^T)=\left(\begin{matrix} 0.02+0.02X_2 & (0.01+0.0015X_2)X_1\\ (0.01+0.0015X_2)X_1 & 0.02+0.0015X_1^2\end{matrix}\right) \]

The stresses can be obtained using the plane stress constitutive relation (Equation 10):

    \[ \left(\begin{array}{c}\sigma_{11}\\\sigma_{22}\\\sigma_{12}\end{array}\right)=\left(\begin{array}{c}33.33+0.714X_1^2+23.81X_2\\33.33+1.786X_1^2+9.524X_2\\X_1(7.143+1.07X_2)\end{array}\right)MPa \]

Thus, the strain component \varepsilon_{33} is:

    \[ \varepsilon_{33}=-\frac{\nu}{E}(\sigma_{11}+\sigma_{22})=-0.027-0.001X_1^2-0.013X_2 \]

The change in thickness in mm can be obtained using the formula:

    \[ \Delta = \varepsilon_{33}t \]

The required contour plot is shown below:
Example7
View Mathematica Code

X={X1,X2,X3};
u={0.02X1+0.02X1*X2,0.0015X2*X1^2+0.02X2};
VectorPlot[u,{X1,0,5},{X2,0,1},AspectRatio->Automatic]  
gradu=Table[D[u[[i]],X[[j]]],{i,1,2},{j,1,2}];
gradu//MatrixForm  
esmall=1/2*(gradu+Transpose[gradu]);
esmall//MatrixForm  
strvector={esmall[[1,1]],esmall[[2,2]],2*esmall[[1,2]]};
strvector//MatrixForm  
Ee=1000;Nu=0.4;Cc=Ee/(1+Nu)*{{1/(1-Nu),Nu/(1-Nu),0},{Nu/(1-Nu),1/(1-Nu),0},{0,0,1/2}};
stressvector=Cc.strvector;
stressvector//MatrixForm  
strain=Table[0,{i,1,3},{j,1,3}];
strain[[1;;2,1;;2]]=esmall;
strain[[3,3]]=-Nu/Ee*(stressvector[[1]]+stressvector[[2]]);
strain//MatrixForm  
thickness=1;
delta=FullSimplify[strain[[3,3]]*thickness]  
ContourPlot[delta,{X1,0,5},{X2,0,1},ContourLabels->All,AspectRatio->Automatic,PlotLabel->”delta Thickness mm”]
 

Example 8:

A cylindrical pressure vessel is made out of a steel material. The material’s Young modulus and Poisson’s ratio are 210 GPa and 0.3 respectively. The pressure vessel is used to hold a gas with an internal pressure P of 3 MPa. The vessel’s average diameter is 2 m, and its thickness is 10mm. Find the von Mises stress, the change in diameter, and the change in thickness of the cylinder in the following two conditions:

  • Condition 1: Capped ends.
  • Condition 2: Restrained ends.

Example8a

Example8b

Solution:

As shown in the sketch, the orthonormal basis vectors are chosen in this example such that e_1 is aligned with the longitudinal axis of the pipe and e_2 is vertical. The analyzed part of the pressure wall is taken to be
perpendicular to e_3.
For condition 1: The normal stress component in the longitudinal direction is:

    \[ \sigma_{11}=\frac{P\pi r^2}{2 \pi r t}=\frac{Pr}{2t}=150 MPa \]

The normal stress component in the direction of the vector e_2 is:

    \[ \sigma_{22}=\frac{Pr}{t}=300 MPa \]

All the remaining components of the stress matrix are zero. Thus, the stress matrix has the following form:

    \[ \sigma=\left(\begin{matrix}150 & 0&0\\0&300&0\\0&0&0\end{matrix}\right)MPa \]

The von Mises stress is:

    \[ \sigma_{vM}=150\sqrt{3}=259.81MPa \]

The strain-stress relationship (Equation 6) can be used to find the vector representation of the strain:

    \[ \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)=\left(\begin{array}{c} 0.000286\\0.001214\\-0.00064\\0\\0\\0\end{array}\right) \]

The strain component \varepsilon_{22} gives the relative change of diameter as follows:

    \[ \varepsilon_{22}=\frac{2\pi (r+\Delta r)-2 \pi r}{2\pi r}=\frac{\Delta r}{r}=\frac{\Delta D}{D}=0.001214 \]

Therefore the change in diameter is:

    \[ \Delta D=0.001214D=0.001214\times 2=0.0024 m. = 2.4mm \]

The component \varepsilon_{33} provides the relative change in the thickness of the pressure vessel. Therefore:

    \[ \Delta t=\varepsilon_{33}t=-0.00064 (10)=-0.0064mm \]

View Mathematica code for Condition 1:
P=3;r=1000;t=10;Ee=210000;Nu=0.3;
s=Table[0,{i,1,3},{j,1,3}];
s[[1,1]]=P*r/2/t  
s[[2,2]]=P*r/t  
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-
sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
a=VonMises[s]  
G=Ee/2/(1+Nu);
Ss={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
stressvector={s[[1,1]],s[[2,2]],s[[3,3]],s[[1,2]],s[[1,3]],s[[2,3]]}  
strainvector=Ss.stressvector;
strainvector//MatrixForm  
deltadiameter=strainvector[[2]]*2*r  
deltathickness=strainvector[[3]]*t  

For condition 2, the normal stress in the longitudinal direction (\sigma_{11}) is not known since the vessel is restrained in the longitudinal direction. The circumferential stress (normal stress in the direction of the vector e_2) is given by:

    \[ \sigma_{22}=\frac{Pr}{t}=300MPa \]

Thus, the stress matrix has the form:

    \[ \sigma=\left(\begin{matrix}\sigma_{11} & 0&0\\0&300&0\\0&0&0\end{matrix}\right)MPa \]

Since for condition 2, the pipe is restrained from expanding or contracting, therefore, \varepsilon_{11}=0. Therefore:

    \[ 0=\varepsilon_{11}=\frac{\sigma_{11}}{E}-\frac{\nu}{E}\left(\sigma_{22}+\sigma_{33}\right)\Rightarrow \sigma_{11}=0.3\times 300 = 90 MPa \]

The von Mises stress is thus given by:

    \[ \sigma_{vM}=353.7 MPa \]

The strain-stress relationship (Equation 6) can be used to find the vector representation of the strain:

    \[ \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)=\left(\begin{array}{c} 0\\0.0013\\-0.00056\\0\\0\\0\end{array}\right) \]

The strain component \varepsilon_{22} gives the relative change of diameter as follows:

    \[ \varepsilon_{22}=\frac{2\pi (r+\Delta r)-2 \pi r}{2\pi r}=\frac{\Delta r}{r}=\frac{\Delta D}{D}=0.0013 \]

Therefore the change in diameter is:

    \[ \Delta D=0.0013D=0.0013\times 2=0.0026 m. = 2.6mm \]

The component \varepsilon_{33} provides the relative change in the thickness of the pressure vessel. Therefore:

    \[ \Delta t=\varepsilon_{33}t=-0.00056 (10)=-0.0056mm \]

View Mathematica code for Condition 2:
P=3;r=1000;t=10;Ee=210000;Nu=0.3;
s=Table[0,{i,1,3},{j,1,3}];
s[[2,2]]=P*r/t  
s[[1,1]]=Nu*s[[2,2]]  
VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];
a=VonMises[s]  
G=Ee/2/(1+Nu);
Ss={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};
stressvector={s[[1,1]],s[[2,2]],s[[3,3]],s[[1,2]],s[[1,3]],s[[2,3]]}  
strainvector=Ss.stressvector;
strainvector//MatrixForm  
deltadiameter=strainvector[[2]]*2*r  
deltathickness=strainvector[[3]]*t  

 

Problems:

  1. What does a negative Poisson’s ratio mean? Conduct a literature search and write a small paragraph supporting or opposing that materials with negative Poisson’s ratio exist (cite any references used).
  2. What does a Poisson’s ratio>0.5 mean? Conduct a literature search and write a small paragraph supporting or opposing that materials with Poisson’s ratio>0.5 exist (cite any references used).
  3. Assuming a linear elastic material such that the relationship between the stress and the strain is described by the following matrix:

        \[ \left(\begin{array}{c} \varepsilon_{11}\\\varepsilon_{22}\\\varepsilon_{33}\\2\varepsilon_{12}\\2\varepsilon_{13}\\2\varepsilon_{23} \end{array}\right)= \left( \begin{matrix} \frac{1}{E} & -\frac{\nu}{E} & -\frac{\nu}{E} & 0 &0 & 0\\ -\frac{\nu}{E}}& \frac{1}{E} & -\frac{\nu}{E} & 0&0& 0 \\ -\frac{\nu}{E}& -\frac{\nu}{E}& \frac{1}{E} & 0&0& 0 \\ 0 &0 &0 & \frac{1}{G}&0& 0 \\ 0 & 0 & 0& 0 & \frac{1}{G}& 0 \\ 0 & 0 &0 & 0 & 0 & \frac{1}{G} \end{matrix} \right) \left(\begin{array}{c} \sigma_{11}\\\sigma_{22}\\\sigma_{33}\\\sigma_{12}\\\sigma_{13}\\\sigma_{23}\end{array}\right) \]

    Assuming that the material is isotropic (i.e., the relationship above applies in any coordinate system) and starting from a state of stress described by the matrix:

        \[ \sigma=\left(\begin{matrix}\sigma_{11}&0&0\\0&0&0\\0&0&0\end{matrix}\right) \]

    show that

        \[ G=\frac{E}{2(1+\nu)} \]

    (Hint: use a procedure similar to the one shown below the definition of the shear modulus).

  4. Show the following relationships for linear elastic isotropic materials:

        \[ \begin{split} S_{ij}=2\mu\varepsilon'_{ij}\\ \sum_{k=1}^3\frac{\sigma_{kk}}{3}=(\frac{2}{3}\mu+\lambda)\sum_{m=1}^3\varepsilon_{mm}\\ \varepsilon_{ij}=\frac{\sigma_{ij}}{2\mu}-\frac{\lambda}{2\mu(2\mu+3\lambda)}\sum_{k=1}^3\sigma_{kk}\delta_{ij} \end{split} \]

    where:
    S_{ij} are the components of the deviatoric stress tensor  
    \varepsilon'_{ij} are the components of the deviatoric strain tensor \varepsilon'_{ij}=\varepsilon_{ij}-\sum_{k=1}^3\frac{\varepsilon_{kk}}{3}\delta_{ij}  
    \lambda and \mu are Lamé’s constants.

  5. A state of stress of a linear isotropic material is represented by a positive definite symmetric stress matrix. Find the volumetric strain in terms of the principal stresses, Young’s modulus and Poisson’s ratio. Comment on whether the volume increases, decreases, stays constant or whether the information provided is not enough to answer this question.
  6. In a confined compression test of a linear elastic isotropic material, the testing apparatus was such that the confining pressure was equal to 50% of the applied vertical compressive stress P. Find the relationship between P and the volumetric strain in terms of Young’s modulus and Poisson’s ratio. Also, find the maximum shear stress, and the von Mises stress in terms of P.
  7. In a confined compression test of a linear elastic isotropic material, the testing apparatus was such that the specimen was restrained from expanding (or contracting) laterally. Find the relationship between the vertical compressive stress and the vertical strain in terms of Young’s modulus and Poisson’s ratio.
  8. Find the value of Poisson’s ratio at which both the shear modulus and the bulk modulus are equal to each other in value.
  9. A cylindrical sample is subjected to a confined compression test in which the horizontal stress is kept constant at –200 MPa, while the magnitude of the compressive vertical stress is increased to reach the value of –600 MPa. Find the volumetric strain if the shear and bulk moduli are 22.72 GPa and 20.833 GPa, respectively. (Answer: –1.6%) .
  10. A linear elastic isotropic material with Young’s modulus of 210GPa and Poisson’s ratio of 0.3 is in a state of strain described by the infinitesimal strain matrix:

        \[ \varepsilon=\left(\begin{matrix}0.03 & 0.002 & 0\\ 0.002 & -0.015 & 0\\0&0&0\end{matrix}\right) \]

    • Find the principal strains and their directions.
    • Find the principal stresses and their directions.
    • Are the principal stresses and the principal strains aligned?
  11. Consider a linear elastic orthotropic material in a state of strain described by the infinitesimal strain matrix as in the previous problem. Assume the following material properties:
    E_{11}=200GPa, E_{22}=100GPa, E_{33}=50GPa, \nu_{12}=\nu_{13}=\nu_{23}=0.3, G_{12}=G_{13}=G_{23}=50GPa.

    • Find the principal stresses and their directions.
    • Are the principal stresses and the principal strains aligned?
  12. Using sound mathematical arguments show that the principal stresses and the principal strains are always aligned for linear elastic isotropic materials. (Hint: use the fact that Equations 6 and 7 are valid in any coordinate system)
  13. The in-plane displacement function of the two dimensional plate shown below is given by:

        \[ u_1=(-0.001 X_1^2+0.005X_1)X_2 \]

        \[ u_2=-0.0016X_1^3+0.0006X_1^2 \]

    FigureExample
    Assuming that the plate is of a linear elastic isotropic material with Young’s modulus E=210 GPa and Poisson’s ratio =0.2, draw the contour plot of the stress components \sigma_{11}, \sigma_{22}, \sigma_{12}, \sigma_{33}, \sigma_{vM}, and the hydrostatic stress for each of the following two cases:

    • the plate is in a state of plane stress.
    • the plate is in a state of plane strain.
  14. A cylindrical pressure vessel with an average diameter of 1m and a wall thickiness of 10mm is made of a material with Young’s modulus of 2.1 GPa and Poisson’s ratio 0.3. A constant pressure P of 1.3MPa is maintained inside the vessel
    • Assuming that the length of the pressure vessel is 5m, find the increase in length. (Answer: 31mm)
    • Find the increase in the diameter of the pressure vessel. (Answer: 26.3mm)
    • Find the change in thickness of the pressure vessel. (Answer: –0.14mm)
    • Find the volumetric strain throughout the pressure vessel wall. (Answer: 1.9%)
  15. An elastic cylindrical pressure vessel with average diameter d, thickness t, length L, Young’s modulus E, Poisson’s ratio  \nu is subjected to an internal pressure P. Assume that the thickness is very small compared to the diameter. Find the circumferential stress, the longitudinal stress, the circumferential strain, and the longitudinal strain in the following two conditions:
    • Condition 1: The pressure vessel has two caps on its ends.
    • Condition 2: The pressure vessel is infinitely long and, thus, cannot expand or contract in the longitudinal direction.

Derivation of the Reduced Number of Independent Matrix Coefficients Due to Material Symmetries:

The following Mathematica code is used to reduce the number of independent coefficients in the fourth-order tensor. In the first part, the symmetry of the stress tensor and the symmetry of the strain tensor are utilized to reduce the number of independent coefficients from 81 to 36. Then, the preservation of energy concept is utilized to reduce the number of independent coefficients from 36 to 21. Next, reflections along the directions of e_1, e_2, and e_3 are applied to find the number of coefficients required for an orthotropic material, with the planes of symmetry being the planes perpendicular to the three basis vectors. In this case, the number of independent constants is shown to be 9. Then, the number of independent constants in a transversely isotropic material is found by utilizing the orthotropic material fourth-order tensor and assuming that the behaviour of the material is independent of a 45^\circ rotation in the plane of e_1 and e_2. The resulting number of independent material constants is shown to be 5. Finally, another 45^\circ rotation around e_1 is used to show that the number of independent material constants for an isotropic material is just 2. Copy and paste the code into your Mathematica browser, and follow it step by step. You should also show that the results for the transversely isotropic and isotropic materials are independent of the angle of rotation used in the derivation, as long as the angle is not a multiple of \frac{\pi}{2}.

View Mathematica Code
(*General Ctensor (fourth order tensor)*)
Ctensor=Table[Row[{c,i,j,k,l}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“A general C matrix relating 9 stress tensor components to 9 strain tensor components:”] 
Ctensor//MatrixForm
(*Applying stress symmetry*)
RuleStressSymmetry={};
Do[AppendTo[RuleStressSymmetry,Ctensor[[i,j,k,l]]->Ctensor[[j,i,k,l]]],{i,2,3},{j,1,i-1},{k,1,3},{l,1,3}];
Print[“Stress tensor symmetry leads to the following relations:”,RuleStressSymmetry] 
CtensorSt=Ctensor/.RuleStressSymmetry;
Print[“The resulting C matrix from stress tensor symmetry is:”] 
CtensorSt//MatrixForm
(*Applying Strain symmetry*)
RuleStrainSymmetry={};
Do[AppendTo[RuleStrainSymmetry,CtensorSt[[i,j,k,l]]->CtensorSt[[i,j,l,k]]],{i,1,3},{j,1,3},{k,2,3},{l,1,k-1}];
Print[“Strain tensor symmetry leads to the following relations:”,RuleStrainSymmetry] 
CtensorStStr=CtensorSt/.RuleStrainSymmetry;
Print[“The resulting C matrix from both stress and strain tensors symmetry is:”] 
CtensorStStr//MatrixForm
(*Using symmetry due to Energy Preservation*)
RuleEnergySymmetry={};
Do[AppendTo[RuleEnergySymmetry,CtensorStStr[[i,j,k,l]]->CtensorStStr[[k,l,i,j]]],{i,1,3},{j,i,3},{k,1,i},{l,k,If[k<i,3,j]}];
Print[“Symmetry due to the existence of an energy function leads to the following relations:”,RuleEnergySymmetry] 
CtensorStStrEn=CtensorStStr/.RuleEnergySymmetry;
Print[“The resulting C matrix after using stress, strain and energy symmetries is:”] 
CtensorStStrEn//MatrixForm
Print[“The original C tensor and the fully symmetric Ctensor:”] 
{Ctensor//MatrixForm,CtensorStStrEn//MatrixForm} 
Print[“The difference (minus) between C tensor and the fully symmetric Ctensor:”] 
Ctensor-CtensorStStrEn//MatrixForm
(*Adopting the new Ctensor after symmetries*)
Ctensor=CtensorStStrEn;
(*Order of the stress vector*)
order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};
(*Converting Ctensor to Cmatrix*)
Cmatrix=Table[ToExpression[ToString[Row[{“Ctensor”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];
Print[“The Following is the resulting C matrix after adopting all the symmetries:”] 
Cmatrix//MatrixForm
(*Orthotropic material*)
Print[“Deriving C tensor and C matrix for Orthotropic Materials”] 
(*The following are three reflections around the three planes of symmetry perpendicular to e1=ex and e2=ey and e3=ez*)
Qx=ReflectionMatrix[{1,0,0}];
Qy=ReflectionMatrix[{0,1,0}];
Qz=ReflectionMatrix[{0,0,1}];
(*Using plane perpendicular to x as a symmetry plane*)
Cdashtensor=Table[Sum[Qx[[i,m]]*Qx[[j,n]]*Qx[[k,o]]*Qx[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“The following is the C tensor after reflection around the plane perpendicular to e1=ex:”] 
Cdashtensor//MatrixForm
Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e1=ex and the original symmetric C tensor:”] 
(Ctensor-Cdashtensor)//MatrixForm
zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e1=ex:”] 
ax=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]] 
Ctensorx=Ctensor/.ax[[1]];
Print[“Following is the resulting C tensor after adopting the previous relations:”] 
Ctensorx//MatrixForm
(*Using plane perpendicular to y as a symmetry plane*)
Cdashtensor=Table[Sum[Qy[[i,m]]*Qy[[j,n]]*Qy[[k,o]]*Qy[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“Following is the C tensor after reflection around the plane perpendicular to e2=ey:”] 
Cdashtensor//MatrixForm
Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e2=ey and the original symmetric C tensor:”] 
(Ctensor-Cdashtensor)//MatrixForm
Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e2=ey:”] 
ay=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]] 
Ctensory=Ctensor/.ay[[1]];
Print[“Following is the resulting C tensor after adopting the previous relations:”] 
Ctensory//MatrixForm
(*Using plane perpendicular to z as a symmetry plane*)
Cdashtensor=Table[Sum[Qz[[i,m]]*Qz[[j,n]]*Qz[[k,o]]*Qz[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“Following is the C tensor after reflection around the plane perpendicular to e3=ez:”] 
Cdashtensor//MatrixForm
Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e3=ez and the original symmetric C tensor:”] 
(Ctensor-Cdashtensor)//MatrixForm
Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e3=ez:”] 
az=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]] 
Ctensorz=Ctensor/.az[[1]];
Print[“Following is the resulting C tensor after adopting the previous relations:”] 
Ctensorz//MatrixForm
CtensorOrthotropic=Ctensor/.Flatten[{az[[1]],ay[[1]],ax[[1]]}];
Print[“This is the resulting C tensor after adopting the orthotropic material symmetries:”] 
CtensorOrthotropic//MatrixForm
(*Order of the stress vector*)
order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};
(*Converting Ctensor to Cmatrix*)
CmatrixOrthotropic=Table[ToExpression[ToString[Row[{“CtensorOrthotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];
Print[“Following is the resulting C matrix after adopting the orthotropic material symmetries:”] 
CmatrixOrthotropic//MatrixForm
(*Transverse Isotropy*)
Print[“Deriving C tensor of a Transversely Isotropic material from C matrix of Orthotropic Materials assumine that e1=ex and e2=ey form the plane of isotropy”] 
(*Rotation around z with 45 Degrees*)
Q=RotationMatrix[45Degree,{0,0,1}];
Cdashtensor=Table[Sum[Q[[i,m]]*Q[[j,n]]*Q[[k,o]]*Q[[l,p]]*CtensorOrthotropic[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“Following is the resulting C tensor after rotating the orthotropic material C matrix with some angle in the plane of isotropy:”] 
Cdashtensor=FullSimplify[Cdashtensor];
Cdashtensor//MatrixForm
Print[“Following are the differences between the previous C tensor and the Orthotropic material C tensor:”] 
(CtensorOrthotropic-Cdashtensor)//MatrixForm
zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Vars={};
Do[If[CtensorOrthotropic[[i,j,k,l]]==0,,,AppendTo[Vars,CtensorOrthotropic[[i,j,k,l]]]],{i,3,1,-1},{j,3,1,-1},{k,3,1,-1},{l,3,1,-1}];
Vars;
Print[“Following are the relationships after equating C tensor after rotating in the plane of isotropy to the Ctensor of orthotropic material:”] 
a=Solve[(CtensorOrthotropic-Cdashtensor)==zero,Vars] CtensorTransverseIsotropic=CtensorOrthotropic/.a[[1]];
Print[“Following is the C tensor for a Transversely Isotropic Material:”] 
CtensorTransverseIsotropic=FullSimplify[CtensorTransverseIsotropic];
CtensorTransverseIsotropic//MatrixForm
(*Order of the stress vector*)
order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};
(*Converting Ctensor to Cmatrix*)
CmatrixTransverseIsotropic=Table[ToExpression[ToString[Row[{“CtensorTransverseIsotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];
CmatrixTransverseIsotropic=FullSimplify[CmatrixTransverseIsotropic];
Print[“Following is the C matrix for a Transversely Isotropic Material:”] 
CmatrixTransverseIsotropic//MatrixForm
(*Full Isotropy*)
Print[“Deriving C tensor of an Isotropic material from the previous Transversely Isotropic Materials C tensor”] 
(*Rotation around x with 45 Degrees*)
Q=RotationMatrix[45Degree,{1,0,0}];
Cdashtensor=Table[Sum[Q[[i,m]]*Q[[j,n]]*Q[[k,o]]*Q[[l,p]]*CtensorTransverseIsotropic[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“Following is the resulting C tensor after rotating the Transversely Isotropic material C matrix with some angle around x:”] 
Cdashtensor=FullSimplify[Cdashtensor];
Cdashtensor//MatrixForm
Print[“Following are the differences between the previous C tensor and the Transversely Isotropic material C tensor:”] 
(CtensorTransverseIsotropic-Cdashtensor)//MatrixForm
zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];
Print[“Following are the relationships after equating C tensor after rotating to the Ctensor of transversely isotropic material:”] 
a=Solve[(CtensorTransverseIsotropic-Cdashtensor)==zero,Vars] 
CtensorIsotropic=CtensorTransverseIsotropic/.a[[1]];
Print[“Following is the C tensor for an Isotropic Material:”] 
CtensorIsotropic=FullSimplify[CtensorIsotropic];
CtensorIsotropic//MatrixForm
(*Order of the stress vector*)
order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};
(*Converting Ctensor to Cmatrix*)
CmatrixIsotropic=Table[ToExpression[ToString[Row[{“CtensorIsotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];
CmatrixIsotropic=FullSimplify[CmatrixIsotropic];
Print[“Following is the C matrix for an Isotropic Material:”] 
CmatrixIsotropic//MatrixForm

Page Comments

  1. Zibo says:

    Hi, Samer. Is there a reason why in example 8 there is no sigma33? in other words, no stress in e3 direction?
    Thanks!

    1. Samer Adeeb says:

      The reason is that this is a thin walled pressure pipe, with large D/t. This causes the normal stress in the direction perpendicular to the wall to be negligible compared to the stresses in the plane tangent to the wall

Leave a Reply

Your email address will not be published. Required fields are marked *