Skip to content

Commit

Permalink
Merge pull request #64 from JuliaFEM/memory
Browse files Browse the repository at this point in the history 
Memory
  • Loading branch information
@jvaara
jvaara authored May 13, 2020
2 parents f5b50aa + 7deec74 commit da04332
Show file tree
Hide file tree
Showing 4 changed files with 354 additions and 0 deletions.
10 changes: 10 additions & 0 deletions src/Materials.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,12 +43,22 @@ function isotropic_elasticity_tensor(lambda, mu)
return jacobian
end

function isotropic_compliance_tensor(lambda, mu)
delta(i,j) = i==j ? 1.0 : 0.0
g(i,j,k,l) = -lambda/(2.0*mu*(3.0*lambda + 2.0*mu))*delta(i,j)*delta(k,l) + 1.0/(4.0*mu)*(delta(i,k)*delta(j,l)+delta(i,l)*delta(j,k))
compliance = SymmetricTensor{4, 3, Float64}(g)
return compliance
end

include("idealplastic.jl")
export IdealPlastic, IdealPlasticDriverState, IdealPlasticParameterState, IdealPlasticVariableState

include("chaboche.jl")
export Chaboche, ChabocheDriverState, ChabocheParameterState, ChabocheVariableState

include("memory.jl")
export Memory, MemoryDriverState, MemoryParameterState, MemoryVariableState

# include("viscoplastic.jl")
# export ViscoPlastic

Expand Down
203 changes: 203 additions & 0 deletions src/memory.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

@with_kw mutable struct MemoryDriverState <: AbstractMaterialState
time :: Float64 = zero(Float64)
strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
end

@with_kw struct MemoryParameterState <: AbstractMaterialState
E :: Float64 = 0.0
nu :: Float64 = 0.0
R0 :: Float64 = 0.0
Kn :: Float64 = 0.0
nn :: Float64 = 0.0
C1 :: Float64 = 0.0
D1 :: Float64 = 0.0
C2 :: Float64 = 0.0
D2 :: Float64 = 0.0
Q0 :: Float64 = 0.0
QM :: Float64 = 0.0
mu :: Float64 = 0.0
b :: Float64 = 0.0
eta :: Float64 = 0.0
m :: Float64 = 0.0
pt :: Float64 = 0.0
xi :: Float64 = 0.0
end

@with_kw struct MemoryVariableState <: AbstractMaterialState
stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
X1 :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
X2 :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
plastic_strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
cumeq :: Float64 = zero(Float64)
R :: Float64 = zero(Float64)
q :: Float64 = zero(Float64)
zeta :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64})
jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64})
end

@with_kw mutable struct Memory <: AbstractMaterial
drivers :: MemoryDriverState = MemoryDriverState()
ddrivers :: MemoryDriverState = MemoryDriverState()
variables :: MemoryVariableState = MemoryVariableState()
variables_new :: MemoryVariableState = MemoryVariableState()
parameters :: MemoryParameterState = MemoryParameterState()
dparameters :: MemoryParameterState = MemoryParameterState()
end

function integrate_material!(material::Memory)
p = material.parameters
v = material.variables
dd = material.ddrivers
d = material.drivers
@unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q0, QM, mu, b, eta, m, pt, xi = p
mu_ = E/(2.0*(1.0+nu))
lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

@unpack strain, time = d
dstrain = dd.strain
dtime = dd.time
@unpack stress, X1, X2, plastic_strain, cumeq, R, q, zeta, jacobian = v


# Elastic trial
jacobian = isotropic_elasticity_tensor(lambda, mu_)
stress += dcontract(jacobian, dstrain)
seff = stress - X1 - X2
seff_dev = dev(seff)
f = sqrt(1.5)*norm(seff_dev) - (R0 + R)
if f > 0.0
g! = create_nonlinear_system_of_equations(material)
x0 = [tovoigt(stress); R; tovoigt(X1); tovoigt(X2)]
F = similar(x0)

res = nlsolve(g!, x0; autodiff = :forward) # Explicit update to memory-surface
res.f_converged || error("Nonlinear system of equations with explicit surface did not converge!")
x = res.zero
stress = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[1:6])
R = x[7]
X1 = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[8:13])
X2 = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[14:19])

seff = stress - X1 - X2
seff_dev = dev(seff)
f = sqrt(1.5)*norm(seff_dev) - (R0 + R)
dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
dp = dotp*dtime
n = sqrt(1.5)*seff_dev/norm(seff_dev)

# Update plastic strain
plastic_strain += dp*n
cumeq += dp

# Strain memory - explicit update
JF = sqrt(1.5)*norm(dev(plastic_strain - zeta))
FF = 2.0/3.0*JF - q
if FF > 0.0
nF = 1.5*dev(plastic_strain - zeta)/JF
nnF = dcontract(n, nF)
if nnF>0
q += 2.0/3.0*eta*nnF*dp
zeta += 2.0/3.0*(1.0 - eta)*nnF*nF*dp
end
else
# Memory evanescence term
if cumeq>=pt
q += -xi*q^m*dp
end
end

