TMM4175 Polymer Composites

Home About Python Links Next Previous Table of Contents

Hooke's Law

Simply described for the scope of this course: Hooke's law represents the linear relations between stresses and strains.

Tensor notation

Using the tensor notation, or index notation, we can state the general Hooke's law as

\begin{equation} \varepsilon_{ij} = S_{ijkl}\sigma_{kl} \tag{1} \end{equation}

From the Einstein summation convention this is equivalent to \begin{equation} \varepsilon_{ij} = \sum_{k=1}^{3}\sum_{l=1}^{3}S_{ijkl}\sigma_{kl} \tag{2} \end{equation}

where $S_{ijkl}$ is called the compliance tensor. The stiffness tensor $C_{ijkl}$ is found in the inverse relation,

\begin{equation} \sigma_{ij} = C_{ijkl}\varepsilon_{kl} \tag{3} \end{equation}

Both the stiffness tensor and the compliance tensor are fourth-order tensors that can be represented by arrays of dimensions (3 x 3 x 3 x 3)

Note that for one component of strain, there are 9 terms in the summation:

\begin{equation} \varepsilon_{ij} = \sum_{k=1}^{3}\sum_{l=1}^{3}S_{ijkl}\sigma_{kl} = S_{ij11}\sigma_{11}+S_{ij12}\sigma_{12}+S_{ij13}\sigma_{13}+S_{ij21}\sigma_{21}+S_{ij22}\sigma_{22}+S_{ij23}\sigma_{23}+S_{ij31}\sigma_{31}S_{ij32}\sigma_{32}+S_{ij33}\sigma_{33} \tag{4} \end{equation}

Hence, there are 9 x 9 = 81 numerical values in both the compliance and in the stiffness tensor.

Most of the values are however zero or redundant due to the inherent symmetry required for satisfying the physics of the relation. Implementation of Hooke's law using this notation is straigth forward but not efficient since the data structure for both stresses, strains and compliance or stiffness must allocate recources that do not contribute to anything.

Matrix notation

Recall from the sections about stress and strain that both the stress tensor and the strain tensor is symmetric. Therfore, there are only 6 independent components of stress and 6 independet components of strain, and we may express the tensors in a vector form. The relation between the stress and the strain will now be a 6 x 6 matrix. For a general anisotropic material we will express Hooke's law by engineering components of stress and strain as:

