# Linear elasticity and notations

## TL;DR

In engineering mechanics there are various way how to present equations and theorems. In this post, I'll go through some linear elasticity using index, tensor and matrix notation. All example are made with Julia programming language. Most of the formulas used here can be found from [here](https://en.wikipedia.org/wiki/Hooke%27s_law).

## Linear elasticity

Introduction to linear elasticity from [Wikipedia](https://en.wikipedia.org/wiki/Linear_elasticity):

`Linear elasticity is the mathematical study of how solid objects deform and become internally stressed due to prescribed loading conditions. Linear elasticity models materials as continua. Linear elasticity is a simplification of the more general nonlinear theory of elasticity and is a branch of continuum mechanics. The fundamental "linearizing" assumptions of linear elasticity are: infinitesimal strains or "small" deformations (or strains) and linear relationships between the components of stress and strain. In addition linear elasticity is valid only for stress states that do not produce yielding. These assumptions are reasonable for many engineering materials and engineering design scenarios. Linear elasticity is therefore used extensively in structural analysis and engineering design, often with the aid of finite element analysis.`

Applying force to a continuum leads to deformation of the medium. This results in changes in shape and size. Hooke's Law is equations, used to descibe this relation:

\begin{equation}
F = -kX
\end{equation}

in which F is force, k is the slope and X is the elongation. The analogous of Hooke's spring law for continuous media is then

\begin{equation}
\sigma = − c \epsilon
\end{equation}

in which $\sigma$ is second order stress tensor, c is a fourth-order stiffness tensor end $\epsilon$ is a strain second order strain tensor. In a Cartesian coordinate system, the stress and strain tensors can be represented as:

<div align="center">
$
\sigma = \begin{bmatrix}
       \sigma_{11} & \sigma_{12} & \sigma_{13} \\[0.3em]
       \sigma_{21} & \sigma_{22} & \sigma_{23} \\[0.3em]
       \sigma_{31} & \sigma_{32} & \sigma_{33}
     \end{bmatrix}
$, 

$
\epsilon = \begin{bmatrix}
       \epsilon_{11} & \epsilon_{12} & \epsilon_{13} \\[0.3em]
       \epsilon_{21} & \epsilon_{22} & \epsilon_{23} \\[0.3em]
       \epsilon_{31} & \epsilon_{32} & \epsilon_{33}
     \end{bmatrix}
$
</div>

stiffness tensor $c$ can be represented as a 3x3x3x3 matrix, which has a total of 81 individual elements. In this post, we'll be only looking at isotropic media (which have the same physical properties in any direction). When isotropic media is used, $c$ can be constructed from two numbers: $K$ (the bulk modulus ) and $G$ (the shear modulus). [Wikipedia](https://en.wikipedia.org/wiki/Linear_elasticity)

## 1D example using Hooke's Law

Before dive into tensors and matrices, we'll calculate a small 1D example with Hooke's law. Equation for stress:

\begin{equation}
\sigma = E\epsilon
\end{equation}

in which E is the Young's modulus. Strain can be calculated using using engineering strain( https://en.wikipedia.org/wiki/Deformation_(mechanics) ):
\begin{equation}
\epsilon = \frac{\Delta L}{L}
\end{equation}



```julia
# 1D
using Plots

L = 1000.0               # Initial length [mm]
ΔL = linspace(0, 2, 100) # Change of length [mm] 

ϵ = ΔL / L               # Strain 

E = 200000.0             # Young's modulus for steel [MPa]
σ = E * ϵ                # Stress [MPa]

plot(ϵ, σ)
xlabel!("Strain")
ylabel!("Stress [MPa]")
```

    [Plots.jl] Initializing backend: pyplot









This is the plot, which we want to archive in each example. In the following examples we'll be using the same $L$ and $\Delta L$ also. But, looking good. Now we can just into the real deal.

## Index notation

Index notation (or in relativity theory, Einstein notation) is a hand way to deal with arrays, vector and matrices. Using the isotropic material (media), stress can be defined as:

\begin{equation}
\sigma_{ij} = \lambda \epsilon_{kk}\delta_{ij} + 2\mu \epsilon_{ij}
\end{equation}

in which $\delta$ is the [Kronecker's delta](https://en.wikipedia.org/wiki/Kronecker_delta) and both $\lambda$ and $\mu$ are the [Lamé's constant](# https://en.wikipedia.org/wiki/Lam%C3%A9_parameters), which can be calculated with:
\begin{equation}
\lambda = K - \frac{2}{3}G,\hspace{1cm}\mu = G
\end{equation}

Lastly, we'll define the bulk modulus ($K$) and shear modulus ($G$):

\begin{equation}
K = \frac{E}{3 (1 - 2\nu)},\hspace{1cm}G = \frac{E}{2 (1 + \nu)}
\end{equation}

in which $\nu$ is a poisson's ratio and $E$ is Young's modulus.



```julia
n = 100                      # Number of steps
E = 200000.0                 # Young's modulus for steel [MPa]
ν = 0.3                      # Poisson's ratio for steel

K = E / (3 *(1 - 2*ν))       # Bulk modulus
G = E / (2 * (1 + ν))        # Shear modulus 

λ = K - 2/3 * G              # Lamé constant
μ = G                        # Lamé constant
L = 1000.0                   # Initial length [mm]
ΔL = linspace(0, 2, n)       # Change of length [mm] 

ϵ_total = ΔL / L             # Total strain

σ = zeros(n, 3, 3)           # Empty stress history

δ = (i, j) -> i == j ? 1 : 0 # Kronecker delta

for t=1:n
    ϵ = zeros(3,3)
    
    # Collect strain and add poisson effect
    ϵ[1, 1] = ϵ_total[t] 
    ϵ[2, 2] = -ϵ_total[t] * ν
    ϵ[3, 3] = -ϵ_total[t] * ν
    
    # Calculate stress
    for i=1:3
        for j=1:3
            ϵ_kk = sum([ϵ[k, k] for k=1:3])
            σ[t, i, j] += λ * ϵ_kk * δ(i, j) + 2 * μ  * ϵ[i, j]
        end
    end
    
end

σ_11 = vec(σ[:, 1, 1]);
plot(ϵ_total, σ_11)
xlabel!("Strain")
ylabel!("Stress [MPa]")
```








Plot remains the same, as it should be.

## Tensor notation

Now that we're comfortable using indeces, we'll proceed to tensor notation. Here's a quotation from [wikipedia]() from tensors:

`Tensors are geometric objects that describe linear relations between geometric vectors, scalars, and other tensors. Elementary examples of such relations include the dot product, the cross product, and linear maps. Euclidean vectors, often used in physics and engineering applications, and scalars themselves are also tensors. A more sophisticated example is the Cauchy stress tensor T, which takes a direction v as input and produces the stress T(v) on the surface normal to this vector for output, thus expressing a relationship between these two vectors.
`

Sound intimidating, but tensor are marely a generalization of scalars and vectors. As said in the quotation, scalars are rank 0 tensors and vectors are rank 1 tensors. One of the nice properties of tensors is that they are independent of any chosen frame of reference. Stress can now be defined in tensor notation:

\begin{equation}
\sigma = \lambda \rm{tr}(\epsilon) \bf{I} + \mu \epsilon
\end{equation}

in which $\lambda$ and $\mu$ are Lamé constants, tr is [trace operator](https://en.wikipedia.org/wiki/Trace_%28linear_algebra%29), $\epsilon$ is strain and $\bf{I}$ is an identity tensor:

\begin{equation}
\bf{I}\rm = \begin{bmatrix}
       1 & 0 & 0 \\[0.3em]
       0 & 1 & 0 \\[0.3em]
       0 & 0 & 1
     \end{bmatrix}
\end{equation}     


```julia
n = 100                  # Number of steps
E = 200000.0             # Young's modulus [GPa]
ν = 0.3                  # Poisson's ratio

K = E / (3 *(1 - 2*ν))   # Bulk modulus
G = E / (2 * (1 + ν))    # Shear modulus 

λ = K - 2/3 * G          # Lamé constant
μ = G                    # Lamé constant

L = 1000.0               # Initial length [mm]
ΔL = linspace(0, 2, n)   # Change of length [mm] 

ϵ_total = ΔL / L
I = eye(3)
σ = zeros(n, 3, 3)       

for i=1:n
    ϵ = zeros(3,3)
    ϵ[1, 1] = ϵ_total[i] 
    ϵ[2, 2] = -ϵ_total[i] * ν
    ϵ[3, 3] = -ϵ_total[i] * ν
    
    # The most general form of Hooke's law for isotropic materials
    σ[i, :, :] = λ * trace(ϵ) * I + 2 * μ * ϵ
end

σ_11 = vec(σ[:, 1, 1])
plot(ϵ_total, σ_11)
xlabel!("Strain")
ylabel!("Stress [MPa]")
```








Usually (at least what I've studied during these years) the stress equation is presented as:

\begin{equation}
\sigma = C : \epsilon
\end{equation}

in which : is [double dot product](https://en.wikipedia.org/wiki/Dyadics) operation, C is fourth order stiffness tensor and $\epsilon$ is second order strain tensor. Stiffness tensor can be presented as:

\begin{equation}
C = \lambda \bf{I} \otimes \bf{I} + 2\mu \bf{|}
\end{equation}

in which $\bf{I}$ is second order identity tensor, $\otimes$ is an [outer product operation](https://en.wikipedia.org/wiki/Outer_product) and $\bf{|}$ is [fourth order identity tensor](http://www.ce.berkeley.edu/~sanjay/ce231mse211/symidentity.pdf). Equation for fourth order identity tensor:
 
\begin{equation}
\bf{|_{\rm ijkl}}\rm = \frac{1}{2}(\delta_{i,k} \delta_{j,l} + \delta_{i,l}\delta_{j,k})
\end{equation}


```julia
δ = (i, j) -> i == j ? 1 : 0 # Kronecker delta function

""" Double dot product operation. This implementation works
only for (3,3) matrices """
function double_dot_product(E, ϵ, i, j)
    stress = 0.0
    for k=1:3
        for l=1:3
            stress += E[i, j, k, l] * ϵ[k, l]
        end
    end
    stress
end

""" Function for creating fourth order identity tensor """
function create_4th_symmetric_identity_tensor()
    II = (i,j,k,l) -> 1/2 * (δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
    sym_fourth_eye = zeros(3,3,3,3)
    for i=1:3
        for j=1:3
            for k=1:3
                for l=1:3
                    sym_fourth_eye[i,j,k,l] = II(i,j,k,l)
                end
            end
        end
    end    
    sym_fourth_eye
end

""" Outer product operation. This implementation works only 
for (3,3) matrices """
function outer_product(A, B)
    tensor = zeros(3,3,3,3)
    for i=1:3
        for j=1:3
            for k=1:3
                for l=1:3
                    tensor[i,j,k,l] = A[i, j] * B[k, l]
                end
            end
        end
    end
    tensor
end


n = 100                  # Number of steps
E = 200000.0             # Young's modulus [GPa]
ν = 0.3                  # Poisson's ratio

K = E / (3 *(1 - 2*ν))   # Bulk modulus
G = E / (2 * (1 + ν))    # Shear modulus 

λ = K - 2/3 * G          # Lamé constant
μ = G                    # Lamé constant
L = 1000.0               # Initial length [mm]
ΔL = linspace(0, 2, n)   # Change of length [mm] 

ϵ_total = ΔL / L        
I = eye(3)
σ = zeros(n, 3, 3)

I = eye(3)
sym_fourth_eye = create_4th_symmetric_identity_tensor()

C = λ * outer_product(I, I) + 2*μ*sym_fourth_eye 

for t=1:n
    ϵ = zeros(3,3)
    ϵ[1, 1] = ϵ_total[t] 
    ϵ[2, 2] = -ϵ_total[t] * ν
    ϵ[3, 3] = -ϵ_total[t] * ν
    
    # Calculating each element for stress
    for i=1:3
        for j=1:3
            σ[t, i, j] = double_dot_product(C, ϵ, i, j)
        end
    end
    
end

σ_11 = vec(σ[:, 1, 1])
plot(ϵ_total, σ_11)
xlabel!("Strain")
ylabel!("Stress [MPa]")
```








## Matrix notation

There are justified reasons presenting equations with tensor notation but the
practical implementation of equations is faster using matrix notation. In matrix notation rank of the tensors is so to speak, dropped. Due to symmetricity, $C$ can be presented as a 2D matrix:

\begin{equation}
C = \frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix}
       1-\nu & \nu & \nu & 0 & 0 & 0 \\[0.3em]
       \nu & 1-\nu & \nu & 0 & 0 & 0 \\[0.3em]
       \nu & \nu & 1-\nu & 0 & 0 & 0 \\[0.3em]
       0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\[0.3em]
       0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\[0.3em]
       0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2}
     \end{bmatrix}
\end{equation}

and the strain tensor can be reduced into vector:

\begin{equation}
\epsilon = \begin{bmatrix}
       \epsilon_{11} \\[0.3em]
       \epsilon_{22} \\[0.3em]
       \epsilon_{33} \\[0.3em]
       2\epsilon_{12} \\[0.3em]
       2\epsilon_{23} \\[0.3em]
       2\epsilon_{13}
     \end{bmatrix}
\end{equation}

The double dot product, used in tensor notation can be now changed into dot product and the stress can be defined as:

\begin{equation}
\sigma = C \cdot \epsilon
\end{equation}


```julia
n = 100                  # Number of steps
E = 200000.0             # Young's modulus [GPa]
ν = 0.3                  # Poisson's ratio

L = 1000.0               # Initial length [mm]
ΔL = linspace(0, 2, n) # Change of length [mm] 

ϵ_total = ΔL / L            # Strain in x-direction
I = eye(3)
σ = zeros(n, 6)

C = E / ((1 + ν)*(1 - 2*ν)) * [
    1-ν   ν   ν         0         0         0
      ν 1-ν   ν         0         0         0
      ν   ν 1-ν         0         0         0
      0   0   0 (1-2*ν)/2         0         0
      0   0   0         0 (1-2*ν)/2         0
      0   0   0         0         0 (1-2*ν)/2]

for i=1:n
    ϵ = zeros(6)
    ϵ[1] = ϵ_total[i] 
    ϵ[2] = -ϵ_total[i] * ν
    ϵ[3] = -ϵ_total[i] * ν
    
    # The most general form of Hooke's law for isotropic materials
    σ[i, :, :] = C * ϵ
end

σ_11 = vec(σ[:, 1])
plot(ϵ_total, σ_11)
xlabel!("Strain")
ylabel!("Stress [MPa]")
```








## Discussion

Juggling with tensors, matrices and indeces can be troublesome but in the end, each of them has their own place. Tensors are much handier when presenting the equations and formulas, matrix formulation is easier (at least in my opinion) to implement.

b813d01f5806a38928677fc9df440022d3d4d929
ipynb
Jupyter Notebook
_posts/unfinished/linear_elasticity/Linear elasticity.ipynb
ovainola/ovainola.github.com
08b75f3822b388fa3db86792a025101baa5f8fbf
CC-BY-3.0
_posts/unfinished/linear_elasticity/Linear elasticity.ipynb
ovainola/ovainola.github.com
08b75f3822b388fa3db86792a025101baa5f8fbf
CC-BY-3.0
_posts/unfinished/linear_elasticity/Linear elasticity.ipynb
ovainola/ovainola.github.com
08b75f3822b388fa3db86792a025101baa5f8fbf
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES