-
Notifications
You must be signed in to change notification settings - Fork 9
/
DSA.jl
309 lines (262 loc) · 10.7 KB
/
DSA.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
module DSAModule
using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters
import ..AbstractMaterial, ..AbstractMaterialState
import ..Utilities: Symm2, Symm4, isotropic_elasticity_tensor, lame, debang
import ..integrate_material! # for method extension
# parametrically polymorphic for any type representing ℝ
export GenericDSA, GenericDSADriverState, GenericDSAParameterState, GenericDSAVariableState
# specialization for Float64
export DSA, DSADriverState, DSAParameterState, DSAVariableState
@with_kw mutable struct GenericDSADriverState{T <: Real} <: AbstractMaterialState
time::T = zero(T)
strain::Symm2{T} = zero(Symm2{T})
end
# TODO: complete this docstring
"""Parameter state for DSA (dynamic strain aging) material.
This is similar to the Chaboche model, but with additional static recovery terms.
`E`: Young's modulus
`nu`: Poisson's ratio
`R0`: initial yield strength
`Kn`: plasticity multiplier divisor (drag stress)
`nn`: plasticity multiplier exponent
`C1`, `D1`: parameters governing behavior of backstress X1
`C2`, `D2`: parameters governing behavior of backstress X2
`Q`: shift parameter for yield strength evolution
`b`: multiplier for yield strength evolution
`w`: ???
`P1`, `P2`: ???
`m`: ???
`m1`, `m2`: ???
`M1`, `M2`: ???
`ba`: ???
`xi`: ???
"""
@with_kw struct GenericDSAParameterState{T <: Real} <: AbstractMaterialState
E::T = 0.0
nu::T = 0.0
R0::T = 0.0
Kn::T = 0.0
nn::T = 0.0
C1::T = 0.0
D1::T = 0.0
C2::T = 0.0
D2::T = 0.0
Q::T = 0.0
b::T = 0.0
w::T = 0.0
P1::T = 0.0
P2::T = 0.0
m::T = 0.0
m1::T = 0.0
m2::T = 0.0
M1::T = 0.0
M2::T = 0.0
ba::T = 0.0
xi::T = 0.0
end
# TODO: complete this docstring
"""Problem state for Chaboche material.
`stress`: stress tensor
`X1`: backstress 1
`X2`: backstress 2
`plastic_strain`: plastic part of strain tensor
`cumeq`: cumulative equivalent plastic strain (scalar, ≥ 0)
`R`: yield strength
`jacobian`: ∂σij/∂εkl
`ta`: ???
`Ra`: ???
"""
@with_kw struct GenericDSAVariableState{T <: Real} <: AbstractMaterialState
stress::Symm2{T} = zero(Symm2{T})
X1::Symm2{T} = zero(Symm2{T})
X2::Symm2{T} = zero(Symm2{T})
plastic_strain::Symm2{T} = zero(Symm2{T})
cumeq::T = zero(T)
R::T = zero(T)
jacobian::Symm4{T} = zero(Symm4{T})
ta::T = zero(T)
Ra::T = zero(T)
end
# TODO: Does this eventually need a {T}?
@with_kw struct DSAOptions <: AbstractMaterialState
nlsolve_method::Symbol = :trust_region
end
@with_kw mutable struct GenericDSA{T <: Real} <: AbstractMaterial
drivers::GenericDSADriverState{T} = GenericDSADriverState{T}()
ddrivers::GenericDSADriverState{T} = GenericDSADriverState{T}()
variables::GenericDSAVariableState{T} = GenericDSAVariableState{T}()
variables_new::GenericDSAVariableState{T} = GenericDSAVariableState{T}()
parameters::GenericDSAParameterState{T} = GenericDSAParameterState{T}()
dparameters::GenericDSAParameterState{T} = GenericDSAParameterState{T}()
options::DSAOptions = DSAOptions()
end
DSADriverState = GenericDSADriverState{Float64}
DSAParameterState = GenericDSAParameterState{Float64}
DSAVariableState = GenericDSAVariableState{Float64}
DSA = GenericDSA{Float64}
"""
state_to_vector(sigma::U, R::T, X1::U, X2::U, ta::T, Ra::T) where U <: Symm2{T} where T <: Real
Adaptor for `nlsolve`. Marshal the problem state into a `Vector`.
"""
function state_to_vector(sigma::U, R::T, X1::U, X2::U, ta::T, Ra::T) where U <: Symm2{T} where T <: Real
return [tovoigt(sigma); R; tovoigt(X1); tovoigt(X2); ta; Ra]::Vector{T}
end
"""
state_from_vector(x::AbstractVector{<:Real})
Adaptor for `nlsolve`. Unmarshal the problem state from a `Vector`.
"""
function state_from_vector(x::AbstractVector{T}) where T <: Real
sigma::Symm2{T} = fromvoigt(Symm2{T}, @view x[1:6])
R::T = x[7]
X1::Symm2{T} = fromvoigt(Symm2{T}, @view x[8:13])
X2::Symm2{T} = fromvoigt(Symm2{T}, @view x[14:19])
ta::T = x[20]
Ra::T = x[21]
return sigma, R, X1, X2, ta, Ra
end
"""
integrate_material!(material::GenericDSA{T}) where T <: Real
Material model with dynamic strain aging (DSA).
This is similar to the Chaboche material with two backstresses, with both
kinematic and isotropic hardening, but this model also features static recovery
terms.
"""
function integrate_material!(material::GenericDSA{T}) where T <: Real
p = material.parameters
v = material.variables
dd = material.ddrivers
d = material.drivers
@unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b, w, P1, P2, m, m1, m2, M1, M2, ba, xi = p
lambda, mu = lame(E, nu)
@unpack strain, time = d
dstrain = dd.strain
dtime = dd.time
@unpack stress, X1, X2, plastic_strain, cumeq, R, jacobian, ta, Ra = v
# elastic part
jacobian = isotropic_elasticity_tensor(lambda, mu)
stress += dcontract(jacobian, dstrain)
# resulting deviatoric plastic stress (accounting for backstresses Xm)
seff_dev = dev(stress - X1 - X2)
# von Mises yield function
f = sqrt(1.5)*norm(seff_dev) - (R0 + R + (1 - xi) * Ra) # using elastic trial problem state
if f > 0.0
g! = create_nonlinear_system_of_equations(material)
x0 = state_to_vector(stress, R, X1, X2, ta + dtime, Ra)
res = nlsolve(g!, x0; method=material.options.nlsolve_method, autodiff = :forward)
converged(res) || error("Nonlinear system of equations did not converge!")
x = res.zero
stress, R, X1, X2, ta, Ra = state_from_vector(x)
# using the new problem state
seff_dev = dev(stress - X1 - X2)
f = sqrt(1.5)*norm(seff_dev) - (R0 + R + (1 - xi) * Ra)
dotp = ((f >= 0.0 ? f : 0.0) / (Kn + xi * Ra))^nn
dp = dotp*dtime
n = sqrt(1.5)*seff_dev/norm(seff_dev)
plastic_strain += dp*n
cumeq += dp
# Compute the new Jacobian, accounting for the plastic contribution.
drdx = ForwardDiff.jacobian(debang(g!), x)
drde = zeros((length(x), 6))
drde[1:6, 1:6] = -tovoigt(jacobian) # elastic Jacobian. Follows from the defn. of g!.
jacobian = fromvoigt(Symm4, (drdx\drde)[1:6, 1:6])
else
ta += dtime
end
variables_new = GenericDSAVariableState{T}(stress = stress,
X1 = X1,
X2 = X2,
R = R,
plastic_strain = plastic_strain,
cumeq = cumeq,
jacobian = jacobian,
ta = ta,
Ra = Ra)
material.variables_new = variables_new
return nothing
end
"""
create_nonlinear_system_of_equations(material::GenericDSA{T}) where T <: Real
Create and return an instance of the equation system for the incremental form of
the evolution equations of the DSA material.
Used internally for computing the plastic contribution in `integrate_material!`.
The input `material` represents the problem state at the end of the previous
timestep. The created equation system will hold its own copy of that state.
The equation system is represented as a mutating function `g!` that computes the
residual:
```julia
g!(F::V, x::V) where V <: AbstractVector{<:Real}
```
Both `F` (output) and `x` (input) are length-21 vectors containing
[sigma, R, X1, X2, ta, Ra], in that order. The tensor quantities
sigma, X1, X2 are encoded in Voigt format.
The function `g!` is intended to be handed over to `nlsolve`.
"""
function create_nonlinear_system_of_equations(material::GenericDSA{T}) where T <: Real
p = material.parameters
v = material.variables
dd = material.ddrivers
d = material.drivers
@unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b, w, P1, P2, m, m1, m2, M1, M2, ba, xi = p
lambda, mu = lame(E, nu)
# Old problem state (i.e. the problem state at the time when this equation
# system instance was created).
#
# Note this does not include the elastic trial; this is the state at the
# end of the previous timestep.
@unpack strain, time = d
dstrain = dd.strain
dtime = dd.time
@unpack stress, X1, X2, plastic_strain, cumeq, R, ta, Ra = v
jacobian = isotropic_elasticity_tensor(lambda, mu)
# Compute the residual. F is output, x is filled by NLsolve.
# The solution is x = x* such that g(x*) = 0.
function g!(F::V, x::V) where V <: AbstractVector{<:Real}
stress_new, R_new, X1_new, X2_new, ta_new, Ra_new = state_from_vector(x) # tentative new values from nlsolve
seff_dev = dev(stress_new - X1_new - X2_new)
f = sqrt(1.5)*norm(seff_dev) - (R0 + R_new + (1 - xi) * Ra_new)
dotp = ((f >= 0.0 ? f : 0.0) / (Kn + xi * Ra_new))^nn
dp = dotp*dtime
n = sqrt(1.5)*seff_dev/norm(seff_dev)
# The equations are written in an incremental form.
# TODO: multiply the equations by -1 to make them easier to understand in the context of the rest of the model.
dstrain_plastic = dp*n
dstrain_elastic = dstrain - dstrain_plastic
tovoigt!(view(F, 1:6), stress - stress_new + dcontract(jacobian, dstrain_elastic))
F[7] = R - R_new + b*(Q - R_new)*dp
# TODO: figure out what's going on here. This doesn't reduce consistently at the limit C1 -> 0.
if isapprox(C1, 0.0)
tovoigt!(view(F, 8:13), X1 - X1_new)
else
ndX1_new = norm(dev(X1_new)) # Static recovery component
if iszero(ndX1_new)
JX1_new = 0.0
else
JX1_new = sqrt(1.5) * ndX1_new
end
# JX1_new = sqrt(1.5)*norm(dev(X1_new))
sr1_new = (JX1_new^(m1 - 1) * X1_new) / (M1^m1) # static recovery
tovoigt!(view(F, 8:13), X1 - X1_new + dp*(2.0/3.0*C1*n - D1*X1_new) - dtime*sr1_new)
end
if isapprox(C2, 0.0)
tovoigt!(view(F, 14:19), X2 - X2_new)
else
ndX2_new = norm(dev(X2_new)) # Static recovery component
if iszero(ndX2_new)
JX2_new = 0.0
else
JX2_new = sqrt(1.5) * ndX2_new
end
# JX2_new = sqrt(1.5)*norm(dev(X2_new))
sr2_new = (JX2_new^(m2 - 1) * X2_new) / (M2^m2) # static recovery
tovoigt!(view(F, 14:19), X2 - X2_new + dp*(2.0/3.0*C2*n - D2*X2_new) - dtime*sr2_new)
end
Ras = P1 * (1.0 - exp(-P2 * ta_new^m))
F[20] = ta - ta_new + dtime - (ta_new / w)*dp
F[21] = Ra - Ra_new + ba*(Ras - Ra_new)*dp
return nothing
end
return g!
end
end