Plasticity
Plasticity Theory
The constitutive theory is about relating stress $\sigma$ and strain $\varepsilon$.
The total strain $\varepsilon$ can be decomposed into
The constitutive relation must characterize the relation between the stress $\varepsilon$ for both $\varepsilon^p$ and $\varepsilon^e$ . The constitutive law for $\varepsilon^e$ is linear, i.e.,
The constitutive law for $\varepsilon^p$ is described by an ordinary differential equation
Here $h_{ij}$ may arise from a potential function $g(\sigma, \xi)$, where $\xi$ is called internal variables
and $\phi$ is a scalar function of the form
where $f$ is the yield surface.
The Tresca yield surface is given by
In $J_2$ plasticity, we have
Thus we have
In the rate-indpendent plasticity, there exists $\lambda$ such that
Since both sides of (1) has derivatives with respect to time, it is time scale independent. That's why we call it rate-independent. Then the flow rule requires that
If $f<0$ (the yield surface is not hit), the inelasticity strain is not "active"; once $f=0$ is hit, the material reacts by increasing plasticity strain (since $\dot \lambda$ can be nonzero).
The deformation theory studies how $g$ is related to $\sigma$ and $\xi$. If $h$ directs along the outward normal of the yield surface, we can the consequent $g$ the normality rule. For example, $J_2$ plasticity can be formulated is subjected to normality rule. A particular case where normality rule holds is $g=f$, in which case we call (1) the associated flow rule with the yield surface.
Finally, the dynamics of internal variable is given by
The three equations (1), (2) and (3) fully characterizes the constitutive relation of $\varepsilon^p$.
Numerical Example
Description | Displacement Field | Vertical Displacement |
---|---|---|
Plasticity | ![]() | ![]() |
Elasticity | ![]() | ![]() |
using Revise
using AdFem
using SparseArrays
using LinearAlgebra
using PyPlot
αm = 0.0
αf = 0.0
β2 = 0.5*(1 - αm + αf)^2
γ = 0.5 - αm + αf
m = 40
n = 20
h = 0.01
NT = 100
Δt = 1/NT
bdedge = []
for j = 1:n
push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
end
bdedge = vcat(bdedge...)
bdnode = Int64[]
for j = 1:n+1
push!(bdnode, (j-1)*(m+1)+1)
end
M = compute_fem_mass_matrix1(m, n, h)
S = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = [M S;S M]
H = diagm(0=>[1,1,0.5])
K = 0.1
σY = 0.03
# σY = 1000.
state = zeros(2(m+1)*(n+1),NT+1)
velo = zeros(2(m+1)*(n+1),NT+1)
acce = zeros(2(m+1)*(n+1),NT+1)
stress = zeros(NT+1, 4*m*n, 3)
internal_variable = zeros(NT+1, 4*m*n)
t = 0.0
for i = 1:NT
@info i
##### STEP 1: Computes the external force ####
T = eval_f_on_boundary_edge((x,y)->0.02*sin(2π*i*Δt), bdedge, m, n, h)
# T = eval_f_on_boundary_edge((x,y)->0.0, bdedge, m, n, h)
T = [zeros(length(T)) -T]
T = compute_fem_traction_term(T, bdedge, m, n, h)
f1 = eval_f_on_gauss_pts((x,y)->0., m, n, h)
f2 = eval_f_on_gauss_pts((x,y)->0., m, n, h)
# f2 = eval_f_on_gauss_pts((x,y)->0.1, m, n, h)
F = compute_fem_source_term(f1, f2, m, n, h)
fext = F+T
##### STEP 2: Extract Variables ####
u = state[:,i]
∂∂u = acce[:,i]
∂u = velo[:,i]
ε0 = eval_strain_on_gauss_pts(u, m, n, h)
σ0 = stress[i,:,:]
α0 = internal_variable[i,:]
##### STEP 3: Newton Iteration ####
global t += (1 - αf)*Δt
∂∂up = ∂∂u[:]
iter = 0
while true
iter += 1
# @info iter
up = (1 - αf)*(u + Δt*∂u + 0.5 * Δt^2 * ((1 - β2)*∂∂u + β2*∂∂up)) + αf*u
global fint, stiff, α, σ = compute_planestressplasticity_stress_and_stiffness_matrix(
up, ε0, σ0, α0, K, σY, H, m, n, h
)
res = M * (∂∂up *(1 - αm) + αm*∂∂u) + fint - fext
A = M*(1 - αm) + (1 - αf) * 0.5 * β2 * Δt^2 * stiff
A, _ = fem_impose_Dirichlet_boundary_condition(A, bdnode, m, n, h)
res[[bdnode; bdnode .+ (m+1)*(n+1)]] .= 0.0
Δ∂∂u = A\res
∂∂up -= Δ∂∂u
err = norm(res)
# @info err
if err<1e-8
break
end
end
global t += αf*Δt
##### STEP 3: Update State Variables ####
u += Δt * ∂u + Δt^2/2 * ((1 - β2) * ∂∂u + β2 * ∂∂up)
∂u += Δt * ((1 - γ) * ∂∂u + γ * ∂∂up)
stress[i+1,:,:] = σ
internal_variable[i+1,:] = α
state[:,i+1] = u
acce[:,i+1] = ∂∂up
velo[:,i+1] = ∂u
end
x = []
y = []
for j= 1:n+1
for i = 1:m+1
push!(x, (i-1)*h)
push!(y, (j-1)*h)
end
end
for i = 1:5:NT+1
close("all")
scatter(x+state[1:(m+1)*(n+1), i], y+state[(m+1)*(n+1)+1:end, i])
scatter(x[m+1]+state[m+1, i],
y[m+1]+state[(m+1)*(n+1)+m+1, i], color="red")
xlabel("x")
ylabel("y")
k = string(i)
k = repeat("0", 3-length(k))*k
title("t = $k")
ylim(-0.05,0.25)
xlim(-0.01, 0.45)
gca().invert_yaxis()
savefig("u$k.png")
close("all");
plot(1:i, -state[(m+1)*(n+1)+m+1, 1:i])
xlim(0, NT+2)
ylim(0, 0.03)
scatter(i, -state[(m+1)*(n+1)+m+1, i], color="red")
savefig("du$k.png")
end
run(`convert -delay 10 -loop 0 u*.png plasticity_u.gif`)
run(`convert -delay 10 -loop 0 du*.png plasticity_du.gif`)