Skip to content

Commit ada3e36

Browse files
KnutAMKnutAM
authored
SimplePlastic (#8)
1 parent de6091e commit ada3e36

4 files changed

Lines changed: 221 additions & 6 deletions

File tree

‎src/MechanicalMaterialModels.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,9 @@ include("ExtraOutputs.jl")
4343
include("PlasticDifferentiate.jl")
4444
export Plastic
4545

46+
include("specialized_models/SimplePlastic.jl")
47+
export SimplePlastic
48+
4649
include("CrystalPlasticity.jl")
4750
export CrystalPlasticity
4851

‎src/specialized_models/SimplePlastic.jl‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
@kwdef struct SimplePlastic{T} <: AbstractMaterial
2+
G::T
3+
K::T
4+
Y0::T
5+
Hiso::T
6+
κ∞::T
7+
Hkin::T
8+
β∞::T
9+
end
10+
11+
function SimplePlastic(mp::Plastic{Ela, Yld, Iso, Kin, Rat}) where {
12+
Ela <: LinearElastic{<:Any, :isotropic},
13+
Yld <: VonMises,
14+
Iso <: Tuple{<:Voce},
15+
Kin <: Tuple{<:ArmstrongFrederick},
16+
Rat <: RateIndependent
17+
}
18+
19+
E, ν = mp.elastic.p
20+
return SimplePlastic(;
21+
G = E/(2 * (1 + ν)), K = E / (3 * (1 - 2 * ν)), Y0 = mp.yield.Y0,
22+
Hiso = only(mp.isotropic).Hiso, κ∞ = only(mp.isotropic).κ∞,
23+
Hkin = only(mp.kinematic).Hkin, β∞ = only(mp.kinematic).β∞)
24+
end
25+
26+
struct SimplePlasticState{T} <: AbstractMaterialState
27+
ϵp::SymmetricTensor{2,3,T,6}
28+
β::SymmetricTensor{2,3,T,6}
29+
κ::T
30+
end
31+
32+
function MMB.initial_material_state(::SimplePlastic)
33+
return SimplePlasticState(zero(SymmetricTensor{2,3}), zero(SymmetricTensor{2,3}), 0.0)
34+
end
35+
36+
function MMB.material_response(m::SimplePlastic, ϵ::SymmetricTensor{2,3}, old::SimplePlasticState, Δt, args...)
37+
σ_trial = 2 * m.G * dev(ϵ - old.ϵp) + 3 * m.K * vol(ϵ)
38+
ϕ_trial = vonmises(σ_trial - old.β) - (m.Y0 + old.κ)
39+
if ϕ_trial ≤ 0
40+
I2 = one(ϵ)
41+
IxI = I2 ⊗ I2
42+
I4 = one(SymmetricTensor{4,3})
43+
E4 = 2 * m.G * (I4 - IxI/3) + m.K * IxI
44+
return σ_trial, E4, old
45+
else
46+
rf(x) = residual(x, m, ϵ, old)
47+
Δλ, ∂r∂x, converged = newtonsolve(rf, 0.0)
48+
#=
49+
# Using bisection
50+
Δλ1, ∂r∂x1, converged = newtonsolve(rf, 0.0)
51+
# Workaround to avoid Δλ and ∂r∂x as `Core.Box`:ed
52+
Δλ, ∂r∂x = if !converged || Δλ1 < 0
53+
Δλ2 = bisect_solve(rf, ϕ_trial / m.G)
54+
∂r∂x2 = ForwardDiff.derivative(rf, Δλ2)
55+
(Δλ2, ∂r∂x2)
56+
else
57+
(Δλ1, ∂r∂x1)
58+
end
59+
converged = true
60+
=#
61+
62+
if converged
63+
∂σ∂ϵ = gradient(e -> calculate_sigma(Δλ, m, e, old), ϵ)
64+
∂σ∂x, σ = gradient(x -> calculate_sigma(x, m, ϵ, old), Δλ, :all)
65+
∂r∂ϵ = gradient(e -> residual(Δλ, m, e, old), ϵ)
66+
# drdϵ = 0 = ∂r∂ϵ + ∂r∂x * dxdϵ => dxdϵ = - inv(∂r∂x) * ∂r∂ϵ
67+
dσdϵ = ∂σ∂ϵ - (∂σ∂x / ∂r∂x) ⊗ ∂r∂ϵ
68+
69+
σdev, β, κ = calculate_evolution(Δλ, m, ϵ, old)
70+
ϵp = old.ϵp + Δλ * (3/2) * (σdev - β) / (m.Y0 + κ)
71+
72+
return σ, dσdϵ, SimplePlasticState(ϵp, β, κ)
73+
74+
else
75+
throw(MMB.NoLocalConvergence("$(typeof(m)): newtonsolve! didn't converge, ϵ = ", ϵ))
76+
end
77+
end
78+
end
79+
80+
function calculate_κ(Δλ, m::SimplePlastic, old::SimplePlasticState)
81+
return m.κ∞ * (old.κ + Δλ * m.Hiso) / (m.κ∞ + Δλ * m.Hiso)
82+
end
83+
84+
function dβ_fun(Δλ, m::SimplePlastic, κ)
85+
return m.Y0 + κ + Δλ * m.Hkin * (m.β∞ + m.Y0 + κ) / m.β∞
86+
end
87+
88+
function calculate_sigma(Δλ, m::SimplePlastic, ϵ::SymmetricTensor, old::SimplePlasticState)
89+
σdev, _, _ = calculate_evolution(Δλ, m, ϵ, old)
90+
return σdev + 3 * m.K * vol(ϵ)
91+
end
92+
93+
function calculate_evolution(Δλ, m::SimplePlastic, ϵ::SymmetricTensor, old::SimplePlasticState)
94+
κ = calculate_κ(Δλ, m, old)
95+
dβ = dβ_fun(Δλ, m, κ)
96+
Y = m.Y0 + κ
97+
σdev = (2 * m.G * dev(ϵ) - 2 * m.G * (old.ϵp - (3 * Δλ / (2 * Y)) * (old.β * Y / dβ))) / (
98+
1 + (3 * m.G * Δλ / Y) * (1 - Δλ * m.Hkin / dβ) )
99+
β = (old.β * Y + Δλ * m.Hkin * σdev) / dβ
100+
return σdev, β, κ
101+
end
102+
103+
function residual(Δλ, m::SimplePlastic, ϵ::SymmetricTensor, old::SimplePlasticState)
104+
σdev, β, κ = calculate_evolution(Δλ, m, ϵ, old)
105+
Φ = vonmises(σdev - β) - (m.Y0 + κ)
106+
return Φ
107+
end
108+
109+
#=
110+
# Extra stability, but doesn't seem to be hit in normal cases.
111+
# Therefore, excluded as it is not tested.
112+
function bisect_solve(rf, x0::Number; maxiter = 200, tol = 1e-6)
113+
x_left = zero(x0)
114+
x_right = x0
115+
r_left = rf(x_left)
116+
@show r_left
117+
@assert r_left < 0 # No need for premature generalization...
118+
119+
# Move right untill r_right is above zero
120+
r_right = rf(x0)
121+
iter = 0
122+
while r_right < 0
123+
x_left = x_right
124+
x_right = 2*x_right
125+
r_right = rf(x_right)
126+
println(iter, ": ", (r_right, x_right), " (move right)")
127+
iter += 1
128+
iter > maxiter && error("Didn't converge (could not find positive value)")
129+
end
130+
131+
# Bisection
132+
for i in iter:maxiter
133+
# 0 = r_center ≈ r_left + (r_right - r_left) * (x_center - x_left) / (x_right - x_left)
134+
x_center = x_left - r_left * (x_right - x_left) / (r_right - r_left)
135+
r_center = rf(x_center)
136+
println(iter, ": ", (r_center, x_center), " (bisect)")
137+
abs(r_center) < tol && return x_center
138+
if r_center < 0
139+
x_left = x_center
140+
r_left = r_center
141+
else
142+
x_right = x_center
143+
r_right = r_center
144+
end
145+
end
146+
error("Did not find sufficiently small residual")
147+
end
148+
=#

‎test/test_plastic.jl‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,3 +67,63 @@
6767
@test contains(show_str, "ArmstrongFrederick")
6868
@test contains(show_str, "Norton overstress")
6969
end
70+
71+
@testset "SimplePlastic" begin
72+
# Run an example, compare with Plastic
73+
E = 210.e3; ν = 0.3
74+
Y0 = 150.0
75+
Hkin = 1 * 10e3; β∞ = 25.0
76+
Hiso = 1 * 13e3; κ∞ = 35.0
77+
G = E/(2 * (1 + ν))
78+
K = E / (3 * (1 - 2 * ν))
79+
elastic = LinearElastic(;E, ν)
80+
kinematic = ArmstrongFrederick(;Hkin, β∞)
81+
isotropic = Voce(;Hiso, κ∞)
82+
m0 = Plastic(;elastic, yield = Y0, kinematic, isotropic)
83+
m0_simple = SimplePlastic(; G, K, Y0, Hiso, κ∞, Hkin, β∞)
84+
m0_simple_from_plastic = SimplePlastic(m0)
85+
for key in fieldnames(typeof(m0_simple))
86+
@test getproperty(m0_simple, key) ≈ getproperty(m0_simple_from_plastic, key)
87+
end
88+
89+
Δϵ21 = 0.01
90+
N = 100
91+
s21, σ, dσdϵ, state, ϵ0 = run_shear(m0, Δϵ21, N)
92+
93+
ϵp_old = state.ϵp
94+
κold = state.κ[1]
95+
βold = state.β[1]
96+
old = MechanicalMaterialModels.SimplePlasticState(ϵp_old, βold, κold)
97+
98+
ϵt = SymmetricTensor{2,3}((i,j) -> i==2 && j==1 ? ϵ0[i, j] + 1e-4 : 0.0)
99+
100+
extras = allocate_differentiation_output(m0)
101+
σ1, dσdϵ1, state1 = material_response(m0, ϵt, state, 0.1, allocate_material_cache(m0), extras)
102+
103+
σ_s, dσdϵ_s, state_s = material_response(m0_simple, ϵt, old, 0.1)
104+
105+
Δλ = extras.X.Δλ
106+
# κ = MechanicalMaterialModels.calculate_κ(Δλ, m0_simple, old)
107+
σdev, β, κ = MechanicalMaterialModels.calculate_evolution(Δλ, m0_simple, ϵt, old)
108+
109+
@test state1.κ[1] ≈ state_s.κ ≈ κ
110+
@test state1.β[1] ≈ state_s.β ≈ β
111+
@test state1.ϵp ≈ state_s.ϵp
112+
@test dev(σ1) ≈ dev(σ_s) ≈ σdev
113+
@test σ1 ≈ σ_s
114+
115+
Δϵ21 = 0.01
116+
N = 100
117+
s21, σ, dσdϵ, state, ϵ = run_shear(m0, Δϵ21, N)
118+
s21_s, σ_s, dσdϵ_s, state_s, ϵ_s = run_shear(m0_simple, Δϵ21, N)
119+
# Using tol=1e-10 in newtonsolve inside material_response makes them pass standard ≈ test...
120+
@test isapprox(s21, s21_s; atol = abs(s21[end]) * 1e-6)
121+
@test σ ≈ σ_s
122+
@test dσdϵ ≈ dσdϵ_s
123+
@test ϵ ≈ ϵ_s
124+
125+
# Test stability
126+
γv = [0.0, 0.5, -0.5, 0.0, -1.0, 1.0]
127+
@test_throws NoLocalConvergence run_shear(m0, γv)
128+
run_shear(m0_simple, γv)
129+
end

‎test/utilities4testing.jl‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,16 @@
1-
function run_shear(m, γmax, numsteps, Δt = nothing; TT=SymmetricTensor)
1+
function run_shear(m, γmax::Number, numsteps::Int, Δt = nothing; TT=SymmetricTensor)
2+
γv = range(0, γmax, numsteps + 1)
3+
return run_shear(m, γv, Δt; TT)
4+
end
5+
6+
function run_shear(m, γv::AbstractVector, Δt = nothing; TT = SymmetricTensor)
27
state = initial_material_state(m)
3-
Δϵ = TT{2,3}((i,j)-> i==2 && j==1 ? γmax/numsteps : zero(γmax))
4-
s21 = zeros(numsteps+1)
8+
s21 = zeros(length(γv))
59
local σ, dσdϵ, ϵ
6-
for i = 1:numsteps
7-
ϵ = i*Δϵ
10+
for (k, ϵ21) in enumerate(γv[2:end])
11+
ϵ = TT{2,3}((i,j)-> i==2 && j==1 ? ϵ21 : zero(ϵ21))
812
σ, dσdϵ, state = material_response(m, ϵ, state, Δt)
9-
s21[i+1] = σ[2,1]
13+
s21[k+1] = σ[2,1]
1014
end
1115
return s21, σ, dσdϵ, state, ϵ
1216
end

0 commit comments

Comments
 (0)