\begin{equation} \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \end{bmatrix} = \begin{bmatrix} S_{11} & S_{12} & S_{13} & S_{14} & S_{15} & S_{16} \\ S_{12} & S_{22} & S_{23} & S_{24} & S_{25} & S_{26} \\ S_{13} & S_{23} & S_{33} & S_{34} & S_{35} & S_{36} \\ S_{14} & S_{24} & S_{34} & S_{44} & S_{45} & S_{46} \\ S_{15} & S_{25} & S_{35} & S_{45} & S_{55} & S_{56} \\ S_{16} & S_{26} & S_{36} & S_{45} & S_{56} & S_{66} \end{bmatrix} \begin{bmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \tau_{23} \\ \tau_{13} \\ \tau_{12} \end{bmatrix} \tag{5} \end{equation}

The short notation is:

\begin{equation} \boldsymbol{\varepsilon}=\mathbf{S}\boldsymbol{\sigma} \tag{6} \end{equation}

Observe that, due to symmetry of the compliance matrix $\mathbf{S}$ (i.e.,$S_{ji}=S_{ij}$), there are only 21 independent elastic constants.

The inverse relation is written

\begin{equation} \begin{bmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \tau_{23} \\ \tau_{13} \\ \tau_{12} \end{bmatrix}= \begin{bmatrix} C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16} \\ C_{12} & C_{22} & C_{23} & C_{24} & C_{25} & C_{26} \\ C_{13} & C_{23} & C_{33} & C_{34} & C_{35} & C_{36} \\ C_{14} & C_{24} & C_{34} & C_{44} & C_{45} & C_{46} \\ C_{15} & C_{25} & C_{35} & C_{45} & C_{55} & C_{56} \\ C_{16} & C_{26} & C_{36} & C_{45} & C_{56} & C_{66} \end{bmatrix} \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \end{bmatrix} \tag{7} \end{equation}

with the short notation \begin{equation}\boldsymbol{\sigma} =\mathbf{C}\boldsymbol{\varepsilon} \tag{8} \end{equation}

The sequence of stresses and strains used in the course is by the incices 1, 2, 3, 23, 13, 12. This is one of two common notations, the other one being the sequence 1, 2, 3, 12, 13, 23.

Orthotropic materials

Orthotropic materials have 3 mutually perpendicular symmetry planes.

Figure-1: Orthotropic and transversely isotropic symmetry

Due to this symmetry, there are no coupling between normal stresses and shear strains, between shear stresses and normal strains, or between a shear stresses and a shear strains on different planes. Hence, the relation takes the form:

\begin{equation} \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \gamma_{23} \\ \gamma_{13} \\ \gamma_{12} \end{bmatrix} = \begin{bmatrix} S_{11} & S_{12} & S_{13} & 0 & 0 & 0 \\ S_{12} & S_{22} & S_{23} & 0 & 0 & 0 \\ S_{13} & S_{23} & S_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & S_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & S_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & S_{66} \end{bmatrix} \begin{bmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \tau_{23} \\ \tau_{13} \\ \tau_{12} \end{bmatrix} \tag{9} \end{equation}

Orthotropic materials have 9 independent elastic constants where the components of the compliance matrix expressed by the engineering constants are:

\begin{equation} \begin{aligned} &S_{11}=\frac{1}{E_1} && S_{22}=\frac{1}{E_2} && S_{33}=\frac{1}{E_3} \\ &S_{12}=-\frac{\nu_{12}}{E_1} && S_{13}=-\frac{\nu_{13}}{E_1} && S_{23}=-\frac{\nu_{23}}{E_2} \\ &S_{44}=\frac{1}{G_{23}} && S_{55}=\frac{1}{G_{13}} && S_{66}=\frac{1}{G_{12}} \end{aligned} \tag{10} \end{equation}

where $E_i$ are elastic moduli, $\nu_{ij}$ are Poisson's ratios, and $G_{ij}$ are shear moduli.

From the symmetry of the compliance matric:

\begin{equation} \frac{\nu_{ij}}{E_i} = \frac{\nu_{ji}}{E_j} \tag{11} \end{equation}

Example

A fictive orthotropic material has the following engineering constants:

In [1]:
import matlib

m1=matlib.newMaterial(E1=150000, E2=20000, E3=10000, v12=0.3, v13=0.4, v23=0.5, G12=5000, G13=4000, G23=3000)

See Material properties and Material Library for explanation of the function newMaterial().

The compliance matrix can conveniently be computed by a function taking the material as an argument:

In [2]:
import numpy as np

def S3D(m):
    return np.array(
    [[        1/m['E1'],-m['v12']/m['E1'],-m['v13']/m['E1'],         0,    0,     0],
     [-m['v12']/m['E1'],        1/m['E2'],-m['v23']/m['E2'],         0,    0,     0],
     [-m['v13']/m['E1'],-m['v23']/m['E2'],        1/m['E3'],         0,    0,     0],
     [                0,                0,                0,1/m['G23'],    0,     0],
     [                0,                0,                0,    0,1/m['G13'],     0],
     [                0,                0,                0,    0,     0,1/m['G12']] ],
    float)

Obtaining the compliance matrix and print formatted output by numpy.array2string():

In [3]:
S=S3D(m1)

print(np.array2string(S, precision=6, suppress_small=True, separator='  ', floatmode='maxprec') )
[[ 0.000007  -0.000002  -0.000003   0.         0.         0.      ]
 [-0.000002   0.00005   -0.000025   0.         0.         0.      ]
 [-0.000003  -0.000025   0.0001     0.         0.         0.      ]
 [ 0.         0.         0.         0.000333   0.         0.      ]
 [ 0.         0.         0.         0.         0.00025    0.      ]
 [ 0.         0.         0.         0.         0.         0.0002  ]]

The stiffness matrix $\mathbf{C}$ is now the inverse of the compliance matrix $\mathbf{S}$:

\begin{equation} \mathbf{C}=\mathbf{S}^{-1} \tag{12} \end{equation}

Another function taking the material as the only argument:

In [4]:
def C3D(m):
    S=S3D(m)
    return np.linalg.inv(S)

C=C3D(m1)
print(np.array2string(C, precision=0, suppress_small=True, separator='  ', floatmode='maxprec_equal') )
[[155448.    9475.    6514.       0.       0.       0.]
 [  9475.   23435.    6111.       0.       0.       0.]
 [  6514.    6111.   11702.       0.       0.       0.]
 [     0.       0.       0.    3000.       0.       0.]
 [     0.       0.       0.       0.    4000.       0.]
 [     0.       0.       0.       0.       0.    5000.]]

As previously stated; normal strains cause only normal stress. For example:

In [5]:
strains=[0.003, 0.002, 0.001, 0.0, 0.0, 0.0]

stresses = np.dot(C,strains)

print(np.array2string(stresses, precision=1, suppress_small=True, separator='  ', floatmode='maxprec_equal') )
[491.8   81.4   43.5    0.0    0.0    0.0]
What if some of the stresses are known and some of the strains are known and we are supposed to find the unknown stresses and strains? The following function does the magic of organizing knowns and unknowns.
In [6]:
def solveHookes(C,**kwargs):
    # C: the stiffness matrix
    # kwargs: prescribed stresses or strains
    stress=np.zeros((6))
    strain=np.zeros((6))
    stressKeys=('s1','s2','s3','s23','s13','s12')   # valid stress keys
    strainKeys=('e1','e2','e3','e23','e13','e12')   # valid strain keys
    knownKeys= ['s1','s2','s3','s23','s13','s12']   #initially assume known stresses
    knownVals=np.array([0,0,0,0,0,0],float)      
    for key in kwargs:
        if key in stressKeys:
            i=stressKeys.index(key)
            knownKeys[i]=key
            knownVals[i]=kwargs[key]
        if key in strainKeys:
            i=strainKeys.index(key)
            knownKeys[i]=key
            knownVals[i]=kwargs[key]
    M1=-np.copy(C)
    M2= np.copy(C)
    ID= np.identity(6)
    for k in range(0,6):
        if knownKeys[k] in stressKeys:
            M2[:,k]=-ID[:,k]
        else:
            M1[:,k]=ID[:,k]
    sol=np.dot(  np.linalg.inv(M1),   np.dot(M2,knownVals))
    for i in range(0,6):
        if knownKeys[i] == stressKeys[i]:
            stress[i]=knownVals[i]
            strain[i]=sol[i]
        if knownKeys[i] == strainKeys[i]:
            strain[i]=knownVals[i]
            stress[i]=sol[i]
    return stress,strain

Example, comparing the previous example having $\sigma_1, \varepsilon_2$ and $\sigma_3$ as known nonzero values while assuming that other stresses are zero:

In [7]:
stress, strain = solveHookes(C=C, s1=491.8, e2=0.002, s3=43.46 )

print('stress=', np.array2string(stress, precision=1, suppress_small=True, separator='  ', floatmode='maxprec_equal') )
print('strain=', np.array2string(strain, precision=4, suppress_small=True, separator='  ', floatmode='maxprec_equal') )
stress= [491.8   81.4   43.5    0.0    0.0    0.0]
strain= [0.003  0.002  0.001  0.000  0.000  0.000]

Transversely isotropic materials

A transversely isotropic material is one with properties that are symmetric about an axis that is normal to a plane of isotropy. This transverse plane has infinite planes of symmetry and thus, within this plane, the material properties are the same in all directions.

Assuming that the plane $2$-$3$ is a plane of isotropy as illustrated in Figure-1, the elastic constants are related by:

\begin{equation} E_3=E_2, \quad \nu_{13}=\nu_{12}, \quad G_{13} = G_{12}, \quad G_{23}=\frac{E_2}{2(1+\nu_{23})} \tag{13} \end{equation}

A transversely isotropic material has therefore 5 independent elastic constants.

Unidirectional fiber composites are usually considered to be transversely isotropic.

The proof of the relations (13) is left as an exercise

Isotropic materials

Isotropic materials have only 2 independent elastic constants:

\begin{equation} E_{ij}=E, \quad \nu_{ij}=\nu, \quad G_{ij}=G = \frac{E}{2(1+\nu)} \tag{14} \end{equation}

Constraints on the engineering constants

Both the stiffness matrix and the compliance matrix must be positive definite.

For orthotropic, transversely isotropic and isotropic materials, this requirement implies that every subdeterminat of the main diagonal is positive. For the compliance matrix that is:

\begin{equation} \text{det} \begin{bmatrix} S_{11} & S_{12} & S_{13} & 0 & 0 & 0 \\ S_{12} & S_{22} & S_{23} & 0 & 0 & 0 \\ S_{13} & S_{23} & S_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & S_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & S_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & S_{66} \end{bmatrix} >0 \tag{15} \end{equation}\begin{equation} \text{det} \begin{bmatrix} S_{11} & S_{12} & S_{13} \\ S_{12} & S_{22} & S_{23} \\ S_{13} & S_{23} & S_{33} \end{bmatrix} >0 \tag{16} \end{equation}\begin{equation} \text{det} \begin{bmatrix} S_{11} & S_{12} \\ S_{12} & S_{22} \end{bmatrix} >0, \quad \text{det} \begin{bmatrix} S_{11} & S_{13} \\ S_{13} & S_{33} \end{bmatrix} >0, \quad \text{det} \begin{bmatrix} S_{22} & S_{23} \\ S_{23} & S_{33} \end{bmatrix} >0 \tag{17} \end{equation}\begin{equation} S_{11}>0,\quad S_{22}>0,\quad S_{33}>0,\quad S_{44}>0,\quad S_{55}>0,\quad S_{66}>0 \tag{18} \end{equation}

Example

In [8]:
S=S3D(m1)

print('Equation 15: det =', np.linalg.det(S)   )
print('Equation 16: det =', np.linalg.det(S[0:3,0:3])   )
Equation 15: det = 4.690740740740756e-25
Equation 16: det = 2.814444444444443e-14

Python tips & tricks

The functions S3D() and C3D() as well as other common functions can be added to a separate Python file for later use. For example, consider the following content in a file compositelib.py:

In [9]:
import numpy as np

def S3D(m):
    '''
    Input:
    - m: dictionary of material properties
    Output:
    - returns the compliance matrix as a 6x6 array
    '''
    return np.array(
    [[        1/m['E1'],-m['v12']/m['E1'],-m['v13']/m['E1'],         0,    0,     0],
     [-m['v12']/m['E1'],        1/m['E2'],-m['v23']/m['E2'],         0,    0,     0],
     [-m['v13']/m['E1'],-m['v23']/m['E2'],        1/m['E3'],         0,    0,     0],
     [                0,                0,                0,1/m['G23'],    0,     0],
     [                0,                0,                0,    0,1/m['G13'],     0],
     [                0,                0,                0,    0,     0,1/m['G12']] ],
    float)

def C3D(m):
    '''
    Input:
    - m: dictionary of material properties
    Output:
    - returns the stiffness matrix as a 6x6 array
    '''
    return np.linalg.inv(S3D(m))

If this file exists in the current directory (i.e., where your main Python code is executed) or in any of the search path (dependent on your system), the functions can be accessed by import:

In [10]:
import compositelib

C=compositelib.C3D(m1)

Alternatively:

In [11]:
from compositelib import C3D

C=C3D(m1)

The comments between the two ''' are available by help():

In [12]:
help(compositelib.C3D)
Help on function C3D in module compositelib:

C3D(m)
    Input:
    - m: dictionary of material properties
    Output:
    - returns the stiffness matrix as a 6x6 array

References and further readings

  1. Hibbeler, R.C., and Kai Beng Yap. Mechanics of Materials. Harlow: Pearson, 2018.
  2. Herakovich, Carl T. Mechanics of Fibrous Composites. New York: Wiley, 1998.
  3. Daniel, Isaac M., and Ori Ishai. Engineering Mechanics of Composite Materials. 2nd ed. New York: Oxford University Press, 2006.
  4. Kollár, Lázló P., and George S. Springer. Mechanics of Composite Structures. Cambridge: Cambridge University Press, 2003.

Disclaimer:This site is about polymer composites, designed for educational purposes. Consumption and use of any sort & kind is solely at your own risk.
Fair use: I spent some time making all the pages, and even the figures and illustrations are my own creations. Obviously, you may steal whatever you find useful here, but please show decency and give some acknowledgment if or when copying. Thanks! Contact me: nils.p.vedvik@ntnu.no www.ntnu.edu/employees/nils.p.vedvik

Copyright 2021, All right reserved, I guess.