Linear elasticity and notations
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.
Linear elasticity
Introduction to linear elasticity from Wikipedia:
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 leads into changes in shape and size. Hooke’s Law is equations, which is used to describe this relation:
\begin{equation} F = -kX \end{equation}
in which F is force, k is the slope and X is the elongation. The analogous equation of Hooke’s spring law for continuous media is then:
\begin{equation} \sigma = − c \epsilon \end{equation}
in which is second order stress tensor, c is a fourth order stiffness tensor end is a strain second order strain tensor. In a Cartesian coordinate system, the stress and strain tensors can be represented as:
Stiffness tensor 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. When isotropic media is used, can be constructed just from two numbers: (the bulk modulus ) and (the shear modulus). Wikipedia
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: \begin{equation} \epsilon = \frac{\Delta L}{L} \end{equation}
# 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 and . But, looking good. Let’s move one.
Index notation
Index notation (or in relativity theory, Einstein notation) is a handy way to deal with arrays, vectors 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 is the Kronecker’s delta and both and are the Lamé’s constants, 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 () and shear modulus ():
\begin{equation} K = \frac{E}{3 (1 - 2\nu)},\hspace{1cm}G = \frac{E}{2 (1 + \nu)} \end{equation}
in which is a poisson’s ratio and is Young’s modulus.
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. Using the tensor notation, stress can now be defined as:
\begin{equation} \sigma = \lambda \rm{tr}(\epsilon) \bf{I} + \mu \epsilon \end{equation}
in which and are Lamé constants, tr is trace operator, is strain and is an identity tensor:
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 operation, C is fourth order stiffness tensor and is second order strain tensor. Stiffness tensor can be presented as:
\begin{equation} C = \lambda \bf{I} \otimes \bf{I} + 2\mu \bf{l} \end{equation}
in which is second order identity tensor, is an outer product operation and is fourth order identity tensor. Equation for fourth order identity tensor:
\begin{equation} \bf{l_{\rm ijkl}}\rm = \frac{1}{2}(\delta_{ik} \delta_{jl} + \delta_{il}\delta_{jk}) \end{equation}
δ = (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, can be presented as a 2D matrix:
and the strain tensor can be reduced into vector:
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}
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
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]")
The end
In these four examples, we’ve calculated stress in 1D, using index, tensor and matrix notations and achieved the same result every time. Juggling with these notations can be troublesome sometimes but, in the end, each of them have their own place. Tensors are much handier when presenting the equations and formulas, matrix formulation is easier (at least in my opinion) for implementation.