diff --git a/examples/abstractmaterial_full_workflow.jl b/examples/abstractmaterial_full_workflow.jl index c6214fc..0f6f04f 100644 --- a/examples/abstractmaterial_full_workflow.jl +++ b/examples/abstractmaterial_full_workflow.jl @@ -14,7 +14,7 @@ end @with_kw mutable struct ChabocheDriverState <: AbstractMaterialState time :: Float64 = zero(Float64) - strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) + strain :: Symm2 = zero(Symm2{Float64}) end @with_kw struct ChabocheParameterState <: AbstractMaterialState @@ -32,13 +32,13 @@ end end @with_kw struct ChabocheVariableState <: 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}) + stress :: Symm2 = zero(Symm2{Float64}) + X1 :: Symm2 = zero(Symm2{Float64}) + X2 :: Symm2 = zero(Symm2{Float64}) + plastic_strain :: Symm2 = zero(Symm2{Float64}) cumeq :: Float64 = zero(Float64) R :: Float64 = zero(Float64) - jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64}) + jacobian :: Symm4 = zero(Symm4{Float64}) end @with_kw mutable struct Chaboche <: AbstractMaterial @@ -50,20 +50,20 @@ end dparameters :: ChabocheParameterState = ChabocheParameterState() end -@with_kw mutable struct IdealPlasticDriverState <: AbstractMaterialState +@with_kw mutable struct PerfectPlasticDriverState <: AbstractMaterialState time :: Float64 = zero(Float64) - strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) + strain :: Symm2 = zero(Symm2{Float64}) end -@with_kw struct IdealPlasticParameterState <: AbstractMaterialState +@with_kw struct PerfectPlasticParameterState <: AbstractMaterialState youngs_modulus :: Float64 = zero(Float64) poissons_ratio :: Float64 = zero(Float64) yield_stress :: Float64 = zero(Float64) end -@with_kw struct IdealPlasticVariableState <: AbstractMaterialState - stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) - plastic_strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) +@with_kw struct PerfectPlasticVariableState <: AbstractMaterialState + stress :: Symm2 = zero(Symm2{Float64}) + plastic_strain :: Symm2 = zero(Symm2{Float64}) cumeq :: Float64 = zero(Float64) end @@ -92,23 +92,23 @@ end # dparameters :: ParameterState{M} = ParameterState{M}() # end -@with_kw mutable struct IdealPlastic <: AbstractMaterial - drivers :: IdealPlasticDriverState = IdealPlasticDriverState() - ddrivers :: IdealPlasticDriverState = IdealPlasticDriverState() - variables :: IdealPlasticVariableState = IdealPlasticVariableState() - variables_new :: IdealPlasticVariableState = IdealPlasticVariableState() - parameters :: IdealPlasticParameterState = IdealPlasticParameterState() - dparameters :: IdealPlasticParameterState = IdealPlasticParameterState() +@with_kw mutable struct PerfectPlastic <: AbstractMaterial + drivers :: PerfectPlasticDriverState = PerfectPlasticDriverState() + ddrivers :: PerfectPlasticDriverState = PerfectPlasticDriverState() + variables :: PerfectPlasticVariableState = PerfectPlasticVariableState() + variables_new :: PerfectPlasticVariableState = PerfectPlasticVariableState() + parameters :: PerfectPlasticParameterState = PerfectPlasticParameterState() + dparameters :: PerfectPlasticParameterState = PerfectPlasticParameterState() end mat = Chaboche() -mat2 = IdealPlastic() +mat2 = PerfectPlastic() function isotropic_elasticity_tensor(lambda, mu) delta(i,j) = i==j ? 1.0 : 0.0 g(i,j,k,l) = lambda*delta(i,j)*delta(k,l) + mu*(delta(i,k)*delta(j,l)+delta(i,l)*delta(j,k)) - jacobian = SymmetricTensor{4, 3, Float64}(g) + jacobian = Symm4{Float64}(g) return jacobian end @@ -140,10 +140,10 @@ function integrate_material!(material::Chaboche) x = res.zero res.f_converged || error("Nonlinear system of equations did not converge!") - stress = fromvoigt(SymmetricTensor{2,3,Float64}, @view x[1:6]) + stress = fromvoigt(Symm2{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]) + X1 = fromvoigt(Symm2{Float64}, @view x[8:13]) + X2 = fromvoigt(Symm2{Float64}, @view x[14:19]) seff = stress - X1 - X2 seff_dev = dev(seff) f = sqrt(1.5)*norm(seff_dev) - (R0 + R) @@ -162,7 +162,7 @@ function integrate_material!(material::Chaboche) 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]) + jacobian = fromvoigt(Symm4, (drdx\drde)[1:6, 1:6]) end variables_new = ChabocheVariableState(stress = stress, X1 = X1, @@ -190,10 +190,10 @@ function create_nonlinear_system_of_equations(material::Chaboche) function g!(F, x::Vector{T}) where {T} # System of non-linear equations jacobian = isotropic_elasticity_tensor(lambda, mu) - stress_ = fromvoigt(SymmetricTensor{2,3,T}, @view x[1:6]) + stress_ = fromvoigt(Symm2{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]) + X1_ = fromvoigt(Symm2{T}, @view x[8:13]) + X2_ = fromvoigt(Symm2{T}, @view x[14:19]) seff = stress_ - X1_ - X2_ seff_dev = dev(seff) @@ -232,7 +232,7 @@ function simple_integration_test() Q = 50.0, b = 0.1) - dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) + dstrain_dtime = fromvoigt(Symm2{Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime) chabmat = Chaboche(parameters = parameters, ddrivers = ddrivers) @info "time = $(chabmat.drivers.time), stress = $(chabmat.variables.stress)" @@ -286,7 +286,7 @@ function test_chaboche() s33s = [chabmat.variables.stress[3,3]] for i=2:length(ts) dtime = ts[i]-ts[i-1] - dstrain = fromvoigt(SymmetricTensor{2,3,Float64}, strains[i]-strains[i-1]; offdiagscale=2.0) + dstrain = fromvoigt(Symm2{Float64}, strains[i]-strains[i-1]; offdiagscale=2.0) chabmat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain) integrate_material!(chabmat) update!(chabmat) @@ -314,7 +314,7 @@ function simple_integration_test_fd_tangent() Q = 50.0, b = 0.1) - dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) + dstrain_dtime = fromvoigt(Symm2{Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime) chabmat = Chaboche(parameters = parameters, ddrivers = ddrivers) @@ -394,7 +394,7 @@ function simple_integration_test_fd_tangent2() Q = 50.0, b = 0.1) - dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) + dstrain_dtime = fromvoigt(Symm2{Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) ddrivers = ChabocheDriverState(time = 0.25, strain = 0.25*dstrain_dtime) chabmat = Chaboche(parameters = parameters, ddrivers = ddrivers) integrate_material!(chabmat) diff --git a/examples/abstractmaterial_testing.jl b/examples/abstractmaterial_testing.jl index 004b378..7d5bbe4 100644 --- a/examples/abstractmaterial_testing.jl +++ b/examples/abstractmaterial_testing.jl @@ -67,9 +67,9 @@ function bench_tensor() end function bench_symtensor() # Random walk test - var = Variable(SymmetricTensor{2, 3}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])) + var = Variable(Symm2([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])) for i in 1:N - var.dvalue += randn(SymmetricTensor{2,3,Float64}) + var.dvalue += randn(Symm2{Float64}) update!(var) end end @@ -106,21 +106,21 @@ function update!(state::T) where {T<:AbstractVariableState} end function bench_chaboche_style_variablestate() - stress = Variable(zero(SymmetricTensor{2,3})) - strain = Variable(zero(SymmetricTensor{2,3})) - backstress1 = Variable(zero(SymmetricTensor{2,3})) - backstress2 = Variable(zero(SymmetricTensor{2,3})) - plastic_strain = Variable(zero(SymmetricTensor{2,3})) + stress = Variable(zero(Symm2)) + strain = Variable(zero(Symm2)) + backstress1 = Variable(zero(Symm2)) + backstress2 = Variable(zero(Symm2)) + plastic_strain = Variable(zero(Symm2)) cumeq = Variable(0.0) R = Variable(0.0) state = VariableState(stress, strain, backstress1, backstress2, plastic_strain, cumeq, R) for i in 1:N - state.stress.dvalue = randn(SymmetricTensor{2,3}) - state.strain.dvalue = randn(SymmetricTensor{2,3}) - state.backstress1.dvalue = randn(SymmetricTensor{2,3}) - state.backstress2.dvalue = randn(SymmetricTensor{2,3}) - state.plastic_strain.dvalue = randn(SymmetricTensor{2,3}) + state.stress.dvalue = randn(Symm2) + state.strain.dvalue = randn(Symm2) + state.backstress1.dvalue = randn(Symm2) + state.backstress2.dvalue = randn(Symm2) + state.plastic_strain.dvalue = randn(Symm2) state.cumeq.dvalue = norm(state.plastic_strain.dvalue) state.R.dvalue = randn() update!(state) diff --git a/examples/abstractmaterial_testing_states.jl b/examples/abstractmaterial_testing_states.jl index 959ee3b..473c4e0 100644 --- a/examples/abstractmaterial_testing_states.jl +++ b/examples/abstractmaterial_testing_states.jl @@ -9,23 +9,23 @@ abstract type AbstractMaterialState end end struct SomeState <: AbstractMaterialState - stress::SymmetricTensor{2,3,Float64} + stress::Symm2{Float64} end -state = SomeState(SymmetricTensor{2,3,Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])) +state = SomeState(Symm2{Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])) N = 1000 function bench_state(N) - state = SomeState(SymmetricTensor{2,3,Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])) + state = SomeState(Symm2{Float64}([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])) for i in 1:N - dstate = SomeState(randn(SymmetricTensor{2,3,Float64})) + dstate = SomeState(randn(Symm2{Float64})) state = state + dstate end return state end -# println("Benchmark State{SymmetricTensor{2,3,Float64}}") +# println("Benchmark State{Symm2{Float64}}") # @btime bench_state(N) struct AnotherState <: AbstractMaterialState @@ -38,21 +38,21 @@ struct AnotherState <: AbstractMaterialState R::Float64 end function bench_chaboche_style_state(N) - stress = zero(SymmetricTensor{2,3}) - strain = zero(SymmetricTensor{2,3}) - backstress1 = zero(SymmetricTensor{2,3}) - backstress2 = zero(SymmetricTensor{2,3}) - plastic_strain = zero(SymmetricTensor{2,3}) + stress = zero(Symm2) + strain = zero(Symm2) + backstress1 = zero(Symm2) + backstress2 = zero(Symm2) + plastic_strain = zero(Symm2) cumeq = 0.0 R = 0.0 state = AnotherState(stress, strain, backstress1, backstress2, plastic_strain, cumeq, R) for i in 1:N - dstress = randn(SymmetricTensor{2,3}) - dstrain = randn(SymmetricTensor{2,3}) - dbackstress1 = randn(SymmetricTensor{2,3}) - dbackstress2 = randn(SymmetricTensor{2,3}) - dplastic_strain = randn(SymmetricTensor{2,3}) + dstress = randn(Symm2) + dstrain = randn(Symm2) + dbackstress1 = randn(Symm2) + dbackstress2 = randn(Symm2) + dplastic_strain = randn(Symm2) dcumeq = norm(dplastic_strain) dR = randn() dstate = AnotherState(dstress, dstrain, dbackstress1, diff --git a/examples/continuum.jl b/examples/continuum.jl index 499af8f..1e833ae 100644 --- a/examples/continuum.jl +++ b/examples/continuum.jl @@ -5,7 +5,7 @@ mutable struct Continuum3D <: FieldProblem material_model :: Symbol end -Continuum3D() = Continuum3D(:IdealPlastic) +Continuum3D() = Continuum3D(:PerfectPlastic) FEMBase.get_unknown_field_name(::Continuum3D) = "displacement" function FEMBase.assemble_elements!(problem::Problem{Continuum3D}, diff --git a/examples/one_element_ideal_plastic.jl b/examples/one_element_ideal_plastic.jl index 88020a4..15fd1cb 100644 --- a/examples/one_element_ideal_plastic.jl +++ b/examples/one_element_ideal_plastic.jl @@ -5,7 +5,7 @@ using Materials, FEMBase, LinearAlgebra # Standard simulation of ideal plastic material model -analysis, problem, element, bc_elements, ip = get_material_analysis(:IdealPlastic) +analysis, problem, element, bc_elements, ip = get_material_analysis(:PerfectPlastic) update!(element, "youngs modulus", 200.0e3) update!(element, "poissons ratio", 0.3) update!(element, "yield stress", 100.0) @@ -44,7 +44,7 @@ t = 0.0:0.1:3.0 plot(e11.(t), s11.(t), label="\$\\sigma_{11}\$") plot!(e22.(t), s22.(t), label="\$\\sigma_{22}\$") plot!(e33.(t), s33.(t), label="\$\\sigma_{33}\$") -title!("Stress-strain curve of idealplastic material model, uniaxial strain") +title!("Stress-strain curve of perfect plastic material model, uniaxial strain") ylabel!("Stress [MPa]") xlabel!("Strain [str]") savefig(joinpath("one_element_ideal_plastic/uniaxial_strain.svg")) diff --git a/src/DSA.jl b/src/DSA.jl index 3ed8745..3844d62 100644 --- a/src/DSA.jl +++ b/src/DSA.jl @@ -1,174 +1,306 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE -@with_kw mutable struct DSADriverState <: AbstractMaterialState - time::Float64 = zero(Float64) - strain::SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) +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 -@with_kw struct DSAParameterState <: 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 - Q::Float64 = 0.0 - b::Float64 = 0.0 - w::Float64 = 0.0 - P1::Float64 = 0.0 - P2::Float64 = 0.0 - m::Float64 = 0.0 - m1::Float64 = 0.0 - m2::Float64 = 0.0 - M1::Float64 = 0.0 - M2::Float64 = 0.0 - ba::Float64 = 0.0 - xi::Float64 = 0.0 +# 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 -@with_kw struct DSAVariableState <: 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) - jacobian::SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64}) - ta::Float64 = zero(Float64) - Ra::Float64 = zero(Float64) +# TODO: Does this eventually need a {T}? +@with_kw struct DSAOptions <: AbstractMaterialState + nlsolve_method::Symbol = :trust_region end -@with_kw mutable struct DSA <: AbstractMaterial - drivers::DSADriverState = DSADriverState() - ddrivers::DSADriverState = DSADriverState() - variables::DSAVariableState = DSAVariableState() - variables_new::DSAVariableState = DSAVariableState() - parameters::DSAParameterState = DSAParameterState() - dparameters::DSAParameterState = DSAParameterState() +@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 -function integrate_material!(material::DSA) +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 - mu = E / ( 2.0 * (1.0 + nu) ) - lambda = E * nu / ( (1.0 + nu) * (1.0 - 2.0 * nu) ) + 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) - seff = stress - X1 - X2 - seff_dev = dev(seff) - f = sqrt(1.5) * norm(seff_dev) - (R0 + R + (1 - xi) * Ra) + # 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 = [tovoigt(stress); R; tovoigt(X1); tovoigt(X2); ta + dtime; Ra] - F = similar(x0) - res = nlsolve(g!, x0; autodiff = :forward) + 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 - res.f_converged || error("Nonlinear system of equations did not converge!") - - 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]) - ta = x[20] - Ra = x[21] - - seff = stress - X1 - X2 - seff_dev = dev(seff) - 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 + 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 - function residuals(x) # Compute Jacobian - F = similar(x) - g!(F, x) - return F - end - drdx = ForwardDiff.jacobian(residuals, x) + # 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) - jacobian = fromvoigt(SymmetricTensor{4,3}, (drdx \ drde)[1:6, 1: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 = DSAVariableState(stress = stress, X1 = X1, X2 = X2, R = R, plastic_strain = plastic_strain, cumeq = cumeq, jacobian = jacobian, ta = ta, Ra = Ra) + 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 -function create_nonlinear_system_of_equations(material::DSA) +""" + 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 - mu = E / ( 2.0 * (1.0 + nu) ) - lambda = E * nu / ( (1.0 + nu) * (1.0 - 2.0 * nu) ) + 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 - function g!(F, x::Vector{T}) where {T} # System of non-linear equations - 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]) - ta_ = x[20] - Ra_ = x[21] - - jacobian = isotropic_elasticity_tensor(lambda, mu) - seff = stress_ - X1_ - X2_ - seff_dev = dev(seff) - 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) - Ras = P1 * (1.0 - exp(-P2 * ta_^m)) - dstrain_plastic = dp * n - - tovoigt!(view(F, 1:6), stress - stress_ + dcontract(jacobian, dstrain - dstrain_plastic)) - F[7] = R - R_ + b * ( Q - R_ ) * dp - if isapprox(C1, 0.0) - tovoigt!(view(F, 8:13), X1 - X1_) + 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 + + # HACK: The zero special case is needed here to make ForwardDiff happy. + # + # Otherwise, when ndX1_new = 0, the components 2:end of the automatic + # derivative of JX1_new will be NaN, which causes the calculation of the + # material jacobian to silently fail. This usually manifests itself as a + # mysterious convergence failure, when this model is used in the strain + # optimizer. + ndX1_new = norm(dev(X1_new)) + if iszero(ndX1_new) + JX1_new = 0.0 else - ndX1_ = norm(dev(X1_)) # Static recovery component - if iszero(ndX1_) - JX1_ = 0.0 - else - JX1_ = sqrt(1.5) * ndX1_ - end - sr1_ = ( JX1_^(m1 - 1) * X1_ ) / (M1^m1) - tovoigt!(view(F, 8:13), X1 - X1_ + 2.0 / 3.0 * C1 * dp * ( n - 1.5 * D1 / C1 * X1_ ) - dtime * sr1_) + JX1_new = sqrt(1.5) * ndX1_new end - if isapprox(C2, 0.0) - tovoigt!(view(F, 14:19), X2 - X2_) + sr1_new = (JX1_new^(m1 - 1) * X1_new) / (M1^m1) # static recovery term + tovoigt!(view(F, 8:13), X1 - X1_new + dp*(2.0/3.0*C1*n - D1*X1_new) - dtime*sr1_new) + + ndX2_new = norm(dev(X2_new)) + if iszero(ndX2_new) + JX2_new = 0.0 else - ndX2_ = norm(dev(X2_)) # Static recovery component - if iszero(ndX2_) - JX2_ = 0.0 - else - JX2_ = sqrt(1.5) * ndX2_ - end - sr2_ = ( JX2_^(m2 - 1) * X2_ ) / (M2^m2) - tovoigt!(view(F, 14:19), X2 - X2_ + 2.0 / 3.0 * C2 * dp * ( n - 1.5 * D2 / C2 * X2_ ) - dtime * sr2_) + JX2_new = sqrt(1.5) * ndX2_new end - F[20] = ta - ta_ + dtime - (ta_ / w) * dp - F[21] = Ra - Ra_ + ba * (Ras - Ra_) * dp + sr2_new = (JX2_new^(m2 - 1) * X2_new) / (M2^m2) # static recovery term + tovoigt!(view(F, 14:19), X2 - X2_new + dp*(2.0/3.0*C2*n - D2*X2_new) - dtime*sr2_new) + + 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 diff --git a/src/Materials.jl b/src/Materials.jl index f864dca..d7150be 100644 --- a/src/Materials.jl +++ b/src/Materials.jl @@ -3,23 +3,53 @@ module Materials -using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters - abstract type AbstractMaterial end abstract type AbstractMaterialState end -@generated function Base.:+(state::T, dstate::T) where {T <: AbstractMaterialState} +export AbstractMaterial, AbstractMaterialState +export integrate_material!, update_material!, reset_material! + +""" + :+(state::T, dstate::T) where T <: AbstractMaterialState + +Addition for material states. + +Given two material states `state` and `dstate` of type `T`, add each field of +`dstate` into the corresponding field of `state`. Return the resulting material +state. +""" +@generated function Base.:+(state::T, dstate::T) where T <: AbstractMaterialState expr = [:(state.$p + dstate.$p) for p in fieldnames(T)] return :(T($(expr...))) end -export AbstractMaterial, AbstractMaterialState +""" + integrate_material!(material::AbstractMaterial) + +Integrate one timestep. The input `material.variables` represents the old +problem state. + +Abstract method. Must be implemented for each material type. When integration is +done, the method **must** write the new state into `material.variables_new`. -function integrate_material!(material::M) where {M<:AbstractMaterial} - error("One needs to define how to integrate material $(M)!") +**Do not** write into `material.variables`; actually committing the timestep +(i.e. accepting that one step of time evolution and applying it permanently) +is the job of `update_material!`. +""" +function integrate_material!(material::M) where M <: AbstractMaterial + error("One needs to define how to integrate material $M!") end -function update_material!(material::M) where {M <: AbstractMaterial} +""" + update_material!(material::AbstractMaterial) + +Commit the result of `integrate_material!`. + +In `material`, we add `ddrivers` into `drivers`, `dparameters` into +`parameters`, and replace `variables` by `variables_new`. Then we +automatically invoke `reset_material!`. +""" +function update_material!(material::AbstractMaterial) material.drivers += material.ddrivers material.parameters += material.dparameters material.variables = material.variables_new @@ -27,51 +57,46 @@ function update_material!(material::M) where {M <: AbstractMaterial} return nothing end -function reset_material!(material::M) where {M <: AbstractMaterial} +""" + reset_material!(material::AbstractMaterial) + +In `material`, we zero out `ddrivers`, `dparameters` and `variables_new`. This +clears out the tentative state produced when a timestep has been computed, but +has not yet been committed. + +Used internally by `update_material!`. +""" +function reset_material!(material::AbstractMaterial) material.ddrivers = typeof(material.ddrivers)() material.dparameters = typeof(material.dparameters)() material.variables_new = typeof(material.variables_new)() return nothing end -export integrate_material!, update_material!, reset_material! - -function isotropic_elasticity_tensor(lambda, mu) - delta(i,j) = i==j ? 1.0 : 0.0 - g(i,j,k,l) = lambda*delta(i,j)*delta(k,l) + mu*(delta(i,k)*delta(j,l)+delta(i,l)*delta(j,k)) - jacobian = SymmetricTensor{4, 3, Float64}(g) - 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("utilities.jl") +using .Utilities +export Symm2, Symm4 +export delta, II, IT, IS, IA, IV, ID, isotropic_elasticity_tensor, isotropic_compliance_tensor +export lame, delame, debang, find_root -include("idealplastic.jl") -export IdealPlastic, IdealPlasticDriverState, IdealPlasticParameterState, IdealPlasticVariableState +include("perfectplastic.jl") +using .PerfectPlasticModule +export PerfectPlastic, PerfectPlasticDriverState, PerfectPlasticParameterState, PerfectPlasticVariableState include("chaboche.jl") +using .ChabocheModule export Chaboche, ChabocheDriverState, ChabocheParameterState, ChabocheVariableState include("memory.jl") +using .MemoryModule export Memory, MemoryDriverState, MemoryParameterState, MemoryVariableState -# include("viscoplastic.jl") -# export ViscoPlastic - -include("uniaxial_increment.jl") -export uniaxial_increment! - -include("biaxial_increment.jl") -export biaxial_increment! - -include("stress_driven_uniaxial_increment.jl") -export stress_driven_uniaxial_increment! - include("DSA.jl") +using .DSAModule export DSA, DSADriverState, DSAParameterState, DSAVariableState +include("increments.jl") +using .Increments +export uniaxial_increment!, biaxial_increment!, stress_driven_uniaxial_increment! + end diff --git a/src/biaxial_increment.jl b/src/biaxial_increment.jl deleted file mode 100644 index 3539162..0000000 --- a/src/biaxial_increment.jl +++ /dev/null @@ -1,23 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - - function biaxial_increment!(material, dstrain11, dstrain12, - dt; - dstrain=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0, 0, dstrain12], # Voigt notation - max_iter=50, norm_acc=1e-9) - converged = false - stress0 = tovoigt(material.variables.stress) - for i=1:max_iter - material.ddrivers.time = dt - material.ddrivers.strain = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0) - integrate_material!(material) - stress = tovoigt(material.variables_new.stress) - dstress = stress - stress0 - D = tovoigt(material.variables_new.jacobian) - dstr = -D[2:end-1,2:end-1]\dstress[2:end-1] - dstrain[2:end-1] .+= dstr - norm(dstr) < norm_acc && (converged = true; break) - end - converged || error("No convergence in strain increment") - return nothing -end diff --git a/src/chaboche.jl b/src/chaboche.jl index 4f62114..8136550 100644 --- a/src/chaboche.jl +++ b/src/chaboche.jl @@ -1,147 +1,267 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE -@with_kw mutable struct ChabocheDriverState <: AbstractMaterialState - time :: Float64 = zero(Float64) - strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) +module ChabocheModule + +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 GenericChaboche, GenericChabocheDriverState, GenericChabocheParameterState, GenericChabocheVariableState + +# specialization for Float64 +export Chaboche, ChabocheDriverState, ChabocheParameterState, ChabocheVariableState + +@with_kw mutable struct GenericChabocheDriverState{T <: Real} <: AbstractMaterialState + time::T = zero(T) + strain::Symm2{T} = zero(Symm2{T}) +end + +# TODO: complete this docstring +"""Parameter state for Chaboche material. + +The classical viscoplastic material is a special case of this model with `C1 = C2 = 0`. + +`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 +""" +@with_kw struct GenericChabocheParameterState{T <: Real} <: AbstractMaterialState + E::T = 0 + nu::T = 0 + R0::T = 0 + Kn::T = 0 + nn::T = 0 + C1::T = 0 + D1::T = 0 + C2::T = 0 + D2::T = 0 + Q::T = 0 + b::T = 0 +end + +"""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 +""" +@with_kw struct GenericChabocheVariableState{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}) end -@with_kw struct ChabocheParameterState <: 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 - Q :: Float64 = 0.0 - b :: Float64 = 0.0 +# TODO: Does this eventually need a {T}? +@with_kw struct ChabocheOptions <: AbstractMaterialState + nlsolve_method::Symbol = :trust_region end -@with_kw struct ChabocheVariableState <: 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) - jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64}) +@with_kw mutable struct GenericChaboche{T <: Real} <: AbstractMaterial + drivers::GenericChabocheDriverState{T} = GenericChabocheDriverState{T}() + ddrivers::GenericChabocheDriverState{T} = GenericChabocheDriverState{T}() + variables::GenericChabocheVariableState{T} = GenericChabocheVariableState{T}() + variables_new::GenericChabocheVariableState{T} = GenericChabocheVariableState{T}() + parameters::GenericChabocheParameterState{T} = GenericChabocheParameterState{T}() + dparameters::GenericChabocheParameterState{T} = GenericChabocheParameterState{T}() + options::ChabocheOptions = ChabocheOptions() +end + +ChabocheDriverState = GenericChabocheDriverState{Float64} +ChabocheParameterState = GenericChabocheParameterState{Float64} +ChabocheVariableState = GenericChabocheVariableState{Float64} +Chaboche = GenericChaboche{Float64} + +""" + state_to_vector(sigma::U, R::T, X1::U, X2::U) 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) where U <: Symm2{T} where T <: Real + return [tovoigt(sigma); R; tovoigt(X1); tovoigt(X2)]::Vector{T} end -@with_kw mutable struct Chaboche <: AbstractMaterial - drivers :: ChabocheDriverState = ChabocheDriverState() - ddrivers :: ChabocheDriverState = ChabocheDriverState() - variables :: ChabocheVariableState = ChabocheVariableState() - variables_new :: ChabocheVariableState = ChabocheVariableState() - parameters :: ChabocheParameterState = ChabocheParameterState() - dparameters :: ChabocheParameterState = ChabocheParameterState() +""" + 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]) + return sigma, R, X1, X2 end -function integrate_material!(material::Chaboche) +""" + integrate_material!(material::GenericChaboche{T}) where T <: Real + +Chaboche material with two backstresses. Both kinematic and isotropic hardening. +""" +function integrate_material!(material::GenericChaboche{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 = p - mu = E/(2.0*(1.0+nu)) - lambda = E*nu/((1.0+nu)*(1.0-2.0*nu)) + lambda, mu = lame(E, nu) @unpack strain, time = d dstrain = dd.strain dtime = dd.time - @unpack stress, X1, X2, plastic_strain, cumeq, R, jacobian = v + @unpack stress, X1, X2, plastic_strain, cumeq, R = v - jacobian = isotropic_elasticity_tensor(lambda, mu) + # elastic part + jacobian = isotropic_elasticity_tensor(lambda, mu) # dσ/dε, i.e. ∂σij/∂εkl + stress += dcontract(jacobian, dstrain) # add the elastic stress increment, get the elastic trial stress - stress += dcontract(jacobian, dstrain) - seff = stress - X1 - X2 - seff_dev = dev(seff) - f = sqrt(1.5)*norm(seff_dev) - (R0 + R) + # 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) # using elastic trial problem state 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) + x0 = state_to_vector(stress, R, X1, X2) + res = nlsolve(g!, x0; method=material.options.nlsolve_method, autodiff=:forward) # user manual: https://github.com/JuliaNLSolvers/NLsolve.jl + converged(res) || error("Nonlinear system of equations did not converge!") x = res.zero - res.f_converged || error("Nonlinear system of equations did not converge!") - - 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) + stress, R, X1, X2 = state_from_vector(x) + + # using the new problem state + seff_dev = dev(stress - X1 - X2) 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) + + dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn # plasticity multiplier, see equations (3) and (4) in Chaboche 2013 + dp = dotp*dtime # |dε_p|, using backward Euler (dotp is ∂ε_p/∂t at the end of the timestep) + n = sqrt(1.5)*seff_dev/norm(seff_dev) # Chaboche: a (tensorial) unit direction, s.t. 2/3 * (n : n) = 1; also n = ∂f/∂σ. plastic_strain += dp*n - cumeq += dp - # 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]) + cumeq += dp # cumulative equivalent plastic strain (note dp ≥ 0) + + # Compute the new Jacobian, accounting for the plastic contribution. Because + # x ≡ [σ R X1 X2] (vector of length 19, with tensors encoded in Voigt format) + # we have + # dσ/dε = (dx/dε)[1:6,1:6] + # for which we can compute the LHS as follows: + # dx/dε = dx/dr dr/dε = inv(dr/dx) dr/dε ≡ (dr/dx) \ (dr/dε) + # where r = r(x) is the residual, given by the function g!. AD can get us dr/dx automatically, + # the other factor we will have to supply manually. + drdx = ForwardDiff.jacobian(debang(g!), x) # Array{19, 19} + drde = zeros((length(x),6)) # Array{19, 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]) end - variables_new = ChabocheVariableState(stress = stress, - X1 = X1, - X2 = X2, - R = R, - plastic_strain = plastic_strain, - cumeq = cumeq, - jacobian = jacobian) + variables_new = GenericChabocheVariableState{T}(stress = stress, + X1 = X1, + X2 = X2, + R = R, + plastic_strain = plastic_strain, + cumeq = cumeq, + jacobian = jacobian) material.variables_new = variables_new + return nothing end -function create_nonlinear_system_of_equations(material::Chaboche) +""" + create_nonlinear_system_of_equations(material::GenericChaboche{T}) where T <: Real + +Create and return an instance of the equation system for the incremental form of +the evolution equations of the Chaboche 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-19 vectors containing +[sigma, R, X1, X2], 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::GenericChaboche{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 = p - mu = E/(2.0*(1.0+nu)) - lambda = E*nu/((1.0+nu)*(1.0-2.0*nu)) + 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 = v - function g!(F, x::Vector{T}) where {T} # System of non-linear equations - 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]) + 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 = state_from_vector(x) # tentative new values from nlsolve - seff = stress_ - X1_ - X2_ - seff_dev = dev(seff) - f = sqrt(1.5)*norm(seff_dev) - (R0 + R_) + seff_dev = dev(stress_new - X1_new - X2_new) + f = sqrt(1.5)*norm(seff_dev) - (R0 + R_new) dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn dp = dotp*dtime n = sqrt(1.5)*seff_dev/norm(seff_dev) + + # The equations are written in an incremental form: + # + # Δσ = (∂σ/∂ε)_e : dε_e = (∂σ/∂ε)_e : (dε - dε_p) (components 1:6) + # ΔR = b (Q - R_new) |dε_p| (component 7) + # ΔX1 = (2/3) C1 |dε_p| (n - (3/2) (D1/C1) X1_new) (components 8:13) + # ΔX2 = (2/3) C2 |dε_p| (n - (3/2) (D2/C2) X2_new) (components 14:19) + # + # where + # + # Δ(...) = (...)_new - (...)_old + # + # Then move the terms on the RHS to the LHS to get the standard form, (stuff) = 0. + # Also, below we avoid the multiplication and division that cancel each other + # in the last terms of the equations for ΔX1 and ΔX2. + # dstrain_plastic = dp*n - tovoigt!(view(F, 1:6), stress - stress_ + dcontract(jacobian, dstrain - dstrain_plastic)) - F[7] = R - R_ + b*(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 + dstrain_elastic = dstrain - dstrain_plastic + tovoigt!(view(F, 1:6), stress_new - stress - dcontract(jacobian, dstrain_elastic)) + + F[7] = R_new - R - b*(Q - R_new)*dp + + tovoigt!(view(F, 8:13), X1_new - X1 - dp*(2.0/3.0*C1*n - D1*X1_new)) + tovoigt!(view(F, 14:19), X2_new - X2 - dp*(2.0/3.0*C2*n - D2*X2_new)) + return nothing end return g! end + +end diff --git a/src/idealplastic.jl b/src/idealplastic.jl deleted file mode 100644 index 50db7e5..0000000 --- a/src/idealplastic.jl +++ /dev/null @@ -1,73 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -@with_kw mutable struct IdealPlasticDriverState <: AbstractMaterialState - time :: Float64 = zero(Float64) - strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) -end - -@with_kw struct IdealPlasticParameterState <: AbstractMaterialState - youngs_modulus :: Float64 = zero(Float64) - poissons_ratio :: Float64 = zero(Float64) - yield_stress :: Float64 = zero(Float64) -end - -@with_kw struct IdealPlasticVariableState <: AbstractMaterialState - stress :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) - plastic_strain :: SymmetricTensor{2,3} = zero(SymmetricTensor{2,3,Float64}) - cumeq :: Float64 = zero(Float64) - jacobian :: SymmetricTensor{4,3} = zero(SymmetricTensor{4,3,Float64}) -end - -@with_kw mutable struct IdealPlastic <: AbstractMaterial - drivers :: IdealPlasticDriverState = IdealPlasticDriverState() - ddrivers :: IdealPlasticDriverState = IdealPlasticDriverState() - variables :: IdealPlasticVariableState = IdealPlasticVariableState() - variables_new :: IdealPlasticVariableState = IdealPlasticVariableState() - parameters :: IdealPlasticParameterState = IdealPlasticParameterState() - dparameters :: IdealPlasticParameterState = IdealPlasticParameterState() -end - -function integrate_material!(material::IdealPlastic) - p = material.parameters - v = material.variables - dd = material.ddrivers - d = material.drivers - - E = p.youngs_modulus - nu = p.poissons_ratio - mu = E/(2.0*(1.0+nu)) - lambda = E*nu/((1.0+nu)*(1.0-2.0*nu)) - R0 = p.yield_stress - - # @unpack strain, time = d - dstrain = dd.strain - dtime = dd.time - @unpack stress, plastic_strain, cumeq, jacobian = v - - jacobian = isotropic_elasticity_tensor(lambda, mu) - stress += dcontract(jacobian, dstrain) - seff_dev = dev(stress) - stress_v = sqrt(1.5)*norm(seff_dev) - f = stress_v - R0 - - if f>0.0 - n = 1.5*seff_dev/stress_v - dp = 1.0/(3.0*mu)*(stress_v - R0) - plastic_strain += dp*n - cumeq += dp - - stress -= dcontract(jacobian, dp*n) - delta(i,j) = i==j ? 1.0 : 0.0 - II = SymmetricTensor{4,3}((i,j,k,l) -> 0.5*(delta(i,k)*delta(j,l)+delta(i,l)*delta(j,k))) - P = II - 1.0/3.0*SymmetricTensor{4,3}((i,j,k,l) -> delta(i,j)*delta(k,l)) - EE = II + dp*dcontract(jacobian, 1.5*P/R0 - otimes(n,n)/R0) - ED = dcontract(inv(EE),jacobian) - jacobian = ED - otimes(dcontract(ED, n), dcontract(n, ED))/dcontract(dcontract(n, ED), n) - end - variables_new = IdealPlasticVariableState(stress=stress, - plastic_strain=plastic_strain, - cumeq=cumeq, - jacobian=jacobian) - material.variables_new = variables_new -end diff --git a/src/increments.jl b/src/increments.jl new file mode 100644 index 0000000..487ff69 --- /dev/null +++ b/src/increments.jl @@ -0,0 +1,165 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + +module Increments + +import LinearAlgebra: norm +import Tensors: tovoigt, fromvoigt + +import ..AbstractMaterial, ..integrate_material! +import ..Utilities: Symm2 + +export find_dstrain!, uniaxial_increment!, biaxial_increment!, stress_driven_uniaxial_increment! + +# The skeleton of the optimizer is always the same, so we provide it as a +# higher-order function. The individual specific optimizer functions +# (`update_dstrain!)` only need to define the "meat" of how to update `dstrain`. +""" + find_dstrain!(material::AbstractMaterial, dstrain::AbstractVector{<:Real}, + dt::Real, update_dstrain!::Function; + max_iter::Integer=50, tol::Real=1e-9) + +Find a compatible strain increment for `material`. + +The functions in this module are made to be able to easily simulate stress +states produced by some of the most common test machines. Take for example the +function `uniaxial_increment!`. In a push-pull machine (with a smooth specimen), +we know that the stress state is uniaxial (in the measuring volume), so given +the strain increment in the direction where the stress is nonzero, we find +a strain increment that produces zero stress in the other directions. + +The `dstrain` supplied to this routine is the initial guess for the +optimization. At each iteration, it must be updated by the user-defined +corrector `update_dstrain!`, whose call signature is expected to be: + + update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T}) + where V <: AbstractVector{T} where T <: Real + -> err::Real + +`dstrain` is the current value of the strain increment, in Voigt format. +Conversion to tensor format uses `offdiagscale=2.0`. The function must update +the Voigt format `dstrain` in-place. + +`dstress = stress - stress0`, where `stress` is the stress state predicted by +integrating the material for one timestep of length `dt`, using the current +value of `dstrain` as a driving strain increment, and `stress0` is the stress +state stored in `materials.variables.stress`. + +`jacobian` is ∂σij/∂εkl (`material.variables_new.jacobian`), as computed by the +material implementation. In many cases, the dstrain optimization can actually be +performed by a Newton-Raphson root finder, so we pass the jacobian to facilitate +writing the update formula for such a root finder. + +The return value `err` must be an error measure (Real, >= 0). + +The update is iterated at most `max_iter` times, until `err` falls below `tol`. + +If `max_iter` is reached and the error measure is still `tol` or greater, +`ErrorException` is thrown. + +Note the timestep is **not** committed; we call `integrate_material!`, but not +`update_material!`. Only `material.variables_new` is updated. +""" +function find_dstrain!(material::AbstractMaterial, dstrain::AbstractVector{<:Real}, + dt::Real, update_dstrain!::Function; + max_iter::Integer=50, tol::Real=1e-9) + stress0 = tovoigt(material.variables.stress) # observed + T = typeof(dstrain[1]) + # @debug "---START---" + for i=1:max_iter + # @debug "$i, $dstrain, $stress0, $(material.variables.stress)" + material.ddrivers.time = dt + material.ddrivers.strain = fromvoigt(Symm2{T}, dstrain; offdiagscale=2.0) + integrate_material!(material) + stress = tovoigt(material.variables_new.stress) # predicted + dstress = stress - stress0 + jacobian = tovoigt(material.variables_new.jacobian) + e = update_dstrain!(dstrain, dstress, jacobian) + if e < tol + return nothing + end + end + error("No convergence in strain increment") +end + +""" + uniaxial_increment!(material::AbstractMaterial, dstrain11::Real, dt::Real; + dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0], + max_iter::Integer=50, norm_acc::Real=1e-9) + +Find a compatible strain increment for `material`. + +The material state (`material.variables`) and the component 11 of the *strain* +increment are taken as prescribed. This routine computes the other components of +the strain increment such that the predicted stress state matches the stored +one. + +See `find_dstrain!`. +""" +function uniaxial_increment!(material::AbstractMaterial, dstrain11::Real, dt::Real; + dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0], + max_iter::Integer=50, norm_acc::Real=1e-9) + function update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T}) where V <: AbstractVector{T} where T <: Real + dstrain_correction = -jacobian[2:end,2:end] \ dstress[2:end] + dstrain[2:end] .+= dstrain_correction + return norm(dstrain_correction) + end + find_dstrain!(material, dstrain, dt, update_dstrain!, max_iter=max_iter, tol=norm_acc) + return nothing +end + +""" + biaxial_increment!(material::AbstractMaterial, dstrain11::Real, dstrain12::Real, dt::Real; + dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0, 0, dstrain12], + max_iter::Integer=50, norm_acc::Real=1e-9) + +Find a compatible strain increment for `material`. + +The material state (`material.variables`) and the components 11 and 12 of the +*strain* increment are taken as prescribed. This routine computes the other +components of the strain increment such that the predicted stress state matches +the stored one. + +See `find_dstrain!`. +""" +function biaxial_increment!(material::AbstractMaterial, dstrain11::Real, dstrain12::Real, dt::Real; + dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0, 0, dstrain12], + max_iter::Integer=50, norm_acc::Real=1e-9) + function update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T}) where V <: AbstractVector{T} where T <: Real + dstrain_correction = -jacobian[2:end-1,2:end-1] \ dstress[2:end-1] + dstrain[2:end-1] .+= dstrain_correction + return norm(dstrain_correction) + end + find_dstrain!(material, dstrain, dt, update_dstrain!, max_iter=max_iter, tol=norm_acc) + return nothing +end + +""" + stress_driven_uniaxial_increment!(material::AbstractMaterial, dstress11::Real, dt::Real; + dstrain::AbstractVector{<:Real}=[dstress11/200e3, -0.3*dstress11/200e3, -0.3*dstress11/200e3, 0.0, 0.0, 0.0], + max_iter::Integer=50, norm_acc::Real=1e-9) + +Find a compatible strain increment for `material`. + +The material state (`material.variables`) and the component 11 of the *stress* +increment are taken as prescribed. This routine computes a strain increment such +that the predicted stress state matches the stored one. + +See `find_dstrain!`. +""" +function stress_driven_uniaxial_increment!(material::AbstractMaterial, dstress11::Real, dt::Real; + dstrain::AbstractVector{<:Real}=[dstress11/200e3, -0.3*dstress11/200e3, -0.3*dstress11/200e3, 0.0, 0.0, 0.0], + max_iter::Integer=50, norm_acc::Real=1e-9) + function update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T}) where V <: AbstractVector{T} where T <: Real + # Mutation of `dstress` doesn't matter, since `dstress` is freshly generated at each iteration. + # The lexical closure property gives us access to `dstress11` in this scope. + dstress[1] -= dstress11 + dstrain_correction = -jacobian \ dstress + dstrain .+= dstrain_correction + return norm(dstrain_correction) + end + find_dstrain!(material, dstrain, dt, update_dstrain!, max_iter=max_iter, tol=norm_acc) + return nothing +end + +end diff --git a/src/memory.jl b/src/memory.jl index f9eec79..08586bc 100644 --- a/src/memory.jl +++ b/src/memory.jl @@ -1,94 +1,139 @@ # 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}) +# TODO: write docstrings for all public functions + +module MemoryModule + +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 GenericMemory, GenericMemoryDriverState, GenericMemoryParameterState, GenericMemoryVariableState + +# specialization for Float64 +export Memory, MemoryDriverState, MemoryParameterState, MemoryVariableState + +@with_kw mutable struct GenericMemoryDriverState{T <: Real} <: AbstractMaterialState + time::T = zero(T) + strain::Symm2{T} = zero(Symm2{T}) +end + +@with_kw struct GenericMemoryParameterState{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 + Q0::T = 0.0 + QM::T = 0.0 + mu::T = 0.0 + b::T = 0.0 + eta::T = 0.0 + m::T = 0.0 + pt::T = 0.0 + xi::T = 0.0 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 +@with_kw struct GenericMemoryVariableState{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) + q::T = zero(T) + zeta::Symm2{T} = zero(Symm2{T}) + jacobian::Symm4{T} = zero(Symm4{T}) 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}) +# TODO: Does this eventually need a {T}? +@with_kw struct MemoryOptions <: AbstractMaterialState + nlsolve_method::Symbol = :trust_region 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() +@with_kw mutable struct GenericMemory{T <: Real} <: AbstractMaterial + drivers::GenericMemoryDriverState{T} = GenericMemoryDriverState{T}() + ddrivers::GenericMemoryDriverState{T} = GenericMemoryDriverState{T}() + variables::GenericMemoryVariableState{T} = GenericMemoryVariableState{T}() + variables_new::GenericMemoryVariableState{T} = GenericMemoryVariableState{T}() + parameters::GenericMemoryParameterState{T} = GenericMemoryParameterState{T}() + dparameters::GenericMemoryParameterState{T} = GenericMemoryParameterState{T}() + options::MemoryOptions = MemoryOptions() end -function integrate_material!(material::Memory) +MemoryDriverState = GenericMemoryDriverState{Float64} +MemoryParameterState = GenericMemoryParameterState{Float64} +MemoryVariableState = GenericMemoryVariableState{Float64} +Memory = GenericMemory{Float64} + +""" + state_to_vector(sigma::U, R::T, X1::U, X2::U) 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) where U <: Symm2{T} where T <: Real + return [tovoigt(sigma); R; tovoigt(X1); tovoigt(X2)]::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]) + return sigma, R, X1, X2 +end + +function integrate_material!(material::GenericMemory{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, 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)) + lambda, elastic_mu = lame(E, 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_) + # elastic part + jacobian = isotropic_elasticity_tensor(lambda, elastic_mu) stress += dcontract(jacobian, dstrain) - seff = stress - X1 - X2 - seff_dev = dev(seff) - f = sqrt(1.5)*norm(seff_dev) - (R0 + R) + + # 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) # using elastic trial problem state if f > 0.0 + # Explicit update to memory-surface 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!") + x0 = state_to_vector(stress, R, X1, X2) + 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 = 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]) + stress, R, X1, X2 = state_from_vector(x) - seff = stress - X1 - X2 - seff_dev = dev(seff) + # using the new problem state + seff_dev = dev(stress - X1 - X2) 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 @@ -98,106 +143,110 @@ function integrate_material!(material::Memory) if FF > 0.0 nF = 1.5*dev(plastic_strain - zeta)/JF nnF = dcontract(n, nF) - if nnF>0 + 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 + 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) + # 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) - jacobian = fromvoigt(SymmetricTensor{4,3}, (drdx\drde)[1:6, 1:6]) + jacobian = fromvoigt(Symm4, (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) + variables_new = GenericMemoryVariableState{T}(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 + return nothing end -function create_nonlinear_system_of_equations(material::Memory) +function create_nonlinear_system_of_equations(material::GenericMemory{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, 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)) + lambda, elastic_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, 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]) + jacobian = isotropic_elasticity_tensor(lambda, elastic_mu) + + # Explicit update of memory surface. + # + # 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 = state_from_vector(x) # tentative new values from nlsolve - seff = stress_ - X1_ - X2_ - seff_dev = dev(seff) - f = sqrt(1.5)*norm(seff_dev) - (R0 + R_) + seff_dev = dev(stress_new - X1_new - X2_new) + f = sqrt(1.5)*norm(seff_dev) - (R0 + R_new) 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)) + plastic_strain_new = plastic_strain + dstrain_plastic + JF = sqrt(1.5)*norm(dev(plastic_strain_new - zeta)) FF = 2.0/3.0*JF - q if FF > 0.0 - nF = 1.5*dev(plastic_strain_ - zeta)/JF + nF = 1.5*dev(plastic_strain_new - 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 + if nnF > 0 + q_new = q + 2.0/3.0*eta*nnF*dp + zeta_new = zeta + 2.0/3.0*(1.0 - eta)*nnF*nF*dp else - q_ = q - zeta_ = zeta + q_new = q + zeta_new = zeta end else # Memory evanescence term - p_ = cumeq + dp - if p_>pt - q_ = q - xi*q^m*dp + p_new = cumeq + dp + if p_new > pt + q_new = q - xi*q^m*dp else - q_ = q + q_new = q end - zeta_ = zeta + zeta_new = 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 + # 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_elastic = dstrain - dstrain_plastic + tovoigt!(view(F, 1:6), stress - stress_new + dcontract(jacobian, dstrain_elastic)) + + F[7] = R - R_new + b*((QM + (Q0 - QM)*exp(-2.0*mu*q_new)) - R_new)*dp + + tovoigt!(view(F, 8:13), X1 - X1_new + dp*(2.0/3.0*C1*n - D1*X1_new)) + tovoigt!(view(F, 14:19), X2 - X2_new + dp*(2.0/3.0*C2*n - D2*X2_new)) + return nothing end return g! end + +end diff --git a/src/perfectplastic.jl b/src/perfectplastic.jl new file mode 100644 index 0000000..0ced933 --- /dev/null +++ b/src/perfectplastic.jl @@ -0,0 +1,103 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + +module PerfectPlasticModule + +using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters + +import ..AbstractMaterial, ..AbstractMaterialState +import ..Utilities: Symm2, Symm4, isotropic_elasticity_tensor, IS, ID, lame +import ..integrate_material! # for method extension + +# parametrically polymorphic for any type representing ℝ +export GenericPerfectPlastic, GenericPerfectPlasticDriverState, GenericPerfectPlasticParameterState, GenericPerfectPlasticVariableState + +# specialization for Float64 +export PerfectPlastic, PerfectPlasticDriverState, PerfectPlasticParameterState, PerfectPlasticVariableState + +@with_kw mutable struct GenericPerfectPlasticDriverState{T <: Real} <: AbstractMaterialState + time::T = zero(T) + strain::Symm2{T} = zero(Symm2{T}) +end + +@with_kw struct GenericPerfectPlasticParameterState{T <: Real} <: AbstractMaterialState + youngs_modulus::T = zero(T) + poissons_ratio::T = zero(T) + yield_stress::T = zero(T) +end + +@with_kw struct GenericPerfectPlasticVariableState{T <: Real} <: AbstractMaterialState + stress::Symm2{T} = zero(Symm2{T}) + plastic_strain::Symm2{T} = zero(Symm2{T}) + cumeq::T = zero(T) + jacobian::Symm4{T} = zero(Symm4{T}) +end + +@with_kw mutable struct GenericPerfectPlastic{T <: Real} <: AbstractMaterial + drivers::GenericPerfectPlasticDriverState{T} = GenericPerfectPlasticDriverState{T}() + ddrivers::GenericPerfectPlasticDriverState{T} = GenericPerfectPlasticDriverState{T}() + variables::GenericPerfectPlasticVariableState{T} = GenericPerfectPlasticVariableState{T}() + variables_new::GenericPerfectPlasticVariableState{T} = GenericPerfectPlasticVariableState{T}() + parameters::GenericPerfectPlasticParameterState{T} = GenericPerfectPlasticParameterState{T}() + dparameters::GenericPerfectPlasticParameterState{T} = GenericPerfectPlasticParameterState{T}() +end + +PerfectPlastic = GenericPerfectPlastic{Float64} +PerfectPlasticDriverState = GenericPerfectPlasticDriverState{Float64} +PerfectPlasticParameterState = GenericPerfectPlasticParameterState{Float64} +PerfectPlasticVariableState = GenericPerfectPlasticVariableState{Float64} + +""" + integrate_material!(material::GenericPerfectPlastic) + +Ideal plastic material: no hardening. The elastic region remains centered on the +origin, and retains its original size. +""" +function integrate_material!(material::GenericPerfectPlastic{T}) where T <: Real + p = material.parameters + v = material.variables + dd = material.ddrivers + d = material.drivers + + E = p.youngs_modulus + nu = p.poissons_ratio + lambda, mu = lame(E, nu) + R0 = p.yield_stress + + # @unpack strain, time = d # not needed for this material + dstrain = dd.strain + dtime = dd.time + @unpack stress, plastic_strain, cumeq, jacobian = v + + jacobian = isotropic_elasticity_tensor(lambda, mu) # dσ/dε, i.e. ∂σij/∂εkl + stress += dcontract(jacobian, dstrain) # add the elastic stress increment, get the elastic trial stress + seff_dev = dev(stress) + f = sqrt(1.5)*norm(seff_dev) - R0 # von Mises yield function; f := J(seff_dev) - Y + + if f > 0.0 + dp = 1.0/(3.0*mu) * f + n = sqrt(1.5)*seff_dev/norm(seff_dev) # a (tensorial) unit direction, s.t. 2/3 * (n : n) = 1 + + plastic_strain += dp*n + cumeq += dp # cumulative equivalent plastic strain (note dp ≥ 0) + + # Ideal plastic material: the stress state cannot be outside the yield surface. + # Project it back to the yield surface. + stress -= dcontract(jacobian, dp*n) + + # Compute ∂σij/∂εkl, accounting for the plastic contribution. + # EE = IS + dp/R0 * (∂σ/∂ε)_e : ((3/2) ID - n ⊗ n) + EE = IS(T) + dp/R0 * dcontract(jacobian, 1.5*ID(T) - otimes(n,n)) # using the elastic jacobian + ED = dcontract(inv(EE), jacobian) + # J = ED - (ED : n) ⊗ (n : ED) / (n : ED : n) + jacobian = ED - otimes(dcontract(ED, n), dcontract(n, ED)) / dcontract(dcontract(n, ED), n) + end + variables_new = GenericPerfectPlasticVariableState(stress=stress, + plastic_strain=plastic_strain, + cumeq=cumeq, + jacobian=jacobian) + material.variables_new = variables_new + return nothing +end + +end diff --git a/src/stress_driven_uniaxial_increment.jl b/src/stress_driven_uniaxial_increment.jl deleted file mode 100644 index 8efa960..0000000 --- a/src/stress_driven_uniaxial_increment.jl +++ /dev/null @@ -1,23 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -function stress_driven_uniaxial_increment!(material, dstress11, - dt; dstrain=[dstress11/200e3, -0.3*dstress11/200e3, -0.3*dstress11/200e3, 0.0, 0.0, 0.0], max_iter=50, norm_acc=1e-9) - converged = false - stress0 = tovoigt(material.variables.stress) - for i=1:max_iter - material.ddrivers.time = dt - material.ddrivers.strain = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0) - integrate_material!(material) - stress = tovoigt(material.variables_new.stress) - dstress = stress - stress0 - r = dstress - r[1] -= dstress11 - D = tovoigt(material.variables_new.jacobian) - dstr = -D\r - dstrain .+= dstr - norm(dstr) < norm_acc && (converged = true; break) - end - converged || error("No convergence in strain increment") - return nothing -end diff --git a/src/uniaxial_increment.jl b/src/uniaxial_increment.jl deleted file mode 100644 index 4479385..0000000 --- a/src/uniaxial_increment.jl +++ /dev/null @@ -1,22 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - - function uniaxial_increment!(material, dstrain11, - dt; dstrain=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0], - max_iter=50, norm_acc=1e-9) - converged = false - stress0 = tovoigt(material.variables.stress) - for i=1:max_iter - material.ddrivers.time = dt - material.ddrivers.strain = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0) - integrate_material!(material) - stress = tovoigt(material.variables_new.stress) - dstress = stress - stress0 - D = tovoigt(material.variables_new.jacobian) - dstr = -D[2:end,2:end]\dstress[2:end] - dstrain[2:end] .+= dstr - norm(dstr) < norm_acc && (converged = true; break) - end - converged || error("No convergence in strain increment") - return nothing -end diff --git a/src/utilities.jl b/src/utilities.jl new file mode 100644 index 0000000..493dc08 --- /dev/null +++ b/src/utilities.jl @@ -0,0 +1,201 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + +module Utilities + +using Tensors, ForwardDiff + +export Symm2, Symm4 +export delta, II, IT, IS, IA, IV, ID, isotropic_elasticity_tensor, isotropic_compliance_tensor +export lame, delame, debang, find_root + +"""Symm2{T} is an alias for SymmetricTensor{2,3,T}.""" +const Symm2{T} = SymmetricTensor{2,3,T} + +"""Symm4{T} is an alias for SymmetricTensor{4,3,T}.""" +const Symm4{T} = SymmetricTensor{4,3,T} + +""" + delta(i::Integer, j::Integer) + +Kronecker delta, defined by delta(i, j) = (i == j) ? 1 : 0. +""" +delta(i::T, j::T) where T <: Integer = (i == j) ? one(T) : zero(T) + +# TODO: We could probably remove the type argument, and just let the results be +# inferred as Symm4{Int64}, Symm4{Rational{Int64}} and similar. Simpler to use, +# and those behave correctly in calculations with types involving other reals +# such as Float64. +# Performance implications? Is the Julia compiler smart enough to optimize? +""" + II(T::Type=Float64) + +Rank-4 unit tensor, defined by II : A = A for any rank-2 tensor A. +""" +II(T::Type=Float64) = Symm4{T}((i,j,k,l) -> delta(i,k)*delta(j,l)) + +""" + IT(T::Type=Float64) + +Rank-4 unit tensor, defined by IT : A = transpose(A) for any rank-2 tensor A. +""" +IT(T::Type=Float64) = Symm4{T}((i,j,k,l) -> delta(i,l)*delta(j,k)) + +""" + IS(T::Type=Float64) + +Rank-4 unit tensor, symmetric. IS ≡ (1/2) (II + IT). +""" +IS(T::Type=Float64) = 1//2 * (II(T) + IT(T)) + +""" + IA(T::Type=Float64) + +Rank-4 unit tensor, screw-symmetric. IA ≡ (1/2) (II - IT). +""" +IA(T::Type=Float64) = 1//2 * (II(T) - IT(T)) + +""" + IV(T::Type=Float64) + +Rank-4 unit tensor, volumetric. IS ≡ (1/3) I ⊗ I, where I is the rank-2 unit tensor. +""" +IV(T::Type=Float64) = Symm4{T}((i,j,k,l) -> 1//3 * delta(i,j)*delta(k,l)) + +""" + ID(T::Type=Float64) + +Rank-4 unit tensor, deviatoric. ID ≡ IS - IV. +""" +ID(T::Type=Float64) = IS(T) - IV(T) + +""" + isotropic_elasticity_tensor(lambda::T, mu::T) where T <: Real + +Compute the elasticity tensor C(i,j,k,l) (rank 4, symmetric) for an isotropic +material having the Lamé parameters `lambda` and `mu`. + +If you have (E, nu) instead, use `lame` to get (lambda, mu). +""" +isotropic_elasticity_tensor(lambda::T, mu::T) where T <: Real = 3 * lambda * IV(T) + 2 * mu * IS(T) + +# TODO: check: original expr from upstream/master: +# 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)) +""" + isotropic_compliance_tensor(lambda::T, mu::T) where T <: Real + +Compute the compliance tensor S(i,j,k,l) (rank 4, symmetric) for an isotropic +material having the Lamé parameters `lambda` and `mu`. + +If you have (E, nu) instead, use `lame` to get (lambda, mu). +""" +isotropic_compliance_tensor(lambda::T, mu::T) where T <: Real = -3 * lambda / (2*mu * (3*lambda + 2*mu)) * IV(T) + 1 / (2*mu) * IS(T) + +""" + lame(E::Real, nu::Real) + +Convert elastic parameters (E, nu) of an isotropic material to Lamé parameters (lambda, mu). + +See: + https://en.wikipedia.org/wiki/Template:Elastic_moduli +""" +function lame(E::Real, nu::Real) + lambda = E * nu / ((1 + nu) * (1 - 2 * nu)) + mu = E / (2 * (1 + nu)) + return lambda, mu +end + +""" + delame(lambda::Real, mu::Real) + +Convert Lamé parameters (lambda, mu) of an isotropic material to elastic parameters (E, nu). + +See: + https://en.wikipedia.org/wiki/Template:Elastic_moduli +""" +function delame(lambda::Real, mu::Real) + E = mu * (3 * lambda + 2 * mu) / (lambda + mu) + nu = lambda / (2 * (lambda + mu)) + return E, nu +end + +""" + debang(f!::Function, ex=nothing) + +Convert a mutating function into non-mutating form. + +`f!` must be a two-argument mutating function, which writes the result into its +first argument. The result of `debang` is then `f`, a single-argument +non-mutating function that allocates and returns the result. Schematically, +`f!(out, x)` becomes `f(x)`. + +When the type, size and shape of `out` is the same as those of `x`, it is enough +to supply just `f!`. When `f` is called, output will be allocated as `similar(x)`. + +When the type, size and/or shape of `out` are different from those of `x`, then +an example instance of the correct type with the correct size and shape for the +output must be supplied, as debang's `ex` argument. When `f` is called, output +will be allocated as `similar(ex)`. The `ex` instance will be automatically kept +alive by the lexical closure of `f`. + +# Note + +While the type of `out` is known at compile time, the size and shape are +typically runtime properties, not encoded into the type. For example, arrays +have the number of dimensions as part of the type, but the length of each +dimension is only defined at run time, when an instance is created. This is why +the `ex` argument is needed. + +# Etymology + +By convention, mutating functions are marked with an exclamation mark, a.k.a. +bang. This function takes away the bang. +""" +function debang(f!::Function, ex=nothing) + if ex === nothing + function f(x) + out = similar(x) + f!(out, x) + return out + end + return f + else + # We use a different name to make incremental compilation happy. + function f_with_ex(x) + out = similar(ex) + f!(out, x) + return out + end + return f_with_ex + end +end + +# This comes from the old viscoplastic.jl, and is currently unused. +# The original wording of the error message suggested this was planned to be used for "radial return". +"""A simple Newton solver for the vector x* such that f(x*) = 0. + +The input `x` is the initial guess. + +The default `dfdx=nothing` uses `ForwardDiff.jacobian` to compute the jacobian +automatically. In this case the output of `f` must be an `AbstractArray`. + +`tol` is measured in the vector norm of the change in `x` between successive +iterations. +""" +function find_root(f::Function, x::AbstractVector{<:Real}, + dfdx::Union{Function, Nothing}=nothing; + max_iter::Integer=50, tol::Real=1e-9) + if dfdx === nothing + dfdx = (x) -> ForwardDiff.jacobian(f, x) + end + for i=1:max_iter + dx = -dfdx(x) \ f(x) + x += dx + if norm(dx) < tol + return x + end + end + error("No convergence!") +end + +end diff --git a/src/viscoplastic.jl b/src/viscoplastic.jl deleted file mode 100644 index 125f0a9..0000000 --- a/src/viscoplastic.jl +++ /dev/null @@ -1,190 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE -using LinearAlgebra -# using ForwardDiff -using NLsolve - -# http://www.solid.iei.liu.se/Education/TMHL55/TMHL55_lp1_2010/lecture_notes/plasticity_flow_rule_isotropic_hardening.pdf -# http://mms.ensmp.fr/msi_paris/transparents/Georges_Cailletaud/2013-GC-plas3D.pdf - -##################################### -# Viscoplastic material definitions # -##################################### -mutable struct ViscoPlastic <: AbstractMaterial - # Material parameters - youngs_modulus :: Float64 - poissons_ratio :: Float64 - yield_stress :: Float64 - potential :: Symbol - params :: Array{Float64, 1} - - # Internal state variables - plastic_strain :: Vector{Float64} - dplastic_strain :: Vector{Float64} - plastic_multiplier :: Float64 - dplastic_multiplier :: Float64 - -end - -POTENTIAL_FUNCTIONS = [:norton] - -function ViscoPlastic(potential::Symbol, params :: Array{Float64,1}) - if !any(POTENTIAL_FUNCTIONS .== potential) - error("Potential $potential not found!") - end - youngs_modulus = 0.0 - poissons_ratio = 0.0 - yield_stress = 0.0 - plastic_strain = zeros(6) - dplastic_strain = zeros(6) - plastic_multiplier = 0.0 - dplastic_multiplier = 0.0 - dt = 0.0 - return ViscoPlastic(youngs_modulus, poissons_ratio, yield_stress, potential, params, - plastic_strain, dplastic_strain, plastic_multiplier, - dplastic_multiplier) -end - -function ViscoPlastic() - youngs_modulus = 0.0 - poissons_ratio = 0.0 - yield_stress = 0.0 - plastic_strain = zeros(6) - dplastic_strain = zeros(6) - plastic_multiplier = 0.0 - dplastic_multiplier = 0.0 - dt = 0.0 - potential = :norton - params = [180.0e3, 0.92] - return ViscoPlastic(youngs_modulus, poissons_ratio, yield_stress, potential, params, - plastic_strain, dplastic_strain, plastic_multiplier, - dplastic_multiplier) -end - -""" Double contraction -""" -function double_contraction(x,y) - return sum(x.*y) -end - -""" Deviatoric stress tensor -""" -function deviatoric_stress(stress) - stress_tensor = [ stress[1] stress[4] stress[6]; - stress[4] stress[2] stress[5]; - stress[6] stress[5] stress[3]] - stress_dev = stress_tensor - 1/3 * tr(stress_tensor) * Matrix(1.0I, 3, 3) - return stress_dev -end - -function J_2_stress(stress) - s = deviatoric_stress(stress) - return 1/2 * double_contraction(s,s) -end - -""" Equivalent stress -""" -function equivalent_stress(stress) - J_2 = J_2_stress(stress) - return sqrt(3 * J_2) -end - -""" Von mises yield function -""" -function von_mises_yield(stress, stress_y) - return equivalent_stress(stress) - stress_y -end - -""" Norton rule -""" -function norton_plastic_potential(stress, params, f) - K = params[1] - n = params[2] - J = f(stress) - return K/(n+1) * (J / K) ^ (n + 1) -end - -function bingham_plastic_potential(stress, params, f) - eta = params[1] - return 0.5 * (f(stress) / eta) ^2 -end - -""" Analytically derivated norton rule -""" -function dNortondStress(stress, params, f) - K = params[1] - n = params[2] - f_ = f(stress) - stress_v = equivalent_stress(stress) - stress_dev_vec = deviatoric_stress(stress) - sol_mat = (f_ / K) ^ n * 3 / 2 * stress_dev_vec / stress_v - return [sol_mat[1,1], sol_mat[2,2], sol_mat[3,3], - 2*sol_mat[1,2], 2*sol_mat[2,3], 2*sol_mat[1,3]] -end - -function find_root(f, df, x; max_iter=50, norm_acc=1e-9) - converged = false - for i=1:max_iter - dx = -df(x) \ f(x) - x += dx - norm(dx) < norm_acc && (converged = true; break) - end - converged || error("No convergence in radial return!") - return x -end - -function integrate_material!(material::Material{ViscoPlastic}) - mat = material.properties - - E = mat.youngs_modulus - nu = mat.poissons_ratio - potential = mat.potential - params = mat.params - mu = E/(2.0*(1.0+nu)) - lambda = E*nu/((1.0+nu)*(1.0-2.0*nu)) - - stress = material.stress - strain = material.strain - dstress = material.dstress - dstrain = material.dstrain - dt = material.dtime - D = material.jacobian - - dedt = dstrain ./ dt - - fill!(D, 0.0) - D[1,1] = D[2,2] = D[3,3] = 2.0*mu + lambda - D[4,4] = D[5,5] = D[6,6] = mu - D[1,2] = D[2,1] = D[2,3] = D[3,2] = D[1,3] = D[3,1] = lambda - - dstress[:] .= D * dedt .* dt - stress_tr = stress + dstress - f = x -> von_mises_yield(x, mat.yield_stress) - yield = f(stress_tr) - - if yield <= 0 - fill!(mat.dplastic_strain, 0.0) - return nothing - else - if potential == :norton - x0 = vec([dstress; 0.0]) - - function g!(F, x) - dsigma = x[1:6] - F[1:6] = dsigma - D * (dedt - dNortondStress(stress + dsigma, params, f)) .* dt - F[end] = von_mises_yield(stress + dsigma, mat.yield_stress) - end - - res = nlsolve(g!, x0) - dsigma = res.zero[1:6] - dstrain_pl = dNortondStress(stress + dsigma, params, f) - - end - mat.dplastic_strain = dstrain_pl - dstress[:] .= D*(dedt - dstrain_pl) .* dt - return nothing - end - - return nothing - -end diff --git a/test/run_lint.jl b/test/run_lint.jl deleted file mode 100644 index 6e7f672..0000000 --- a/test/run_lint.jl +++ /dev/null @@ -1,17 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -using Lint -using Base.Test - -results = lintpkg("Materials") -if !isempty(results) - info("Lint.jl is a tool that uses static analysis to assist in the development process by detecting common bugs and potential issues.") - info("For this package, Lint.jl report is following:") - display(results) - info("For more information, see https://lintjl.readthedocs.io/en/stable/") - warn("Package syntax test has failed.") - @test isempty(results) -else - info("Lint.jl: syntax check pass.") -end diff --git a/test/runtests.jl b/test/runtests.jl index 731a9e6..dd717a2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,37 +4,34 @@ using Materials, Test @testset "Test Materials.jl" begin - # @testset "test abstract material" begin - # include("test_abstract_material.jl") - # end - @testset "test ideal plastic material model" begin - include("test_idealplastic.jl") + @testset "test utilities" begin + include("test_utilities.jl") end - @testset "test ideal plastic pure shear" begin - include("test_idealplastic_shear.jl") + @testset "test perfect plastic uniaxial stress" begin + include("test_perfectplastic.jl") end - @testset "test uniaxial increment" begin - include("test_uniaxial_increment.jl") + @testset "test perfect plastic pure shear" begin + include("test_perfectplastic_shear.jl") + end + @testset "test chaboche uniaxial stress" begin + include("test_chaboche.jl") end @testset "test chaboche pure shear" begin include("test_chaboche_shear.jl") end - @testset "test chaboche uniaxial stress" begin - include("test_chaboche.jl") + @testset "test memory material model" begin + include("test_memory.jl") + end + @testset "test DSA material model" begin + include("test_DSA.jl") + end + @testset "test uniaxial increment" begin + include("test_uniaxial_increment.jl") end @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 - @testset "test DSA material model" begin - include("test_DSA.jl") - end - # @testset "test viscoplastic material model" begin - # include("test_viscoplastic.jl") - # end end diff --git a/test/test_abstract_material.jl b/test/test_abstract_material.jl deleted file mode 100644 index ba1af89..0000000 --- a/test/test_abstract_material.jl +++ /dev/null @@ -1,15 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -using Test - -struct MyMaterial <: AbstractMaterial - E :: Float64 -end - -function MyMaterial() - error("MyMaterial needs at least `E` defined.") -end - -mat = Material(MyMaterial, (E=210.0,)) -@test mat.properties.E == 210.0 diff --git a/test/test_biaxial_increment.jl b/test/test_biaxial_increment.jl index 7e9ec89..3d1f10b 100644 --- a/test/test_biaxial_increment.jl +++ b/test/test_biaxial_increment.jl @@ -1,39 +1,39 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + using Test, Tensors -dtime = 0.25 -parameters = ChabocheParameterState(E = 200.0e3, - nu = 0.3, - R0 = 100.0, - Kn = 100.0, - nn = 10.0, - C1 = 10000.0, - D1 = 100.0, - C2 = 50000.0, - D2 = 1000.0, - Q = 50.0, - b = 0.1) -mat = Chaboche(parameters = parameters) -times = [mat.drivers.time] -stresses = [copy(tovoigt(mat.variables.stress))] -dstrain11 = 1e-3*dtime -dstrain12 = 1e-3*dtime -dtimes = [dtime, dtime, dtime, dtime, 1.0] -dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11] -dstrains12 = [dstrain12, dstrain12, dstrain12, -dstrain12, -4*dstrain12] -plasticity_test = zeros(length(dstrains11)-1) -for i in 1:length(dtimes) - dstrain11 = dstrains11[i] - dstrain12 = dstrains12[i] - dtime = dtimes[i] - biaxial_increment!(mat, dstrain11, dstrain12, dtime) - update_material!(mat) - push!(stresses, copy(tovoigt(mat.variables.stress))) - if i > 1 - plasticity_test[i-1] = mat.variables.cumeq > 0.0 +let dtime = 0.25, + parameters = ChabocheParameterState(E=200.0e3, + nu=0.3, + R0=100.0, + Kn=100.0, + nn=10.0, + C1=10000.0, + D1=100.0, + C2=50000.0, + D2=1000.0, + Q=50.0, + b=0.1), + mat = Chaboche(parameters = parameters), + dstrain11 = 1e-3*dtime, + dstrain12 = 1e-3*dtime, + dtimes = [dtime, dtime, dtime, dtime, 1.0], + dstrains11 = dstrain11*[1.0, 1.0, 1.0, -1.0, -4.0], + dstrains12 = dstrain12*[1.0, 1.0, 1.0, -1.0, -4.0] + + plastic_flow_occurred = zeros(Bool, length(dstrains11) - 1) + for i in 1:length(dtimes) + dstrain11 = dstrains11[i] + dstrain12 = dstrains12[i] + dtime = dtimes[i] + biaxial_increment!(mat, dstrain11, dstrain12, dtime) + update_material!(mat) + if i > 1 + plastic_flow_occurred[i-1] = (mat.variables.cumeq > 0.0) + end + @test !iszero(mat.variables.stress[1,1]) && !iszero(mat.variables.stress[1,2]) + @test isapprox(tovoigt(mat.variables.stress)[2:5], zeros(4); atol=1e-8) end - @test !iszero(tovoigt(mat.variables.stress)[1]) && !iszero(tovoigt(mat.variables.stress)[end]) - @test isapprox(tovoigt(mat.variables.stress; offdiagscale=2.0)[2:5],zeros(4); atol=1e-5) + @test any(plastic_flow_occurred) end -@test any(x->x==1,plasticity_test) diff --git a/test/test_chaboche.jl b/test/test_chaboche.jl index 610a45a..5d84548 100644 --- a/test/test_chaboche.jl +++ b/test/test_chaboche.jl @@ -4,81 +4,43 @@ using Test, Tensors using DelimitedFiles -path = joinpath("test_chaboche", "unitelement_results.rpt") -data = readdlm(path, Float64; skipstart=4) -ts = data[:,1] -s11_ = data[:,2] -s12_ = data[:,3] -s13_ = data[:,4] -s22_ = data[:,5] -s23_ = data[:,6] -s33_ = data[:,7] -e11_ = data[:,8] -e12_ = data[:,9] -e13_ = data[:,10] -e22_ = data[:,11] -e23_ = data[:,12] -e33_ = data[:,13] +let path = joinpath("test_chaboche", "unitelement_results.rpt"), + data = readdlm(path, Float64; skipstart=4), + ts = data[:,1], + s11_ = data[:,2], + s12_ = data[:,3], + s13_ = data[:,4], + s22_ = data[:,5], + s23_ = data[:,6], + s33_ = data[:,7], + e11_ = data[:,8], + e12_ = data[:,9], + e13_ = data[:,10], + e22_ = data[:,11], + e23_ = data[:,12], + e33_ = data[:,13], + strains = [[e11_[i], e22_[i], e33_[i], e23_[i], e13_[i], e12_[i]] for i in 1:length(ts)], + parameters = ChabocheParameterState(E=200.0e3, + nu=0.3, + R0=100.0, + Kn=100.0, + nn=10.0, + C1=10000.0, + D1=100.0, + C2=50000.0, + D2=1000.0, + Q=50.0, + b=0.1), + mat = Chaboche(parameters=parameters) -strains = [[e11_[i], e22_[i], e33_[i], e23_[i], e13_[i], e12_[i]] for i in 1:length(ts)] - -parameters = ChabocheParameterState(E = 200.0e3, - nu = 0.3, - R0 = 100.0, - Kn = 100.0, - nn = 10.0, - C1 = 10000.0, - D1 = 100.0, - C2 = 50000.0, - D2 = 1000.0, - Q = 50.0, - b = 0.1) -chabmat = Chaboche(parameters = parameters) -s33s = [chabmat.variables.stress[3,3]] -for i=2:length(ts) - dtime = ts[i]-ts[i-1] - dstrain = fromvoigt(SymmetricTensor{2,3,Float64}, strains[i]-strains[i-1]; offdiagscale=2.0) - chabmat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain) - integrate_material!(chabmat) - update_material!(chabmat) - push!(s33s, chabmat.variables.stress[3,3]) + s33s = [mat.variables.stress[3,3]] + for i=2:length(ts) + dtime = ts[i] - ts[i-1] + dstrain = fromvoigt(Symm2{Float64}, strains[i] - strains[i-1]; offdiagscale=2.0) + mat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain) + integrate_material!(mat) + update_material!(mat) + push!(s33s, mat.variables.stress[3,3]) + end + @test isapprox(s33s, s33_; rtol=0.05) end -@test isapprox(s33s, s33_; rtol=0.05) - -# function test_chab() -# mat = Material(Chaboche, ()) -# mat.properties.youngs_modulus = 200.0e3 -# mat.properties.poissons_ratio = 0.3 -# mat.properties.yield_stress = 100.0 -# mat.properties.K_n = 100.0 -# mat.properties.n_n = 10.0 -# mat.properties.C_1 = 10000.0 -# mat.properties.D_1 = 100.0 -# mat.properties.C_2 = 50000.0 -# mat.properties.D_2 = 1000.0 -# mat.properties.Q = 50.0 -# mat.properties.b = 0.1 -# -# mat.stress = zeros(6) -# s33s = [0.0] -# for i=2:length(ts) -# dtime = ts[i]-ts[i-1] -# dstrain = strains[i]-strains[i-1] -# mat.dtime = dtime -# mat.dstrain = dstrain -# integrate_material!(mat) -# mat.time += mat.dtime -# mat.strain .+= mat.dstrain -# mat.stress .+= mat.dstress -# mat.properties.plastic_strain .+= mat.properties.dplastic_strain -# mat.properties.backstress1 .+= mat.properties.dbackstress1 -# mat.properties.backstress2 .+= mat.properties.dbackstress2 -# mat.properties.R += mat.properties.dR -# mat.properties.cumulative_equivalent_plastic_strain += mat.properties.dcumulative_equivalent_plastic_strain -# push!(s33s, copy(mat.stress[3])) -# end -# # @test isapprox(s33s, s33_; rtol=0.001) -# end - -# using BenchmarkTools -# @btime test_chab() diff --git a/test/test_chaboche_shear.jl b/test/test_chaboche_shear.jl index d60a854..785057d 100644 --- a/test/test_chaboche_shear.jl +++ b/test/test_chaboche_shear.jl @@ -3,85 +3,68 @@ using Test, Tensors -# mat = Material(Chaboche, tuple()) -# props = mat.properties -# props.youngs_modulus = 200.0e3 -# props.poissons_ratio = 0.3 -# props.yield_stress = 100.0 -# props.K_n = 100.0 -# props.n_n = 3.0 -# props.C_1 = 0.0 -# props.D_1 = 100.0 -# props.C_2 = 0.0 -# props.D_2 = 1000.0 -# props.Q = 0.0 -# props.b = 0.1 -E = 200.0e3 -nu = 0.3 -parameters = ChabocheParameterState(E = E, - nu = nu, - R0 = 100.0, - Kn = 100.0, - nn = 3.0, - C1 = 0.0, - D1 = 100.0, - C2 = 0.0, - D2 = 1000.0, - Q = 0.0, - b = 0.1) -mat = Chaboche(parameters = parameters) -times = [0.0] -loads = [0.0] -dt = 0.5 -G = 0.5*E/(1+nu) -syield = 100.0 -# vm = sqrt(3)*G*ga | ea = ga -ea = 2*syield/(sqrt(3)*G) -# Go to elastic border -push!(times, times[end]+dt) -push!(loads, loads[end] + ea*dt) - # Proceed to plastic flow -push!(times, times[end]+dt) -push!(loads, loads[end] + ea*dt) - # Reverse direction -push!(times, times[end]+dt) -push!(loads, loads[end] - ea*dt) - # Continue and pass yield criterion -push!(times, times[end]+dt) -push!(loads, loads[end] - ea*dt) -push!(times, times[end]+dt) -push!(loads, loads[end] - ea*dt) +let E = 200.0e3, + nu = 0.3, + parameters = ChabocheParameterState(E=E, + nu=nu, + R0=100.0, # yield in shear = R0 / sqrt(3) + Kn=100.0, + nn=3.0, + C1=0.0, + D1=100.0, + C2=0.0, + D2=1000.0, + Q=0.0, + b=0.1), + mat = Chaboche(parameters=parameters), + times = [0.0], + loads = [0.0], + dt = 1.0, + G = 0.5*E/(1+nu), + yield_strength = 100.0, + # vonMises = sqrt(3 J_2) = sqrt(3/2 tr(s^2)) = sqrt(3) |tau| = sqrt(3)*G*|gamma| + # gamma = 2 e12 + # set vonMises = Y + gamma_yield = yield_strength/(sqrt(3)*G) -eeqs = [mat.variables.cumeq] -stresses = [copy(tovoigt(mat.variables.stress))] -for i=2:length(times) - dtime = times[i]-times[i-1] - dstrain31 = loads[i]-loads[i-1] - dstrain = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain31] - # mat.dtime = dtime - # mat.dstrain = dstrain - dstrain_ = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0) - mat.ddrivers = ChabocheDriverState(time = dtime, strain = dstrain_) - integrate_material!(mat) - # mat.time += mat.dtime - # mat.strain .+= mat.dstrain - # mat.stress .+= mat.dstress - # mat.properties.plastic_strain .+= mat.properties.dplastic_strain - # mat.properties.backstress1 .+= mat.properties.dbackstress1 - # mat.properties.backstress2 .+= mat.properties.dbackstress2 - # mat.properties.R += mat.properties.dR - update_material!(mat) - push!(stresses, copy(tovoigt(mat.variables.stress))) - push!(eeqs, mat.variables.cumeq) - # @info "time = $(mat.time), stress = $(mat.stress), cumeq = $(mat.properties.cumulative_equivalent_plastic_strain))" -end + # Go to elastic border + push!(times, times[end] + dt) + push!(loads, loads[end] + gamma_yield*dt) + # Proceed to plastic flow + push!(times, times[end] + dt) + push!(loads, loads[end] + gamma_yield*dt) + # Reverse direction + push!(times, times[end] + dt) + push!(loads, loads[end] - gamma_yield*dt) + # Continue and pass yield criterion + push!(times, times[end] + dt) + push!(loads, loads[end] - gamma_yield*dt) + push!(times, times[end] + dt) + push!(loads, loads[end] - gamma_yield*dt) -for i in 1:length(times) - @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6) -end -s31 = [s[6] for s in stresses] + eeqs = [mat.variables.cumeq] + stresses = [copy(tovoigt(mat.variables.stress))] + for i=2:length(times) + dtime = times[i] - times[i-1] + dstrain12 = loads[i] - loads[i-1] + dstrain_voigt = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain12] + dstrain_tensor = fromvoigt(Symm2{Float64}, dstrain_voigt; offdiagscale=2.0) + mat.ddrivers = ChabocheDriverState(time=dtime, strain=dstrain_tensor) + integrate_material!(mat) + # @info "$i, $gamma_yield, $(mat.variables_new.stress[1,2]), $(2.0*mat.variables_new.plastic_strain[1,2])\n" + update_material!(mat) + push!(stresses, copy(tovoigt(mat.variables.stress))) + push!(eeqs, mat.variables.cumeq) + # @info "time = $(mat.time), stress = $(mat.stress), cumeq = $(mat.properties.cumulative_equivalent_plastic_strain))" + end -@test isapprox(s31[2], syield/sqrt(3.0)) -@test isapprox(s31[3]*sqrt(3.0), syield + 100.0*((eeqs[3]-eeqs[2])/dt)^(1.0/3.0); rtol=1e-2) -@test isapprox(s31[4], s31[3]-G*ea*dt) -@test isapprox(s31[6]*sqrt(3.0), -(syield + 100.0*((eeqs[6]-eeqs[5])/dt)^(1.0/3.0)); rtol=1e-2) + for i in 1:length(times) + @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6) + end + s31 = [s[6] for s in stresses] + + @test isapprox(s31[2], yield_strength/sqrt(3.0)) + @test isapprox(s31[3]*sqrt(3.0), yield_strength + 100.0*((eeqs[3] - eeqs[2])/dt)^(1.0/3.0); rtol=1e-2) + @test isapprox(s31[4], s31[3] - G*gamma_yield*dt) + @test isapprox(s31[6]*sqrt(3.0), -(yield_strength + 100.0*((eeqs[6] - eeqs[5])/dt)^(1.0/3.0)); rtol=1e-2) +end diff --git a/test/test_idealplastic.jl b/test/test_idealplastic.jl deleted file mode 100644 index 9fb32f5..0000000 --- a/test/test_idealplastic.jl +++ /dev/null @@ -1,76 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -using Test, Tensors -# time = [...] -# strain = [...] -# stress = run_simulation(mat, time, strain) -# @test isapprox(stress, stress_expected) - -# -# mat = Material(IdealPlastic, tuple()) -# mat.properties.youngs_modulus = 200.0e3 -# mat.properties.poissons_ratio = 0.3 -# mat.properties.yield_stress = 100.0 - -parameters = IdealPlasticParameterState(youngs_modulus = 200.0e3, - poissons_ratio = 0.3, - yield_stress = 100.0) - -dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64},1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) -ddrivers = IdealPlasticDriverState(time = 0.25, strain = 0.25*dstrain_dtime) -mat = IdealPlastic(parameters=parameters, ddrivers=ddrivers) -# mat.dtime = 0.25 -# mat.dstrain .= 1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0]*mat.dtime -integrate_material!(mat) -update_material!(mat) -@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [50.0, 0.0, 0.0, 0.0, 0.0, 0.0])) - -# mat.time += mat.dtime -# mat.strain .+= mat.dstrain -# mat.stress .+= mat.dstress -mat.ddrivers = ddrivers -integrate_material!(mat) -update_material!(mat) -@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [100.0, 0.0, 0.0, 0.0, 0.0, 0.0])) -@test isapprox(mat.variables.cumeq, 0.0; atol=1.0e-12) - -# mat.time += mat.dtime -# mat.strain .+= mat.dstrain -# mat.stress .+= mat.dstress -# mat.dstrain[:] .= 1.0e-3*[-0.5, -0.5, 1.0, 0.0, 0.0, 0.0]*mat.dtime -dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]; offdiagscale=2.0) -ddrivers = IdealPlasticDriverState(time = 0.25, strain = 0.25*dstrain_dtime) -mat.ddrivers = ddrivers -integrate_material!(mat) -update_material!(mat) -@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [100.0, 0.0, 0.0, 0.0, 0.0, 0.0]); atol=1.0e-12) -@test isapprox(mat.variables.cumeq, 0.25*1.0e-3) - -# mat.time += mat.dtime -# mat.strain .+= mat.dstrain -# mat.stress .+= mat.dstress -# mat.dstrain[:] .= -1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0]*mat.dtime -dstrain_dtime = fromvoigt(SymmetricTensor{2,3,Float64}, -1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) -ddrivers = IdealPlasticDriverState(time = 0.25, strain = 0.25*dstrain_dtime) -mat.ddrivers = ddrivers -integrate_material!(mat) -update_material!(mat) -@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [50.0, 0.0, 0.0, 0.0, 0.0, 0.0]); atol=1.0e-12) - -# mat.time += mat.dtime -# mat.strain .+= mat.dstrain -# mat.stress .+= mat.dstress -# mat.dtime = 1.0 -# m1 = 1.0e-3*[-0.3, -0.3, 1.0, 0.0, 0.0, 0.0] -# m2 = 1.0e-3*[-0.5, -0.5, 1.0, 0.0, 0.0, 0.0] -# mat.dstrain[:] .= -m1*0.75 - m2*0.25 -dstrain_dtime = (-0.75*fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]; offdiagscale=2.0) - -0.25*fromvoigt(SymmetricTensor{2,3,Float64}, 1e-3*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]; offdiagscale=2.0)) -ddrivers = IdealPlasticDriverState(time = 1.0, strain = dstrain_dtime) -mat.ddrivers = ddrivers -integrate_material!(mat) -integrate_material!(mat) -update_material!(mat) -@test isapprox(mat.variables.stress, fromvoigt(SymmetricTensor{2,3}, [-100.0, 0.0, 0.0, 0.0, 0.0, 0.0])) -# @test isapprox(mat.properties.dplastic_multiplier, 0.25*1.0e-3) diff --git a/test/test_idealplastic_shear.jl b/test/test_idealplastic_shear.jl deleted file mode 100644 index 2702b43..0000000 --- a/test/test_idealplastic_shear.jl +++ /dev/null @@ -1,54 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -using Test, Tensors - -E = 200.0e3 -nu = 0.3 -syield = 100.0 -parameters = IdealPlasticParameterState(youngs_modulus = E, - poissons_ratio = nu, - yield_stress = syield) - -mat = IdealPlastic(parameters=parameters) - -times = [0.0] -loads = [0.0] -dt = 0.5 -G = 0.5*E/(1+nu) -# vm = sqrt(3)*G*ga | ea = ga -ea = 2*syield/(sqrt(3)*G) -# Go to elastic border -push!(times, times[end]+dt) -push!(loads, loads[end] + ea*dt) - # Proceed to plastic flow -push!(times, times[end]+dt) -push!(loads, loads[end] + ea*dt) - # Reverse direction -push!(times, times[end]+dt) -push!(loads, loads[end] - ea*dt) - # Continue and pass yield criterion -push!(times, times[end]+dt) -push!(loads, loads[end] - 2*ea*dt) -stresses = [copy(tovoigt(mat.variables.stress))] -for i=2:length(times) - dtime = times[i]-times[i-1] - dstrain31 = loads[i]-loads[i-1] - dstrain = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain31] - dstrain_ = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0) - ddrivers = IdealPlasticDriverState(time = dtime, strain = dstrain_) - #mat.dtime = dtime - #mat.dstrain = dstrain - mat.ddrivers = ddrivers - integrate_material!(mat) - update_material!(mat) - push!(stresses, copy(tovoigt(mat.variables.stress))) -end - -for i in 1:length(times) - @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6) -end -s31 = [s[6] for s in stresses] - -s31_expected = [0.0, syield/sqrt(3.0), syield/sqrt(3.0), 0.0, -syield/sqrt(3.0)] -@test isapprox(s31, s31_expected; rtol=1.0e-2) diff --git a/test/test_perfectplastic.jl b/test/test_perfectplastic.jl new file mode 100644 index 0000000..1fa8cca --- /dev/null +++ b/test/test_perfectplastic.jl @@ -0,0 +1,52 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + +using Test, Tensors + +let parameters = PerfectPlasticParameterState(youngs_modulus=200.0e3, + poissons_ratio=0.3, + yield_stress=100.0), + epsilon=1e-3, + mat, # scope the name to this level; actual definition follows later + tostrain(vec) = fromvoigt(Symm2, vec; offdiagscale=2.0), + tostress(vec) = fromvoigt(Symm2, vec), + uniaxial_stress(sigma) = tostress([sigma, 0, 0, 0, 0, 0]) + let dtime=0.25 + dstrain_dtime = tostrain(epsilon*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]) + ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_dtime*dtime) + mat = PerfectPlastic(parameters=parameters, ddrivers=ddrivers) + integrate_material!(mat) + update_material!(mat) + @test isapprox(mat.variables.stress, uniaxial_stress(50.0)) + + mat.ddrivers = ddrivers + integrate_material!(mat) + update_material!(mat) + @test isapprox(mat.variables.stress, uniaxial_stress(100.0)) + @test isapprox(mat.variables.cumeq, 0.0; atol=1.0e-12) + + dstrain_dtime = tostrain(epsilon*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]) + ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_dtime*dtime) + mat.ddrivers = ddrivers + integrate_material!(mat) + update_material!(mat) + @test isapprox(mat.variables.stress, uniaxial_stress(100.0); atol=1.0e-12) + @test isapprox(mat.variables.cumeq, dtime*epsilon) + + dstrain_dtime = tostrain(-epsilon*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]) + ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_dtime*dtime) + mat.ddrivers = ddrivers + integrate_material!(mat) + update_material!(mat) + @test isapprox(mat.variables.stress, uniaxial_stress(50.0); atol=1.0e-12) + end + + dstrain_dtime = (-0.75*tostrain(epsilon*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]) + -0.25*tostrain(epsilon*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0])) + ddrivers = PerfectPlasticDriverState(time=1.0, strain=dstrain_dtime*1.0) + mat.ddrivers = ddrivers + integrate_material!(mat) + integrate_material!(mat) + update_material!(mat) + @test isapprox(mat.variables.stress, uniaxial_stress(-100.0)) +end diff --git a/test/test_perfectplastic_shear.jl b/test/test_perfectplastic_shear.jl new file mode 100644 index 0000000..1d7f896 --- /dev/null +++ b/test/test_perfectplastic_shear.jl @@ -0,0 +1,53 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + +using Test, Tensors + +let E = 200.0e3, + nu = 0.3, + yield_stress = 100.0, + parameters = PerfectPlasticParameterState(youngs_modulus=E, + poissons_ratio=nu, + yield_stress=yield_stress), + mat = PerfectPlastic(parameters=parameters), + times = [0.0], + loads = [0.0], + dt = 0.5, + G = 0.5*E/(1+nu), + # vm = sqrt(3)*G*ga | ea = ga + ea = 2*yield_stress/(sqrt(3)*G) + + # Go to elastic border + push!(times, times[end] + dt) + push!(loads, loads[end] + ea*dt) + # Proceed to plastic flow + push!(times, times[end] + dt) + push!(loads, loads[end] + ea*dt) + # Reverse direction + push!(times, times[end] + dt) + push!(loads, loads[end] - ea*dt) + # Continue and pass yield criterion + push!(times, times[end] + dt) + push!(loads, loads[end] - 2*ea*dt) + + stresses = [copy(tovoigt(mat.variables.stress))] + for i=2:length(times) + dtime = times[i] - times[i-1] + dstrain31 = loads[i] - loads[i-1] + dstrain = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain31] + dstrain_ = fromvoigt(Symm2{Float64}, dstrain; offdiagscale=2.0) + ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_) + mat.ddrivers = ddrivers + integrate_material!(mat) + update_material!(mat) + push!(stresses, copy(tovoigt(mat.variables.stress))) + end + + for i in 1:length(times) + @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6) + end + s31 = [s[6] for s in stresses] + + s31_expected = [0.0, yield_stress/sqrt(3.0), yield_stress/sqrt(3.0), 0.0, -yield_stress/sqrt(3.0)] + @test isapprox(s31, s31_expected; rtol=1.0e-2) +end diff --git a/test/test_stress_driven_uniaxial_increment.jl b/test/test_stress_driven_uniaxial_increment.jl index a0d359d..0f5efc1 100644 --- a/test/test_stress_driven_uniaxial_increment.jl +++ b/test/test_stress_driven_uniaxial_increment.jl @@ -1,52 +1,54 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE -using Test, Tensors, Materials -parameters = ChabocheParameterState(E = 200.0e3, - nu = 0.3, - R0 = 100.0, - Kn = 100.0, - nn = 10.0, - C1 = 10000.0, - D1 = 100.0, - C2 = 50000.0, - D2 = 1000.0, - Q = 0.0, - b = 0.1) -material = Chaboche(parameters = parameters) +using Test, Tensors -dtime = 0.25 -times = [material.drivers.time] -stresses = [copy(tovoigt(material.variables.stress))] -strains = [copy(tovoigt(material.drivers.strain; offdiagscale=2.0))] -cumeqs = [copy(material.variables.cumeq)] +let dtime = 0.25, + parameters = ChabocheParameterState(E=200.0e3, + nu=0.3, + R0=100.0, + Kn=100.0, + nn=10.0, + C1=10000.0, + D1=100.0, + C2=50000.0, + D2=1000.0, + Q=0.0, + b=0.1), + material = Chaboche(parameters=parameters), + times = [material.drivers.time], + stresses = [copy(tovoigt(material.variables.stress))], + strains = [copy(tovoigt(material.drivers.strain; offdiagscale=2.0))], + cumeqs = [copy(material.variables.cumeq)], + tostrain(vec) = fromvoigt(Symm2, vec; offdiagscale=2.0), + tostress(vec) = fromvoigt(Symm2, vec), + uniaxial_stress(sigma) = tostress([sigma, 0, 0, 0, 0, 0]), + stresses_expected = [uniaxial_stress(50.0), + uniaxial_stress(100.0), + uniaxial_stress(150.0), + uniaxial_stress(150.0), + uniaxial_stress(100.0), + uniaxial_stress(-100.0)], + dstress = 50.0, + dstresses11 = dstress*[1.0, 1.0, 1.0, 0.0, -1.0, -4.0] + dtimes = [dtime, dtime, dtime, 1e3, dtime, 1e3] -stresses_expected = [[ 50.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [ 100.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [ 150.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [ 150.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [ 100.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [ -100.0, 0.0, 0.0, 0.0, 0.0, 0.0]] -dstress = 50.0 + for i in 1:length(dtimes) + dstress11 = dstresses11[i] + dtime = dtimes[i] + stress_driven_uniaxial_increment!(material, dstress11, dtime) + update_material!(material) + push!(times, material.drivers.time) + push!(stresses, copy(tovoigt(material.variables.stress))) + push!(strains, copy(tovoigt(material.drivers.strain; offdiagscale=2.0))) + push!(cumeqs, copy(material.variables.cumeq)) + @test isapprox(material.variables.stress, stresses_expected[i]; atol=1e-4) + end -dtimes = [dtime, dtime, dtime, 1e3, dtime, 1e3] -dstresses11 = [dstress, dstress, dstress, 0.0, -dstress, -4*dstress] + dstrain_creep = strains[5] - strains[4] + @test isapprox(dstrain_creep[2], -dstrain_creep[1]*0.5; atol=1e-4) + @test isapprox(dstrain_creep[3], -dstrain_creep[1]*0.5; atol=1e-4) -for i in 1:length(dtimes) - dstress11 = dstresses11[i] - dtime = dtimes[i] - stress_driven_uniaxial_increment!(material, dstress11, dtime) - update_material!(material) - push!(times, material.drivers.time) - push!(stresses, copy(tovoigt(material.variables.stress))) - push!(strains, copy(tovoigt(material.drivers.strain; offdiagscale=2.0))) - push!(cumeqs, copy(material.variables.cumeq)) - @test isapprox(tovoigt(material.variables.stress), stresses_expected[i]; atol=1e-4) + dcumeq = cumeqs[end] - cumeqs[end-1] + @test dcumeq > 0 end - -dstrain_creep = strains[5] - strains[4] -@test isapprox(dstrain_creep[2], -dstrain_creep[1]*0.5; atol=1e-4) -@test isapprox(dstrain_creep[3], -dstrain_creep[1]*0.5; atol=1e-4) - -dcumeq = cumeqs[end] - cumeqs[end-1] -@test dcumeq > 0 \ No newline at end of file diff --git a/test/test_uniaxial_increment.jl b/test/test_uniaxial_increment.jl index 0886e4c..c8012b5 100644 --- a/test/test_uniaxial_increment.jl +++ b/test/test_uniaxial_increment.jl @@ -1,42 +1,36 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + using Test, Tensors -dtime = 0.25 -# mat = Material(IdealPlastic, tuple()) -# mat.properties.youngs_modulus = 200.0e3 -# mat.properties.poissons_ratio = 0.3 -# mat.properties.yield_stress = 100.0 -parameters = IdealPlasticParameterState(youngs_modulus = 200.0e3, - poissons_ratio = 0.3, - yield_stress = 100.0) -mat = IdealPlastic(parameters=parameters) -times = [mat.drivers.time] -stresses = [copy(tovoigt(mat.variables.stress))] -stresses_expected = [[50.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [100.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [100.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [50.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [-100.0, 0.0, 0.0, 0.0, 0.0, 0.0]] -dstrain11 = 1e-3*dtime -strains_expected = [[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0], - [2*dstrain11, -0.3*dstrain11*2, -0.3*dstrain11*2, 0.0, 0.0, 0.0], - [3*dstrain11, -0.3*dstrain11*2 - 0.5*dstrain11, -0.3*dstrain11*2 - 0.5*dstrain11, 0.0, 0.0, 0.0], - [2*dstrain11, -0.3*dstrain11 - 0.5*dstrain11, -0.3*dstrain11 - 0.5*dstrain11, 0.0, 0.0, 0.0], - [-2*dstrain11, 0.3*dstrain11*2, 0.3*dstrain11*2, 0.0, 0.0, 0.0]] -dtimes = [dtime, dtime, dtime, dtime, 1.0] -dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11] -for i in 1:length(dtimes) - dstrain11 = dstrains11[i] - dtime = dtimes[i] - uniaxial_increment!(mat, dstrain11, dtime) - # uniaxial_increment!(mat, dstrain11, dtime; dstrain = copy(tovoigt(mat.ddrivers.strain))*dstrain11/mat.ddrivers.strain[1,1]*dtime/mat.ddrivers.time) - # mat.time += mat.dtime - # mat.strain .+= mat.dstrain - # mat.stress .+= mat.dstress - update_material!(mat) - # push!(times, mat.drivers.time) - # push!(stresses, copy(tovoigt(mat.variables.stress))) - #@info(tovoigt(mat.variables.stress), stresses_expected[i]) - @test isapprox(tovoigt(mat.variables.stress), stresses_expected[i]) - @test isapprox(tovoigt(mat.drivers.strain; offdiagscale=2.0), strains_expected[i]) + +let dtime = 0.25, + parameters = PerfectPlasticParameterState(youngs_modulus=200.0e3, + poissons_ratio=0.3, + yield_stress=100.0), + mat = PerfectPlastic(parameters=parameters), + tostrain(vec) = fromvoigt(Symm2, vec; offdiagscale=2.0), + tostress(vec) = fromvoigt(Symm2, vec), + uniaxial_stress(sigma) = tostress([sigma, 0, 0, 0, 0, 0]), + stresses_expected = [uniaxial_stress(50.0), + uniaxial_stress(100.0), + uniaxial_stress(100.0), + uniaxial_stress(50.0), + uniaxial_stress(-100.0)], + dstrain11 = 1e-3*dtime, + strains_expected = [tostrain(dstrain11*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]), + tostrain(dstrain11*[2, -0.3*2, -0.3*2, 0.0, 0.0, 0.0]), + tostrain(dstrain11*[3, -0.3*2 - 0.5, -0.3*2 - 0.5, 0.0, 0.0, 0.0]), + tostrain(dstrain11*[2, -0.3 - 0.5, -0.3 - 0.5, 0.0, 0.0, 0.0]), + tostrain(dstrain11*[-2, 0.3*2, 0.3*2, 0.0, 0.0, 0.0])], + dtimes = [dtime, dtime, dtime, dtime, 1.0], + dstrains11 = dstrain11*[1.0, 1.0, 1.0, -1.0, -4.0] + + for i in 1:length(dtimes) + dstrain11 = dstrains11[i] + dtime = dtimes[i] + uniaxial_increment!(mat, dstrain11, dtime) + update_material!(mat) + @test isapprox(mat.variables.stress, stresses_expected[i]) + @test isapprox(mat.drivers.strain, strains_expected[i]) + end end diff --git a/test/test_utilities.jl b/test/test_utilities.jl new file mode 100644 index 0000000..fcc284f --- /dev/null +++ b/test/test_utilities.jl @@ -0,0 +1,72 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE + +using Test, Tensors, LinearAlgebra + +# Kronecker delta +@test delta(1, 1) == 1 +@test delta(1, 2) == 0 +@test_throws MethodError delta(1.0, 2.0) # indices must be integers +@test_throws MethodError delta(1, BigInt(2)) # both must have the same type +@test delta(1, 2) isa Int # the output type matches the input +@test delta(BigInt(1), BigInt(2)) isa BigInt + +# Various tensors +let Z3 = zeros(3, 3), + O3 = ones(3, 3), + I3 = I(3) + @test isapprox(tovoigt(II()), I(6)) + @test isapprox(tovoigt(IT()), [I3 Z3; + Z3 Z3]) + + @test isapprox(tovoigt(IS()), [I3 Z3; + Z3 1//2*I3]) + @test isapprox(tovoigt(IA()), [Z3 Z3; + Z3 1//2*I3]) + + @test isapprox(tovoigt(IV()), [1//3*O3 Z3; + Z3 Z3]) + @test isapprox(tovoigt(ID()), [(I3 - 1//3*O3) Z3; + Z3 1//2*I3]) + + @test let lambda = 10.0, + mu = 1.0 + isapprox(tovoigt(isotropic_elasticity_tensor(lambda, mu)), [(lambda*O3 + 2*mu*I3) Z3; + Z3 mu*I3]) + end +end + +# Lamé parameters for isotropic solids +@test all(isapprox(result, expected) + for (result, expected) in zip(lame(1e11, 0.3), (5.769230769230769e10, 3.846153846153846e10))) +@test all(isapprox(result, expected) + for (result, expected) in zip(delame(lame(1e11, 0.3)...), (1e11, 0.3))) + +# Mutating function to non-mutating function conversion +let # introduce a local scope so the name `f!` is only defined locally for this test. + function f!(out, x) + out[:] = [sin(elt) for elt in x] + return nothing + end + + let + out = [0.0] + @test all([f!(out, [pi/4]) == nothing, + isapprox(out, [1/sqrt(2)])]) + end + + let + out = [0.0] + f = debang(f!) + @test f isa Function + @test all([isapprox(f([pi/4]), [1/sqrt(2)]), + out == [0.0]]) + end +end + +# Newton root finder +let g(x) = [(1 - x[1]^2) + x[2]], + x0 = [0.8, 0.2] + @test !isapprox(g(x0), [0.0], atol=1e-15) # find_root should have to actually do something + @test isapprox(g(find_root(g, x0)), [0.0], atol=1e-15) +end diff --git a/test/test_viscoplastic.jl b/test/test_viscoplastic.jl deleted file mode 100644 index f2822fb..0000000 --- a/test/test_viscoplastic.jl +++ /dev/null @@ -1,61 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE - -using Materials, Test, ForwardDiff - - -@testset "Calculate plastic strain" begin - # time = [...] - # strain = [...] - # stress = run_simulation(mat, time, strain) - # @test isapprox(stress, stress_expected) - - n = 1.0 - K = 100.0 - - mat = Material(ViscoPlastic, (:norton, [K, n])) - mat.properties.youngs_modulus = 200.0e3 - mat.properties.poissons_ratio = 0.3 - mat.properties.yield_stress = 100.0 - - mat.dtime = 0.25 - mat.dstrain .= 1.0e-3*[1.0, -0.3, -0.3, 0.0, 0.0, 0.0]*mat.dtime - integrate_material!(mat) - @test isapprox(mat.stress+mat.dstress, [50.0, 0.0, 00.0, 0.0, 0.0, 0.0]) - mat.time += mat.dtime - mat.strain .+= mat.dstrain - mat.stress .+= mat.dstress - mat.properties.plastic_strain .+= mat.properties.dplastic_strain - mat.properties.plastic_multiplier += mat.properties.dplastic_multiplier - - - integrate_material!(mat) - @test isapprox(mat.stress+mat.dstress, [100.0, 0.0, 0.0, 0.0, 0.0, 0.0]) - @test isapprox(mat.properties.dplastic_multiplier, 0.0; atol=1.0e-12) - - mat.time += mat.dtime - mat.strain .+= mat.dstrain - mat.stress .+= mat.dstress - mat.properties.plastic_strain .+= mat.properties.dplastic_strain - mat.properties.plastic_multiplier += mat.properties.dplastic_multiplier - - integrate_material!(mat) - @test mat.properties.dplastic_strain[1] > 0 - @test mat.properties.dplastic_strain[2] < 0 - @test mat.properties.dplastic_strain[3] < 0 - - mat.time += mat.dtime - mat.strain .+= mat.dstrain - mat.stress .+= mat.dstress - mat.properties.plastic_strain .+= mat.properties.dplastic_strain - mat.properties.plastic_multiplier += mat.properties.dplastic_multiplier - - stress = mat.stress - dstress = mat.dstress - - @test Materials.von_mises_yield(stress + dstress, 200) <= 0 -end - -@testset "Invalid potential is not defined" begin - @test_throws ErrorException Material(ViscoPlastic, (:notNorton, [0.0, 0.0])) -end diff --git a/test/test_viscoplastic/plot_material.jl b/test/test_viscoplastic/plot_material.jl deleted file mode 100644 index 1be750f..0000000 --- a/test/test_viscoplastic/plot_material.jl +++ /dev/null @@ -1,107 +0,0 @@ -using Materials, Test -using Materials: Simulator -using PyPlot - -num = 100 # Number of steps -E = 200000.0 # Young's modulus [GPa] -poissons_ratio = 0.3 # Poisson's ratio -stress_y = 200.0 # Yield stress -dt = 0.1 -n = 0.92 -K = 180.0e3 - - -# Simulation 1 -mat = Material(ViscoPlastic, (:norton, [K, n])) -mat.properties.youngs_modulus = E -mat.properties.poissons_ratio = poissons_ratio -mat.properties.yield_stress = stress_y -times = [0.0] -strains = [zeros(6)] - -L = 1000.0 # Initial length [mm] -dL = range(0, stop=3.1, length=num) # Change of length [mm] - -strain = collect(dL / L) -# strain_total = vec([strain; reverse(strain[2:end-1]); -strain; reverse(-strain[2:end-1])]) -strain_total = strain -nsteps = length(strain_total) -dt = 0.03 - -for i=1:nsteps - str = zeros(6) - str[1] = strain_total[i] - str[2] = -strain_total[i] * poissons_ratio - str[3] = -strain_total[i] * poissons_ratio - - push!(times, times[end] + dt) - push!(strains, str) -end - -println(strains[2][1] / dt) -sim = Simulator(mat) -Materials.initialize!(sim, strains, times) -Materials.run!(sim) - -s11s = [s[1] for s in sim.stresses] -e11s = [e[1] for e in strains] - -# Simulation 2 -mat = Material(ViscoPlastic, (:norton, [K, n])) -mat.properties.youngs_modulus = E -mat.properties.poissons_ratio = poissons_ratio -mat.properties.yield_stress = stress_y -times = [0.0] -strains = [zeros(6)] - -dt = 0.08 - -for i=1:nsteps - str = zeros(6) - str[1] = strain_total[i] - str[2] = -strain_total[i] * poissons_ratio - str[3] = -strain_total[i] * poissons_ratio - - push!(times, times[end] + dt) - push!(strains, str) -end -println(strains[3][1] / dt) -sim = Simulator(mat) -Materials.initialize!(sim, strains, times) -Materials.run!(sim) - -s11s2 = [s[1] for s in sim.stresses] -e11s2 = [e[1] for e in strains] - -# Simulation 3 -mat = Material(ViscoPlastic, (:norton, [K, n])) -mat.properties.youngs_modulus = E -mat.properties.poissons_ratio = poissons_ratio -mat.properties.yield_stress = stress_y -times = [0.0] -strains = [zeros(6)] - -dt = 0.5 - -for i=1:nsteps - str = zeros(6) - str[1] = strain_total[i] - str[2] = -strain_total[i] * poissons_ratio - str[3] = -strain_total[i] * poissons_ratio - - push!(times, times[end] + dt) - push!(strains, str) -end -println(strains[3][1] / dt) -sim = Simulator(mat) -Materials.initialize!(sim, strains, times) -Materials.run!(sim) - -s11s3 = [s[1] for s in sim.stresses] -e11s3 = [e[1] for e in strains] -legend(bbox_to_anchor=[1.05,1],loc=2,borderaxespad=0) -plot(e11s, s11s, label="de/dt=0.001346") -plot(e11s2, s11s2, label="de/dt=0.00081") -plot(e11s3, s11s3, label="de/dt=0.000008") -grid() -show() \ No newline at end of file