# Compute Jacobian
function residuals(x)
F = similar(x)
g!(F, x)
return F
end
drdx = ForwardDiff.jacobian(residuals, x)
drde = zeros((length(x),6))
drde[1:6, 1:6] = -tovoigt(jacobian)
jacobian = fromvoigt(SymmetricTensor{4,3}, (drdx\drde)[1:6, 1:6])
end
variables_new = MemoryVariableState(stress = stress,
X1 = X1,
X2 = X2,
R = R,
plastic_strain = plastic_strain,
cumeq = cumeq,
q = q,
zeta = zeta,
jacobian = jacobian)
material.variables_new = variables_new
end

function create_nonlinear_system_of_equations(material::Memory)
p = material.parameters
v = material.variables
dd = material.ddrivers
d = material.drivers
@unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q0, QM, mu, b, eta, m, pt, xi = p
mu_ = E/(2.0*(1.0+nu))
lambda = E*nu/((1.0+nu)*(1.0-2.0*nu))

@unpack strain, time = d
dstrain = dd.strain
dtime = dd.time
@unpack stress, X1, X2, plastic_strain, cumeq, R, q, zeta, jacobian = v

function g!(F, x::Vector{T}) where {T} # Explicit update of memory surface
jacobian = isotropic_elasticity_tensor(lambda, mu_)
stress_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[1:6])
R_ = x[7]
X1_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[8:13])
X2_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[14:19])

seff = stress_ - X1_ - X2_
seff_dev = dev(seff)
f = sqrt(1.5)*norm(seff_dev) - (R0 + R_)

dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
dp = dotp*dtime
n = sqrt(1.5)*seff_dev/norm(seff_dev)
dstrain_plastic = dp*n
plastic_strain_ = plastic_strain + dstrain_plastic
# Strain memory - explicit update
JF = sqrt(1.5)*norm(dev(plastic_strain_ - zeta))
FF = 2.0/3.0*JF - q
if FF > 0.0
nF = 1.5*dev(plastic_strain_ - zeta)/JF
nnF = dcontract(n, nF)
if nnF>0
q_ = q + 2.0/3.0*eta*nnF*dp
zeta_ = zeta + 2.0/3.0*(1.0 - eta)*nnF*nF*dp
else
q_ = q
zeta_ = zeta
end
else
# Memory evanescence term
p_ = cumeq + dp
if p_>pt
q_ = q - xi*q^m*dp
else
q_ = q
end
zeta_ = zeta
end

tovoigt!(view(F, 1:6), stress - stress_ + dcontract(jacobian, dstrain - dstrain_plastic))
F[7] = R - R_ + b*((QM + (Q0 - QM)*exp(-2.0*mu*q_))-R_)*dp
if isapprox(C1, 0.0)
tovoigt!(view(F, 8:13), X1 - X1_)
else
tovoigt!(view(F, 8:13), X1 - X1_ + 2.0/3.0*C1*dp*(n - 1.5*D1/C1*X1_))
end
if isapprox(C2, 0.0)
tovoigt!(view(F, 14:19), X2 - X2_)
else
tovoigt!(view(F, 14:19), X2 - X2_ + 2.0/3.0*C2*dp*(n - 1.5*D2/C2*X2_))
end
end
return g!
end
3 changes: 3 additions & 0 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,9 @@ using Materials, Test
@testset "test biaxial increment" begin
include("test_biaxial_increment.jl")
end
@testset "test memory" begin
include("test_memory.jl")
end
@testset "test stress-driven uniaxial increment" begin
include("test_stress_driven_uniaxial_increment.jl")
end
Expand Down
138 changes: 138 additions & 0 deletions test/test_memory.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
using Test, Tensors

parameters = MemoryParameterState(E = 200.0e3,
nu = 0.3,
R0 = 100.0,
Kn = 20.0,
nn = 3.0,
C1 = 10000.0,
D1 = 100.0,
C2 = 50000.0,
D2 = 1000.0,
Q0 = 100.0,
#QM = 300.0,
QM = 500.0,
mu = 100.0,
b = 30.0,
#eta = 0.0,
# eta = 0.05,
eta = 0.5,
m = 0.5,
pt = 0.0,
xi = 0.3)
# xi = 0.0)
mat = Memory(parameters=parameters)

times = [copy(mat.drivers.time)]
stresses = [copy(tovoigt(mat.variables.stress))]
strains = [copy(tovoigt(mat.drivers.strain; offdiagscale=2.0))]
plastic_strains = [copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0))]
cumeqs = [copy(mat.variables.cumeq)]
qs = [copy(mat.variables.q)]
Rs = [copy(mat.variables.R)]
zetas = [copy(tovoigt(mat.variables.zeta; offdiagscale=2.0))]

