Viscoplasticity
TL;DR
After practicing with elastic material model, we’ll continue to nonlinear material modelling: viscoplasticity. In this post we’ll explore Norton and Bingham viscoplastic potential equations, von Mises yield function and automatic derivation using ForwardDiff. Most of the theory can be found from here and here.
Note! Equations are written in tensor notation, but the implemented code is in matrix notation!
Recap to linear material modelling
Here’s a small recap example from elastic material modelling with matrix notation. Do note that we’re using now cyclic loading.
\begin{equation} \sigma = C : \epsilon \end{equation}
using PyPlot
using ForwardDiff
num = 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, 5, num) # Change of length [mm]
ϵ_ = collect(ΔL / L)
ϵ_total = vec([ϵ_; reverse(ϵ_[2:end-1]); -ϵ_; reverse(-ϵ_[2:end-1])])
Δt = 0.1
accumulator = 0.0
I = eye(3)
total_num = length(ϵ_total)
time_ = zeros(total_num)
σ = zeros(total_num, 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:(total_num)
ϵ = zeros(6)
ϵ[1] = ϵ_total[i]
ϵ[2] = -ϵ_total[i] * ν
ϵ[3] = -ϵ_total[i] * ν
# The most general form of Hooke's law for isotropic materials
σ[i, :] = C * ϵ
accumulator += Δt
time_[i] = accumulator
end
σ_11 = vec(σ[:, 1])
fig = figure()
grid()
plot(ϵ_total, σ_11)
ylabel("Stress [MPa]")
xlabel("Strain")
title("Stress-strain cyclic curve")
fig = figure()
grid()
title("Strain history")
plot(time_, ϵ_total)
xlabel("Time [s]")
ylabel("Strain")
PyObject <matplotlib.text.Text object at 0x7f38d9d57950>
Nonlinear material
Stress is always related only to the elastic part of the the strain. The theory of plasticity is used to define the non-linear region of the material behaviour. Moving from linear to nonlinear modelling, we’ll have to introduce a new variable: plastic strain \(\epsilon^p\). Here’s a definition for plasticity from wiki:
In physics and materials science, plasticity describes the deformation of a (solid) material undergoing non-reversible changes of shape in response to applied forces. For example, a solid piece of metal being bent or pounded into a new > shape displays plasticity as permanent changes occur within the material itself. In engineering, the transition from > elastic behavior to plastic behavior is called yield.
Decomposition of strain
Theory presented here can be found from here. Strain increment \(d\epsilon\) is decomposed into the elastic part \(d\epsilon^e\) and the viscoplastic part \(d\epsilon^{vp}\) using kinematics. Additive decomposition is used for small strains and multiplicative decomposition for large strains. In elastic-plastic materials the strain increment is decomposed into an elastic and a viscoplastic part:
\begin{equation} d\epsilon = d\epsilon^e + d\epsilon^{vp} \end{equation}
Dividing both sides with the differential time increment \(dt\) produces the rate form:
\begin{equation} \dot\epsilon = \dot\epsilon^e + \dot\epsilon^{vp} \end{equation}
Substituting the elastic strain rate into Hooke’s law produces the stress-rate relation:
\begin{equation} \dot\sigma = C : \dot\epsilon^e = C : (\dot\epsilon - \dot\epsilon^{vp}) \end{equation}
In the non-linear regime, the stress-strain relation is defined as:
\begin{equation} \dot\sigma = C : \dot\epsilon^e = C^{tan} : \dot\epsilon \end{equation}
Stress can now be expressed as: \begin{equation} \sigma = \int_0^t \dot\sigma dt’ \end{equation}
but usually, stress is more conveniently written as:
\begin{equation} \sigma_{n+1} = \sigma_n + C : (\dot\epsilon - \dot\epsilon^{vp})\,\Delta t \end{equation}
Yield function
Yield function defines the transition between elastic and plastic zones. One of the most common yield functions is von Mises yield function, which is defined as:
\begin{equation} f(\sigma) = \sigma_{eq}(\sigma) - \sigma_y \end{equation}
in which \(\sigma_y\) is yield stress, which defines the limit, when plastic strain starts and \(\sigma_{eq}\) is equivalent stress, which is expressed as:
\begin{equation} \sigma_{eq}(\sigma) = \sqrt{3J_2(\sigma)} \end{equation}
in which \(J_2\) is second deviatoric stress invariant:
\begin{equation} J_2(\sigma) = \frac{1}{2}s:s \end{equation}
and \(s\) is deviatoric stress, which is defined as:
\begin{equation} s = \sigma - \frac{1}{3}\rm{tr}(\sigma)I \end{equation}
in which tr is a trace operator and \(\textbf{I}\) is second order identity tensor. What the yield function basically tells us: if the return value of given stress tensor is negative, we are in the elastic region and if the returned value is positive, we have exceeded the yield limit and ventured into plastic region. In order to take this into account in the implementation, we’ll make a simple comparison operation with:
\begin{equation} f(\sigma) > 0 \end{equation}
which will define if we apply viscoplastic strain or not to our stress.
""" Deviatoric stress tensor
"""
function deviatoric_stress(σ)
σ_tensor = [σ[1] σ[4] σ[6];
σ[4] σ[2] σ[5];
σ[6] σ[5] σ[3]]
σ_dev = σ_tensor - 1/3 * trace(σ_tensor) * eye(3)
return [σ_dev[1, 1], σ_dev[2, 2], σ_dev[3, 3], σ_dev[1, 2], σ_dev[2, 3], σ_dev[1, 3]]
end
""" Second deviatoric stress invariant
"""
function J_2_stress(σ)
s = deviatoric_stress(σ)
s_vec = vec([s; s[4]; s[5]; s[6]]) # Adding missing shear elements for double dot product
return 1/2 * dot(s_vec, s_vec) # equivalent to 1/2 * s : s, in tensor notation
end
""" Equivalent stress
"""
function equivalent_stress(σ)
J_2 = J_2_stress(σ)
return sqrt(3 * J_2)
end
""" Von mises yield function
"""
function von_mises_yield(σ, σ_y)
return equivalent_stress(σ) - σ_y
end
von_mises_yield (generic function with 1 method)
Here are presented von Mises yield surfaces both in 2D and 3D. Used yield stress is 200 MPa. As you can see, the von Mises yield surface is an ellipse in the 2D principal-stress plane. Inside the ellipse represents the elastic region and outside the plastic region. Do note that the x- and y-axes are principal stresses and not \(\sigma_{x}\) and \(\sigma_{y}\).
The 3D presentation reveals the whole truth and we can see that the elliptical shape, which we saw in 2D plane, is merely a cross-section of the cylinder, which is the von Mises yield surface. The axis of the cylinder resides at the hydrostatic pressure axis.
Viscoplastic potential
Next we’ll need to define the viscoplastic strain rate \(\dot\epsilon^{vp}\):
\begin{equation} \dot\epsilon^{vp} = \frac{\partial\Phi}{\partial\sigma} \end{equation} in which \(\Phi\) is a plastic potential, which is a function of the stress and the material parameters. In the post, we’ll look at both Norton rule and Bingham model. But first, the Norton rule:
\begin{equation} \Phi = \frac{K}{n + 1} \left(\frac{J}{K}\right)^{n+1}, \hspace{1cm} J = f(\sigma) \end{equation}
in which \(K\) and \(n\) are material parameters and \(f\) is a yield function. When we use von Mises yield function, we can differentiate this function with respect to \(\sigma\):
\begin{equation} \frac{\partial\Phi}{\partial\sigma} = \left(\frac{J}{K}\right)^{n}\frac{3}{2}\frac{s}{\sigma_{eq}} \end{equation}
""" Norton rule
"""
function norton_plastic_potential(σ, K, n, f)
f_ = f(σ)
return K/(n+1) * (f_ / K) ^ (n + 1)
end
""" Analytically differentiated norton rule
"""
function analytical_dΦdσ(σ, K, n, f)
f_ = f(σ)
σ_v = equivalent_stress(σ)
σ_dev_vec = deviatoric_stress(σ)
voigt_factor = [1.0, 1.0, 1.0, 2.0, 2.0, 2.0]
return (f_ / K) ^ n * 3/2 * (σ_dev_vec .* voigt_factor) / σ_v
end
analytical_dΦdσ (generic function with 1 method)
Summary of needed functions
Here’s a summary of the main functions, which we need:
\begin{equation} \sigma_{n+1} = \sigma_n + C : (\dot\epsilon - \dot\epsilon^{vp})\,\Delta t \hspace{1.5cm}\dot\epsilon^{vp} = \left(\frac{J}{K}\right)^{n}\frac{3}{2}\frac{s}{\sigma_{eq}}\hspace{1.5cm} J = f(\sigma) = \sigma_{eq}(\sigma) - \sigma_y \end{equation}
Now, we can proceed to the simulation:
# Create an empty history vectors
σ_hist = zeros(total_num, 6)
ϵ_hist = zeros(total_num, 6)
# Initialize arrays
ϵ_last = zeros(6)
σ = zeros(6)
ϵ = zeros(6)
# We'll define our yield limit to be 200 MPa
σ_y = 200.0
# Define Δt, K and n
Δt = 0.1
n = 0.92
K = 180.0e3
# Let's make a curried function to make it pretty!
f = x->von_mises_yield(x, 200.)
# Ok, going through time steps
for i=2:total_num
ϵ = zeros(6)
ϵ[1] = ϵ_total[i]
ϵ[2] = -ϵ_total[i] * ν
ϵ[3] = -ϵ_total[i] * ν
# Trial stress
dϵdt = (ϵ - ϵ_last) / Δt
dσdt = C * dϵdt
σ_tria = σ + dσdt * Δt
# Checking if yielded
yield_ = f(σ_tria)
if yield_ > 0
# If yielded, we'll add the viscoplastic strain to stress
dϵ_vp = analytical_dΦdσ(σ_tria, K, n, f)
σ += C * (dϵdt - dϵ_vp) * Δt
else
σ[:] = σ_tria
end
# Collect data
σ_hist[i, :, :] = σ
ϵ_last[:] = ϵ
ϵ_hist[i, :] = ϵ
end
# Plot
plot(ϵ_hist[:, 1], σ_hist[:, 1])
grid()
xlabel("Strain")
ylabel("Stress x [MPa]")
PyObject <matplotlib.text.Text object at 0x7ff74adb2610>
As you can see, plot has now changed. At 200 MPa our stress exceeds the yield limit and viscoplastic strain is added. After reaching strain level of 0.005, we’ll start to go back to zero and return to the elastic strain region. At ~20 MPa we exceed the yield limit again and start adding viscoplastic strain. I also wanted to show the power of automatic derivation using ForwardDiff in the same calculation.
σ_hist = zeros(total_num, 6)
ϵ_hist = zeros(total_num, 6)
ϵ_last = zeros(6)
ϵ = zeros(6)
σ = zeros(6)
σ_y = 200.0
Δt = 0.1
n = 0.92
K = 180.0e3
f = x->von_mises_yield(x, σ_y)
for i=2:total_num
ϵ = zeros(6)
ϵ[1] = ϵ_total[i]
ϵ[2] = -ϵ_total[i] * ν
ϵ[3] = -ϵ_total[i] * ν
dϵdt = (ϵ - ϵ_last) / Δt
dσdt = C * dϵdt
σ_tria = σ + dσdt * Δt
yield_ = f(σ_tria)
if yield_ > 0
wrap = x-> norton_plastic_potential(x, K, n, f)
# This part is different. As you can see, we'll take an automatic derivative
# Way more convenient compared to hand-differentiated potential function
dϵ_vp = ForwardDiff.gradient(wrap, σ_tria)
σ += C * (dϵdt - dϵ_vp) * Δt
else
σ[:] = σ_tria
end
σ_hist[i, :, :] = σ
ϵ_last[:] = ϵ
ϵ_hist[i, :] = ϵ
end
plot(ϵ_hist[:, 1], σ_hist[:, 1])
grid()
Using the ForwardDiff now enables me to differentiate quite complex looking functions and also saves time from debugging hand-differentiated functions. Lastly, we’ll try another viscoplastic potential (Bingham model) for our plastic strain rate, which is defined as:
\begin{equation} \Phi = \frac{1}{2} \left(\frac{J}{\eta}\right)^2, \hspace{1cm} J=f(\sigma) \end{equation}
function bingham_plastic_potential(σ, η, f)
return 0.5 * (f(σ) / η) ^2
end
bingham_plastic_potential (generic function with 1 method)
Now that we’ve discovered the power of ForwardDiff, we’ll also use automatic derivation in this simulation.
ϵ_last = zeros(6)
ϵ_hist = zeros(total_num, 6)
σ_hist = zeros(total_num, 6)
σ = zeros(6)
σ_y = 200.0
Δt = 0.01
ϵ = zeros(6)
# Parameter for Bingham model
η=190.
f = x->von_mises_yield(x, σ_y)
for i=2:total_num
ϵ = zeros(6)
ϵ[1] = ϵ_total[i]
ϵ[2] = -ϵ_total[i] * ν
ϵ[3] = -ϵ_total[i] * ν
dϵdt = (ϵ - ϵ_last) / Δt
dσdt = C * dϵdt
σ_tria = σ + dσdt * Δt
yield_ = f(σ_tria)
if yield_ > 0
wrap = x->bingham_plastic_potential(x, η, f)
# Again, I'm a bit lazy and just use ForwardDiff
dϵ_vp = ForwardDiff.gradient(wrap, σ_tria)
σ += C * (dϵdt - dϵ_vp) * Δt
else
σ[:] = σ_tria
end
σ_hist[i, :, :] = σ
ϵ_last[:] = ϵ
ϵ_hist[i, :] = ϵ
end
plot(ϵ_hist[:, 1], σ_hist[:, 1])
grid()
The End
In this post we have seen how to implement non-linear viscoplastic material model. Varying yield function or viscoplastic potential we can create different material responses.