Plasticity

Plasticity Theory

The constitutive theory is about relating stress $\sigma$ and strain $\varepsilon$.

The total strain $\varepsilon$ can be decomposed into

\[\varepsilon = \varepsilon^p + \varepsilon^e\]

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.,

\[\sigma = C\varepsilon^e\Leftrightarrow\varepsilon^e = C^{-1}\sigma\]

The constitutive law for $\varepsilon^p$ is described by an ordinary differential equation

\[\boxed{\dot\varepsilon^p_{ij} = \phi h_{ij}}\tag{1}\]

Here $h_{ij}$ may arise from a potential function $g(\sigma, \xi)$, where $\xi$ is called internal variables

\[h_{ij} = \frac{\partial g}{\partial \sigma_{ij}}\]

and $\phi$ is a scalar function of the form

\[\phi = \eta(\sigma,\xi) f(\sigma, \xi)_+ := \eta(\sigma,\xi)\max(0, f(\sigma, \xi))\]

where $f$ is the yield surface.

In the rate-indpendent plasticity, there exists $\lambda$ such that

\[\phi = \dot\lambda\]

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

\[\boxed{\dot \lambda f = 0\quad \lambda\geq 0\quad f\leq 0}\tag{2}\]

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

\[\boxed{\dot\xi_\alpha = \tilde h_\alpha(\sigma, \xi)}\tag{3}\]

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`)