n_cycles = 30
ppc = 40
t = range(0.0; stop=Float64(n_cycles), length=n_cycles*ppc+1)
dtime = t[end]/(length(t)-1)

# Amplitude 1
ea = 0.003
strains11 = ea*sin.(2*pi*t)
for dstrain11 in diff(strains11)
uniaxial_increment!(mat, dstrain11, dtime)
update_material!(mat)
push!(times, mat.drivers.time)
push!(stresses, copy(tovoigt(mat.variables.stress)))
push!(strains, copy(tovoigt(mat.drivers.strain; offdiagscale=2.0)))
push!(plastic_strains, copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0)))
push!(cumeqs, copy(mat.variables.cumeq))
push!(qs, copy(mat.variables.q))
push!(Rs, copy(mat.variables.R))
push!(zetas, copy(tovoigt(mat.variables.zeta; offdiagscale=2.0)))
end
R1 = copy(Rs[end])

# Amplitude 2
ea = 0.005
strains11 = ea*sin.(2*pi*t)
for dstrain11 in diff(strains11)
uniaxial_increment!(mat, dstrain11, dtime)
update_material!(mat)
push!(times, mat.drivers.time)
push!(stresses, copy(tovoigt(mat.variables.stress)))
push!(strains, copy(tovoigt(mat.drivers.strain; offdiagscale=2.0)))
push!(plastic_strains, copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0)))
push!(cumeqs, copy(mat.variables.cumeq))
push!(qs, copy(mat.variables.q))
push!(Rs, copy(mat.variables.R))
push!(zetas, copy(tovoigt(mat.variables.zeta; offdiagscale=2.0)))
end
R2 = copy(Rs[end])

# Amplitude 3
ea = 0.007
strains11 = ea*sin.(2*pi*t)
for dstrain11 in diff(strains11)
uniaxial_increment!(mat, dstrain11, dtime)
update_material!(mat)
push!(times, mat.drivers.time)
push!(stresses, copy(tovoigt(mat.variables.stress)))
push!(strains, copy(tovoigt(mat.drivers.strain; offdiagscale=2.0)))
push!(plastic_strains, copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0)))
push!(cumeqs, copy(mat.variables.cumeq))
push!(qs, copy(mat.variables.q))
push!(Rs, copy(mat.variables.R))
push!(zetas, copy(tovoigt(mat.variables.zeta; offdiagscale=2.0)))
end
R3 = copy(Rs[end])

# Amplitude 4 - Test evanescence
ea = 0.003
strains11 = ea*sin.(2*pi*t)
for dstrain11 in diff(strains11)
uniaxial_increment!(mat, dstrain11, dtime)
update_material!(mat)
push!(times, mat.drivers.time)
push!(stresses, copy(tovoigt(mat.variables.stress)))
push!(strains, copy(tovoigt(mat.drivers.strain; offdiagscale=2.0)))
push!(plastic_strains, copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0)))
push!(cumeqs, copy(mat.variables.cumeq))
push!(qs, copy(mat.variables.q))
push!(Rs, copy(mat.variables.R))
push!(zetas, copy(tovoigt(mat.variables.zeta; offdiagscale=2.0)))
end

for dstrain11 in diff(strains11)
uniaxial_increment!(mat, dstrain11, dtime)
update_material!(mat)
push!(times, mat.drivers.time)
push!(stresses, copy(tovoigt(mat.variables.stress)))
push!(strains, copy(tovoigt(mat.drivers.strain; offdiagscale=2.0)))
push!(plastic_strains, copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0)))
push!(cumeqs, copy(mat.variables.cumeq))
push!(qs, copy(mat.variables.q))
push!(Rs, copy(mat.variables.R))
push!(zetas, copy(tovoigt(mat.variables.zeta; offdiagscale=2.0)))
end

for dstrain11 in diff(strains11)
uniaxial_increment!(mat, dstrain11, dtime)
update_material!(mat)
push!(times, mat.drivers.time)
push!(stresses, copy(tovoigt(mat.variables.stress)))
push!(strains, copy(tovoigt(mat.drivers.strain; offdiagscale=2.0)))
push!(plastic_strains, copy(tovoigt(mat.variables.plastic_strain; offdiagscale=2.0)))
push!(cumeqs, copy(mat.variables.cumeq))
push!(qs, copy(mat.variables.q))
push!(Rs, copy(mat.variables.R))
push!(zetas, copy(tovoigt(mat.variables.zeta; offdiagscale=2.0)))
end
R4 = copy(Rs[end])
@test R2>R1
@test R3>R2
@test R4<R3
@test isapprox(R1, R4; atol=1.0)

0 comments on commit da04332

Please sign in to comment.