From a180c00a4bff4e695c56b26fc7b1a8d47bacdade Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 7 Jun 2021 18:24:16 +0200 Subject: [PATCH 01/64] Draft of modular small strain plasticity model --- src/MaterialModels.jl | 11 + src/SmallStrainPlasticity/Elasticity.jl | 32 +++ .../IsotropicHardening.jl | 20 ++ .../KinematicHardening.jl | 37 +++ src/SmallStrainPlasticity/Residual.jl | 47 ++++ .../SmallStrainPlasticity.jl | 240 ++++++++++++++++++ src/SmallStrainPlasticity/UtilityFunctions.jl | 14 + 7 files changed, 401 insertions(+) create mode 100644 src/SmallStrainPlasticity/Elasticity.jl create mode 100644 src/SmallStrainPlasticity/IsotropicHardening.jl create mode 100644 src/SmallStrainPlasticity/KinematicHardening.jl create mode 100644 src/SmallStrainPlasticity/Residual.jl create mode 100644 src/SmallStrainPlasticity/SmallStrainPlasticity.jl create mode 100644 src/SmallStrainPlasticity/UtilityFunctions.jl diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 088ba5f..23f4e6e 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -8,6 +8,11 @@ import ForwardDiff import DiffResults using TimerOutputs using StaticArrays +using LinearAlgebra + + + + # Write your package code here. """ @@ -71,6 +76,7 @@ include("Plastic.jl") include("CrystalViscoPlastic/slipsystems.jl") include("CrystalViscoPlastic/CrystalViscoPlastic.jl") include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl") +include("SmallStrainPlasticity/SmallStrainPlasticity.jl") include("nonlinear_solver.jl") @@ -81,4 +87,9 @@ export AbstractMaterial export LinearElastic, Plastic export LinearElasticState, PlasticState +export Chaboche +export Elastic +export Iso_Voce, Iso_Swift +export Kin_AF, Kin_DB, Kin_OW + end diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl new file mode 100644 index 0000000..777cc86 --- /dev/null +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -0,0 +1,32 @@ +using Tensors + +# Elasticity +struct Elastic{T} + G::T # Shear modulus + K::T # Bulk modulus +end +# Overload initialization method to use more common input parameters +# E: Young's modulus, ν: Poissons ratio +function Elastic(E::Number, ν::Number) + T = promote_type(typeof(E), typeof(ν)) + G = E / 2(1 + ν) + K = E / 3(1 - 2ν) + return Elastic{T}(G, K) +end + +# Elastic material +function material_model(cache, material::Elastic, ϵ::T2, state_old, Δt::AbstractFloat) where{T2<:SymmetricTensor{2,3,T}} where T + ν = (3material.K - 2material.G)/(2*(3material.K+material.G)) # Calculate poissons ratio + + σ = 2 * material.G*dev(ϵ) + 3 * material.K*vol(ϵ) # Calculate stress + + # Create stiffness matrix + δ(i,j) = i == j ? 1.0 : 0.0 # helper function + Dfun(i,j,k,l) = 2.0*material.G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l)) + 𝔻 = SymmetricTensor{4, 3}(Dfun) + + # Return updated values + converged = true + state = state_old + return σ, 𝔻, state, converged +end \ No newline at end of file diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl new file mode 100644 index 0000000..6a913af --- /dev/null +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -0,0 +1,20 @@ +# Isotropic hardening +abstract type AbstractIsoHard{T} end + +struct Iso_Voce{T} <:AbstractIsoHard{T} + Hiso::T # Initial hardening modulus + κ∞::T # Saturation stress +end +function IsotropicHardening(param::Iso_Voce, λ::Number) + param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) +end + +struct Iso_Swift{T} <:AbstractIsoHard{T} + K::T + λ0::T + n::T +end +function IsotropicHardening(param::Iso_Swift, λ::Number) + param.K * (param.λ0 + λ)^n +end + diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl new file mode 100644 index 0000000..bd72e26 --- /dev/null +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -0,0 +1,37 @@ +using Tensors + +# Kinematic hardening +abstract type AbstractKinHard{T} end +struct Kin_AF{T} <: AbstractKinHard{T} + Hkin::T # Initial hardening modulus + β∞::T # Saturation stress +end +function KinematicEvolution(param::Kin_AF{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} + param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) +end + +struct Kin_DB{T} <: AbstractKinHard{T} + Hkin::T # Initial hardening modulus + β∞::T # Saturation stress + δ::T # Amount of Armstrong-Frederick hardening +end +function KinematicEvolution(param::Kin_DB{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} + AF_Term = (param.δ/param.β∞) * 𝛃ᵢ # Armstrong Frederick term + BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡𝛃ᵢ)/param.β∞)*ν # Burlet Cailletaud term + return param.Hkin * ((2.0/3.0) * 𝛎 - AF_Term - BC_Term) # Complete evolution +end + +struct Kin_OW{T} <: AbstractKinHard{T} + Hkin::T # Initial hardening modulus + β∞::T # Saturation stress + mexp::T # Ohno Wang exponent +end +function KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} + β_vm = vonMisesDev(𝛃ᵢ) + if β_vm < param.β∞ * eps(T) + return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ + end + mac_term = (macaulay(𝛎⊡𝛃ᵢ) /param.β∞) + exp_term = (β_vm/param.β∞)^param.mexp + return param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ * mac_term * exp_term / β_vm ) +end \ No newline at end of file diff --git a/src/SmallStrainPlasticity/Residual.jl b/src/SmallStrainPlasticity/Residual.jl new file mode 100644 index 0000000..dba60e2 --- /dev/null +++ b/src/SmallStrainPlasticity/Residual.jl @@ -0,0 +1,47 @@ +# General residual function +function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::ChabocheState, ϵ) where{NKin_R} + Δλ = X.λ - old.λ + + σ_vm = vonMisesDev(X.σ_red_dev) + ν = X.σ_red_dev * ((3/2)*σ_vm) + + Φ = yieldCriterion(material, X.σ_red_dev, X.λ) + + σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, Δλ) + + if NKin_R > 0 + β_hat0 = σ_dev - X.σ_red_dev - sum(X.β1) + β_hat1 = X.β1 + β1 = ntuple(i->old.β[i+1] + Δλ * KinematicEvolution(material.kinematic[i+1], ν, β_hat1[i]), Val{NKin_R}()) + else + β_hat0 = σ_dev - X.σ_red_dev + end + + β0 = old.β[1] + Δλ * KinematicEvolution(material.kinematic[1], ν, β_hat0) + + if NKin_R > 0 + σ_red_dev = σ_dev - β0 - sum(β1) + R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev, ntuple(i->β1[i]-β_hat1[i], Val{NKin_R}())) + else + σ_red_dev = σ_dev - β0 + R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) + end + return R +end + +# Specialized for only one backstress (NKin_R=0) +function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState, ϵ) + + σ_vm = vonMisesDev(X.σ_red_dev) + ν = X.σ_red_dev * ((3/2)*σ_vm) + Φ = yieldCriterion(material, σ_vm, X.λ) + + σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, X.λ - old.λ) + β_hat0 = σ_dev - X.σ_red_dev + + β0 = old.β[1] + (X.λ - old.λ) * KinematicEvolution(material.kinematic[1], ν, β_hat0) + + σ_red_dev = σ_dev - β0 + + return ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) +end \ No newline at end of file diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl new file mode 100644 index 0000000..0226930 --- /dev/null +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -0,0 +1,240 @@ +include("UtilityFunctions.jl") +include("Elasticity.jl") +include("IsotropicHardening.jl") +include("KinematicHardening.jl") + + +# Definition of material properties +struct Chaboche{T,ElasticType,IsoType,KinType} + elastic::ElasticType # Elastic definition + σ_y0::T # Initial yield limit + isotropic::IsoType # Tuple of isotropic hardening definitions + kinematic::KinType # Tuple of kinematic hardening definitions +end + +# Definition of material state +struct ChabocheState{Nkin,T,N} + ϵₚ::SymmetricTensor{2,3,T,N} + λ::T + β::NTuple{Nkin, SymmetricTensor{2,3,T,N}} +end + +struct ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens} + λ::Tλ + σ_red_dev::SymmetricTensor{2,3,Tσ,N_tens} + β1::NTuple{NKin_R, SymmetricTensor{2,3,Tβ,N_tens}} +end +# Specialize for only one backstress +function ChabocheResidual(λ::Tλ,σ_red_dev::SymmetricTensor{2,3,Tσ,N_tens}) where {Tλ,Tσ,N_tens} + ChabocheResidual{0,Tλ,Tσ,Float64,N_tens}(λ,σ_red_dev,()) +end + +Tensors.get_base(::Type{ChabocheResidual}) = ChabocheResidual # needed for frommandel + +function Tensors.tomandel!(v::AbstractVector{T}, r::ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens}) where {T,NKin_R,Tλ,Tσ,Tβ,N_tens} + v[1] = r.λ + tomandel!(v, r.σ_red_dev, offset=1) + for i=1:NKin_R + tomandel!(v, r.β1[i], offset=1+N_tens*6) + end + return v +end + +function Tensors.frommandel(::Type{ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens}}, v::AbstractVector{Tv}) where {Tv,NKin_R,Tλ,Tσ,Tβ,N_tens} + λ = v[1] + σ_red_dev = frommandel(SymmetricTensor{2,3,Tv}, v, offset=1) + if NKin_R > 0 + β1 = ntuple(i->frommandel(SymmetricTensor{2,3,Tv}, v, offset=1+N_tens*i), NKin_R) + return ChabocheResidual(λ,σ_red_dev,β1) + else + return ChabocheResidual(λ,σ_red_dev) + end +end + +function initial_material_state(material::Chaboche{T}) where {T} + ChabocheState(zero(SymmetricTensor{2,3,T}), 0.0, ntuple(i->zero(SymmetricTensor{2,3,T}), Val{length(material.kinematic)}())) +end + +# Definition of material cache +##= +struct ChabocheCache{T, nx, nx6, NL_TF, NL_TDF, NL_TX, NL_NC} + # General purpose + vnx6::MMatrix{nx,6,T,nx6} + v6xn::MMatrix{6,nx,T,nx6} + v6x6::MMatrix{6,6,T,36} + v6::MVector{6,T} + # For solving R(X)=0 + #X0::MVector{nx,T} + R_X_oncediff::OnceDifferentiable{NL_TF, NL_TDF, NL_TX} + R_X_newton::NLsolve.NewtonCache{NL_NC} +end + +function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType,IsoType,KinType} + nx = 1 + 6*length(material.kinematic) + + # Construct residual function and create OnceDifferentiable object + state_tmp = initial_material_state(material) + σ_trial_dev = zero(SymmetricTensor{2,3,T}) + X_tensor = initial_guess(material, state_tmp, zero(SymmetricTensor{2,3,T})) + rf_tens(X_tensor) = residual(X_tensor, material, state_tmp, σ_trial_dev) + rf!(R, X) = vector_residual!(rf_tens, R, X, typeof(X_tensor)) + X0 = MVector{nx}(zeros(T, nx)) + # X0 only for shape and type information here: + R_X_oncediff = OnceDifferentiable(rf!, X0, X0; autodiff = :forward) + R_X_newton = NLsolve.NewtonCache(R_X_oncediff) + return ChabocheCache(MMatrix{nx,6}(zeros(nx,6)), MMatrix{6,nx}(zeros(6,nx)), MMatrix{6,6}(zeros(6,6)), MVector{6}(zeros(6)), + R_X_oncediff, R_X_newton) +end + +# Material model +function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} + + σ_trial, dσdϵ_elastic, _, _ = material_model(cache, material.elastic, ϵ-state_old.ϵₚ, nothing, Δt) + + Φ_trial = yieldCriterion(material, σ_trial-sum(state_old.β), state_old.λ) + if Φ_trial < 0 + return σ_trial, dσdϵ_elastic, state_old, true + else + num_blas_threads = LinearAlgebra.BLAS.get_num_threads() # (~ 2ns) + LinearAlgebra.BLAS.set_num_threads(1) # Big performance benefit, takes ~8ns + converged = solve_local_problem!(cache, material, state_old, ϵ, options) + if converged + σ, dσdϵ, state = get_plastic_output(cache, material, state_old, σ_trial, dσdϵ_elastic, options) + else + σ, dσdϵ, state = (σ_trial, dσdϵ_elastic, state_old) + println("Did not converge!") + end + LinearAlgebra.BLAS.set_num_threads(num_blas_threads) + return σ, dσdϵ, state, converged + end + +end + +include("Residual.jl") + +function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}, options::Dict{Symbol, Any}) where {NKin,Ts,N} + # σ = σ(X(ϵ), ϵ) yields + # dσ/dϵ = ∂σ/∂ϵ + ∂σ/∂X : dX/dϵ [1] + # R = R(X(ϵ), ϵ) yields + # dR/dϵ = 0 = ∂R/∂ϵ + ∂R/∂X : dXdϵ [2] + # Solve [2] for dX/dϵ + # dX/dϵ = - [∂R/∂X]^-1 : ∂R∂ϵ [3] + # Insert [3] in [1] + # dσ/dϵ = ∂σ/∂ϵ - ∂σ/∂X : [∂R/∂X]^-1 : ∂R/∂ϵ [4] + + X_vec = cache.R_X_oncediff.x_f + dRdX = cache.R_X_oncediff.DF + X_tensor = frommandel(ChabocheResidual{NKin-1,Ts,Ts,Ts,N}, X_vec) + + # ∂σ∂X + # - Stress function + σ_X(X_arg) = get_sigma(material, state_old, X_arg, ϵ) + # - Stress (vector) function: + σ_X_vec!(σv_arg, Xv_arg) = vector_residual!(σ_X, σv_arg, Xv_arg, typeof(X_tensor)) + σ_vec = cache.v6 + # - Preallocate (should have been done beforehand, problem with Dual Tag values?) + cfg = ForwardDiff.JacobianConfig(σ_X_vec!, σ_vec, X_vec, ForwardDiff.Chunk{length(X_vec)}()) + # - Create DiffResult (this should be non-allocating) + ∂σ∂X = cache.v6xn + diff_result = DiffResults.MutableDiffResult(σ_vec, (∂σ∂X,)) + # - Calculate σ and ∂σ∂X + ForwardDiff.jacobian!(diff_result, σ_X_vec!, σ_vec, X_vec, cfg) + σ = frommandel(SymmetricTensor{2,3}, diff_result.value) + + # ∂R/∂ϵ + # - Specialized residual (tensor) function: + R_ϵ(ϵ_arg) = residual(X_tensor, material, state_old, ϵ_arg) + # - Specialized residual (vector) function: + R_ϵ_vec!(Rv_arg, ϵv_arg) = vector_residual!(R_ϵ, Rv_arg, ϵv_arg, SymmetricTensor{2,3}) + ϵ_vec = cache.v6 # Use cache value (give name that makes more sense) + tomandel!(ϵ_vec, ϵ) + # - Preallocate (should have been done beforehand, problem with Dual Tag values?) + R_vec = X_vec # Use as cache (ok as X_vec is not used anymore) + cfg = ForwardDiff.JacobianConfig(R_ϵ_vec!, R_vec, ϵ_vec, ForwardDiff.Chunk{length(6)}()) + # - Calculate ∂R∂ϵ + ∂R∂ϵ = cache.vnx6 + ForwardDiff.jacobian!(∂R∂ϵ, R_ϵ_vec!, R_vec, ϵ_vec, cfg) + + # Calculate full tangent stiffness + dσdϵ = dσdϵ_elastic - frommandel(SymmetricTensor{4,3}, ∂σ∂X*(dRdX\∂R∂ϵ)) + + λ = X_tensor.λ + Δλ = λ-state_old.λ + σ_red_dev = X_tensor.σ_red_dev + ν = (3.0/2.0) * σ_red_dev / vonMisesDev(σ_red_dev) + + if NKin > 1 + β = ntuple(i-> i==1 ? dev(σ) - σ_red_dev - sum(X_tensor.β1) : X_tensor.β1[i-1], NKin) + else + β = (dev(σ) - σ_red_dev,) + end + + ϵₚ = state_old.ϵₚ + Δλ * ν + state = ChabocheState(ϵₚ, λ, β) + + return σ, dσdϵ, state +end + +function initial_guess(material::Chaboche, state_old::ChabocheState{Nkin, T, N}, ϵ) where {T, Nkin, N} + if Nkin<1 + error("Nkin < 1 is not supported") + end + # Becomes trial by setting Δλ=0 + σ_trial_dev = calc_sigma_dev(material.elastic, state_old, ϵ, ϵ, 0.0) + λ = state_old.λ + σ_red_trial = σ_trial_dev - sum(state_old.β) + + if Nkin > 1 + return ChabocheResidual(λ,σ_red_trial,ntuple(i->state_old.β[i], Val{Nkin-1}())) + else + return ChabocheResidual(λ,σ_red_trial) + end +end + + +function solve_local_problem!(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{Nkin,Ts,N}, ϵ::SymmetricTensor{2,3}, options::Dict{Symbol, Any}) where {Ts,Nkin,N} + + X_tensor = initial_guess(material, state_old, ϵ) + rf_tens(X_tensor_arg) = residual(X_tensor_arg, material, state_old, ϵ) + rf!(R, X) = vector_residual!(rf_tens, R, X, typeof(X_tensor)) + update_cache!(cache.R_X_oncediff, rf!) + + tomandel!(cache.R_X_oncediff.x_f, X_tensor) + # Should this be centrally managed? I.e. process_options or similar? + nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton)) + haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method + nlsolve_options[:method] == :newton || merge!(nlsolve_options, Dict(:cache=>cache.R_X_newton)) + # Need to call newton directly to allow newton caching... + #result = NLsolve.newton(; ftol=1.e-6, cache=cache.R_X_newton) + my_newton(df, x0; xtol=0.0, ftol=1.e-8, iterations=100, store_trace=false, show_trace=false, extended_trace=false, linesearch=NLsolve.LineSearches.Static(),cache=NewtonCache(df)) = NLsolve.newton(df, x0, xtol, ftol, iterations, store_trace, show_trace, extended_trace, linesearch,cache) + result = my_newton(cache.R_X_oncediff, cache.R_X_oncediff.x_f, cache=cache.R_X_newton) + # Is this necessary? + cache.R_X_oncediff.x_f = result.zero + + return result.f_converged +end + +function get_sigma(material::Chaboche, state_old::ChabocheState, X::ChabocheResidual, ϵ::SymmetricTensor{2,3}) + Δλ = X.λ - state_old.λ + σ_red_dev = X.σ_red_dev + ν = (3.0/2.0) * σ_red_dev / vonMises(σ_red_dev) + σ = calc_sigma(material.elastic, state_old, ϵ, ν, Δλ) + return σ +end + +function calc_sigma(material::Elastic, state_old::ChabocheState, ϵ, ν, Δλ) + return 3 * material.K*vol(ϵ) + calc_sigma_dev(material, state_old, ϵ, Δλ, ν) +end + +function calc_sigma_dev(material::Elastic, state_old::ChabocheState, ϵ, ν, Δλ) + return 2 * material.G * (dev(ϵ - state_old.ϵₚ) - Δλ*ν) +end + +function yieldCriterion(material::Chaboche{Tp,ElType,IsoType,KinType}, σ_vm_red::Number, λ) where {Tp,ElType,IsoType<:NTuple{Niso,Any},KinType} where {Niso} + κ = sum(ntuple(i->IsotropicHardening(material.isotropic[i], λ), Val{Niso}())) + return σ_vm_red - (κ + material.σ_y0) +end + +function yieldCriterion(material::Chaboche, σ_red_dev::AbstractTensor, λ) + return yieldCriterion(material, vonMises(σ_red_dev), λ) +end diff --git a/src/SmallStrainPlasticity/UtilityFunctions.jl b/src/SmallStrainPlasticity/UtilityFunctions.jl new file mode 100644 index 0000000..3873491 --- /dev/null +++ b/src/SmallStrainPlasticity/UtilityFunctions.jl @@ -0,0 +1,14 @@ +# Should be fixed in ForwardDiff? +DiffResults.DiffResult(value::MArray, derivs::Tuple{Vararg{MArray}}) = DiffResults.MutableDiffResult(value, derivs) + +# Generic functions, should be defined elsewhere? +function vonMises(𝛔::SymmetricTensor{2,3}) + 𝛔_dev = dev(𝛔) + return sqrt((3.0/2.0) * (𝛔_dev ⊡ 𝛔_dev)) +end + +function vonMisesDev(𝛔_dev::SymmetricTensor{2,3}) + return sqrt((3.0/2.0) * (𝛔_dev ⊡ 𝛔_dev)) +end + +macaulay(x::T) where {T} = x > 0.0 ? x : zero(T) \ No newline at end of file From 230a9c40a002c3320f4d42f5b806f8a5a96b0b81 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 17:39:02 +0200 Subject: [PATCH 02/64] Added keyword argument constructor and corrected naming for elastic model --- src/SmallStrainPlasticity/Elasticity.jl | 7 +++++-- src/SmallStrainPlasticity/IsotropicHardening.jl | 6 ++++++ src/SmallStrainPlasticity/KinematicHardening.jl | 9 +++++++++ src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 2 +- 4 files changed, 21 insertions(+), 3 deletions(-) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index 777cc86..3b405f9 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -1,6 +1,6 @@ using Tensors -# Elasticity +# Linear isotropic elasticity struct Elastic{T} G::T # Shear modulus K::T # Bulk modulus @@ -13,9 +13,12 @@ function Elastic(E::Number, ν::Number) K = E / 3(1 - 2ν) return Elastic{T}(G, K) end +Elastic(;E, ν) = Elastic(E, ν) # Keyword argument constructor # Elastic material -function material_model(cache, material::Elastic, ϵ::T2, state_old, Δt::AbstractFloat) where{T2<:SymmetricTensor{2,3,T}} where T +get_cache(::Elastic) = nothing + +function material_response(material::Elastic, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} ν = (3material.K - 2material.G)/(2*(3material.K+material.G)) # Calculate poissons ratio σ = 2 * material.G*dev(ϵ) + 3 * material.K*vol(ϵ) # Calculate stress diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl index 6a913af..092cdc7 100644 --- a/src/SmallStrainPlasticity/IsotropicHardening.jl +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -1,19 +1,25 @@ # Isotropic hardening abstract type AbstractIsoHard{T} end +# Voce type of isotropic hardening (exponentially saturating) struct Iso_Voce{T} <:AbstractIsoHard{T} Hiso::T # Initial hardening modulus κ∞::T # Saturation stress end +Iso_Voce(;Hiso, κ∞) = Iso_Voce(Hiso, κ∞) # Keyword argument constructor + function IsotropicHardening(param::Iso_Voce, λ::Number) param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) end +# Swift type of kinematic hardening (power law) struct Iso_Swift{T} <:AbstractIsoHard{T} K::T λ0::T n::T end +Iso_Swift(;K, λ0, n) = Iso_Swift(K, λ0, n) # Keyword argument constructor + function IsotropicHardening(param::Iso_Swift, λ::Number) param.K * (param.λ0 + λ)^n end diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index bd72e26..3100708 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -2,30 +2,39 @@ using Tensors # Kinematic hardening abstract type AbstractKinHard{T} end + +# Armstrong-Frederick struct Kin_AF{T} <: AbstractKinHard{T} Hkin::T # Initial hardening modulus β∞::T # Saturation stress end +Kin_AF(;Hkin, β∞) = Kin_AF(Hkin, β∞) # Keyword argument constructor + function KinematicEvolution(param::Kin_AF{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) end +# Delobelle (Combination of Armstrong-Frederick and Burlet-Cailletaud) struct Kin_DB{T} <: AbstractKinHard{T} Hkin::T # Initial hardening modulus β∞::T # Saturation stress δ::T # Amount of Armstrong-Frederick hardening end +Kin_DB(;Hkin, β∞, δ) = Kin_DB(Hkin, β∞, δ) # Keyword argument constructor function KinematicEvolution(param::Kin_DB{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} AF_Term = (param.δ/param.β∞) * 𝛃ᵢ # Armstrong Frederick term BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡𝛃ᵢ)/param.β∞)*ν # Burlet Cailletaud term return param.Hkin * ((2.0/3.0) * 𝛎 - AF_Term - BC_Term) # Complete evolution end +# Ohno-Wang struct Kin_OW{T} <: AbstractKinHard{T} Hkin::T # Initial hardening modulus β∞::T # Saturation stress mexp::T # Ohno Wang exponent end +Kin_OW(;Hkin, β∞, mexp) = Kin_OW(Hkin, β∞, mexp) # Keyword argument constructor + function KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} β_vm = vonMisesDev(𝛃ᵢ) if β_vm < param.β∞ * eps(T) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 0226930..35f5852 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -89,7 +89,7 @@ end # Material model function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} - σ_trial, dσdϵ_elastic, _, _ = material_model(cache, material.elastic, ϵ-state_old.ϵₚ, nothing, Δt) + σ_trial, dσdϵ_elastic, _, _ = material_response(material.elastic, ϵ-state_old.ϵₚ, nothing, Δt) Φ_trial = yieldCriterion(material, σ_trial-sum(state_old.β), state_old.λ) if Φ_trial < 0 From 6126f7813507b3fe7db46b029f1f21ba48cc2ae2 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 20:24:47 +0200 Subject: [PATCH 03/64] Fixed bugs caused by merging to package functions --- Project.toml | 1 + src/SmallStrainPlasticity/Elasticity.jl | 2 +- .../SmallStrainPlasticity.jl | 19 ++++++++++--------- src/nonlinear_solver.jl | 2 +- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 4a0ae9d..01a96d8 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index 3b405f9..44830ab 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -18,7 +18,7 @@ Elastic(;E, ν) = Elastic(E, ν) # Keyword argument constructor # Elastic material get_cache(::Elastic) = nothing -function material_response(material::Elastic, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} +function material_response(material::Elastic, ϵ::SymmetricTensor{2,3}, state_old, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) ν = (3material.K - 2material.G)/(2*(3material.K+material.G)) # Calculate poissons ratio σ = 2 * material.G*dev(ϵ) + 3 * material.K*vol(ϵ) # Calculate stress diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 35f5852..2b6e305 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -11,6 +11,7 @@ struct Chaboche{T,ElasticType,IsoType,KinType} isotropic::IsoType # Tuple of isotropic hardening definitions kinematic::KinType # Tuple of kinematic hardening definitions end +Chaboche(;elastic, σ_y0, isotropic, kinematic) = Chaboche(elastic, σ_y0, isotropic, kinematic) # Definition of material state struct ChabocheState{Nkin,T,N} @@ -29,22 +30,22 @@ function ChabocheResidual(λ::Tλ,σ_red_dev::SymmetricTensor{2,3,Tσ,N_tens}) w ChabocheResidual{0,Tλ,Tσ,Float64,N_tens}(λ,σ_red_dev,()) end -Tensors.get_base(::Type{ChabocheResidual}) = ChabocheResidual # needed for frommandel +Tensors.get_base(::Type{<:ChabocheResidual{NKin_R}}) where{NKin_R} = ChabocheResidual{NKin_R} # needed for frommandel function Tensors.tomandel!(v::AbstractVector{T}, r::ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens}) where {T,NKin_R,Tλ,Tσ,Tβ,N_tens} v[1] = r.λ tomandel!(v, r.σ_red_dev, offset=1) for i=1:NKin_R - tomandel!(v, r.β1[i], offset=1+N_tens*6) + tomandel!(v, r.β1[i], offset=1+N_tens*i) end return v end -function Tensors.frommandel(::Type{ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens}}, v::AbstractVector{Tv}) where {Tv,NKin_R,Tλ,Tσ,Tβ,N_tens} +function Tensors.frommandel(::Type{<:ChabocheResidual{NKin_R}}, v::AbstractVector{Tv}) where {Tv,NKin_R} λ = v[1] - σ_red_dev = frommandel(SymmetricTensor{2,3,Tv}, v, offset=1) + σ_red_dev = frommandel(SymmetricTensor{2,3}, v, offset=1) if NKin_R > 0 - β1 = ntuple(i->frommandel(SymmetricTensor{2,3,Tv}, v, offset=1+N_tens*i), NKin_R) + β1 = ntuple(i->frommandel(SymmetricTensor{2,3,Tv}, v, offset=1+6*i), NKin_R) return ChabocheResidual(λ,σ_red_dev,β1) else return ChabocheResidual(λ,σ_red_dev) @@ -77,7 +78,7 @@ function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType σ_trial_dev = zero(SymmetricTensor{2,3,T}) X_tensor = initial_guess(material, state_tmp, zero(SymmetricTensor{2,3,T})) rf_tens(X_tensor) = residual(X_tensor, material, state_tmp, σ_trial_dev) - rf!(R, X) = vector_residual!(rf_tens, R, X, typeof(X_tensor)) + rf!(R, X) = vector_residual!(rf_tens, R, X, X_tensor) X0 = MVector{nx}(zeros(T, nx)) # X0 only for shape and type information here: R_X_oncediff = OnceDifferentiable(rf!, X0, X0; autodiff = :forward) @@ -130,7 +131,7 @@ function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old: # - Stress function σ_X(X_arg) = get_sigma(material, state_old, X_arg, ϵ) # - Stress (vector) function: - σ_X_vec!(σv_arg, Xv_arg) = vector_residual!(σ_X, σv_arg, Xv_arg, typeof(X_tensor)) + σ_X_vec!(σv_arg, Xv_arg) = vector_residual!(σ_X, σv_arg, Xv_arg, X_tensor) σ_vec = cache.v6 # - Preallocate (should have been done beforehand, problem with Dual Tag values?) cfg = ForwardDiff.JacobianConfig(σ_X_vec!, σ_vec, X_vec, ForwardDiff.Chunk{length(X_vec)}()) @@ -145,7 +146,7 @@ function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old: # - Specialized residual (tensor) function: R_ϵ(ϵ_arg) = residual(X_tensor, material, state_old, ϵ_arg) # - Specialized residual (vector) function: - R_ϵ_vec!(Rv_arg, ϵv_arg) = vector_residual!(R_ϵ, Rv_arg, ϵv_arg, SymmetricTensor{2,3}) + R_ϵ_vec!(Rv_arg, ϵv_arg) = vector_residual!(R_ϵ, Rv_arg, ϵv_arg, ϵ) ϵ_vec = cache.v6 # Use cache value (give name that makes more sense) tomandel!(ϵ_vec, ϵ) # - Preallocate (should have been done beforehand, problem with Dual Tag values?) @@ -196,7 +197,7 @@ function solve_local_problem!(cache::ChabocheCache, material::Chaboche, state_ol X_tensor = initial_guess(material, state_old, ϵ) rf_tens(X_tensor_arg) = residual(X_tensor_arg, material, state_old, ϵ) - rf!(R, X) = vector_residual!(rf_tens, R, X, typeof(X_tensor)) + rf!(R, X) = vector_residual!(rf_tens, R, X, X_tensor) update_cache!(cache.R_X_oncediff, rf!) tomandel!(cache.R_X_oncediff.x_f, X_tensor) diff --git a/src/nonlinear_solver.jl b/src/nonlinear_solver.jl index d7da852..991188c 100644 --- a/src/nonlinear_solver.jl +++ b/src/nonlinear_solver.jl @@ -18,7 +18,7 @@ # return start + 1 # end -function vector_residual!(R::Function, r_vector::Vector{T}, x_vector::Vector{T}, m) where T +function vector_residual!(R::Function, r_vector::AbstractVector{T}, x_vector::AbstractVector{T}, m) where T # construct residuals with type T x_tensor = frommandel(Tensors.get_base(typeof(m)), x_vector) r_tensor = R(x_tensor) From e033d31d7a726337753f8a5eb7ca66441a9a14dd Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 20:25:22 +0200 Subject: [PATCH 04/64] Added a basic test showing a simple example --- test/runtests.jl | 1 + test/test_small_strain_plasticity.jl | 28 ++++++++++++++++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 test/test_small_strain_plasticity.jl diff --git a/test/runtests.jl b/test/runtests.jl index e4ad04d..141bc59 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,3 +6,4 @@ using Rotations include("test_utils.jl") include("test_linear_elastic.jl") include("test_plastic.jl") +include("test_small_strain_plasticity.jl") diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl new file mode 100644 index 0000000..d77b869 --- /dev/null +++ b/test/test_small_strain_plasticity.jl @@ -0,0 +1,28 @@ +@testset "SmallStrainPlasticity" begin + # constructor + m = Chaboche(elastic=Elastic(E=210.e3, ν=0.3), + σ_y0=100.0, + isotropic=(Iso_Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) + kinematic=(Kin_AF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple + ) + cache = get_cache(m) + + # initial state (not used here) + state = initial_material_state(m) + + # strain at yield point for uniaxial stress some test case at some plastic strain: + ϵ = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.004 : -0.000782017) : 0.0) + + ϵₚ_old = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.00153388 : -0.000766938) : 0.0) + λ_old = 0.0015338757291717328 + β_old = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 133.268 : -66.6339) : 0.0) + + state_old = MaterialModels.ChabocheState(ϵₚ_old,λ_old, (β_old,)) + + Δt = 1.0 # No influence... + + σ, ∂σ∂ε, temp_state, converged = material_response(m, ϵ, state, Δt; cache=cache) + + @test converged + +end \ No newline at end of file From f1d7ba0085fd03d575bd2dca67712c23630152d4 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 21:04:05 +0200 Subject: [PATCH 05/64] Added some docstrings --- .../KinematicHardening.jl | 45 +++++++++++++++++-- .../SmallStrainPlasticity.jl | 41 +++++++++++++++-- 2 files changed, 78 insertions(+), 8 deletions(-) diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index 3100708..0874c05 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -10,7 +10,16 @@ struct Kin_AF{T} <: AbstractKinHard{T} end Kin_AF(;Hkin, β∞) = Kin_AF(Hkin, β∞) # Keyword argument constructor -function KinematicEvolution(param::Kin_AF{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} +""" + KinematicEvolution(param::Kin_AF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + + Armstrong-Frederick kinematic hardening law + + ```math + g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) + ``` +""" +function KinematicEvolution(param::Kin_AF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) end @@ -21,7 +30,21 @@ struct Kin_DB{T} <: AbstractKinHard{T} δ::T # Amount of Armstrong-Frederick hardening end Kin_DB(;Hkin, β∞, δ) = Kin_DB(Hkin, β∞, δ) # Keyword argument constructor -function KinematicEvolution(param::Kin_DB{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} + +""" + KinematicEvolution(param::Kin_DB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + + Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law + + ```math + g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin \\left[\\frac{2}{3}\\boldsymbol{\\nu} + - \\delta\\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty} + - \\frac{2}{3\\beta_\\infty}\\left[1 - \\delta\\right]\\left[\\boldsymbol{\\nu}:\\boldsymbol{\\beta}_i\right]\\boldsymbol{\\nu} + \\right] + ``` + +""" +function KinematicEvolution(param::Kin_DB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) AF_Term = (param.δ/param.β∞) * 𝛃ᵢ # Armstrong Frederick term BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡𝛃ᵢ)/param.β∞)*ν # Burlet Cailletaud term return param.Hkin * ((2.0/3.0) * 𝛎 - AF_Term - BC_Term) # Complete evolution @@ -35,9 +58,23 @@ struct Kin_OW{T} <: AbstractKinHard{T} end Kin_OW(;Hkin, β∞, mexp) = Kin_OW(Hkin, β∞, mexp) # Keyword argument constructor -function KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor{dim,Tν}, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tν,Tβ,dim} +""" + KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} + + Kinematic hardening law according to Ohno-Wang + + ```math + g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin \\left[\\frac{2}{3}\\boldsymbol{\\nu} + - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty} + \\frac{\\langle \\boldsymbol{\\nu}:\\boldsymbol{\\beta}_i \\rangle}{\\beta_\\infty} + \\left[\\frac{\\beta_\\mathrm{vM}}{\\beta_\\infty}\right]^\\mathrm{mexp} + \\right] + ``` + +""" +function KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} β_vm = vonMisesDev(𝛃ᵢ) - if β_vm < param.β∞ * eps(T) + if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ end mac_term = (macaulay(𝛎⊡𝛃ᵢ) /param.β∞) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 2b6e305..d1993fb 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -87,8 +87,41 @@ function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType R_X_oncediff, R_X_newton) end -# Material model -function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} +""" + material_response(m::Plastic, ϵ::SymmetricTensor{2,3}, state::ChabocheState, Δt; ) + +Return the stress tensor, stress tangent, the new `MaterialState` and boolean if local iterations converged for the given strain ε and previous material state `state`. + +Von Mises yield function: +```math +\\Phi = \\sqrt{\\frac{3}{2}} \\left| \\text{dev} \\left( \\boldsymbol{\\sigma} - \\boldsymbol{\\alpha} \\right) \\right| - \\sigma_y - \\kappa +``` +Associative plastic flow: +```math +\\dot{\\epsilon}_{\\mathrm{p}} = \\dot{\\lambda} \\frac{\\partial \\Phi}{\\boldsymbol{\\sigma}} += \\dot{\\lambda} \\boldsymbol{\\nu} +``` + +Isotropic hardening is formulated as +```math +\\kappa = \\sum_{i=1}^{N_{\\mathrm{iso}}} g_{\\mathrm{iso},i}(\\lambda) +``` +where ``$g_{\\mathrm{iso},i}(\\lambda)$`` is specified by structs of subtype to `AbstractIsoHard` + +Kinematic hardening is formulated as +```math +\\boldsymbol{\\beta}_i = \\dot{\\lambda} g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) +``` +where ``$g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)$`` is specified by structs of subtype to `AbstractKinHard` +and ``i\\in[1,N_\\mathrm{kin}]``. + +``` +# Keyword arguments +- `cache`: Cache for the iterative solver, used by NLsolve.jl. It is strongly recommended to pre-allocate the cache for repeated calls to `material_response`. See [`get_cache`](@ref). +- `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over. +See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-options). By default the Newton solver will be used. +""" +function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} σ_trial, dσdϵ_elastic, _, _ = material_response(material.elastic, ϵ-state_old.ϵₚ, nothing, Δt) @@ -100,7 +133,7 @@ function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_o LinearAlgebra.BLAS.set_num_threads(1) # Big performance benefit, takes ~8ns converged = solve_local_problem!(cache, material, state_old, ϵ, options) if converged - σ, dσdϵ, state = get_plastic_output(cache, material, state_old, σ_trial, dσdϵ_elastic, options) + σ, dσdϵ, state = get_plastic_output(cache, material, state_old, σ_trial, dσdϵ_elastic) else σ, dσdϵ, state = (σ_trial, dσdϵ_elastic, state_old) println("Did not converge!") @@ -113,7 +146,7 @@ end include("Residual.jl") -function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}, options::Dict{Symbol, Any}) where {NKin,Ts,N} +function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}) where {NKin,Ts,N} # σ = σ(X(ϵ), ϵ) yields # dσ/dϵ = ∂σ/∂ϵ + ∂σ/∂X : dX/dϵ [1] # R = R(X(ϵ), ϵ) yields From 1ca196c83c3b0d55674c66e86d504734912dcf83 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 21:13:26 +0200 Subject: [PATCH 06/64] Fixed minor bug in documentation --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index d1993fb..154c920 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -106,13 +106,13 @@ Isotropic hardening is formulated as ```math \\kappa = \\sum_{i=1}^{N_{\\mathrm{iso}}} g_{\\mathrm{iso},i}(\\lambda) ``` -where ``$g_{\\mathrm{iso},i}(\\lambda)$`` is specified by structs of subtype to `AbstractIsoHard` +where ``g_{\\mathrm{iso},i}(\\lambda)`` is specified by structs of subtype to `AbstractIsoHard` Kinematic hardening is formulated as ```math \\boldsymbol{\\beta}_i = \\dot{\\lambda} g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) ``` -where ``$g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)$`` is specified by structs of subtype to `AbstractKinHard` +where ``g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)`` is specified by structs of subtype to `AbstractKinHard` and ``i\\in[1,N_\\mathrm{kin}]``. ``` @@ -256,11 +256,11 @@ function get_sigma(material::Chaboche, state_old::ChabocheState, X::ChabocheResi return σ end -function calc_sigma(material::Elastic, state_old::ChabocheState, ϵ, ν, Δλ) +function calc_sigma(material::LinearIsotropicElasticity, state_old::ChabocheState, ϵ, ν, Δλ) return 3 * material.K*vol(ϵ) + calc_sigma_dev(material, state_old, ϵ, Δλ, ν) end -function calc_sigma_dev(material::Elastic, state_old::ChabocheState, ϵ, ν, Δλ) +function calc_sigma_dev(material::LinearIsotropicElasticity, state_old::ChabocheState, ϵ, ν, Δλ) return 2 * material.G * (dev(ϵ - state_old.ϵₚ) - Δλ*ν) end From 069cf1d471db60d776830e39c0047e57bd56d2f7 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 21:13:53 +0200 Subject: [PATCH 07/64] Updated naming of Elastic to LinearIsotropicElasticity --- src/MaterialModels.jl | 2 +- src/SmallStrainPlasticity/Elasticity.jl | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 23f4e6e..550953b 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -88,7 +88,7 @@ export LinearElastic, Plastic export LinearElasticState, PlasticState export Chaboche -export Elastic +export LinearIsotropicElasticity export Iso_Voce, Iso_Swift export Kin_AF, Kin_DB, Kin_OW diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index 44830ab..f84bf88 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -1,31 +1,31 @@ using Tensors - +abstract type AbstractElasticity{T} end # Linear isotropic elasticity -struct Elastic{T} +struct LinearIsotropicElasticity{T} <:AbstractElasticity{T} G::T # Shear modulus K::T # Bulk modulus end # Overload initialization method to use more common input parameters # E: Young's modulus, ν: Poissons ratio -function Elastic(E::Number, ν::Number) +function LinearIsotropicElasticity(E::Number, ν::Number) T = promote_type(typeof(E), typeof(ν)) G = E / 2(1 + ν) K = E / 3(1 - 2ν) - return Elastic{T}(G, K) + return LinearIsotropicElasticity{T}(G, K) end -Elastic(;E, ν) = Elastic(E, ν) # Keyword argument constructor +LinearIsotropicElasticity(;E, ν) = LinearIsotropicElasticity(E, ν) # Keyword argument constructor # Elastic material -get_cache(::Elastic) = nothing +get_cache(::LinearIsotropicElasticity) = nothing -function material_response(material::Elastic, ϵ::SymmetricTensor{2,3}, state_old, Δt::AbstractFloat; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) - ν = (3material.K - 2material.G)/(2*(3material.K+material.G)) # Calculate poissons ratio +function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt::AbstractFloat; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) + ν = (3*m.K - 2*m.G)/(2*(3*m.K+m.G)) # Calculate poissons ratio - σ = 2 * material.G*dev(ϵ) + 3 * material.K*vol(ϵ) # Calculate stress + σ = 2 * m.G*dev(ϵ) + 3 * m.K*vol(ϵ) # Calculate stress # Create stiffness matrix δ(i,j) = i == j ? 1.0 : 0.0 # helper function - Dfun(i,j,k,l) = 2.0*material.G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l)) + Dfun(i,j,k,l) = 2.0*m.G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l)) 𝔻 = SymmetricTensor{4, 3}(Dfun) # Return updated values From 9954621dd1f2b5b2906630432d612f694647d591 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 9 Jun 2021 21:14:17 +0200 Subject: [PATCH 08/64] Updated test example to reflect changes --- test/test_small_strain_plasticity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl index d77b869..5eb94b5 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_small_strain_plasticity.jl @@ -1,6 +1,6 @@ @testset "SmallStrainPlasticity" begin # constructor - m = Chaboche(elastic=Elastic(E=210.e3, ν=0.3), + m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), σ_y0=100.0, isotropic=(Iso_Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) kinematic=(Kin_AF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple From 37a7e48e7c8d7ef14e5ae55739d580a7edd4ee67 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Thu, 10 Jun 2021 13:47:36 +0200 Subject: [PATCH 09/64] Corrected some docstrings --- .../SmallStrainPlasticity.jl | 31 ++++++++++++++----- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 154c920..2a1ff3b 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -88,31 +88,48 @@ function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType end """ - material_response(m::Plastic, ϵ::SymmetricTensor{2,3}, state::ChabocheState, Δt; ) + material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, state::ChabocheState, Δt; ) -Return the stress tensor, stress tangent, the new `MaterialState` and boolean if local iterations converged for the given strain ε and previous material state `state`. +Return the stress tensor, stress tangent, the new `MaterialState` and a boolean specifying if local iterations converged. +The total strain ε and previous material state `state` are given as input, Δt has no influence as the material is rate independent. + +By specifying different laws in `m.elastic`, `m.isotropic`, and `m.kinematic`, different models can be obtained, +within the general equations specified below. Note that `m.isotropic` and `m.kinematic` are of type Tuple, in which each +element contain a hardening law of type ``AbstractiIsoHard`` and ``AbstractKinHard``, respectively. + +The stress is calculated from the elastic strains, ``\\boldsymbol{\\epsilon}_\\mathrm{e}``, obtained via the +additive decomposition, ``\\boldsymbol{\\epsilon} = \\boldsymbol{\\epsilon}_\\mathrm{e} + \\boldsymbol{\\epsilon}_\\mathrm{p}``. +The elastic law, i.e. `` is specified by `m.elastic` Von Mises yield function: ```math -\\Phi = \\sqrt{\\frac{3}{2}} \\left| \\text{dev} \\left( \\boldsymbol{\\sigma} - \\boldsymbol{\\alpha} \\right) \\right| - \\sigma_y - \\kappa +\\Phi = \\sqrt{\\frac{3}{2}} \\left| \\text{dev} \\left( \\boldsymbol{\\sigma} - \\boldsymbol{\\beta} \\right) \\right| - \\sigma_y - \\kappa ``` -Associative plastic flow: +where ``\\boldsymbol{\\beta} = \\sum_{i=1}^{N_\\mathrm{kin} \\boldsymbol{\\beta}_i`` is the total back-stress. +The evolution of ``\\boldsymbol{\\beta}_i`` is given by the kinematic hardening, specified below and using `m.kinematic` + +Associative plastic flow is used to obtain the plastic strains, ```math \\dot{\\epsilon}_{\\mathrm{p}} = \\dot{\\lambda} \\frac{\\partial \\Phi}{\\boldsymbol{\\sigma}} = \\dot{\\lambda} \\boldsymbol{\\nu} ``` -Isotropic hardening is formulated as +The plastic multiplier, ``\\lambda``, is obtained via the KKT-conditions, +```math +\\dot{\\lambda} \\geq 0, \\quad \\Phi \\leq 0, \\quad \\dot{\\lambda}\\Phi = 0 +``` + +The isotropic hardening is formulated as ```math \\kappa = \\sum_{i=1}^{N_{\\mathrm{iso}}} g_{\\mathrm{iso},i}(\\lambda) ``` -where ``g_{\\mathrm{iso},i}(\\lambda)`` is specified by structs of subtype to `AbstractIsoHard` +where ``g_{\\mathrm{iso},i}(\\lambda)`` is specified by `m.isotropic[i]` Kinematic hardening is formulated as ```math \\boldsymbol{\\beta}_i = \\dot{\\lambda} g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) ``` -where ``g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)`` is specified by structs of subtype to `AbstractKinHard` +where ``g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)`` is specified by `m.kinematic[i]` and ``i\\in[1,N_\\mathrm{kin}]``. ``` From c97f99b2f96df2d867691b31b00589308bcc6480 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Thu, 10 Jun 2021 16:00:29 +0200 Subject: [PATCH 10/64] Added abstract supertypes to match other routines --- src/SmallStrainPlasticity/Elasticity.jl | 2 +- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index f84bf88..0b6e219 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -1,5 +1,5 @@ using Tensors -abstract type AbstractElasticity{T} end +abstract type AbstractElasticity{T} <:AbstractMaterial end # Linear isotropic elasticity struct LinearIsotropicElasticity{T} <:AbstractElasticity{T} G::T # Shear modulus diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 2a1ff3b..c5b99ff 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -5,7 +5,7 @@ include("KinematicHardening.jl") # Definition of material properties -struct Chaboche{T,ElasticType,IsoType,KinType} +struct Chaboche{T,ElasticType,IsoType,KinType} <:AbstractMaterial elastic::ElasticType # Elastic definition σ_y0::T # Initial yield limit isotropic::IsoType # Tuple of isotropic hardening definitions @@ -14,13 +14,13 @@ end Chaboche(;elastic, σ_y0, isotropic, kinematic) = Chaboche(elastic, σ_y0, isotropic, kinematic) # Definition of material state -struct ChabocheState{Nkin,T,N} +struct ChabocheState{Nkin,T,N} <:AbstractMaterialState ϵₚ::SymmetricTensor{2,3,T,N} λ::T β::NTuple{Nkin, SymmetricTensor{2,3,T,N}} end -struct ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens} +struct ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens} <:AbstractResiduals λ::Tλ σ_red_dev::SymmetricTensor{2,3,Tσ,N_tens} β1::NTuple{NKin_R, SymmetricTensor{2,3,Tβ,N_tens}} From 10b152579271497694e7b731a3b5f1b86553f257 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 11 Jun 2021 14:16:01 +0200 Subject: [PATCH 11/64] Moved content of Residual.jl into SmallStrainPlasticity.jl --- src/SmallStrainPlasticity/Residual.jl | 47 ------------------ .../SmallStrainPlasticity.jl | 48 ++++++++++++++++++- 2 files changed, 47 insertions(+), 48 deletions(-) diff --git a/src/SmallStrainPlasticity/Residual.jl b/src/SmallStrainPlasticity/Residual.jl index dba60e2..e69de29 100644 --- a/src/SmallStrainPlasticity/Residual.jl +++ b/src/SmallStrainPlasticity/Residual.jl @@ -1,47 +0,0 @@ -# General residual function -function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::ChabocheState, ϵ) where{NKin_R} - Δλ = X.λ - old.λ - - σ_vm = vonMisesDev(X.σ_red_dev) - ν = X.σ_red_dev * ((3/2)*σ_vm) - - Φ = yieldCriterion(material, X.σ_red_dev, X.λ) - - σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, Δλ) - - if NKin_R > 0 - β_hat0 = σ_dev - X.σ_red_dev - sum(X.β1) - β_hat1 = X.β1 - β1 = ntuple(i->old.β[i+1] + Δλ * KinematicEvolution(material.kinematic[i+1], ν, β_hat1[i]), Val{NKin_R}()) - else - β_hat0 = σ_dev - X.σ_red_dev - end - - β0 = old.β[1] + Δλ * KinematicEvolution(material.kinematic[1], ν, β_hat0) - - if NKin_R > 0 - σ_red_dev = σ_dev - β0 - sum(β1) - R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev, ntuple(i->β1[i]-β_hat1[i], Val{NKin_R}())) - else - σ_red_dev = σ_dev - β0 - R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) - end - return R -end - -# Specialized for only one backstress (NKin_R=0) -function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState, ϵ) - - σ_vm = vonMisesDev(X.σ_red_dev) - ν = X.σ_red_dev * ((3/2)*σ_vm) - Φ = yieldCriterion(material, σ_vm, X.λ) - - σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, X.λ - old.λ) - β_hat0 = σ_dev - X.σ_red_dev - - β0 = old.β[1] + (X.λ - old.λ) * KinematicEvolution(material.kinematic[1], ν, β_hat0) - - σ_red_dev = σ_dev - β0 - - return ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) -end \ No newline at end of file diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index c5b99ff..db4c30c 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -161,7 +161,53 @@ function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_o end -include("Residual.jl") +# General residual function +function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::ChabocheState, ϵ) where{NKin_R} + Δλ = X.λ - old.λ + + σ_vm = vonMisesDev(X.σ_red_dev) + ν = X.σ_red_dev * ((3/2)*σ_vm) + + Φ = yieldCriterion(material, X.σ_red_dev, X.λ) + + σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, Δλ) + + if NKin_R > 0 + β_hat0 = σ_dev - X.σ_red_dev - sum(X.β1) + β_hat1 = X.β1 + β1 = ntuple(i->old.β[i+1] + Δλ * KinematicEvolution(material.kinematic[i+1], ν, β_hat1[i]), Val{NKin_R}()) + else + β_hat0 = σ_dev - X.σ_red_dev + end + + β0 = old.β[1] + Δλ * KinematicEvolution(material.kinematic[1], ν, β_hat0) + + if NKin_R > 0 + σ_red_dev = σ_dev - β0 - sum(β1) + R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev, ntuple(i->β1[i]-β_hat1[i], Val{NKin_R}())) + else + σ_red_dev = σ_dev - β0 + R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) + end + return R +end + +# Specialized for only one backstress (NKin_R=0) +function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState, ϵ) + + σ_vm = vonMisesDev(X.σ_red_dev) + ν = X.σ_red_dev * ((3/2)*σ_vm) + Φ = yieldCriterion(material, σ_vm, X.λ) + + σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, X.λ - old.λ) + β_hat0 = σ_dev - X.σ_red_dev + + β0 = old.β[1] + (X.λ - old.λ) * KinematicEvolution(material.kinematic[1], ν, β_hat0) + + σ_red_dev = σ_dev - β0 + + return ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) +end function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}) where {NKin,Ts,N} # σ = σ(X(ϵ), ϵ) yields From b1186066b656a4bf9e491c3cce1c1527b086fb38 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 11 Jun 2021 14:16:31 +0200 Subject: [PATCH 12/64] Deleted empty Residual.jl --- src/SmallStrainPlasticity/Residual.jl | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/SmallStrainPlasticity/Residual.jl diff --git a/src/SmallStrainPlasticity/Residual.jl b/src/SmallStrainPlasticity/Residual.jl deleted file mode 100644 index e69de29..0000000 From 6a78f2a32a82402f1b0698b1365ff6771bd06c67 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 11 Jun 2021 14:18:57 +0200 Subject: [PATCH 13/64] Removed typespec for abstract type --- src/SmallStrainPlasticity/Elasticity.jl | 4 ++-- src/SmallStrainPlasticity/IsotropicHardening.jl | 6 +++--- src/SmallStrainPlasticity/KinematicHardening.jl | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index 0b6e219..a0898a7 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -1,7 +1,7 @@ using Tensors -abstract type AbstractElasticity{T} <:AbstractMaterial end +abstract type AbstractElasticity <:AbstractMaterial end # Linear isotropic elasticity -struct LinearIsotropicElasticity{T} <:AbstractElasticity{T} +struct LinearIsotropicElasticity{T} <:AbstractElasticity G::T # Shear modulus K::T # Bulk modulus end diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl index 092cdc7..757b75c 100644 --- a/src/SmallStrainPlasticity/IsotropicHardening.jl +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -1,8 +1,8 @@ # Isotropic hardening -abstract type AbstractIsoHard{T} end +abstract type AbstractIsoHard end # Voce type of isotropic hardening (exponentially saturating) -struct Iso_Voce{T} <:AbstractIsoHard{T} +struct Iso_Voce{T} <:AbstractIsoHard Hiso::T # Initial hardening modulus κ∞::T # Saturation stress end @@ -13,7 +13,7 @@ function IsotropicHardening(param::Iso_Voce, λ::Number) end # Swift type of kinematic hardening (power law) -struct Iso_Swift{T} <:AbstractIsoHard{T} +struct Iso_Swift{T} <:AbstractIsoHard K::T λ0::T n::T diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index 0874c05..b40fa3c 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -1,10 +1,10 @@ using Tensors # Kinematic hardening -abstract type AbstractKinHard{T} end +abstract type AbstractKinHard end # Armstrong-Frederick -struct Kin_AF{T} <: AbstractKinHard{T} +struct Kin_AF{T} <: AbstractKinHard Hkin::T # Initial hardening modulus β∞::T # Saturation stress end @@ -24,7 +24,7 @@ function KinematicEvolution(param::Kin_AF, 𝛎::SecondOrderTensor, 𝛃ᵢ::Sec end # Delobelle (Combination of Armstrong-Frederick and Burlet-Cailletaud) -struct Kin_DB{T} <: AbstractKinHard{T} +struct Kin_DB{T} <: AbstractKinHard Hkin::T # Initial hardening modulus β∞::T # Saturation stress δ::T # Amount of Armstrong-Frederick hardening @@ -51,7 +51,7 @@ function KinematicEvolution(param::Kin_DB, 𝛎::SecondOrderTensor, 𝛃ᵢ::Sec end # Ohno-Wang -struct Kin_OW{T} <: AbstractKinHard{T} +struct Kin_OW{T} <: AbstractKinHard Hkin::T # Initial hardening modulus β∞::T # Saturation stress mexp::T # Ohno Wang exponent From 3b2426e69c7cc907c26365efe123fb226b281c72 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 11 Jun 2021 14:31:06 +0200 Subject: [PATCH 14/64] Updated naming of hardening types --- src/MaterialModels.jl | 4 +-- .../IsotropicHardening.jl | 14 +++++----- .../KinematicHardening.jl | 26 +++++++++---------- .../SmallStrainPlasticity.jl | 2 +- test/test_small_strain_plasticity.jl | 4 +-- 5 files changed, 24 insertions(+), 26 deletions(-) diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 550953b..c7baaae 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -89,7 +89,7 @@ export LinearElasticState, PlasticState export Chaboche export LinearIsotropicElasticity -export Iso_Voce, Iso_Swift -export Kin_AF, Kin_DB, Kin_OW +export IsotropicHardeningVoce, IsotropicHardeningSwift +export KinematicHardeningAF, KinematicHardeningDB, KinematicHardeningOW end diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl index 757b75c..ee47227 100644 --- a/src/SmallStrainPlasticity/IsotropicHardening.jl +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -1,26 +1,26 @@ # Isotropic hardening -abstract type AbstractIsoHard end +abstract type AbstractIsotropicHardening end # Voce type of isotropic hardening (exponentially saturating) -struct Iso_Voce{T} <:AbstractIsoHard +struct IsotropicHardeningVoce{T} <:AbstractIsotropicHardening Hiso::T # Initial hardening modulus κ∞::T # Saturation stress end -Iso_Voce(;Hiso, κ∞) = Iso_Voce(Hiso, κ∞) # Keyword argument constructor +IsotropicHardeningVoce(;Hiso, κ∞) = IsotropicHardeningVoce(Hiso, κ∞) # Keyword argument constructor -function IsotropicHardening(param::Iso_Voce, λ::Number) +function IsotropicHardening(param::IsotropicHardeningVoce, λ::Number) param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) end # Swift type of kinematic hardening (power law) -struct Iso_Swift{T} <:AbstractIsoHard +struct IsotropicHardeningSwift{T} <:AbstractIsotropicHardening K::T λ0::T n::T end -Iso_Swift(;K, λ0, n) = Iso_Swift(K, λ0, n) # Keyword argument constructor +IsotropicHardeningSwift(;K, λ0, n) = IsotropicHardeningSwift(K, λ0, n) # Keyword argument constructor -function IsotropicHardening(param::Iso_Swift, λ::Number) +function IsotropicHardening(param::IsotropicHardeningSwift, λ::Number) param.K * (param.λ0 + λ)^n end diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index b40fa3c..1b105d8 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -1,17 +1,15 @@ -using Tensors - # Kinematic hardening -abstract type AbstractKinHard end +abstract type AbstractKinematicHardening end # Armstrong-Frederick -struct Kin_AF{T} <: AbstractKinHard +struct KinematicHardeningAF{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress end -Kin_AF(;Hkin, β∞) = Kin_AF(Hkin, β∞) # Keyword argument constructor +KinematicHardeningAF(;Hkin, β∞) = KinematicHardeningAF(Hkin, β∞) # Keyword argument constructor """ - KinematicEvolution(param::Kin_AF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + KinematicEvolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Armstrong-Frederick kinematic hardening law @@ -19,20 +17,20 @@ Kin_AF(;Hkin, β∞) = Kin_AF(Hkin, β∞) # Keyword argument constructor g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) ``` """ -function KinematicEvolution(param::Kin_AF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) +function KinematicEvolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) end # Delobelle (Combination of Armstrong-Frederick and Burlet-Cailletaud) -struct Kin_DB{T} <: AbstractKinHard +struct KinematicHardeningDB{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress δ::T # Amount of Armstrong-Frederick hardening end -Kin_DB(;Hkin, β∞, δ) = Kin_DB(Hkin, β∞, δ) # Keyword argument constructor +KinematicHardeningDB(;Hkin, β∞, δ) = KinematicHardeningDB(Hkin, β∞, δ) # Keyword argument constructor """ - KinematicEvolution(param::Kin_DB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + KinematicEvolution(param::KinematicHardeningDB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law @@ -51,15 +49,15 @@ function KinematicEvolution(param::Kin_DB, 𝛎::SecondOrderTensor, 𝛃ᵢ::Sec end # Ohno-Wang -struct Kin_OW{T} <: AbstractKinHard +struct KinematicHardeningOW{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress mexp::T # Ohno Wang exponent end -Kin_OW(;Hkin, β∞, mexp) = Kin_OW(Hkin, β∞, mexp) # Keyword argument constructor +KinematicHardeningOW(;Hkin, β∞, mexp) = KinematicHardeningOW(Hkin, β∞, mexp) # Keyword argument constructor """ - KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} + KinematicEvolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} Kinematic hardening law according to Ohno-Wang @@ -72,7 +70,7 @@ Kin_OW(;Hkin, β∞, mexp) = Kin_OW(Hkin, β∞, mexp) # Keyword argument con ``` """ -function KinematicEvolution(param::Kin_OW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} +function KinematicEvolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} β_vm = vonMisesDev(𝛃ᵢ) if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index db4c30c..11b226d 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -95,7 +95,7 @@ The total strain ε and previous material state `state` are given as input, Δt By specifying different laws in `m.elastic`, `m.isotropic`, and `m.kinematic`, different models can be obtained, within the general equations specified below. Note that `m.isotropic` and `m.kinematic` are of type Tuple, in which each -element contain a hardening law of type ``AbstractiIsoHard`` and ``AbstractKinHard``, respectively. +element contain a hardening law of type ``AbstractIsotropicHardening`` and ``AbstractKinematicHardening``, respectively. The stress is calculated from the elastic strains, ``\\boldsymbol{\\epsilon}_\\mathrm{e}``, obtained via the additive decomposition, ``\\boldsymbol{\\epsilon} = \\boldsymbol{\\epsilon}_\\mathrm{e} + \\boldsymbol{\\epsilon}_\\mathrm{p}``. diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl index 5eb94b5..7fe1aa6 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_small_strain_plasticity.jl @@ -2,8 +2,8 @@ # constructor m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), σ_y0=100.0, - isotropic=(Iso_Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) - kinematic=(Kin_AF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple + isotropic=(IsotropicHardeningVoce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) + kinematic=(KinematicHardeningAF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple ) cache = get_cache(m) From ce5d372d34e2222152e8efe68a61f8eeaacf3f04 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 11 Jun 2021 14:35:35 +0200 Subject: [PATCH 15/64] Updated names of hardening and evolution functions to be more general and avoid constructor apperance --- src/SmallStrainPlasticity/IsotropicHardening.jl | 4 ++-- src/SmallStrainPlasticity/KinematicHardening.jl | 12 ++++++------ src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 8 ++++---- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl index ee47227..9dbc950 100644 --- a/src/SmallStrainPlasticity/IsotropicHardening.jl +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -8,7 +8,7 @@ struct IsotropicHardeningVoce{T} <:AbstractIsotropicHardening end IsotropicHardeningVoce(;Hiso, κ∞) = IsotropicHardeningVoce(Hiso, κ∞) # Keyword argument constructor -function IsotropicHardening(param::IsotropicHardeningVoce, λ::Number) +function get_hardening(param::IsotropicHardeningVoce, λ::Number) param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) end @@ -20,7 +20,7 @@ struct IsotropicHardeningSwift{T} <:AbstractIsotropicHardening end IsotropicHardeningSwift(;K, λ0, n) = IsotropicHardeningSwift(K, λ0, n) # Keyword argument constructor -function IsotropicHardening(param::IsotropicHardeningSwift, λ::Number) +function get_hardening(param::IsotropicHardeningSwift, λ::Number) param.K * (param.λ0 + λ)^n end diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index 1b105d8..d04bc1b 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -9,7 +9,7 @@ end KinematicHardeningAF(;Hkin, β∞) = KinematicHardeningAF(Hkin, β∞) # Keyword argument constructor """ - KinematicEvolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + get_evolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Armstrong-Frederick kinematic hardening law @@ -17,7 +17,7 @@ KinematicHardeningAF(;Hkin, β∞) = KinematicHardeningAF(Hkin, β∞) # Keyw g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) ``` """ -function KinematicEvolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) +function get_evolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) end @@ -30,7 +30,7 @@ end KinematicHardeningDB(;Hkin, β∞, δ) = KinematicHardeningDB(Hkin, β∞, δ) # Keyword argument constructor """ - KinematicEvolution(param::KinematicHardeningDB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + get_evolution(param::KinematicHardeningDB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law @@ -42,7 +42,7 @@ KinematicHardeningDB(;Hkin, β∞, δ) = KinematicHardeningDB(Hkin, β∞, δ) ``` """ -function KinematicEvolution(param::Kin_DB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) +function get_evolution(param::KinematicHardeningDB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) AF_Term = (param.δ/param.β∞) * 𝛃ᵢ # Armstrong Frederick term BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡𝛃ᵢ)/param.β∞)*ν # Burlet Cailletaud term return param.Hkin * ((2.0/3.0) * 𝛎 - AF_Term - BC_Term) # Complete evolution @@ -57,7 +57,7 @@ end KinematicHardeningOW(;Hkin, β∞, mexp) = KinematicHardeningOW(Hkin, β∞, mexp) # Keyword argument constructor """ - KinematicEvolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} + get_evolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} Kinematic hardening law according to Ohno-Wang @@ -70,7 +70,7 @@ KinematicHardeningOW(;Hkin, β∞, mexp) = KinematicHardeningOW(Hkin, β∞, mex ``` """ -function KinematicEvolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} +function get_evolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} β_vm = vonMisesDev(𝛃ᵢ) if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 11b226d..cda0e14 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -175,12 +175,12 @@ function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::Chaboche if NKin_R > 0 β_hat0 = σ_dev - X.σ_red_dev - sum(X.β1) β_hat1 = X.β1 - β1 = ntuple(i->old.β[i+1] + Δλ * KinematicEvolution(material.kinematic[i+1], ν, β_hat1[i]), Val{NKin_R}()) + β1 = ntuple(i->old.β[i+1] + Δλ * get_evolution(material.kinematic[i+1], ν, β_hat1[i]), Val{NKin_R}()) else β_hat0 = σ_dev - X.σ_red_dev end - β0 = old.β[1] + Δλ * KinematicEvolution(material.kinematic[1], ν, β_hat0) + β0 = old.β[1] + Δλ * get_evolution(material.kinematic[1], ν, β_hat0) if NKin_R > 0 σ_red_dev = σ_dev - β0 - sum(β1) @@ -202,7 +202,7 @@ function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, X.λ - old.λ) β_hat0 = σ_dev - X.σ_red_dev - β0 = old.β[1] + (X.λ - old.λ) * KinematicEvolution(material.kinematic[1], ν, β_hat0) + β0 = old.β[1] + (X.λ - old.λ) * get_evolution(material.kinematic[1], ν, β_hat0) σ_red_dev = σ_dev - β0 @@ -328,7 +328,7 @@ function calc_sigma_dev(material::LinearIsotropicElasticity, state_old::Chaboche end function yieldCriterion(material::Chaboche{Tp,ElType,IsoType,KinType}, σ_vm_red::Number, λ) where {Tp,ElType,IsoType<:NTuple{Niso,Any},KinType} where {Niso} - κ = sum(ntuple(i->IsotropicHardening(material.isotropic[i], λ), Val{Niso}())) + κ = sum(ntuple(i->get_hardening(material.isotropic[i], λ), Val{Niso}())) return σ_vm_red - (κ + material.σ_y0) end From 05bb83ccc675a57006045759b1fdd21c1f32ebec Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 11 Jun 2021 14:37:27 +0200 Subject: [PATCH 16/64] Removed uncessary using Tensors --- src/SmallStrainPlasticity/Elasticity.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index a0898a7..399cf6e 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -1,4 +1,3 @@ -using Tensors abstract type AbstractElasticity <:AbstractMaterial end # Linear isotropic elasticity struct LinearIsotropicElasticity{T} <:AbstractElasticity From f634ef896be7bada0f8e0f1892202face6da8052 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 10:45:26 +0200 Subject: [PATCH 17/64] Removed custom newton caching --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index cda0e14..b7a5146 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -300,11 +300,10 @@ function solve_local_problem!(cache::ChabocheCache, material::Chaboche, state_ol # Should this be centrally managed? I.e. process_options or similar? nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton)) haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method - nlsolve_options[:method] == :newton || merge!(nlsolve_options, Dict(:cache=>cache.R_X_newton)) - # Need to call newton directly to allow newton caching... - #result = NLsolve.newton(; ftol=1.e-6, cache=cache.R_X_newton) - my_newton(df, x0; xtol=0.0, ftol=1.e-8, iterations=100, store_trace=false, show_trace=false, extended_trace=false, linesearch=NLsolve.LineSearches.Static(),cache=NewtonCache(df)) = NLsolve.newton(df, x0, xtol, ftol, iterations, store_trace, show_trace, extended_trace, linesearch,cache) - result = my_newton(cache.R_X_oncediff, cache.R_X_oncediff.x_f, cache=cache.R_X_newton) + + # Solve local problem: + result = NLsolve.nlsolve(cache, cache.R_X_oncediff.x_f; nlsolve_options...) + # Is this necessary? cache.R_X_oncediff.x_f = result.zero From 6ecbe6bd55ece2e42fc3dba32d174551b82ebeee Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 10:52:41 +0200 Subject: [PATCH 18/64] Cleaned solve_local_problem to (approx) match Plasticity.jl --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index b7a5146..ed4f251 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -289,14 +289,13 @@ function initial_guess(material::Chaboche, state_old::ChabocheState{Nkin, T, N}, end -function solve_local_problem!(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{Nkin,Ts,N}, ϵ::SymmetricTensor{2,3}, options::Dict{Symbol, Any}) where {Ts,Nkin,N} +function solve_local_problem!(cache::ChabocheCache, m::Chaboche, state_old::ChabocheState{Nkin,Ts,N}, ϵ::SymmetricTensor{2,3}, options::Dict{Symbol, Any}) where {Ts,Nkin,N} - X_tensor = initial_guess(material, state_old, ϵ) - rf_tens(X_tensor_arg) = residual(X_tensor_arg, material, state_old, ϵ) - rf!(R, X) = vector_residual!(rf_tens, R, X, X_tensor) + x0 = initial_guess(m, state_old, ϵ) + rf!(r_vector, x_vector) = vector_residual!((x)->residual(x,m,state_old,ϵ), r_vector, x_vector, x0) # Using x0 as template for residual instead of material as this is related to Tensors update_cache!(cache.R_X_oncediff, rf!) - tomandel!(cache.R_X_oncediff.x_f, X_tensor) + tomandel!(cache.R_X_oncediff.x_f, x0) # Should this be centrally managed? I.e. process_options or similar? nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton)) haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method From 9ff99b5a2fd8dab639e9dd522eca6f31e9732d16 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 11:20:25 +0200 Subject: [PATCH 19/64] Commented and simplified get_plastic_output --- .../SmallStrainPlasticity.jl | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index ed4f251..35019d7 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -209,7 +209,7 @@ function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState return ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) end -function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}) where {NKin,Ts,N} +function get_plastic_output(cache::ChabocheCache, m::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}) where {NKin,Ts,N} # σ = σ(X(ϵ), ϵ) yields # dσ/dϵ = ∂σ/∂ϵ + ∂σ/∂X : dX/dϵ [1] # R = R(X(ϵ), ϵ) yields @@ -221,40 +221,37 @@ function get_plastic_output(cache::ChabocheCache, material::Chaboche, state_old: X_vec = cache.R_X_oncediff.x_f dRdX = cache.R_X_oncediff.DF - X_tensor = frommandel(ChabocheResidual{NKin-1,Ts,Ts,Ts,N}, X_vec) + X_tensor = frommandel(ChabocheResidual{NKin-1}, X_vec) # ∂σ∂X - # - Stress function - σ_X(X_arg) = get_sigma(material, state_old, X_arg, ϵ) # - Stress (vector) function: - σ_X_vec!(σv_arg, Xv_arg) = vector_residual!(σ_X, σv_arg, Xv_arg, X_tensor) + σ_X_vec!(σ_vector, x_vector) = vector_residual!((x)->get_sigma(m, state_old, x, ϵ), σ_vector, x_vector, X_tensor) σ_vec = cache.v6 # - Preallocate (should have been done beforehand, problem with Dual Tag values?) cfg = ForwardDiff.JacobianConfig(σ_X_vec!, σ_vec, X_vec, ForwardDiff.Chunk{length(X_vec)}()) # - Create DiffResult (this should be non-allocating) ∂σ∂X = cache.v6xn - diff_result = DiffResults.MutableDiffResult(σ_vec, (∂σ∂X,)) + diff_result = DiffResults.DiffResult(σ_vec, (∂σ∂X,)) # - Calculate σ and ∂σ∂X ForwardDiff.jacobian!(diff_result, σ_X_vec!, σ_vec, X_vec, cfg) σ = frommandel(SymmetricTensor{2,3}, diff_result.value) # ∂R/∂ϵ - # - Specialized residual (tensor) function: - R_ϵ(ϵ_arg) = residual(X_tensor, material, state_old, ϵ_arg) # - Specialized residual (vector) function: - R_ϵ_vec!(Rv_arg, ϵv_arg) = vector_residual!(R_ϵ, Rv_arg, ϵv_arg, ϵ) + R_ϵ_vec!(r_vector, ϵ_vector) = vector_residual!((ϵ_arg)->residual(X_tensor, m, state_old, ϵ_arg), r_vector, ϵ_vector, ϵ) ϵ_vec = cache.v6 # Use cache value (give name that makes more sense) tomandel!(ϵ_vec, ϵ) - # - Preallocate (should have been done beforehand, problem with Dual Tag values?) - R_vec = X_vec # Use as cache (ok as X_vec is not used anymore) + # - Preallocate cfg (should have been done beforehand, problem with Dual Tag values?) + R_vec = cache.R_X_oncediff.x_df cfg = ForwardDiff.JacobianConfig(R_ϵ_vec!, R_vec, ϵ_vec, ForwardDiff.Chunk{length(6)}()) - # - Calculate ∂R∂ϵ + # - Calculate ∂R∂ϵ (only jacobian required, no need for diff results) ∂R∂ϵ = cache.vnx6 ForwardDiff.jacobian!(∂R∂ϵ, R_ϵ_vec!, R_vec, ϵ_vec, cfg) # Calculate full tangent stiffness dσdϵ = dσdϵ_elastic - frommandel(SymmetricTensor{4,3}, ∂σ∂X*(dRdX\∂R∂ϵ)) + # Calculate new state variables λ = X_tensor.λ Δλ = λ-state_old.λ σ_red_dev = X_tensor.σ_red_dev From 21a520b02997d500439cc8decaf8fe916da62d7b Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 11:26:38 +0200 Subject: [PATCH 20/64] Fixed minor bug due to typo --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 35019d7..fbae0d8 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -298,7 +298,7 @@ function solve_local_problem!(cache::ChabocheCache, m::Chaboche, state_old::Chab haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method # Solve local problem: - result = NLsolve.nlsolve(cache, cache.R_X_oncediff.x_f; nlsolve_options...) + result = NLsolve.nlsolve(cache.R_X_oncediff, cache.R_X_oncediff.x_f; nlsolve_options...) # Is this necessary? cache.R_X_oncediff.x_f = result.zero From a5773c6c74ae3739d48ce47bbea9f2b37bbdeedd Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 11:36:52 +0200 Subject: [PATCH 21/64] Updated to follow naming convention of lower case function names --- .../KinematicHardening.jl | 2 +- .../SmallStrainPlasticity.jl | 20 +++++++++---------- src/SmallStrainPlasticity/UtilityFunctions.jl | 6 +++--- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index d04bc1b..5e31cd5 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -71,7 +71,7 @@ KinematicHardeningOW(;Hkin, β∞, mexp) = KinematicHardeningOW(Hkin, β∞, mex """ function get_evolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} - β_vm = vonMisesDev(𝛃ᵢ) + β_vm = vonmises_dev(𝛃ᵢ) if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ end diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index fbae0d8..726516f 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -142,7 +142,7 @@ function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_o σ_trial, dσdϵ_elastic, _, _ = material_response(material.elastic, ϵ-state_old.ϵₚ, nothing, Δt) - Φ_trial = yieldCriterion(material, σ_trial-sum(state_old.β), state_old.λ) + Φ_trial = yield_function(material, σ_trial-sum(state_old.β), state_old.λ) if Φ_trial < 0 return σ_trial, dσdϵ_elastic, state_old, true else @@ -165,10 +165,10 @@ end function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::ChabocheState, ϵ) where{NKin_R} Δλ = X.λ - old.λ - σ_vm = vonMisesDev(X.σ_red_dev) + σ_vm = vonmises_dev(X.σ_red_dev) ν = X.σ_red_dev * ((3/2)*σ_vm) - Φ = yieldCriterion(material, X.σ_red_dev, X.λ) + Φ = yield_function(material, X.σ_red_dev, X.λ) σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, Δλ) @@ -195,9 +195,9 @@ end # Specialized for only one backstress (NKin_R=0) function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState, ϵ) - σ_vm = vonMisesDev(X.σ_red_dev) + σ_vm = vonmises_dev(X.σ_red_dev) ν = X.σ_red_dev * ((3/2)*σ_vm) - Φ = yieldCriterion(material, σ_vm, X.λ) + Φ = yield_function(material, σ_vm, X.λ) σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, X.λ - old.λ) β_hat0 = σ_dev - X.σ_red_dev @@ -255,7 +255,7 @@ function get_plastic_output(cache::ChabocheCache, m::Chaboche, state_old::Chaboc λ = X_tensor.λ Δλ = λ-state_old.λ σ_red_dev = X_tensor.σ_red_dev - ν = (3.0/2.0) * σ_red_dev / vonMisesDev(σ_red_dev) + ν = (3.0/2.0) * σ_red_dev / vonmises_dev(σ_red_dev) if NKin > 1 β = ntuple(i-> i==1 ? dev(σ) - σ_red_dev - sum(X_tensor.β1) : X_tensor.β1[i-1], NKin) @@ -309,7 +309,7 @@ end function get_sigma(material::Chaboche, state_old::ChabocheState, X::ChabocheResidual, ϵ::SymmetricTensor{2,3}) Δλ = X.λ - state_old.λ σ_red_dev = X.σ_red_dev - ν = (3.0/2.0) * σ_red_dev / vonMises(σ_red_dev) + ν = (3.0/2.0) * σ_red_dev / vonmises_dev(σ_red_dev) σ = calc_sigma(material.elastic, state_old, ϵ, ν, Δλ) return σ end @@ -322,11 +322,11 @@ function calc_sigma_dev(material::LinearIsotropicElasticity, state_old::Chaboche return 2 * material.G * (dev(ϵ - state_old.ϵₚ) - Δλ*ν) end -function yieldCriterion(material::Chaboche{Tp,ElType,IsoType,KinType}, σ_vm_red::Number, λ) where {Tp,ElType,IsoType<:NTuple{Niso,Any},KinType} where {Niso} +function yield_function(material::Chaboche{Tp,ElType,IsoType,KinType}, σ_vm_red::Number, λ) where {Tp,ElType,IsoType<:NTuple{Niso,Any},KinType} where {Niso} κ = sum(ntuple(i->get_hardening(material.isotropic[i], λ), Val{Niso}())) return σ_vm_red - (κ + material.σ_y0) end -function yieldCriterion(material::Chaboche, σ_red_dev::AbstractTensor, λ) - return yieldCriterion(material, vonMises(σ_red_dev), λ) +function yield_function(material::Chaboche, σ_red_dev::AbstractTensor, λ) + return yield_function(material, vonmises(σ_red_dev), λ) end diff --git a/src/SmallStrainPlasticity/UtilityFunctions.jl b/src/SmallStrainPlasticity/UtilityFunctions.jl index 3873491..6a5cf78 100644 --- a/src/SmallStrainPlasticity/UtilityFunctions.jl +++ b/src/SmallStrainPlasticity/UtilityFunctions.jl @@ -2,12 +2,12 @@ DiffResults.DiffResult(value::MArray, derivs::Tuple{Vararg{MArray}}) = DiffResults.MutableDiffResult(value, derivs) # Generic functions, should be defined elsewhere? -function vonMises(𝛔::SymmetricTensor{2,3}) +function vonmises(𝛔::SymmetricTensor{2,3}) 𝛔_dev = dev(𝛔) - return sqrt((3.0/2.0) * (𝛔_dev ⊡ 𝛔_dev)) + return vonmises_dev(𝛔_dev) end -function vonMisesDev(𝛔_dev::SymmetricTensor{2,3}) +function vonmises_dev(𝛔_dev::SymmetricTensor{2,3}) return sqrt((3.0/2.0) * (𝛔_dev ⊡ 𝛔_dev)) end From a1a16be6f0420c5c32d0459db2980a36ff891d8f Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 11:48:17 +0200 Subject: [PATCH 22/64] Bug fix --- src/SmallStrainPlasticity/IsotropicHardening.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl index 9dbc950..d6a4464 100644 --- a/src/SmallStrainPlasticity/IsotropicHardening.jl +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -21,6 +21,6 @@ end IsotropicHardeningSwift(;K, λ0, n) = IsotropicHardeningSwift(K, λ0, n) # Keyword argument constructor function get_hardening(param::IsotropicHardeningSwift, λ::Number) - param.K * (param.λ0 + λ)^n + param.K * (param.λ0 + λ)^param.n end From 838702b258c33c6ca90402741b209962e5156093 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 14 Jun 2021 11:48:34 +0200 Subject: [PATCH 23/64] Added another example to show modularity --- test/test_small_strain_plasticity.jl | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl index 7fe1aa6..02bc860 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_small_strain_plasticity.jl @@ -1,9 +1,11 @@ @testset "SmallStrainPlasticity" begin + # Basic setup with Voce isotropic hardening and one back-stress of Armstrong-Frederick time + # constructor m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), - σ_y0=100.0, - isotropic=(IsotropicHardeningVoce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) - kinematic=(KinematicHardeningAF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple + σ_y0=100.0, + isotropic=(IsotropicHardeningVoce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) + kinematic=(KinematicHardeningAF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple ) cache = get_cache(m) @@ -25,4 +27,24 @@ @test converged + # Example with a more advanced material: + # Linear isotropic elasticity + # Two isotropic hardening laws: Voce and Swift + # Two back-stresses, one Armstrong-Frederick and one Ohno-Wang + m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), + σ_y0=100.0, + isotropic=(IsotropicHardeningVoce(Hiso=100000.0, κ∞=100.0), + IsotropicHardeningSwift(K=100.0, λ0=1.0e-2, n=0.5)), + kinematic=(KinematicHardeningAF(Hkin=40.e3, β∞=200.0), + KinematicHardeningOW(Hkin=30.e3, β∞=200.0, mexp=4.0)) + ) + + cache = get_cache(m) + state = initial_material_state(m) + Δt = 1.0 # No influence... + ϵ = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.1/100.0 : -0.3*0.1/100.0) : 0.0) + σ, ∂σ∂ε, temp_state, converged = material_response(m, ϵ, state, Δt; cache=cache) + + @test converged + end \ No newline at end of file From fd57c0a6a7baf29831e02f3fded5c1f556458303 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 16 Jun 2021 19:49:30 +0200 Subject: [PATCH 24/64] Remove print statement --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 726516f..d3a071a 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -153,7 +153,6 @@ function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_o σ, dσdϵ, state = get_plastic_output(cache, material, state_old, σ_trial, dσdϵ_elastic) else σ, dσdϵ, state = (σ_trial, dσdϵ_elastic, state_old) - println("Did not converge!") end LinearAlgebra.BLAS.set_num_threads(num_blas_threads) return σ, dσdϵ, state, converged From ec6be8baf225c87946585e951ed2919e7daf69ad Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 16 Jun 2021 19:50:46 +0200 Subject: [PATCH 25/64] Remove explicit setting of blas threads --- Project.toml | 1 - src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 3 --- 2 files changed, 4 deletions(-) diff --git a/Project.toml b/Project.toml index 01a96d8..4a0ae9d 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index d3a071a..15e5ee9 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -146,15 +146,12 @@ function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_o if Φ_trial < 0 return σ_trial, dσdϵ_elastic, state_old, true else - num_blas_threads = LinearAlgebra.BLAS.get_num_threads() # (~ 2ns) - LinearAlgebra.BLAS.set_num_threads(1) # Big performance benefit, takes ~8ns converged = solve_local_problem!(cache, material, state_old, ϵ, options) if converged σ, dσdϵ, state = get_plastic_output(cache, material, state_old, σ_trial, dσdϵ_elastic) else σ, dσdϵ, state = (σ_trial, dσdϵ_elastic, state_old) end - LinearAlgebra.BLAS.set_num_threads(num_blas_threads) return σ, dσdϵ, state, converged end From cce47cd3f49fc0255649267ef26a18e62edd56ae Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 16 Jun 2021 20:00:05 +0200 Subject: [PATCH 26/64] Updated naming of isotropic/kinematic hardening --- src/MaterialModels.jl | 5 ++-- .../IsotropicHardening.jl | 12 +++++----- .../KinematicHardening.jl | 24 +++++++++---------- test/test_small_strain_plasticity.jl | 12 +++++----- 4 files changed, 26 insertions(+), 27 deletions(-) diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index c7baaae..e42c749 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -8,7 +8,6 @@ import ForwardDiff import DiffResults using TimerOutputs using StaticArrays -using LinearAlgebra @@ -89,7 +88,7 @@ export LinearElasticState, PlasticState export Chaboche export LinearIsotropicElasticity -export IsotropicHardeningVoce, IsotropicHardeningSwift -export KinematicHardeningAF, KinematicHardeningDB, KinematicHardeningOW +export Voce, Swift +export ArmstrongFrederick, Delobelle, OhnoWang end diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/SmallStrainPlasticity/IsotropicHardening.jl index d6a4464..01af6d5 100644 --- a/src/SmallStrainPlasticity/IsotropicHardening.jl +++ b/src/SmallStrainPlasticity/IsotropicHardening.jl @@ -2,25 +2,25 @@ abstract type AbstractIsotropicHardening end # Voce type of isotropic hardening (exponentially saturating) -struct IsotropicHardeningVoce{T} <:AbstractIsotropicHardening +struct Voce{T} <:AbstractIsotropicHardening Hiso::T # Initial hardening modulus κ∞::T # Saturation stress end -IsotropicHardeningVoce(;Hiso, κ∞) = IsotropicHardeningVoce(Hiso, κ∞) # Keyword argument constructor +Voce(;Hiso, κ∞) = Voce(Hiso, κ∞) # Keyword argument constructor -function get_hardening(param::IsotropicHardeningVoce, λ::Number) +function get_hardening(param::Voce, λ::Number) param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) end # Swift type of kinematic hardening (power law) -struct IsotropicHardeningSwift{T} <:AbstractIsotropicHardening +struct Swift{T} <:AbstractIsotropicHardening K::T λ0::T n::T end -IsotropicHardeningSwift(;K, λ0, n) = IsotropicHardeningSwift(K, λ0, n) # Keyword argument constructor +Swift(;K, λ0, n) = Swift(K, λ0, n) # Keyword argument constructor -function get_hardening(param::IsotropicHardeningSwift, λ::Number) +function get_hardening(param::Swift, λ::Number) param.K * (param.λ0 + λ)^param.n end diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/SmallStrainPlasticity/KinematicHardening.jl index 5e31cd5..2c69657 100644 --- a/src/SmallStrainPlasticity/KinematicHardening.jl +++ b/src/SmallStrainPlasticity/KinematicHardening.jl @@ -2,14 +2,14 @@ abstract type AbstractKinematicHardening end # Armstrong-Frederick -struct KinematicHardeningAF{T} <: AbstractKinematicHardening +struct ArmstrongFrederick{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress end -KinematicHardeningAF(;Hkin, β∞) = KinematicHardeningAF(Hkin, β∞) # Keyword argument constructor +ArmstrongFrederick(;Hkin, β∞) = ArmstrongFrederick(Hkin, β∞) # Keyword argument constructor """ - get_evolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + get_evolution(param::ArmstrongFrederick, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Armstrong-Frederick kinematic hardening law @@ -17,20 +17,20 @@ KinematicHardeningAF(;Hkin, β∞) = KinematicHardeningAF(Hkin, β∞) # Keyw g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) ``` """ -function get_evolution(param::KinematicHardeningAF, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) +function get_evolution(param::ArmstrongFrederick, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) end # Delobelle (Combination of Armstrong-Frederick and Burlet-Cailletaud) -struct KinematicHardeningDB{T} <: AbstractKinematicHardening +struct Delobelle{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress δ::T # Amount of Armstrong-Frederick hardening end -KinematicHardeningDB(;Hkin, β∞, δ) = KinematicHardeningDB(Hkin, β∞, δ) # Keyword argument constructor +Delobelle(;Hkin, β∞, δ) = Delobelle(Hkin, β∞, δ) # Keyword argument constructor """ - get_evolution(param::KinematicHardeningDB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + get_evolution(param::Delobelle, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law @@ -42,22 +42,22 @@ KinematicHardeningDB(;Hkin, β∞, δ) = KinematicHardeningDB(Hkin, β∞, δ) ``` """ -function get_evolution(param::KinematicHardeningDB, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) +function get_evolution(param::Delobelle, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) AF_Term = (param.δ/param.β∞) * 𝛃ᵢ # Armstrong Frederick term BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡𝛃ᵢ)/param.β∞)*ν # Burlet Cailletaud term return param.Hkin * ((2.0/3.0) * 𝛎 - AF_Term - BC_Term) # Complete evolution end # Ohno-Wang -struct KinematicHardeningOW{T} <: AbstractKinematicHardening +struct OhnoWang{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress mexp::T # Ohno Wang exponent end -KinematicHardeningOW(;Hkin, β∞, mexp) = KinematicHardeningOW(Hkin, β∞, mexp) # Keyword argument constructor +OhnoWang(;Hkin, β∞, mexp) = OhnoWang(Hkin, β∞, mexp) # Keyword argument constructor """ - get_evolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} + get_evolution(param::OhnoWang{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} Kinematic hardening law according to Ohno-Wang @@ -70,7 +70,7 @@ KinematicHardeningOW(;Hkin, β∞, mexp) = KinematicHardeningOW(Hkin, β∞, mex ``` """ -function get_evolution(param::KinematicHardeningOW{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} +function get_evolution(param::OhnoWang{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} β_vm = vonmises_dev(𝛃ᵢ) if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl index 02bc860..70226ed 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_small_strain_plasticity.jl @@ -4,8 +4,8 @@ # constructor m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), σ_y0=100.0, - isotropic=(IsotropicHardeningVoce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) - kinematic=(KinematicHardeningAF(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple + isotropic=(Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) + kinematic=(ArmstrongFrederick(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple ) cache = get_cache(m) @@ -33,10 +33,10 @@ # Two back-stresses, one Armstrong-Frederick and one Ohno-Wang m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), σ_y0=100.0, - isotropic=(IsotropicHardeningVoce(Hiso=100000.0, κ∞=100.0), - IsotropicHardeningSwift(K=100.0, λ0=1.0e-2, n=0.5)), - kinematic=(KinematicHardeningAF(Hkin=40.e3, β∞=200.0), - KinematicHardeningOW(Hkin=30.e3, β∞=200.0, mexp=4.0)) + isotropic=(Voce(Hiso=100000.0, κ∞=100.0), + Swift(K=100.0, λ0=1.0e-2, n=0.5)), + kinematic=(ArmstrongFrederick(Hkin=40.e3, β∞=200.0), + OhnoWang(Hkin=30.e3, β∞=200.0, mexp=4.0)) ) cache = get_cache(m) From 4ac46ea7b60080f8bd60b24d7028382a5be8cf82 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 16 Jun 2021 20:11:24 +0200 Subject: [PATCH 27/64] Removed chunking from get_plastic_output --- .../SmallStrainPlasticity.jl | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 15e5ee9..9b9da55 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -219,30 +219,21 @@ function get_plastic_output(cache::ChabocheCache, m::Chaboche, state_old::Chaboc dRdX = cache.R_X_oncediff.DF X_tensor = frommandel(ChabocheResidual{NKin-1}, X_vec) - # ∂σ∂X - # - Stress (vector) function: + # Calculate σ and ∂σ∂X using automatic differentiation σ_X_vec!(σ_vector, x_vector) = vector_residual!((x)->get_sigma(m, state_old, x, ϵ), σ_vector, x_vector, X_tensor) - σ_vec = cache.v6 - # - Preallocate (should have been done beforehand, problem with Dual Tag values?) - cfg = ForwardDiff.JacobianConfig(σ_X_vec!, σ_vec, X_vec, ForwardDiff.Chunk{length(X_vec)}()) - # - Create DiffResult (this should be non-allocating) ∂σ∂X = cache.v6xn + σ_vec = cache.v6 diff_result = DiffResults.DiffResult(σ_vec, (∂σ∂X,)) - # - Calculate σ and ∂σ∂X - ForwardDiff.jacobian!(diff_result, σ_X_vec!, σ_vec, X_vec, cfg) + ForwardDiff.jacobian!(diff_result, σ_X_vec!, σ_vec, X_vec) σ = frommandel(SymmetricTensor{2,3}, diff_result.value) - # ∂R/∂ϵ - # - Specialized residual (vector) function: + # Calculate ∂R/∂ϵ using automatic differentiation R_ϵ_vec!(r_vector, ϵ_vector) = vector_residual!((ϵ_arg)->residual(X_tensor, m, state_old, ϵ_arg), r_vector, ϵ_vector, ϵ) ϵ_vec = cache.v6 # Use cache value (give name that makes more sense) tomandel!(ϵ_vec, ϵ) - # - Preallocate cfg (should have been done beforehand, problem with Dual Tag values?) R_vec = cache.R_X_oncediff.x_df - cfg = ForwardDiff.JacobianConfig(R_ϵ_vec!, R_vec, ϵ_vec, ForwardDiff.Chunk{length(6)}()) - # - Calculate ∂R∂ϵ (only jacobian required, no need for diff results) ∂R∂ϵ = cache.vnx6 - ForwardDiff.jacobian!(∂R∂ϵ, R_ϵ_vec!, R_vec, ϵ_vec, cfg) + ForwardDiff.jacobian!(∂R∂ϵ, R_ϵ_vec!, R_vec, ϵ_vec) # Calculate full tangent stiffness dσdϵ = dσdϵ_elastic - frommandel(SymmetricTensor{4,3}, ∂σ∂X*(dRdX\∂R∂ϵ)) From 7115cd6cb1fc94e5d4c54153ccee5759f839b7dc Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 16 Jun 2021 20:37:09 +0200 Subject: [PATCH 28/64] Removed uncessary Val types --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 9b9da55..f5f97a9 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -53,7 +53,7 @@ function Tensors.frommandel(::Type{<:ChabocheResidual{NKin_R}}, v::AbstractVecto end function initial_material_state(material::Chaboche{T}) where {T} - ChabocheState(zero(SymmetricTensor{2,3,T}), 0.0, ntuple(i->zero(SymmetricTensor{2,3,T}), Val{length(material.kinematic)}())) + ChabocheState(zero(SymmetricTensor{2,3,T}), 0.0, ntuple(i->zero(SymmetricTensor{2,3,T}), length(material.kinematic))) end # Definition of material cache @@ -171,7 +171,7 @@ function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::Chaboche if NKin_R > 0 β_hat0 = σ_dev - X.σ_red_dev - sum(X.β1) β_hat1 = X.β1 - β1 = ntuple(i->old.β[i+1] + Δλ * get_evolution(material.kinematic[i+1], ν, β_hat1[i]), Val{NKin_R}()) + β1 = ntuple(i->old.β[i+1] + Δλ * get_evolution(material.kinematic[i+1], ν, β_hat1[i]), NKin_R) else β_hat0 = σ_dev - X.σ_red_dev end @@ -180,7 +180,7 @@ function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::Chaboche if NKin_R > 0 σ_red_dev = σ_dev - β0 - sum(β1) - R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev, ntuple(i->β1[i]-β_hat1[i], Val{NKin_R}())) + R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev, ntuple(i->β1[i]-β_hat1[i], NKin_R)) else σ_red_dev = σ_dev - β0 R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) @@ -266,7 +266,7 @@ function initial_guess(material::Chaboche, state_old::ChabocheState{Nkin, T, N}, σ_red_trial = σ_trial_dev - sum(state_old.β) if Nkin > 1 - return ChabocheResidual(λ,σ_red_trial,ntuple(i->state_old.β[i], Val{Nkin-1}())) + return ChabocheResidual(λ,σ_red_trial,ntuple(i->state_old.β[i], Nkin-1)) else return ChabocheResidual(λ,σ_red_trial) end @@ -310,7 +310,7 @@ function calc_sigma_dev(material::LinearIsotropicElasticity, state_old::Chaboche end function yield_function(material::Chaboche{Tp,ElType,IsoType,KinType}, σ_vm_red::Number, λ) where {Tp,ElType,IsoType<:NTuple{Niso,Any},KinType} where {Niso} - κ = sum(ntuple(i->get_hardening(material.isotropic[i], λ), Val{Niso}())) + κ = sum(ntuple(i->get_hardening(material.isotropic[i], λ), Niso)) return σ_vm_red - (κ + material.σ_y0) end From 81058210646aa30264feafee94b1cb06981ca197 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 16 Jun 2021 22:59:33 +0200 Subject: [PATCH 29/64] Changed to non-reduced residual, equation system increase by 6 unknowns --- .../SmallStrainPlasticity.jl | 273 +++++------------- 1 file changed, 77 insertions(+), 196 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index f5f97a9..af080ee 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -20,36 +20,28 @@ struct ChabocheState{Nkin,T,N} <:AbstractMaterialState β::NTuple{Nkin, SymmetricTensor{2,3,T,N}} end -struct ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens} <:AbstractResiduals +struct ChabocheResidual{NKin,Tλ,Tσ,Tβ,N_tens} <:AbstractResiduals + σ::SymmetricTensor{2,3,Tσ,N_tens} λ::Tλ - σ_red_dev::SymmetricTensor{2,3,Tσ,N_tens} - β1::NTuple{NKin_R, SymmetricTensor{2,3,Tβ,N_tens}} -end -# Specialize for only one backstress -function ChabocheResidual(λ::Tλ,σ_red_dev::SymmetricTensor{2,3,Tσ,N_tens}) where {Tλ,Tσ,N_tens} - ChabocheResidual{0,Tλ,Tσ,Float64,N_tens}(λ,σ_red_dev,()) + β::NTuple{NKin, SymmetricTensor{2,3,Tβ,N_tens}} end -Tensors.get_base(::Type{<:ChabocheResidual{NKin_R}}) where{NKin_R} = ChabocheResidual{NKin_R} # needed for frommandel +Tensors.get_base(::Type{<:ChabocheResidual{NKin}}) where{NKin} = ChabocheResidual{NKin} # needed for frommandel -function Tensors.tomandel!(v::AbstractVector{T}, r::ChabocheResidual{NKin_R,Tλ,Tσ,Tβ,N_tens}) where {T,NKin_R,Tλ,Tσ,Tβ,N_tens} - v[1] = r.λ - tomandel!(v, r.σ_red_dev, offset=1) - for i=1:NKin_R - tomandel!(v, r.β1[i], offset=1+N_tens*i) +function Tensors.tomandel!(v::AbstractVector{T}, r::ChabocheResidual{NKin,Tλ,Tσ,Tβ,N_tens}) where {T,NKin,Tλ,Tσ,Tβ,N_tens} + tomandel!(v, r.σ) + v[7] = r.λ + for i=1:NKin + tomandel!(v, r.β[i], offset=1+N_tens*i) end return v end -function Tensors.frommandel(::Type{<:ChabocheResidual{NKin_R}}, v::AbstractVector{Tv}) where {Tv,NKin_R} - λ = v[1] - σ_red_dev = frommandel(SymmetricTensor{2,3}, v, offset=1) - if NKin_R > 0 - β1 = ntuple(i->frommandel(SymmetricTensor{2,3,Tv}, v, offset=1+6*i), NKin_R) - return ChabocheResidual(λ,σ_red_dev,β1) - else - return ChabocheResidual(λ,σ_red_dev) - end +function Tensors.frommandel(::Type{<:ChabocheResidual{NKin}}, v::AbstractVector{Tv}) where {Tv,NKin} + σ = frommandel(SymmetricTensor{2,3}, v) + λ = v[7] + β = ntuple(i->frommandel(SymmetricTensor{2,3,Tv}, v, offset=1+6*i), NKin) + return ChabocheResidual(σ,λ,β) end function initial_material_state(material::Chaboche{T}) where {T} @@ -57,34 +49,23 @@ function initial_material_state(material::Chaboche{T}) where {T} end # Definition of material cache -##= -struct ChabocheCache{T, nx, nx6, NL_TF, NL_TDF, NL_TX, NL_NC} - # General purpose - vnx6::MMatrix{nx,6,T,nx6} - v6xn::MMatrix{6,nx,T,nx6} - v6x6::MMatrix{6,6,T,36} - v6::MVector{6,T} - # For solving R(X)=0 - #X0::MVector{nx,T} +struct ChabocheCache{NL_TF, NL_TDF, NL_TX} R_X_oncediff::OnceDifferentiable{NL_TF, NL_TDF, NL_TX} - R_X_newton::NLsolve.NewtonCache{NL_NC} end function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType,IsoType,KinType} - nx = 1 + 6*length(material.kinematic) + nx = 7 + 6*length(material.kinematic) # Construct residual function and create OnceDifferentiable object state_tmp = initial_material_state(material) - σ_trial_dev = zero(SymmetricTensor{2,3,T}) + ϵ = zero(SymmetricTensor{2,3,T}) X_tensor = initial_guess(material, state_tmp, zero(SymmetricTensor{2,3,T})) - rf_tens(X_tensor) = residual(X_tensor, material, state_tmp, σ_trial_dev) + rf_tens(X_tensor) = residual(X_tensor, material, state_tmp, ϵ) rf!(R, X) = vector_residual!(rf_tens, R, X, X_tensor) X0 = MVector{nx}(zeros(T, nx)) # X0 only for shape and type information here: R_X_oncediff = OnceDifferentiable(rf!, X0, X0; autodiff = :forward) - R_X_newton = NLsolve.NewtonCache(R_X_oncediff) - return ChabocheCache(MMatrix{nx,6}(zeros(nx,6)), MMatrix{6,nx}(zeros(6,nx)), MMatrix{6,6}(zeros(6,6)), MVector{6}(zeros(6)), - R_X_oncediff, R_X_newton) + return ChabocheCache(R_X_oncediff) end """ @@ -132,188 +113,88 @@ Kinematic hardening is formulated as where ``g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)`` is specified by `m.kinematic[i]` and ``i\\in[1,N_\\mathrm{kin}]``. -``` +# Algorithmic tangent stiffness computations +σ(X(ϵ), ϵ): dσdϵ = ∂σ∂X dXdϵ + ∂σ∂ϵ +R(X(ϵ), ϵ): dRdX=0=∂R∂X dXdϵ + ∂R∂ϵ +dσdϵ = - ∂σ∂X [∂R∂X]^-1 ∂R∂ϵ + ∂σ∂ϵ + +For this specific case, + +- ∂σ∂X is identity in upper left corner, zero the rest +- ∂R∂ϵ is the elastic stiffness in upper left corner, zero the rest +- ∂σ∂ϵ is zero (the entire stress σ is only a function of x.σ) + # Keyword arguments - `cache`: Cache for the iterative solver, used by NLsolve.jl. It is strongly recommended to pre-allocate the cache for repeated calls to `material_response`. See [`get_cache`](@ref). - `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over. See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-options). By default the Newton solver will be used. """ -function material_response(material::Chaboche, ϵ::SymmetricTensor{2,3}, state_old::ChabocheState{Nkin,T,N}, Δt; cache=get_cache(material), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} +function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheState{Nkin,T,N}, Δt; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} - σ_trial, dσdϵ_elastic, _, _ = material_response(material.elastic, ϵ-state_old.ϵₚ, nothing, Δt) + σ_trial, dσdϵ_elastic, _, _ = material_response(m.elastic, ϵ-old.ϵₚ, nothing, Δt) + Φ_trial = vonmises(σ_trial-sum(old.β)) - (m.σ_y0 + sum(get_hardening.(m.isotropic, old.λ))) - Φ_trial = yield_function(material, σ_trial-sum(state_old.β), state_old.λ) if Φ_trial < 0 - return σ_trial, dσdϵ_elastic, state_old, true + return σ_trial, dσdϵ_elastic, old, true else - converged = solve_local_problem!(cache, material, state_old, ϵ, options) - if converged - σ, dσdϵ, state = get_plastic_output(cache, material, state_old, σ_trial, dσdϵ_elastic) + x0 = initial_guess(m, old, ϵ) + rf!(r_vector, x_vector) = vector_residual!((x)->residual(x,m,old,ϵ), r_vector, x_vector, x0) # Using x0 as template for residual instead of material as this is related to Tensors + update_cache!(cache.R_X_oncediff, rf!) + + tomandel!(cache.R_X_oncediff.x_f, x0) + # Should this be centrally managed? I.e. process_options or similar? + nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton)) + haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method + + # Solve local problem: + result = NLsolve.nlsolve(cache.R_X_oncediff, cache.R_X_oncediff.x_f; nlsolve_options...) + + if result.f_converged + x = frommandel(ChabocheResidual{Nkin}, result.zero) + dRdx = cache.R_X_oncediff.DF + inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx)) + dσdϵ = inv_J_σσ ⊡ dσdϵ_elastic + σ_red_dev = dev(x.σ) - sum(x.β) + ϵₚ = calculate_plastic_strain(old, σ_red_dev * ((3/2)*vonmises_dev(σ_red_dev)), x.λ-old.λ) + return x.σ, dσdϵ, ChabocheState(ϵₚ, x.λ, x.β), true else - σ, dσdϵ, state = (σ_trial, dσdϵ_elastic, state_old) + return σ_trial, dσdϵ_elastic, old, false end - return σ, dσdϵ, state, converged end end # General residual function -function residual(X::ChabocheResidual{NKin_R}, material::Chaboche, old::ChabocheState, ϵ) where{NKin_R} - Δλ = X.λ - old.λ - - σ_vm = vonmises_dev(X.σ_red_dev) - ν = X.σ_red_dev * ((3/2)*σ_vm) - - Φ = yield_function(material, X.σ_red_dev, X.λ) - - σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, Δλ) - - if NKin_R > 0 - β_hat0 = σ_dev - X.σ_red_dev - sum(X.β1) - β_hat1 = X.β1 - β1 = ntuple(i->old.β[i+1] + Δλ * get_evolution(material.kinematic[i+1], ν, β_hat1[i]), NKin_R) - else - β_hat0 = σ_dev - X.σ_red_dev - end - - β0 = old.β[1] + Δλ * get_evolution(material.kinematic[1], ν, β_hat0) - - if NKin_R > 0 - σ_red_dev = σ_dev - β0 - sum(β1) - R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev, ntuple(i->β1[i]-β_hat1[i], NKin_R)) - else - σ_red_dev = σ_dev - β0 - R = ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) - end - return R -end +function residual(x::ChabocheResidual{NKin}, m::Chaboche, old::ChabocheState, ϵ) where{NKin} + σ_red_dev = dev(x.σ) - sum(x.β) + σ_vm = vonmises_dev(σ_red_dev) + Δλ = x.λ-old.λ + ν = σ_red_dev * ((3/2)*σ_vm) # Gradient of von mises yield surface + ϵₑ = calculate_elastic_strain(old, ϵ, ν, Δλ) # Using assumption of associative plasticity + κ = sum(get_hardening.(m.isotropic, x.λ)) -# Specialized for only one backstress (NKin_R=0) -function residual(X::ChabocheResidual{0}, material::Chaboche, old::ChabocheState, ϵ) + Rσ = x.σ - calculate_sigma(m.elastic, ϵₑ) # Using the specific elastic law + Rλ = σ_vm - (m.σ_y0 + κ) + Rβ = ntuple((i) -> x.β[i] - old.β[i] - Δλ * get_evolution(m.kinematic[i], ν, x.β[i]), NKin) - σ_vm = vonmises_dev(X.σ_red_dev) - ν = X.σ_red_dev * ((3/2)*σ_vm) - Φ = yield_function(material, σ_vm, X.λ) - - σ_dev = calc_sigma_dev(material.elastic, old, ϵ, ν, X.λ - old.λ) - β_hat0 = σ_dev - X.σ_red_dev - - β0 = old.β[1] + (X.λ - old.λ) * get_evolution(material.kinematic[1], ν, β_hat0) - - σ_red_dev = σ_dev - β0 - - return ChabocheResidual(Φ, X.σ_red_dev-σ_red_dev) -end - -function get_plastic_output(cache::ChabocheCache, m::Chaboche, state_old::ChabocheState{NKin,Ts,N}, ϵ::SymmetricTensor{2,3}, dσdϵ_elastic::SymmetricTensor{4,3}) where {NKin,Ts,N} - # σ = σ(X(ϵ), ϵ) yields - # dσ/dϵ = ∂σ/∂ϵ + ∂σ/∂X : dX/dϵ [1] - # R = R(X(ϵ), ϵ) yields - # dR/dϵ = 0 = ∂R/∂ϵ + ∂R/∂X : dXdϵ [2] - # Solve [2] for dX/dϵ - # dX/dϵ = - [∂R/∂X]^-1 : ∂R∂ϵ [3] - # Insert [3] in [1] - # dσ/dϵ = ∂σ/∂ϵ - ∂σ/∂X : [∂R/∂X]^-1 : ∂R/∂ϵ [4] - - X_vec = cache.R_X_oncediff.x_f - dRdX = cache.R_X_oncediff.DF - X_tensor = frommandel(ChabocheResidual{NKin-1}, X_vec) - - # Calculate σ and ∂σ∂X using automatic differentiation - σ_X_vec!(σ_vector, x_vector) = vector_residual!((x)->get_sigma(m, state_old, x, ϵ), σ_vector, x_vector, X_tensor) - ∂σ∂X = cache.v6xn - σ_vec = cache.v6 - diff_result = DiffResults.DiffResult(σ_vec, (∂σ∂X,)) - ForwardDiff.jacobian!(diff_result, σ_X_vec!, σ_vec, X_vec) - σ = frommandel(SymmetricTensor{2,3}, diff_result.value) - - # Calculate ∂R/∂ϵ using automatic differentiation - R_ϵ_vec!(r_vector, ϵ_vector) = vector_residual!((ϵ_arg)->residual(X_tensor, m, state_old, ϵ_arg), r_vector, ϵ_vector, ϵ) - ϵ_vec = cache.v6 # Use cache value (give name that makes more sense) - tomandel!(ϵ_vec, ϵ) - R_vec = cache.R_X_oncediff.x_df - ∂R∂ϵ = cache.vnx6 - ForwardDiff.jacobian!(∂R∂ϵ, R_ϵ_vec!, R_vec, ϵ_vec) - - # Calculate full tangent stiffness - dσdϵ = dσdϵ_elastic - frommandel(SymmetricTensor{4,3}, ∂σ∂X*(dRdX\∂R∂ϵ)) - - # Calculate new state variables - λ = X_tensor.λ - Δλ = λ-state_old.λ - σ_red_dev = X_tensor.σ_red_dev - ν = (3.0/2.0) * σ_red_dev / vonmises_dev(σ_red_dev) - - if NKin > 1 - β = ntuple(i-> i==1 ? dev(σ) - σ_red_dev - sum(X_tensor.β1) : X_tensor.β1[i-1], NKin) - else - β = (dev(σ) - σ_red_dev,) - end - - ϵₚ = state_old.ϵₚ + Δλ * ν - state = ChabocheState(ϵₚ, λ, β) - - return σ, dσdϵ, state -end - -function initial_guess(material::Chaboche, state_old::ChabocheState{Nkin, T, N}, ϵ) where {T, Nkin, N} - if Nkin<1 - error("Nkin < 1 is not supported") - end - # Becomes trial by setting Δλ=0 - σ_trial_dev = calc_sigma_dev(material.elastic, state_old, ϵ, ϵ, 0.0) - λ = state_old.λ - σ_red_trial = σ_trial_dev - sum(state_old.β) - - if Nkin > 1 - return ChabocheResidual(λ,σ_red_trial,ntuple(i->state_old.β[i], Nkin-1)) - else - return ChabocheResidual(λ,σ_red_trial) - end -end - - -function solve_local_problem!(cache::ChabocheCache, m::Chaboche, state_old::ChabocheState{Nkin,Ts,N}, ϵ::SymmetricTensor{2,3}, options::Dict{Symbol, Any}) where {Ts,Nkin,N} - - x0 = initial_guess(m, state_old, ϵ) - rf!(r_vector, x_vector) = vector_residual!((x)->residual(x,m,state_old,ϵ), r_vector, x_vector, x0) # Using x0 as template for residual instead of material as this is related to Tensors - update_cache!(cache.R_X_oncediff, rf!) - - tomandel!(cache.R_X_oncediff.x_f, x0) - # Should this be centrally managed? I.e. process_options or similar? - nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton)) - haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method - - # Solve local problem: - result = NLsolve.nlsolve(cache.R_X_oncediff, cache.R_X_oncediff.x_f; nlsolve_options...) - - # Is this necessary? - cache.R_X_oncediff.x_f = result.zero - - return result.f_converged -end - -function get_sigma(material::Chaboche, state_old::ChabocheState, X::ChabocheResidual, ϵ::SymmetricTensor{2,3}) - Δλ = X.λ - state_old.λ - σ_red_dev = X.σ_red_dev - ν = (3.0/2.0) * σ_red_dev / vonmises_dev(σ_red_dev) - σ = calc_sigma(material.elastic, state_old, ϵ, ν, Δλ) - return σ + return ChabocheResidual(Rσ, Rλ, Rβ) end -function calc_sigma(material::LinearIsotropicElasticity, state_old::ChabocheState, ϵ, ν, Δλ) - return 3 * material.K*vol(ϵ) + calc_sigma_dev(material, state_old, ϵ, Δλ, ν) +function initial_guess(m::Chaboche, old::ChabocheState{Nkin}, ϵ) where {Nkin} + σ_trial = calculate_sigma(m.elastic, ϵ-old.ϵₚ) + λ = old.λ + β = ntuple(i->old.β[i], Nkin) + return ChabocheResidual(σ_trial,λ,β) end -function calc_sigma_dev(material::LinearIsotropicElasticity, state_old::ChabocheState, ϵ, ν, Δλ) - return 2 * material.G * (dev(ϵ - state_old.ϵₚ) - Δλ*ν) +function calculate_sigma(m::LinearIsotropicElasticity, ϵₑ) + return 3*m.K*vol(ϵₑ) + 2 * m.G * dev(ϵₑ) end -function yield_function(material::Chaboche{Tp,ElType,IsoType,KinType}, σ_vm_red::Number, λ) where {Tp,ElType,IsoType<:NTuple{Niso,Any},KinType} where {Niso} - κ = sum(ntuple(i->get_hardening(material.isotropic[i], λ), Niso)) - return σ_vm_red - (κ + material.σ_y0) +function calculate_elastic_strain(old::ChabocheState, ϵ, ν, Δλ) + return ϵ - calculate_plastic_strain(old, ν, Δλ) end -function yield_function(material::Chaboche, σ_red_dev::AbstractTensor, λ) - return yield_function(material, vonmises(σ_red_dev), λ) +function calculate_plastic_strain(old::ChabocheState, ν, Δλ) + return old.ϵₚ + Δλ*ν end From add46a6b2ed6e385e6a2122860956012c0d4ec78 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 18 Jun 2021 11:46:42 +0200 Subject: [PATCH 30/64] Added get initial material state to elastic --- src/SmallStrainPlasticity/Elasticity.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index 399cf6e..51f3dec 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -17,6 +17,8 @@ LinearIsotropicElasticity(;E, ν) = LinearIsotropicElasticity(E, ν) # Keywor # Elastic material get_cache(::LinearIsotropicElasticity) = nothing +initial_material_state(::LinearIsotropicElasticity) = nothing + function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt::AbstractFloat; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) ν = (3*m.K - 2*m.G)/(2*(3*m.K+m.G)) # Calculate poissons ratio From 3d83cb77cb060907f50885327bdfdbe333971c31 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 18 Jun 2021 11:47:14 +0200 Subject: [PATCH 31/64] Fixed bug in dphi/dsigma calculation --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index af080ee..f58ef37 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -155,7 +155,7 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx)) dσdϵ = inv_J_σσ ⊡ dσdϵ_elastic σ_red_dev = dev(x.σ) - sum(x.β) - ϵₚ = calculate_plastic_strain(old, σ_red_dev * ((3/2)*vonmises_dev(σ_red_dev)), x.λ-old.λ) + ϵₚ = calculate_plastic_strain(old, σ_red_dev * ((3/2)/vonmises_dev(σ_red_dev)), x.λ-old.λ) return x.σ, dσdϵ, ChabocheState(ϵₚ, x.λ, x.β), true else return σ_trial, dσdϵ_elastic, old, false @@ -169,7 +169,7 @@ function residual(x::ChabocheResidual{NKin}, m::Chaboche, old::ChabocheState, ϵ σ_red_dev = dev(x.σ) - sum(x.β) σ_vm = vonmises_dev(σ_red_dev) Δλ = x.λ-old.λ - ν = σ_red_dev * ((3/2)*σ_vm) # Gradient of von mises yield surface + ν = σ_red_dev * ((3/2)/σ_vm) # Gradient of von mises yield surface ϵₑ = calculate_elastic_strain(old, ϵ, ν, Δλ) # Using assumption of associative plasticity κ = sum(get_hardening.(m.isotropic, x.λ)) From e3bf097c095f34edd1427ed2b77a07a9da375b91 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:08:58 +0200 Subject: [PATCH 32/64] Fixed type instability at nlopt output --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index f58ef37..03942db 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -150,7 +150,7 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS result = NLsolve.nlsolve(cache.R_X_oncediff, cache.R_X_oncediff.x_f; nlsolve_options...) if result.f_converged - x = frommandel(ChabocheResidual{Nkin}, result.zero) + x = frommandel(ChabocheResidual{Nkin}, result.zero::MVector{7 + 6*Nkin, T}) dRdx = cache.R_X_oncediff.DF inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx)) dσdϵ = inv_J_σσ ⊡ dσdϵ_elastic From dcf8fafb747dff2eb8fc482760988d8be5d5ae2e Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:29:01 +0200 Subject: [PATCH 33/64] Updated LinearElastic to use new convention of full strain --- Project.toml | 5 +++-- src/LinearElastic.jl | 5 ++--- test/test_linear_elastic.jl | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 856e27a..c2ea052 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" @@ -13,7 +14,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" [compat] -julia = "1" DiffResults = "1.0" ForwardDiff = "0.10" NLsolve = "4.5" @@ -21,10 +21,11 @@ Reexport = "1.0" Rotations = "1.0" StaticArrays = "1.2" Tensors = "1.4" +julia = "1" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "JLD2"] diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index 6fdff06..c2482f1 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -53,8 +53,7 @@ Return the stress tensor, stress tangent and the new `MaterialState` for the giv \\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\Delta \\boldsymbol{\\varepsilon} . ``` """ -function material_response(m::LinearElastic, Δε::SymmetricTensor{2,3}, state::LinearElasticState{3}, Δt=nothing; cache=nothing, options=nothing) - Δσ = m.Eᵉ ⊡ Δε - σ = state.σ + Δσ +function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState{3}, Δt=nothing; cache=nothing, options=nothing) + σ = m.Eᵉ ⊡ ε return σ, m.Eᵉ, LinearElasticState(σ) end \ No newline at end of file diff --git a/test/test_linear_elastic.jl b/test/test_linear_elastic.jl index 0ef46e2..1b8a1d3 100644 --- a/test/test_linear_elastic.jl +++ b/test/test_linear_elastic.jl @@ -8,8 +8,8 @@ @test state.σ == zero(SymmetricTensor{2,3}) # constitutive driver - Δε = rand(SymmetricTensor{2,3}) - σ, ∂σ∂ε, temp_state = material_response(m, Δε, state) + ε = rand(SymmetricTensor{2,3}) + σ, ∂σ∂ε, temp_state = material_response(m, ε, state) @test σ == temp_state.σ @test ∂σ∂ε == m.Eᵉ From e5f3e63c7088373c001458d53e5d2983cc2ba75a Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:33:54 +0200 Subject: [PATCH 34/64] Moved AbstractElasticity to LinearElasticity.jl --- src/LinearElastic.jl | 10 +++++----- src/SmallStrainPlasticity/Elasticity.jl | 1 - 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index c2482f1..6541d28 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -1,4 +1,4 @@ - +abstract type AbstractElasticity <:AbstractMaterial end """ LinearElastic(E, ν) @@ -7,7 +7,7 @@ Isotropic linear elasticity. - `E::Float64`: Young's modulus - `ν::Float64`: Poisson's ratio """ -struct LinearElastic <: AbstractMaterial +struct LinearElastic <: AbstractElasticity # parameters E::Float64 # Young's modulus ν::Float64 # Poisson's ratio @@ -45,12 +45,12 @@ initial_material_state(::LinearElastic) = zero(LinearElasticState{3,Float64,6}) # constitutive drivers generally operate in 3D # (we could specialize for lower dimensions if needed for performance) """ - material_response(m::LinearElastic, Δε::SymmetricTensor{2,3}, state::LinearElasticState{3}) + material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState{3}) -Return the stress tensor, stress tangent and the new `MaterialState` for the given strain step Δε such that +Return the stress tensor, stress tangent and the new `MaterialState` for the given strain ε such that ```math -\\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\Delta \\boldsymbol{\\varepsilon} . +\\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\boldsymbol{\\varepsilon} . ``` """ function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState{3}, Δt=nothing; cache=nothing, options=nothing) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index 51f3dec..dd00993 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -1,4 +1,3 @@ -abstract type AbstractElasticity <:AbstractMaterial end # Linear isotropic elasticity struct LinearIsotropicElasticity{T} <:AbstractElasticity G::T # Shear modulus From b0aec3c5c2ee844dde765b151cc82e4565f0cb2d Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:41:02 +0200 Subject: [PATCH 35/64] Included support for a general elasticity law (removed converged from elastic material response) --- src/SmallStrainPlasticity/Elasticity.jl | 5 ++--- .../SmallStrainPlasticity.jl | 15 ++++++++++----- test/test_small_strain_plasticity.jl | 2 +- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/SmallStrainPlasticity/Elasticity.jl index dd00993..5184a0c 100644 --- a/src/SmallStrainPlasticity/Elasticity.jl +++ b/src/SmallStrainPlasticity/Elasticity.jl @@ -18,7 +18,7 @@ get_cache(::LinearIsotropicElasticity) = nothing initial_material_state(::LinearIsotropicElasticity) = nothing -function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt::AbstractFloat; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) +function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) ν = (3*m.K - 2*m.G)/(2*(3*m.K+m.G)) # Calculate poissons ratio σ = 2 * m.G*dev(ϵ) + 3 * m.K*vol(ϵ) # Calculate stress @@ -29,7 +29,6 @@ function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3 𝔻 = SymmetricTensor{4, 3}(Dfun) # Return updated values - converged = true state = state_old - return σ, 𝔻, state, converged + return σ, 𝔻, state end \ No newline at end of file diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index 03942db..f27af20 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -129,9 +129,9 @@ For this specific case, - `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over. See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-options). By default the Newton solver will be used. """ -function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheState{Nkin,T,N}, Δt; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} +function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheState{Nkin,T,N}, Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} - σ_trial, dσdϵ_elastic, _, _ = material_response(m.elastic, ϵ-old.ϵₚ, nothing, Δt) + σ_trial, dσdϵ_elastic, _ = material_response(m.elastic, ϵ-old.ϵₚ, initial_material_state(m.elastic)) Φ_trial = vonmises(σ_trial-sum(old.β)) - (m.σ_y0 + sum(get_hardening.(m.isotropic, old.λ))) if Φ_trial < 0 @@ -158,7 +158,7 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS ϵₚ = calculate_plastic_strain(old, σ_red_dev * ((3/2)/vonmises_dev(σ_red_dev)), x.λ-old.λ) return x.σ, dσdϵ, ChabocheState(ϵₚ, x.λ, x.β), true else - return σ_trial, dσdϵ_elastic, old, false + return σ_trial, dσdϵ_elastic, old end end @@ -187,10 +187,15 @@ function initial_guess(m::Chaboche, old::ChabocheState{Nkin}, ϵ) where {Nkin} return ChabocheResidual(σ_trial,λ,β) end -function calculate_sigma(m::LinearIsotropicElasticity, ϵₑ) - return 3*m.K*vol(ϵₑ) + 2 * m.G * dev(ϵₑ) +# Functions for different elastic laws (elasticity = path independent, no state required) +function calculate_sigma(m::AbstractElasticity, ϵₑ) + σ, _, _ = material_response(m, ϵₑ, initial_material_state(m)) + return σ end +# Specialized function (faster) +calculate_sigma(m::LinearIsotropicElasticity, ϵₑ) = 3*m.K*vol(ϵₑ) + 2 * m.G * dev(ϵₑ) + function calculate_elastic_strain(old::ChabocheState, ϵ, ν, Δλ) return ϵ - calculate_plastic_strain(old, ν, Δλ) end diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl index 70226ed..ed79a7a 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_small_strain_plasticity.jl @@ -31,7 +31,7 @@ # Linear isotropic elasticity # Two isotropic hardening laws: Voce and Swift # Two back-stresses, one Armstrong-Frederick and one Ohno-Wang - m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), + m = Chaboche(elastic=LinearElastic(E=210.e3, ν=0.3), σ_y0=100.0, isotropic=(Voce(Hiso=100000.0, κ∞=100.0), Swift(K=100.0, λ0=1.0e-2, n=0.5)), From 83c6f4c86dfad0fd6fef143cd6f5cfe20764f69a Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:44:05 +0200 Subject: [PATCH 36/64] Removed converged bool to follow standard material_response interface --- src/SmallStrainPlasticity/SmallStrainPlasticity.jl | 6 +++--- test/test_small_strain_plasticity.jl | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl index f27af20..ade97c2 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/SmallStrainPlasticity/SmallStrainPlasticity.jl @@ -135,7 +135,7 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS Φ_trial = vonmises(σ_trial-sum(old.β)) - (m.σ_y0 + sum(get_hardening.(m.isotropic, old.λ))) if Φ_trial < 0 - return σ_trial, dσdϵ_elastic, old, true + return σ_trial, dσdϵ_elastic, old else x0 = initial_guess(m, old, ϵ) rf!(r_vector, x_vector) = vector_residual!((x)->residual(x,m,old,ϵ), r_vector, x_vector, x0) # Using x0 as template for residual instead of material as this is related to Tensors @@ -156,9 +156,9 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS dσdϵ = inv_J_σσ ⊡ dσdϵ_elastic σ_red_dev = dev(x.σ) - sum(x.β) ϵₚ = calculate_plastic_strain(old, σ_red_dev * ((3/2)/vonmises_dev(σ_red_dev)), x.λ-old.λ) - return x.σ, dσdϵ, ChabocheState(ϵₚ, x.λ, x.β), true + return x.σ, dσdϵ, ChabocheState(ϵₚ, x.λ, x.β) else - return σ_trial, dσdϵ_elastic, old + error("Material model not converged. Could not find material state.") end end diff --git a/test/test_small_strain_plasticity.jl b/test/test_small_strain_plasticity.jl index ed79a7a..cd3fd73 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_small_strain_plasticity.jl @@ -23,9 +23,9 @@ Δt = 1.0 # No influence... - σ, ∂σ∂ε, temp_state, converged = material_response(m, ϵ, state, Δt; cache=cache) + σ, ∂σ∂ε, temp_state = material_response(m, ϵ, state, Δt; cache=cache) - @test converged + @test true # Check that it ran (throws error if not converged) # Example with a more advanced material: # Linear isotropic elasticity @@ -43,8 +43,8 @@ state = initial_material_state(m) Δt = 1.0 # No influence... ϵ = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.1/100.0 : -0.3*0.1/100.0) : 0.0) - σ, ∂σ∂ε, temp_state, converged = material_response(m, ϵ, state, Δt; cache=cache) + σ, ∂σ∂ε, temp_state = material_response(m, ϵ, state, Δt; cache=cache) - @test converged + @test true # Check that it ran (throws error if not converged) end \ No newline at end of file From 430bc111a07d1f31fa88d65207295b81b54240b9 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:54:00 +0200 Subject: [PATCH 37/64] Updated naming from SmallStrainPlasticity to VonMisesPlasticity --- src/MaterialModels.jl | 4 +- .../Elasticity.jl | 0 .../IsotropicHardening.jl | 0 .../KinematicHardening.jl | 0 .../UtilityFunctions.jl | 0 .../VonMisesPlasticity.jl} | 46 +++++++++---------- ...sticity.jl => test_vonmises_plasticity.jl} | 24 +++++----- 7 files changed, 37 insertions(+), 37 deletions(-) rename src/{SmallStrainPlasticity => VonMisesPlasticity}/Elasticity.jl (100%) rename src/{SmallStrainPlasticity => VonMisesPlasticity}/IsotropicHardening.jl (100%) rename src/{SmallStrainPlasticity => VonMisesPlasticity}/KinematicHardening.jl (100%) rename src/{SmallStrainPlasticity => VonMisesPlasticity}/UtilityFunctions.jl (100%) rename src/{SmallStrainPlasticity/SmallStrainPlasticity.jl => VonMisesPlasticity/VonMisesPlasticity.jl} (77%) rename test/{test_small_strain_plasticity.jl => test_vonmises_plasticity.jl} (59%) diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 99b8e22..2ffac1a 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -70,7 +70,7 @@ include("Plastic.jl") include("CrystalViscoPlastic/slipsystems.jl") include("CrystalViscoPlastic/CrystalViscoPlastic.jl") include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl") -include("SmallStrainPlasticity/SmallStrainPlasticity.jl") +include("VonMisesPlasticity/VonMisesPlasticity.jl") include("nonlinear_solver.jl") include("wrappers.jl") @@ -83,7 +83,7 @@ export LinearElastic, Plastic export LinearElasticState, PlasticState export OneD, UniaxialStrain, UniaxialStress, PlaneStrain, PlaneStress -export Chaboche +export VonMisesPlasticity export LinearIsotropicElasticity export Voce, Swift export ArmstrongFrederick, Delobelle, OhnoWang diff --git a/src/SmallStrainPlasticity/Elasticity.jl b/src/VonMisesPlasticity/Elasticity.jl similarity index 100% rename from src/SmallStrainPlasticity/Elasticity.jl rename to src/VonMisesPlasticity/Elasticity.jl diff --git a/src/SmallStrainPlasticity/IsotropicHardening.jl b/src/VonMisesPlasticity/IsotropicHardening.jl similarity index 100% rename from src/SmallStrainPlasticity/IsotropicHardening.jl rename to src/VonMisesPlasticity/IsotropicHardening.jl diff --git a/src/SmallStrainPlasticity/KinematicHardening.jl b/src/VonMisesPlasticity/KinematicHardening.jl similarity index 100% rename from src/SmallStrainPlasticity/KinematicHardening.jl rename to src/VonMisesPlasticity/KinematicHardening.jl diff --git a/src/SmallStrainPlasticity/UtilityFunctions.jl b/src/VonMisesPlasticity/UtilityFunctions.jl similarity index 100% rename from src/SmallStrainPlasticity/UtilityFunctions.jl rename to src/VonMisesPlasticity/UtilityFunctions.jl diff --git a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl similarity index 77% rename from src/SmallStrainPlasticity/SmallStrainPlasticity.jl rename to src/VonMisesPlasticity/VonMisesPlasticity.jl index ade97c2..204d0b5 100644 --- a/src/SmallStrainPlasticity/SmallStrainPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -5,30 +5,30 @@ include("KinematicHardening.jl") # Definition of material properties -struct Chaboche{T,ElasticType,IsoType,KinType} <:AbstractMaterial +struct VonMisesPlasticity{T,ElasticType,IsoType,KinType} <:AbstractMaterial elastic::ElasticType # Elastic definition σ_y0::T # Initial yield limit isotropic::IsoType # Tuple of isotropic hardening definitions kinematic::KinType # Tuple of kinematic hardening definitions end -Chaboche(;elastic, σ_y0, isotropic, kinematic) = Chaboche(elastic, σ_y0, isotropic, kinematic) +VonMisesPlasticity(;elastic, σ_y0, isotropic, kinematic) = VonMisesPlasticity(elastic, σ_y0, isotropic, kinematic) # Definition of material state -struct ChabocheState{Nkin,T,N} <:AbstractMaterialState +struct VonMisesPlasticityState{Nkin,T,N} <:AbstractMaterialState ϵₚ::SymmetricTensor{2,3,T,N} λ::T β::NTuple{Nkin, SymmetricTensor{2,3,T,N}} end -struct ChabocheResidual{NKin,Tλ,Tσ,Tβ,N_tens} <:AbstractResiduals +struct VonMisesPlasticityResidual{NKin,Tλ,Tσ,Tβ,N_tens} <:AbstractResiduals σ::SymmetricTensor{2,3,Tσ,N_tens} λ::Tλ β::NTuple{NKin, SymmetricTensor{2,3,Tβ,N_tens}} end -Tensors.get_base(::Type{<:ChabocheResidual{NKin}}) where{NKin} = ChabocheResidual{NKin} # needed for frommandel +Tensors.get_base(::Type{<:VonMisesPlasticityResidual{NKin}}) where{NKin} = VonMisesPlasticityResidual{NKin} # needed for frommandel -function Tensors.tomandel!(v::AbstractVector{T}, r::ChabocheResidual{NKin,Tλ,Tσ,Tβ,N_tens}) where {T,NKin,Tλ,Tσ,Tβ,N_tens} +function Tensors.tomandel!(v::AbstractVector{T}, r::VonMisesPlasticityResidual{NKin,Tλ,Tσ,Tβ,N_tens}) where {T,NKin,Tλ,Tσ,Tβ,N_tens} tomandel!(v, r.σ) v[7] = r.λ for i=1:NKin @@ -37,23 +37,23 @@ function Tensors.tomandel!(v::AbstractVector{T}, r::ChabocheResidual{NKin,Tλ,T return v end -function Tensors.frommandel(::Type{<:ChabocheResidual{NKin}}, v::AbstractVector{Tv}) where {Tv,NKin} +function Tensors.frommandel(::Type{<:VonMisesPlasticityResidual{NKin}}, v::AbstractVector{Tv}) where {Tv,NKin} σ = frommandel(SymmetricTensor{2,3}, v) λ = v[7] β = ntuple(i->frommandel(SymmetricTensor{2,3,Tv}, v, offset=1+6*i), NKin) - return ChabocheResidual(σ,λ,β) + return VonMisesPlasticityResidual(σ,λ,β) end -function initial_material_state(material::Chaboche{T}) where {T} - ChabocheState(zero(SymmetricTensor{2,3,T}), 0.0, ntuple(i->zero(SymmetricTensor{2,3,T}), length(material.kinematic))) +function initial_material_state(material::VonMisesPlasticity{T}) where {T} + VonMisesPlasticityState(zero(SymmetricTensor{2,3,T}), 0.0, ntuple(i->zero(SymmetricTensor{2,3,T}), length(material.kinematic))) end # Definition of material cache -struct ChabocheCache{NL_TF, NL_TDF, NL_TX} +struct VonMisesPlasticityCache{NL_TF, NL_TDF, NL_TX} R_X_oncediff::OnceDifferentiable{NL_TF, NL_TDF, NL_TX} end -function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType,IsoType,KinType} +function get_cache(material::VonMisesPlasticity{T,ElType,IsoType,KinType}) where {T,ElType,IsoType,KinType} nx = 7 + 6*length(material.kinematic) # Construct residual function and create OnceDifferentiable object @@ -65,11 +65,11 @@ function get_cache(material::Chaboche{T,ElType,IsoType,KinType}) where {T,ElType X0 = MVector{nx}(zeros(T, nx)) # X0 only for shape and type information here: R_X_oncediff = OnceDifferentiable(rf!, X0, X0; autodiff = :forward) - return ChabocheCache(R_X_oncediff) + return VonMisesPlasticityCache(R_X_oncediff) end """ - material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, state::ChabocheState, Δt; ) + material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, state::VonMisesPlasticityState, Δt; ) Return the stress tensor, stress tangent, the new `MaterialState` and a boolean specifying if local iterations converged. The total strain ε and previous material state `state` are given as input, Δt has no influence as the material is rate independent. @@ -129,7 +129,7 @@ For this specific case, - `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over. See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-options). By default the Newton solver will be used. """ -function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheState{Nkin,T,N}, Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} +function material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, old::VonMisesPlasticityState{Nkin,T,N}, Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} σ_trial, dσdϵ_elastic, _ = material_response(m.elastic, ϵ-old.ϵₚ, initial_material_state(m.elastic)) Φ_trial = vonmises(σ_trial-sum(old.β)) - (m.σ_y0 + sum(get_hardening.(m.isotropic, old.λ))) @@ -150,13 +150,13 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS result = NLsolve.nlsolve(cache.R_X_oncediff, cache.R_X_oncediff.x_f; nlsolve_options...) if result.f_converged - x = frommandel(ChabocheResidual{Nkin}, result.zero::MVector{7 + 6*Nkin, T}) + x = frommandel(VonMisesPlasticityResidual{Nkin}, result.zero::MVector{7 + 6*Nkin, T}) dRdx = cache.R_X_oncediff.DF inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx)) dσdϵ = inv_J_σσ ⊡ dσdϵ_elastic σ_red_dev = dev(x.σ) - sum(x.β) ϵₚ = calculate_plastic_strain(old, σ_red_dev * ((3/2)/vonmises_dev(σ_red_dev)), x.λ-old.λ) - return x.σ, dσdϵ, ChabocheState(ϵₚ, x.λ, x.β) + return x.σ, dσdϵ, VonMisesPlasticityState(ϵₚ, x.λ, x.β) else error("Material model not converged. Could not find material state.") end @@ -165,7 +165,7 @@ function material_response(m::Chaboche, ϵ::SymmetricTensor{2,3}, old::ChabocheS end # General residual function -function residual(x::ChabocheResidual{NKin}, m::Chaboche, old::ChabocheState, ϵ) where{NKin} +function residual(x::VonMisesPlasticityResidual{NKin}, m::VonMisesPlasticity, old::VonMisesPlasticityState, ϵ) where{NKin} σ_red_dev = dev(x.σ) - sum(x.β) σ_vm = vonmises_dev(σ_red_dev) Δλ = x.λ-old.λ @@ -177,14 +177,14 @@ function residual(x::ChabocheResidual{NKin}, m::Chaboche, old::ChabocheState, ϵ Rλ = σ_vm - (m.σ_y0 + κ) Rβ = ntuple((i) -> x.β[i] - old.β[i] - Δλ * get_evolution(m.kinematic[i], ν, x.β[i]), NKin) - return ChabocheResidual(Rσ, Rλ, Rβ) + return VonMisesPlasticityResidual(Rσ, Rλ, Rβ) end -function initial_guess(m::Chaboche, old::ChabocheState{Nkin}, ϵ) where {Nkin} +function initial_guess(m::VonMisesPlasticity, old::VonMisesPlasticityState{Nkin}, ϵ) where {Nkin} σ_trial = calculate_sigma(m.elastic, ϵ-old.ϵₚ) λ = old.λ β = ntuple(i->old.β[i], Nkin) - return ChabocheResidual(σ_trial,λ,β) + return VonMisesPlasticityResidual(σ_trial,λ,β) end # Functions for different elastic laws (elasticity = path independent, no state required) @@ -196,10 +196,10 @@ end # Specialized function (faster) calculate_sigma(m::LinearIsotropicElasticity, ϵₑ) = 3*m.K*vol(ϵₑ) + 2 * m.G * dev(ϵₑ) -function calculate_elastic_strain(old::ChabocheState, ϵ, ν, Δλ) +function calculate_elastic_strain(old::VonMisesPlasticityState, ϵ, ν, Δλ) return ϵ - calculate_plastic_strain(old, ν, Δλ) end -function calculate_plastic_strain(old::ChabocheState, ν, Δλ) +function calculate_plastic_strain(old::VonMisesPlasticityState, ν, Δλ) return old.ϵₚ + Δλ*ν end diff --git a/test/test_small_strain_plasticity.jl b/test/test_vonmises_plasticity.jl similarity index 59% rename from test/test_small_strain_plasticity.jl rename to test/test_vonmises_plasticity.jl index cd3fd73..377182e 100644 --- a/test/test_small_strain_plasticity.jl +++ b/test/test_vonmises_plasticity.jl @@ -1,11 +1,11 @@ -@testset "SmallStrainPlasticity" begin +@testset "VonMisesPlasticity" begin # Basic setup with Voce isotropic hardening and one back-stress of Armstrong-Frederick time # constructor - m = Chaboche(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), - σ_y0=100.0, - isotropic=(Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) - kinematic=(ArmstrongFrederick(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple + m = VonMisesPlasticity(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), + σ_y0=100.0, + isotropic=(Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) + kinematic=(ArmstrongFrederick(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple ) cache = get_cache(m) @@ -19,7 +19,7 @@ λ_old = 0.0015338757291717328 β_old = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 133.268 : -66.6339) : 0.0) - state_old = MaterialModels.ChabocheState(ϵₚ_old,λ_old, (β_old,)) + state_old = MaterialModels.VonMisesPlasticityState(ϵₚ_old,λ_old, (β_old,)) Δt = 1.0 # No influence... @@ -31,12 +31,12 @@ # Linear isotropic elasticity # Two isotropic hardening laws: Voce and Swift # Two back-stresses, one Armstrong-Frederick and one Ohno-Wang - m = Chaboche(elastic=LinearElastic(E=210.e3, ν=0.3), - σ_y0=100.0, - isotropic=(Voce(Hiso=100000.0, κ∞=100.0), - Swift(K=100.0, λ0=1.0e-2, n=0.5)), - kinematic=(ArmstrongFrederick(Hkin=40.e3, β∞=200.0), - OhnoWang(Hkin=30.e3, β∞=200.0, mexp=4.0)) + m = VonMisesPlasticity(elastic=LinearElastic(E=210.e3, ν=0.3), + σ_y0=100.0, + isotropic=(Voce(Hiso=100000.0, κ∞=100.0), + Swift(K=100.0, λ0=1.0e-2, n=0.5)), + kinematic=(ArmstrongFrederick(Hkin=40.e3, β∞=200.0), + OhnoWang(Hkin=30.e3, β∞=200.0, mexp=4.0)) ) cache = get_cache(m) From 6bc710a7e22f274c7d1848eb3273072acb8de3a3 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 22:59:36 +0200 Subject: [PATCH 38/64] Fixed wrong include file --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index fc5610b..ca7da1b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ using JLD2 include("test_utils.jl") include("test_linear_elastic.jl") include("test_plastic.jl") -include("test_small_strain_plasticity.jl") +include("test_vonmises_plasticity.jl") include("test_crystal_visco_plastic.jl") include("test_crystal_visco_plastic_red.jl") include("test_wrappers.jl") From 03efb1fc8e02002fc684be08e9b2ea832884ec56 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Mon, 21 Jun 2021 23:00:23 +0200 Subject: [PATCH 39/64] Moved utility function to general routine --- src/MaterialModels.jl | 1 + src/VonMisesPlasticity/VonMisesPlasticity.jl | 1 - .../UtilityFunctions.jl => utility_functions.jl} | 0 3 files changed, 1 insertion(+), 1 deletion(-) rename src/{VonMisesPlasticity/UtilityFunctions.jl => utility_functions.jl} (100%) diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 2ffac1a..1641398 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -72,6 +72,7 @@ include("CrystalViscoPlastic/CrystalViscoPlastic.jl") include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl") include("VonMisesPlasticity/VonMisesPlasticity.jl") +include("utility_functions.jl") include("nonlinear_solver.jl") include("wrappers.jl") diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index 204d0b5..a663a44 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -1,4 +1,3 @@ -include("UtilityFunctions.jl") include("Elasticity.jl") include("IsotropicHardening.jl") include("KinematicHardening.jl") diff --git a/src/VonMisesPlasticity/UtilityFunctions.jl b/src/utility_functions.jl similarity index 100% rename from src/VonMisesPlasticity/UtilityFunctions.jl rename to src/utility_functions.jl From 00fde5ca3dd053c78267db8c47a649f53963d8e2 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 10:20:49 +0200 Subject: [PATCH 40/64] Updated documentation --- src/VonMisesPlasticity/IsotropicHardening.jl | 24 ++++++++++++++++++++ src/VonMisesPlasticity/KinematicHardening.jl | 5 ++-- src/VonMisesPlasticity/VonMisesPlasticity.jl | 13 +---------- 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/src/VonMisesPlasticity/IsotropicHardening.jl b/src/VonMisesPlasticity/IsotropicHardening.jl index 01af6d5..5faaa7a 100644 --- a/src/VonMisesPlasticity/IsotropicHardening.jl +++ b/src/VonMisesPlasticity/IsotropicHardening.jl @@ -8,6 +8,20 @@ struct Voce{T} <:AbstractIsotropicHardening end Voce(;Hiso, κ∞) = Voce(Hiso, κ∞) # Keyword argument constructor +""" + get_hardening(param::Voce, λ::Number) + + Exponentially saturating isotropic hardening + + ```math + \\kappa_i = g_{\\mathrm{iso},i}(\\lambda) = \\kappa_\\infty \\left[1 - \\mathrm{exp}\\left(\\frac{H_\\mathrm{iso}}{\\kappa_\\infty} \\lambda \\right)\\right] + ``` + or alternatively as differential equations + ```math + \\dot{\\kappa_i} = \\dot{\\lambda} H_\\mathrm{iso} \\left[1 - \\frac{\\kappa_i}{\\kappa_\\infty}\\right] + ``` + +""" function get_hardening(param::Voce, λ::Number) param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) end @@ -20,6 +34,16 @@ struct Swift{T} <:AbstractIsotropicHardening end Swift(;K, λ0, n) = Swift(K, λ0, n) # Keyword argument constructor +""" + get_hardening(param::Voce, λ::Number) + + Isotropic hardening by the Swift power law + + ```math + \\kappa_i = g_{\\mathrm{iso},i}(\\lambda) = K \\left[\\lambda_0 + \\lambda \\right]^n + ``` + +""" function get_hardening(param::Swift, λ::Number) param.K * (param.λ0 + λ)^param.n end diff --git a/src/VonMisesPlasticity/KinematicHardening.jl b/src/VonMisesPlasticity/KinematicHardening.jl index 2c69657..36fb4a5 100644 --- a/src/VonMisesPlasticity/KinematicHardening.jl +++ b/src/VonMisesPlasticity/KinematicHardening.jl @@ -11,7 +11,7 @@ ArmstrongFrederick(;Hkin, β∞) = ArmstrongFrederick(Hkin, β∞) # Keyword """ get_evolution(param::ArmstrongFrederick, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) - Armstrong-Frederick kinematic hardening law + Armstrong-Frederick kinematic hardening law (doi: 10.1179/096034007X207589) ```math g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) @@ -33,6 +33,7 @@ Delobelle(;Hkin, β∞, δ) = Delobelle(Hkin, β∞, δ) # Keyword argument c get_evolution(param::Delobelle, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law + (doi: 10.1016/S0749-6419(95)00001-1) ```math g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin \\left[\\frac{2}{3}\\boldsymbol{\\nu} @@ -59,7 +60,7 @@ OhnoWang(;Hkin, β∞, mexp) = OhnoWang(Hkin, β∞, mexp) # Keyword argument """ get_evolution(param::OhnoWang{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} - Kinematic hardening law according to Ohno-Wang + Kinematic hardening law according to Ohno-Wang (doi: 10.1016/0749-6419(93)90042-O) ```math g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin \\left[\\frac{2}{3}\\boldsymbol{\\nu} diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index a663a44..f8c2257 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -79,7 +79,7 @@ element contain a hardening law of type ``AbstractIsotropicHardening`` and ``Abs The stress is calculated from the elastic strains, ``\\boldsymbol{\\epsilon}_\\mathrm{e}``, obtained via the additive decomposition, ``\\boldsymbol{\\epsilon} = \\boldsymbol{\\epsilon}_\\mathrm{e} + \\boldsymbol{\\epsilon}_\\mathrm{p}``. -The elastic law, i.e. `` is specified by `m.elastic` +The elastic law is specified by `m.elastic` and is evaluated by giving it the elastic strain. Von Mises yield function: ```math @@ -112,17 +112,6 @@ Kinematic hardening is formulated as where ``g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)`` is specified by `m.kinematic[i]` and ``i\\in[1,N_\\mathrm{kin}]``. -# Algorithmic tangent stiffness computations -σ(X(ϵ), ϵ): dσdϵ = ∂σ∂X dXdϵ + ∂σ∂ϵ -R(X(ϵ), ϵ): dRdX=0=∂R∂X dXdϵ + ∂R∂ϵ -dσdϵ = - ∂σ∂X [∂R∂X]^-1 ∂R∂ϵ + ∂σ∂ϵ - -For this specific case, - -- ∂σ∂X is identity in upper left corner, zero the rest -- ∂R∂ϵ is the elastic stiffness in upper left corner, zero the rest -- ∂σ∂ϵ is zero (the entire stress σ is only a function of x.σ) - # Keyword arguments - `cache`: Cache for the iterative solver, used by NLsolve.jl. It is strongly recommended to pre-allocate the cache for repeated calls to `material_response`. See [`get_cache`](@ref). - `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over. From a62dce18941b9740c21b19093d7585cf5ed53239 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 10:23:06 +0200 Subject: [PATCH 41/64] Added description of AbstractElasticity --- src/LinearElastic.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index 6541d28..406c9c8 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -1,4 +1,8 @@ +# Definition of materials of type AbstractElasticity: +# It does not rely on previous states to calculate the stress and stiffness +# This allows it to be used inside other materials abstract type AbstractElasticity <:AbstractMaterial end + """ LinearElastic(E, ν) From b4e170b731acafa7dd2efe1218515da355f807e4 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 10:34:22 +0200 Subject: [PATCH 42/64] Added documentation to LinearIsotropicElasticity --- src/VonMisesPlasticity/Elasticity.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/VonMisesPlasticity/Elasticity.jl b/src/VonMisesPlasticity/Elasticity.jl index 5184a0c..8ddd9dc 100644 --- a/src/VonMisesPlasticity/Elasticity.jl +++ b/src/VonMisesPlasticity/Elasticity.jl @@ -18,7 +18,25 @@ get_cache(::LinearIsotropicElasticity) = nothing initial_material_state(::LinearIsotropicElasticity) = nothing -function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) +""" + material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt=nothing; ) + +Return the stress tensor, stress tangent and the new `MaterialState` +The total (elastic) strain ϵ is given as input (state_old and Δt can be supplied, but have no influence) +The stress is calculated as + +```math +\\boldsymbol{\\sigma} = 2G \\boldsymbol{\\epsilon}^\\mathrm{dev} + 3K \\boldsymbol{\\epsilon}^\\mathrm{vol} +``` +where ``\\boldsymbol{\\epsilon}^\\mathrm{vol} = \\boldsymbol{I}:\\boldsymbol{\\epsilon} \\boldsymbol{I}/3.0`` +is the volumetric strain and ``\\boldsymbol{\\epsilon}^\\mathrm{dev} = \\boldsymbol{\\epsilon} - \\boldsymbol{\\epsilon}^\\mathrm{vol}`` +is the deviatoric strain. + +# Keyword arguments +- `cache`: Not used +- `options::Dict{Symbol, Any}`: Not used +""" +function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old=initial_material_state(m), Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) ν = (3*m.K - 2*m.G)/(2*(3*m.K+m.G)) # Calculate poissons ratio σ = 2 * m.G*dev(ϵ) + 3 * m.K*vol(ϵ) # Calculate stress From 210fd6621819cfb8726af4dd7de50dc445f7178a Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 12:15:22 +0200 Subject: [PATCH 43/64] Removed old elastic routine --- src/LinearElastic.jl | 6 ++- src/VonMisesPlasticity/Elasticity.jl | 52 -------------------- src/VonMisesPlasticity/VonMisesPlasticity.jl | 10 ---- test/test_vonmises_plasticity.jl | 2 +- 4 files changed, 5 insertions(+), 65 deletions(-) delete mode 100644 src/VonMisesPlasticity/Elasticity.jl diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index 406c9c8..33438d5 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -58,6 +58,8 @@ Return the stress tensor, stress tangent and the new `MaterialState` for the giv ``` """ function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState{3}, Δt=nothing; cache=nothing, options=nothing) - σ = m.Eᵉ ⊡ ε + σ = calculate_sigma(m, ε) return σ, m.Eᵉ, LinearElasticState(σ) -end \ No newline at end of file +end + +calculate_sigma(m::LinearElastic, ε) = m.Eᵉ ⊡ ε # Useful when calling from other materials \ No newline at end of file diff --git a/src/VonMisesPlasticity/Elasticity.jl b/src/VonMisesPlasticity/Elasticity.jl deleted file mode 100644 index 8ddd9dc..0000000 --- a/src/VonMisesPlasticity/Elasticity.jl +++ /dev/null @@ -1,52 +0,0 @@ -# Linear isotropic elasticity -struct LinearIsotropicElasticity{T} <:AbstractElasticity - G::T # Shear modulus - K::T # Bulk modulus -end -# Overload initialization method to use more common input parameters -# E: Young's modulus, ν: Poissons ratio -function LinearIsotropicElasticity(E::Number, ν::Number) - T = promote_type(typeof(E), typeof(ν)) - G = E / 2(1 + ν) - K = E / 3(1 - 2ν) - return LinearIsotropicElasticity{T}(G, K) -end -LinearIsotropicElasticity(;E, ν) = LinearIsotropicElasticity(E, ν) # Keyword argument constructor - -# Elastic material -get_cache(::LinearIsotropicElasticity) = nothing - -initial_material_state(::LinearIsotropicElasticity) = nothing - -""" - material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old, Δt=nothing; ) - -Return the stress tensor, stress tangent and the new `MaterialState` -The total (elastic) strain ϵ is given as input (state_old and Δt can be supplied, but have no influence) -The stress is calculated as - -```math -\\boldsymbol{\\sigma} = 2G \\boldsymbol{\\epsilon}^\\mathrm{dev} + 3K \\boldsymbol{\\epsilon}^\\mathrm{vol} -``` -where ``\\boldsymbol{\\epsilon}^\\mathrm{vol} = \\boldsymbol{I}:\\boldsymbol{\\epsilon} \\boldsymbol{I}/3.0`` -is the volumetric strain and ``\\boldsymbol{\\epsilon}^\\mathrm{dev} = \\boldsymbol{\\epsilon} - \\boldsymbol{\\epsilon}^\\mathrm{vol}`` -is the deviatoric strain. - -# Keyword arguments -- `cache`: Not used -- `options::Dict{Symbol, Any}`: Not used -""" -function material_response(m::LinearIsotropicElasticity, ϵ::SymmetricTensor{2,3}, state_old=initial_material_state(m), Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) - ν = (3*m.K - 2*m.G)/(2*(3*m.K+m.G)) # Calculate poissons ratio - - σ = 2 * m.G*dev(ϵ) + 3 * m.K*vol(ϵ) # Calculate stress - - # Create stiffness matrix - δ(i,j) = i == j ? 1.0 : 0.0 # helper function - Dfun(i,j,k,l) = 2.0*m.G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l)) - 𝔻 = SymmetricTensor{4, 3}(Dfun) - - # Return updated values - state = state_old - return σ, 𝔻, state -end \ No newline at end of file diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index f8c2257..941c674 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -1,4 +1,3 @@ -include("Elasticity.jl") include("IsotropicHardening.jl") include("KinematicHardening.jl") @@ -175,15 +174,6 @@ function initial_guess(m::VonMisesPlasticity, old::VonMisesPlasticityState{Nkin} return VonMisesPlasticityResidual(σ_trial,λ,β) end -# Functions for different elastic laws (elasticity = path independent, no state required) -function calculate_sigma(m::AbstractElasticity, ϵₑ) - σ, _, _ = material_response(m, ϵₑ, initial_material_state(m)) - return σ -end - -# Specialized function (faster) -calculate_sigma(m::LinearIsotropicElasticity, ϵₑ) = 3*m.K*vol(ϵₑ) + 2 * m.G * dev(ϵₑ) - function calculate_elastic_strain(old::VonMisesPlasticityState, ϵ, ν, Δλ) return ϵ - calculate_plastic_strain(old, ν, Δλ) end diff --git a/test/test_vonmises_plasticity.jl b/test/test_vonmises_plasticity.jl index 377182e..10ddd1b 100644 --- a/test/test_vonmises_plasticity.jl +++ b/test/test_vonmises_plasticity.jl @@ -46,5 +46,5 @@ σ, ∂σ∂ε, temp_state = material_response(m, ϵ, state, Δt; cache=cache) @test true # Check that it ran (throws error if not converged) - + end \ No newline at end of file From 03de2d03336590fc5dfff231419b23d2eb0f99b6 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 12:23:15 +0200 Subject: [PATCH 44/64] Updated documentation and test --- src/VonMisesPlasticity/VonMisesPlasticity.jl | 10 ++++++++++ test/test_vonmises_plasticity.jl | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index 941c674..73ab083 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -115,6 +115,16 @@ and ``i\\in[1,N_\\mathrm{kin}]``. - `cache`: Cache for the iterative solver, used by NLsolve.jl. It is strongly recommended to pre-allocate the cache for repeated calls to `material_response`. See [`get_cache`](@ref). - `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over. See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-options). By default the Newton solver will be used. + +# Example +```julia +m = VonMisesPlasticity(elastic=LinearElastic(E=210.e3, ν=0.3), + σ_y0=100.0, + isotropic=(Voce(Hiso=-100.e3, κ∞=-100.0),Voce(Hiso=10.e3, κ∞=200.0)), + kinematic=(ArmstrongFrederick(Hkin=200.e3, β∞=200.0), + OhnoWang(Hkin=1000.e3, β∞=200.0, mexp=3.0)) +``` + """ function material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, old::VonMisesPlasticityState{Nkin,T,N}, Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {T,N,Nkin} diff --git a/test/test_vonmises_plasticity.jl b/test/test_vonmises_plasticity.jl index 10ddd1b..672b262 100644 --- a/test/test_vonmises_plasticity.jl +++ b/test/test_vonmises_plasticity.jl @@ -2,7 +2,7 @@ # Basic setup with Voce isotropic hardening and one back-stress of Armstrong-Frederick time # constructor - m = VonMisesPlasticity(elastic=LinearIsotropicElasticity(E=210.e3, ν=0.3), + m = VonMisesPlasticity(elastic=LinearElastic(E=210.e3, ν=0.3), σ_y0=100.0, isotropic=(Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) kinematic=(ArmstrongFrederick(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple From 535298df0ee9b3235453f4f79fd4cd5f2b749b10 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 20:16:38 +0200 Subject: [PATCH 45/64] Adopted LinearElastic.jl from branch ka/use_full_strain --- src/LinearElastic.jl | 40 +++++++------------- src/VonMisesPlasticity/VonMisesPlasticity.jl | 4 ++ 2 files changed, 18 insertions(+), 26 deletions(-) diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index 33438d5..ba18d88 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -1,17 +1,11 @@ -# Definition of materials of type AbstractElasticity: -# It does not rely on previous states to calculate the stress and stiffness -# This allows it to be used inside other materials -abstract type AbstractElasticity <:AbstractMaterial end - """ LinearElastic(E, ν) - Isotropic linear elasticity. # Arguments - `E::Float64`: Young's modulus - `ν::Float64`: Poisson's ratio """ -struct LinearElastic <: AbstractElasticity +struct LinearElastic <: AbstractMaterial # parameters E::Float64 # Young's modulus ν::Float64 # Poisson's ratio @@ -35,31 +29,25 @@ function elastic_tangent_3D(E::T, ν::T) where T return C end -# TODO not sure if we could hardcode Float64 here - do we want to differentiate states? -struct LinearElasticState{dim, T, M} <: AbstractMaterialState - σ::SymmetricTensor{2,dim,T,M} # stress -end - -# TODO could these be automatically generated? -Base.zero(::Type{LinearElasticState{dim,T,M}}) where {dim,T,M} = LinearElasticState(zero(SymmetricTensor{2,dim,T,M})) +struct LinearElasticState <: AbstractMaterialState end # define which state belongs to the material -initial_material_state(::LinearElastic) = zero(LinearElasticState{3,Float64,6}) +initial_material_state(::LinearElastic) = LinearElasticState() + +# for restricted stress states to create buffer +get_stress_type(::LinearElasticState) = SymmetricTensor{2,3,Float64,6} # constitutive drivers generally operate in 3D # (we could specialize for lower dimensions if needed for performance) """ - material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState{3}) - -Return the stress tensor, stress tangent and the new `MaterialState` for the given strain ε such that - + material_response(m::LinearElastic, Δε::SymmetricTensor{2,3}) +Return the stress tensor and the stress tangent for the given strain ε such that ```math -\\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\boldsymbol{\\varepsilon} . +\\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\Delta \\boldsymbol{\\varepsilon} . ``` +No `MaterialState` is needed for the stress computation, thus if a state is handed over to `material_response`, the same state is returned. """ -function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState{3}, Δt=nothing; cache=nothing, options=nothing) - σ = calculate_sigma(m, ε) - return σ, m.Eᵉ, LinearElasticState(σ) -end - -calculate_sigma(m::LinearElastic, ε) = m.Eᵉ ⊡ ε # Useful when calling from other materials \ No newline at end of file +function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState=LinearElasticState(), Δt=nothing; cache=nothing, options=nothing) + σ = m.Eᵉ ⊡ ε + return σ, m.Eᵉ, state +end \ No newline at end of file diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index 73ab083..42e63c2 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -191,3 +191,7 @@ end function calculate_plastic_strain(old::VonMisesPlasticityState, ν, Δλ) return old.ϵₚ + Δλ*ν end + +# This function belongs in LinearElastic.jl, but to avoid conflict with ka/use_full_strain +# it is kept here until that branch is merged with main. +calculate_sigma(m::LinearElastic, ε) = m.Eᵉ ⊡ ε \ No newline at end of file From 89a676635094d00e8ef27d63c12ff087752acb16 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 20:21:46 +0200 Subject: [PATCH 46/64] Adopted test_linear_elastic.jl from branch ka/use_full_strain --- test/test_linear_elastic.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_linear_elastic.jl b/test/test_linear_elastic.jl index 1b8a1d3..32377e7 100644 --- a/test/test_linear_elastic.jl +++ b/test/test_linear_elastic.jl @@ -5,12 +5,11 @@ # initial state state = initial_material_state(m) - @test state.σ == zero(SymmetricTensor{2,3}) # constitutive driver ε = rand(SymmetricTensor{2,3}) σ, ∂σ∂ε, temp_state = material_response(m, ε, state) - @test σ == temp_state.σ + @test σ == m.Eᵉ ⊡ ε @test ∂σ∂ε == m.Eᵉ end From b2402dc3caa577c30a208ef9d389e2df424b80cb Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 23:36:24 +0200 Subject: [PATCH 47/64] Removed all boldface --- src/VonMisesPlasticity/KinematicHardening.jl | 28 ++++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/VonMisesPlasticity/KinematicHardening.jl b/src/VonMisesPlasticity/KinematicHardening.jl index 36fb4a5..f78ef25 100644 --- a/src/VonMisesPlasticity/KinematicHardening.jl +++ b/src/VonMisesPlasticity/KinematicHardening.jl @@ -9,7 +9,7 @@ end ArmstrongFrederick(;Hkin, β∞) = ArmstrongFrederick(Hkin, β∞) # Keyword argument constructor """ - get_evolution(param::ArmstrongFrederick, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + get_evolution(param::ArmstrongFrederick, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) Armstrong-Frederick kinematic hardening law (doi: 10.1179/096034007X207589) @@ -17,8 +17,8 @@ ArmstrongFrederick(;Hkin, β∞) = ArmstrongFrederick(Hkin, β∞) # Keyword g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) ``` """ -function get_evolution(param::ArmstrongFrederick, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) - param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ/param.β∞) +function get_evolution(param::ArmstrongFrederick, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) + param.Hkin * ((2.0/3.0) * ν - βᵢ/param.β∞) end # Delobelle (Combination of Armstrong-Frederick and Burlet-Cailletaud) @@ -30,7 +30,7 @@ end Delobelle(;Hkin, β∞, δ) = Delobelle(Hkin, β∞, δ) # Keyword argument constructor """ - get_evolution(param::Delobelle, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) + get_evolution(param::Delobelle, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law (doi: 10.1016/S0749-6419(95)00001-1) @@ -43,10 +43,10 @@ Delobelle(;Hkin, β∞, δ) = Delobelle(Hkin, β∞, δ) # Keyword argument c ``` """ -function get_evolution(param::Delobelle, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor) - AF_Term = (param.δ/param.β∞) * 𝛃ᵢ # Armstrong Frederick term - BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡𝛃ᵢ)/param.β∞)*ν # Burlet Cailletaud term - return param.Hkin * ((2.0/3.0) * 𝛎 - AF_Term - BC_Term) # Complete evolution +function get_evolution(param::Delobelle, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) + AF_Term = (param.δ/param.β∞) * βᵢ # Armstrong Frederick term + BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡βᵢ)/param.β∞)*ν # Burlet Cailletaud term + return param.Hkin * ((2.0/3.0) * ν - AF_Term - BC_Term) # Complete evolution end # Ohno-Wang @@ -58,7 +58,7 @@ end OhnoWang(;Hkin, β∞, mexp) = OhnoWang(Hkin, β∞, mexp) # Keyword argument constructor """ - get_evolution(param::OhnoWang{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} + get_evolution(param::OhnoWang{Tp}, ν::SecondOrderTensor, βᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} Kinematic hardening law according to Ohno-Wang (doi: 10.1016/0749-6419(93)90042-O) @@ -71,12 +71,12 @@ OhnoWang(;Hkin, β∞, mexp) = OhnoWang(Hkin, β∞, mexp) # Keyword argument ``` """ -function get_evolution(param::OhnoWang{Tp}, 𝛎::SecondOrderTensor, 𝛃ᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} - β_vm = vonmises_dev(𝛃ᵢ) +function get_evolution(param::OhnoWang{Tp}, ν::SecondOrderTensor, βᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} + β_vm = vonmises_dev(βᵢ) if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) - return param.Hkin * (2.0/3.0) * 𝛎 + 0*𝛃ᵢ + return param.Hkin * (2.0/3.0) * ν + 0*βᵢ end - mac_term = (macaulay(𝛎⊡𝛃ᵢ) /param.β∞) + mac_term = (macaulay(ν⊡βᵢ) /param.β∞) exp_term = (β_vm/param.β∞)^param.mexp - return param.Hkin * ((2.0/3.0) * 𝛎 - 𝛃ᵢ * mac_term * exp_term / β_vm ) + return param.Hkin * ((2.0/3.0) * ν - βᵢ * mac_term * exp_term / β_vm ) end \ No newline at end of file From c573fd66f33c62bad8332eb1f674fb67ddb870b4 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 23:45:14 +0200 Subject: [PATCH 48/64] Modified cache to be compatible with wrappers.jl --- src/VonMisesPlasticity/VonMisesPlasticity.jl | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index 42e63c2..71cbb8f 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -46,11 +46,6 @@ function initial_material_state(material::VonMisesPlasticity{T}) where {T} VonMisesPlasticityState(zero(SymmetricTensor{2,3,T}), 0.0, ntuple(i->zero(SymmetricTensor{2,3,T}), length(material.kinematic))) end -# Definition of material cache -struct VonMisesPlasticityCache{NL_TF, NL_TDF, NL_TX} - R_X_oncediff::OnceDifferentiable{NL_TF, NL_TDF, NL_TX} -end - function get_cache(material::VonMisesPlasticity{T,ElType,IsoType,KinType}) where {T,ElType,IsoType,KinType} nx = 7 + 6*length(material.kinematic) @@ -62,8 +57,7 @@ function get_cache(material::VonMisesPlasticity{T,ElType,IsoType,KinType}) where rf!(R, X) = vector_residual!(rf_tens, R, X, X_tensor) X0 = MVector{nx}(zeros(T, nx)) # X0 only for shape and type information here: - R_X_oncediff = OnceDifferentiable(rf!, X0, X0; autodiff = :forward) - return VonMisesPlasticityCache(R_X_oncediff) + return OnceDifferentiable(rf!, X0, X0; autodiff = :forward) end """ @@ -136,19 +130,19 @@ function material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, old: else x0 = initial_guess(m, old, ϵ) rf!(r_vector, x_vector) = vector_residual!((x)->residual(x,m,old,ϵ), r_vector, x_vector, x0) # Using x0 as template for residual instead of material as this is related to Tensors - update_cache!(cache.R_X_oncediff, rf!) + update_cache!(cache, rf!) - tomandel!(cache.R_X_oncediff.x_f, x0) + tomandel!(cache.x_f, x0) # Should this be centrally managed? I.e. process_options or similar? nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton)) haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method # Solve local problem: - result = NLsolve.nlsolve(cache.R_X_oncediff, cache.R_X_oncediff.x_f; nlsolve_options...) + result = NLsolve.nlsolve(cache, cache.x_f; nlsolve_options...) if result.f_converged x = frommandel(VonMisesPlasticityResidual{Nkin}, result.zero::MVector{7 + 6*Nkin, T}) - dRdx = cache.R_X_oncediff.DF + dRdx = cache.DF inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx)) dσdϵ = inv_J_σσ ⊡ dσdϵ_elastic σ_red_dev = dev(x.σ) - sum(x.β) From ab5efdd649b635bce1c373a9faef610430e0227d Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Tue, 22 Jun 2021 23:47:00 +0200 Subject: [PATCH 49/64] Added proper testing to VonMisesPlasticity --- test/jld2_files/VonMisesPlasticity1.jld2 | Bin 0 -> 15133 bytes test/test_vonmises_plasticity.jl | 162 ++++++++++++++++++----- 2 files changed, 130 insertions(+), 32 deletions(-) create mode 100644 test/jld2_files/VonMisesPlasticity1.jld2 diff --git a/test/jld2_files/VonMisesPlasticity1.jld2 b/test/jld2_files/VonMisesPlasticity1.jld2 new file mode 100644 index 0000000000000000000000000000000000000000..59a83a7eb4936326b9627bb64cd103bfc0393881 GIT binary patch literal 15133 zcmeHO2~bo=8tys7;ZhKzXgr3jA|BveTycYgB2b<|QB2m23^GOy0>a>m(Tc@+7F!SPWNTLT_>6R|E<_)qN7^77 zq7MlR(P_g>!71Zy+JXJLB3Jv(GEfj76YK5K)9ka?~sA^2$BOg1&b#L zr1j0nvWAdrJCS4|=rBmToBE^Eb5+dlLCw^jn3I`CI^m{felMz0e{5QY#TI5FrY6h} zv#4o70~;&dSpU`xKhc_=ospFtGAwr@W-@F17-~RGS)7)|Sz6K;E186RXPB;G@y_Pu z4tR`a=k7;*D1sz1FL{M1L-}# zdNqYR$>8weUx^D0%~_aN7Xt5zfu1R4K&69Vo@ZqV6+VZ9T{Dg3pMtzzp*Zg`d7R0v zt*Cu2led_RWwn8gOsf59+!!VgF&WT?+Gj9%g~_P4)LwP;A79XUT~sZ>uvf{+=RJdvqhg^&Shq|8R-BV{%re<^b)(n^`bP>__lQ6(Z$(0~#V9wx%G5cDSa zbz9rf_heiQq^<*S_^NrXy<^E6Y==_>C{*h(Yq(xQYCWmAw1m)altXeM!K4X)Wa6MV z#jm|^X-kQ#cN*k{#*Vk9rH&@kCH1;&p`J$;1gjmPdSBcLx44{dTT<>rxP|D26gmwl z5t3#`8ce$yYq0$CHI2T9}UxBCe-#f zoa&8zE-&6bR0jOKqVoksnjCYfUoZS;G((!nDP1O~G);L<@3;T05i;PHd_OIrP@88? zSemrIcDD$h*=pvF^t(9%Jm}Sb+@B-^{!Yq2KeZL-o4uM9J$h0m!p|Svc2mEpc>+A> zuaEB0LI%8b)SRbDfs4%50eg>l?G@p_JbeJ2eUdN0gYG%%l|wS%QyKha4*nDe5Bk2F z%f`xpubTNJ`GW(q&F6QF$bEZ+uQ1y72lIwhT4p_`4?9Q2?;~4yU4j$lOHfiY01ug;<8Fya9=ZQ`8)w0 z^m_&Qa`v4154Rsb6~ENnzbtR;{rw{Rtf;LkR~QQgc+d~`xqO(8{{M5=)f(UZk@jw9 zD)$6Z*))jCB=&F3&E2T}usfA;!Bj2}q4HKwD#P_u<{PLyXQZ-AD3yN-qq1gBx5xH* zDNd+VGgam=JPs}X;3VCpZi4WAm!uuvcuCsvotLB?-+D>f@x7O%JsjyI?GY$K(jJLF z#-+xu*P(c6haSBt?J%H`(heg^mUfuXIB7>H%9M75;Tvzs;fA9b(vAo;SK1MY7E3w| zIT+F?Yyq#Y)-SK1MZDy1D^_|{Hp=7wE%zv# z^D6Fy+taPV>_)heM`4nfw9PGwH_`i`Mt1zwM`5y~k+ywl^HRE##VDlC!h<@}I87%7 zrbnx2iO;B?rTauvtN4tMu;&Z9KQV)a*U8bg^yodo4{;62hWhb&0FNK1vp9$c$>(#{ zkGDgJWEs&M(1RdtVpgzz#B)pT+lAI$ z>z(2N9-q${{umB#4a--%tb$hQ8hL#l^gg+=FXnSqa!_2{9FfnqxqmUdi~MDBI?V4A>(xtRxUa{w3DJ%hyy#wuUMW8u79u1qXKjjJcZBq%N4io_B!()r#OJe=R3{t zGdcW=9Nrr&-<|hMF3Q(Dmw^S%Aat;aK?o;t+= zJU$=H<7W=P1Bdq=%ct95>Me)Q-&aQ|MwMaS#taYI*zIH`jN{s(&M$6 zpWM&)(#NYhti}6ho%iSQXDh`0GeG&bOr692sb%_O7FRrR!@`StasRw-w~PB{HK)V; z3S*9xZBjblKLc5w`&|DPZ;s3iTF&*E zAM)wlCbQJhCvyM$P3YDAl`!uE;{r9)`-tlu^l}y#IUzmEaO$cV;;ZfU5~n)E?S6R9 zW*l-&hyy#wGM2Z3!&~^pYi~8%qY>77>B+J+w&qCCp+7JAe1QKshtK0dzRUN9HkZTa ztda>M-7ueagBS2w z&hqTz`qvB{Ss38J=Sbz{m6k^E$>%%3@XI)S9uM+8b{pbv=M!yF+J9{ZzuyDtp(Imt z!@mSm_$V|w+VKD7Xyfb#`0@^CH^Nt7IJ*fdoSZ!r%7n;=C2|l;Pk7CuLtgR1A>(hX z;URlx+p?_L+1BjQ6Vg%*U3}TNnbT%Tf7hCBqs}7s-8MH&1n%i(|50J~(QuX i==2 && j==1 ? ϵ21_init : zero(typeof(ϵ21_init))) + σ, dσdϵ, state = material_response(m, ϵ, state, t_max/num_steps; cache=cache, options=options) + for k=1:num_steps + ϵ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? get_ramp_value(m, ϵ21_init, ϵ21_max, k, num_steps) : zero(typeof(ϵ21_max))) + σ, dσdϵ, state = material_response(m, ϵ, state, t_max/num_steps; cache=cache, options=options) + end + return σ, dσdϵ, state +end + @testset "VonMisesPlasticity" begin # Basic setup with Voce isotropic hardening and one back-stress of Armstrong-Frederick time # constructor - m = VonMisesPlasticity(elastic=LinearElastic(E=210.e3, ν=0.3), - σ_y0=100.0, + E=210.e3; ν=0.3; σ_y0=100.0 + m = VonMisesPlasticity(elastic=LinearElastic(E=E, ν=ν), + σ_y0=σ_y0, isotropic=(Voce(Hiso=100000.0, κ∞=100.0),), # Can add more dragstresses by more entries in Tuple) kinematic=(ArmstrongFrederick(Hkin=1000000.0, β∞=200.0),) # Can add more backstresses by more entries in Tuple ) - cache = get_cache(m) + t_max = 1.0 + num_steps = 100 + ϵ11_max = 2.0 * σ_y0 / E + σ, dσdϵ, state = uniaxial_loading(m, 0.0, ϵ11_max, num_steps, t_max) + # Check that yield criterion is zero + σ_full = SymmetricTensor{2,3}((i,j)->i==j && i==1 ? σ[1] : 0.0) + @test MaterialModels.vonmises(σ_full-sum(state.β)) ≈ σ_y0 + sum(MaterialModels.get_hardening.(m.isotropic, state.λ)) + + G = E/(2*(1+ν)) + ϵ21_max = 2.0 * (σ_y0/sqrt(3.0))/(2*G) + σ, dσdϵ, state = shear_loading(m, 0.0, ϵ21_max, num_steps, t_max) + # Check that yield criterion is zero + @test MaterialModels.vonmises(σ-sum(state.β)) ≈ σ_y0 + sum(MaterialModels.get_hardening.(m.isotropic, state.λ)) - # initial state (not used here) + # Check that error is thrown if material doesn't converge state = initial_material_state(m) + cache = get_cache(m) + options = Dict{Symbol, Any}(:nlsolve_params=> Dict{Symbol,Any}(:method=>:newton, :ftol=>1.e-999)) + ϵ = SymmetricTensor{2,3}((i,j)->i==j && i==1 ? ϵ11_max : zero(typeof(ϵ11_max))) + @test_throws ErrorException σ, dσdϵ, state = material_response(m, ϵ, state, nothing; cache=cache, options=options) - # strain at yield point for uniaxial stress some test case at some plastic strain: - ϵ = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.004 : -0.000782017) : 0.0) +end + - ϵₚ_old = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.00153388 : -0.000766938) : 0.0) - λ_old = 0.0015338757291717328 - β_old = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 133.268 : -66.6339) : 0.0) +@testset "VonMisesPlasticity vs Plastic" begin + # Setup materials + E = 200.e3; ν=0.3; σ_y0=200.0; Hiso=25.0; κ∞=13.0; Hkin=25.0; β∞=13.0 + plastic = Plastic(E=E, ν=ν, σ_y=σ_y0, H=Hiso+Hkin, r=Hiso/(Hiso+Hkin), κ_∞=κ∞, α_∞=β∞) + vonmises_plasticity = VonMisesPlasticity(elastic=LinearElastic(E=E, ν=ν), σ_y0=σ_y0, + isotropic=(Voce(Hiso=Hiso, κ∞=κ∞),), + kinematic=(ArmstrongFrederick(Hkin=Hkin, β∞=β∞),)) + # Uniaxial loading test + t_max = 1.0 + options = Dict{Symbol, Any}(:nlsolve_params=> Dict{Symbol,Any}(:method=>:newton, :ftol=>1.e-12)) + num_steps = 100 + ϵ11_yld = σ_y0 / E + ϵ11_init = 0.99 * ϵ11_yld + ϵ11_max = 1.5 * ϵ11_yld + σ_plastic, dσdϵ_plastic, _ = uniaxial_loading(plastic, ϵ11_init, ϵ11_max, num_steps, t_max, options) + σ_vmplast, dσdϵ_vmplast, _ = uniaxial_loading(vonmises_plasticity, ϵ11_init, ϵ11_max, num_steps, t_max, options) - state_old = MaterialModels.VonMisesPlasticityState(ϵₚ_old,λ_old, (β_old,)) + @test σ_plastic ≈ σ_vmplast + @test isapprox(dσdϵ_plastic, dσdϵ_vmplast; rtol=1.e-4) # Requires many time steps and low tolerance to avoid small error due to different implementations with default tolerance - Δt = 1.0 # No influence... + # Shear loading test + G = E/(2*(1+ν)) + ϵ21_yld = (σ_y0/sqrt(3.0))/(2*G) + ϵ21_init = 0.9 * ϵ21_yld + ϵ21_max = 2.0 * ϵ21_yld + σ_plastic, dσdϵ_plastic, _ = shear_loading(plastic, ϵ21_init, ϵ21_max, num_steps, t_max) + σ_vmplast, dσdϵ_vmplast, _ = shear_loading(vonmises_plasticity, ϵ21_init, ϵ21_max, num_steps, t_max) + + @test σ_plastic ≈ σ_vmplast + @test dσdϵ_plastic ≈ dσdϵ_vmplast +end + +@testset "VonMisesPlasticity jld2" begin + E = 200.e3; ν=0.3; σ_y0=200.0; Hiso=25.0; κ∞=13.0; Hkin=25.0; β∞=13.0 + vonmises_plasticity = VonMisesPlasticity(elastic=LinearElastic(E=E, ν=ν), σ_y0=σ_y0, + isotropic=(Voce(Hiso=Hiso, κ∞=κ∞),), + kinematic=(ArmstrongFrederick(Hkin=Hkin, β∞=β∞),)) + + loading = get_Plastic_loading() + check_jld2(vonmises_plasticity, loading, "VonMisesPlasticity1"; OVERWRITE_JLD2=false) +end + +get_evolution(af, db, ow, ν, β) = MaterialModels.get_evolution.((af, db, ow), ntuple(i->ν, 3), ntuple(i->β, 3)) + +@testset "KinematicHardening" begin + Hkin=10.e3; β∞=30.0 + af = ArmstrongFrederick(Hkin=Hkin, β∞=β∞) + db = Delobelle(Hkin=Hkin, β∞=β∞, δ=0.5) + ow = OhnoWang(Hkin=Hkin, β∞=β∞, mexp=3.0) + σ_red_dev = dev(rand(SymmetricTensor{2,3})) + ν = (3/2)*σ_red_dev/MaterialModels.vonmises(σ_red_dev) + β = zero(SymmetricTensor{2,3}) - σ, ∂σ∂ε, temp_state = material_response(m, ϵ, state, Δt; cache=cache) - - @test true # Check that it ran (throws error if not converged) - - # Example with a more advanced material: - # Linear isotropic elasticity - # Two isotropic hardening laws: Voce and Swift - # Two back-stresses, one Armstrong-Frederick and one Ohno-Wang - m = VonMisesPlasticity(elastic=LinearElastic(E=210.e3, ν=0.3), - σ_y0=100.0, - isotropic=(Voce(Hiso=100000.0, κ∞=100.0), - Swift(K=100.0, λ0=1.0e-2, n=0.5)), - kinematic=(ArmstrongFrederick(Hkin=40.e3, β∞=200.0), - OhnoWang(Hkin=30.e3, β∞=200.0, mexp=4.0)) - ) + # Check that all give the same initial hardening modulus + af_h0, db_h0, ow_h0 = get_evolution(af, db, ow, ν, β) + @test af_h0 ≈ db_h0 + @test af_h0 ≈ ow_h0 - cache = get_cache(m) - state = initial_material_state(m) - Δt = 1.0 # No influence... - ϵ = SymmetricTensor{2, 3}((i,j) -> i==j ? (i==1 ? 0.1/100.0 : -0.3*0.1/100.0) : 0.0) - σ, ∂σ∂ε, temp_state = material_response(m, ϵ, state, Δt; cache=cache) + # Check that ArmstrongFrederick and Delobelle give the same result when β and ν are aligned + β = (1.0 + rand())*ν # Ensure scaling > 0 + af_h1, db_h1, ow_h1 = get_evolution(af, db, ow, ν, β) + @test af_h1 ≈ db_h1 # Should be equal + @test !(af_h0 ≈ ow_h1) # Should not be equal +end + +@testset "IsotropicHardening" begin + λ = 0.01 + rand() + dλ = 1.e-8 + + # Voce hardening + Hiso = 200.0; κ∞=10.0 + voce = Voce(Hiso=Hiso, κ∞=κ∞) + κ = MaterialModels.get_hardening(voce, λ) + dκ = MaterialModels.get_hardening(voce, λ+dλ) - κ + @test isapprox(dκ/dλ, Hiso*(1 - κ/κ∞), rtol=1.e-3) - @test true # Check that it ran (throws error if not converged) + # Swift hardening + K=10.0 + λ0=1.e-3 + n=2.0 + swift = Swift(K=K, λ0=λ0, n=n) + κ = MaterialModels.get_hardening(swift, λ) + @test κ ≈ K*(λ0 + λ)^n # Same equation, but at least protects against changes due to optimizations end \ No newline at end of file From 4b58faa59286f097c41b359a9a26dd1591e8a8c7 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 23 Jun 2021 08:31:55 +0200 Subject: [PATCH 50/64] Fixed documentation errors for LinearElastic --- docs/src/index.md | 2 +- src/LinearElastic.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 81f20b5..2f68836 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,7 +1,7 @@ # Materials ```@docs LinearElastic -material_response(m::LinearElastic, Δε::SymmetricTensor{2,3}, state::LinearElasticState{3}) +material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState) Plastic material_response(m::Plastic, Δε::SymmetricTensor{2,3,T,6}, state::PlasticState{3}; kwargs...) where T ``` diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index ba18d88..55e7a7c 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -40,7 +40,7 @@ get_stress_type(::LinearElasticState) = SymmetricTensor{2,3,Float64,6} # constitutive drivers generally operate in 3D # (we could specialize for lower dimensions if needed for performance) """ - material_response(m::LinearElastic, Δε::SymmetricTensor{2,3}) + material_response(m::LinearElastic, ε::SymmetricTensor{2,3}) Return the stress tensor and the stress tangent for the given strain ε such that ```math \\boldsymbol{\\sigma} = \\mathbf{E}^\\text{e} : \\Delta \\boldsymbol{\\varepsilon} . From 48d43495083b8436754d8082803fb1b136740c92 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 23 Jun 2021 08:35:02 +0200 Subject: [PATCH 51/64] Added VonMisesPlasticity to documentation --- docs/src/index.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/index.md b/docs/src/index.md index 2f68836..b02275d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,4 +4,6 @@ LinearElastic material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState) Plastic material_response(m::Plastic, Δε::SymmetricTensor{2,3,T,6}, state::PlasticState{3}; kwargs...) where T +VonMisesPlasticity +material_response(m::VonMisesPlasticity, ε::SymmetricTensor{2,3}, state::VonMisesPlasticityState; kwargs...) where T ``` From ceb953d02eaf14bbcc6625bbbbf8351413164221 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 23 Jun 2021 08:55:45 +0200 Subject: [PATCH 52/64] Removed explicit state for VonMisesPlasticity (not exported) --- docs/src/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index b02275d..dd97af2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -5,5 +5,5 @@ material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElast Plastic material_response(m::Plastic, Δε::SymmetricTensor{2,3,T,6}, state::PlasticState{3}; kwargs...) where T VonMisesPlasticity -material_response(m::VonMisesPlasticity, ε::SymmetricTensor{2,3}, state::VonMisesPlasticityState; kwargs...) where T +material_response(m::VonMisesPlasticity, ε::SymmetricTensor{2,3}, state, Δt; kwargs...) where T ``` From 7941f1198d1f7090abd5f6a8a16ca26553cc20ee Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Wed, 23 Jun 2021 08:58:20 +0200 Subject: [PATCH 53/64] Changed to using ForwardDiff when checking isotropic hardening derivative --- test/runtests.jl | 1 + test/test_vonmises_plasticity.jl | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index ca7da1b..e170067 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test using Rotations using StaticArrays using JLD2 +using ForwardDiff include("test_utils.jl") include("test_linear_elastic.jl") diff --git a/test/test_vonmises_plasticity.jl b/test/test_vonmises_plasticity.jl index 0c2c2b0..13c0093 100644 --- a/test/test_vonmises_plasticity.jl +++ b/test/test_vonmises_plasticity.jl @@ -134,9 +134,9 @@ end Hiso = 200.0; κ∞=10.0 voce = Voce(Hiso=Hiso, κ∞=κ∞) κ = MaterialModels.get_hardening(voce, λ) - dκ = MaterialModels.get_hardening(voce, λ+dλ) - κ - @test isapprox(dκ/dλ, Hiso*(1 - κ/κ∞), rtol=1.e-3) - + dκdλ = ForwardDiff.derivative(λarg->MaterialModels.get_hardening(voce, λarg), λ) + @test dκdλ ≈ Hiso*(1 - κ/κ∞) + # Swift hardening K=10.0 λ0=1.e-3 From 0fe9f3c69fd900ce3ebbaaf49a3cbe2b9359b5c1 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 25 Jun 2021 00:20:52 +0200 Subject: [PATCH 54/64] Updated documentation formatting --- docs/Project.toml | 3 +- docs/src/index.md | 17 +++- src/VonMisesPlasticity/IsotropicHardening.jl | 58 ++++++----- src/VonMisesPlasticity/KinematicHardening.jl | 101 +++++++++++-------- src/VonMisesPlasticity/VonMisesPlasticity.jl | 21 ++-- test/test_vonmises_plasticity.jl | 2 +- 6 files changed, 125 insertions(+), 77 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 3879473..1fe7fea 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +MaterialModels = "d7e62f90-dad9-4fee-a9a8-264218a34fd2" [compat] -Documenter = "0.26" \ No newline at end of file +Documenter = "0.26" diff --git a/docs/src/index.md b/docs/src/index.md index dd97af2..f3b22d1 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,9 +1,24 @@ # Materials +## Elastic materials ```@docs LinearElastic material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState) +``` +## Plastic materials +```@docs Plastic material_response(m::Plastic, Δε::SymmetricTensor{2,3,T,6}, state::PlasticState{3}; kwargs...) where T VonMisesPlasticity -material_response(m::VonMisesPlasticity, ε::SymmetricTensor{2,3}, state, Δt; kwargs...) where T +material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, old::MaterialModels.VonMisesPlasticityState; kwargs...) ``` +### Isotropic hardening laws +```@docs +Voce +Swift +``` +### Kinematic hardening laws +```@docs +ArmstrongFrederick +Delobelle +OhnoWang +``` \ No newline at end of file diff --git a/src/VonMisesPlasticity/IsotropicHardening.jl b/src/VonMisesPlasticity/IsotropicHardening.jl index 5faaa7a..6f7777b 100644 --- a/src/VonMisesPlasticity/IsotropicHardening.jl +++ b/src/VonMisesPlasticity/IsotropicHardening.jl @@ -1,32 +1,49 @@ # Isotropic hardening abstract type AbstractIsotropicHardening end -# Voce type of isotropic hardening (exponentially saturating) +""" + Voce(Hiso, κ∞) + +Exponentially saturating isotropic hardening + +```math +\\kappa_i = g_{\\mathrm{iso},i}(\\lambda) = \\kappa_\\infty \\left[1 - \\mathrm{exp}\\left(\\frac{H_\\mathrm{iso}}{\\kappa_\\infty} \\lambda \\right)\\right] +``` +or alternatively as differential equations +```math +\\dot{\\kappa_i} = \\dot{\\lambda} H_\\mathrm{iso} \\left[1 - \\frac{\\kappa_i}{\\kappa_\\infty}\\right] +``` + +# Arguments +- `Hiso`: Isotropic hardening modulus, ``H_\\mathrm{iso}`` +- `κ∞`: Saturation hardening value, ``\\kappa_\\infty`` + +""" struct Voce{T} <:AbstractIsotropicHardening Hiso::T # Initial hardening modulus κ∞::T # Saturation stress end Voce(;Hiso, κ∞) = Voce(Hiso, κ∞) # Keyword argument constructor -""" - get_hardening(param::Voce, λ::Number) - - Exponentially saturating isotropic hardening - - ```math - \\kappa_i = g_{\\mathrm{iso},i}(\\lambda) = \\kappa_\\infty \\left[1 - \\mathrm{exp}\\left(\\frac{H_\\mathrm{iso}}{\\kappa_\\infty} \\lambda \\right)\\right] - ``` - or alternatively as differential equations - ```math - \\dot{\\kappa_i} = \\dot{\\lambda} H_\\mathrm{iso} \\left[1 - \\frac{\\kappa_i}{\\kappa_\\infty}\\right] - ``` - -""" function get_hardening(param::Voce, λ::Number) param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) end -# Swift type of kinematic hardening (power law) +""" + Swift(K, λ0, n) + +Isotropic hardening by the Swift power law + +```math +\\kappa_i = g_{\\mathrm{iso},i}(\\lambda) = K \\left[\\lambda_0 + \\lambda \\right]^n +``` + +# Arguments +- `K`: ``K`` +- `λ0`: ``\\lambda_0`` +- `n`: ``n`` + +""" struct Swift{T} <:AbstractIsotropicHardening K::T λ0::T @@ -34,16 +51,7 @@ struct Swift{T} <:AbstractIsotropicHardening end Swift(;K, λ0, n) = Swift(K, λ0, n) # Keyword argument constructor -""" - get_hardening(param::Voce, λ::Number) - - Isotropic hardening by the Swift power law - ```math - \\kappa_i = g_{\\mathrm{iso},i}(\\lambda) = K \\left[\\lambda_0 + \\lambda \\right]^n - ``` - -""" function get_hardening(param::Swift, λ::Number) param.K * (param.λ0 + λ)^param.n end diff --git a/src/VonMisesPlasticity/KinematicHardening.jl b/src/VonMisesPlasticity/KinematicHardening.jl index f78ef25..16439c2 100644 --- a/src/VonMisesPlasticity/KinematicHardening.jl +++ b/src/VonMisesPlasticity/KinematicHardening.jl @@ -1,27 +1,48 @@ # Kinematic hardening abstract type AbstractKinematicHardening end -# Armstrong-Frederick +""" + ArmstrongFrederick(Hkin, β∞) + +Armstrong-Frederick kinematic hardening law (doi: 10.1179/096034007X207589) + +```math +g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) = H_\\mathrm{kin} (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) +``` + +# Arguments +- `Hkin`: Kinematic hardening modulus, ``H_\\mathrm{kin}`` +- `β∞`: Effective back-stress saturation value, ``\\beta_\\infty`` +""" struct ArmstrongFrederick{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress end ArmstrongFrederick(;Hkin, β∞) = ArmstrongFrederick(Hkin, β∞) # Keyword argument constructor -""" - get_evolution(param::ArmstrongFrederick, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) - - Armstrong-Frederick kinematic hardening law (doi: 10.1179/096034007X207589) - - ```math - g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin (\\frac{2}{3}\\boldsymbol{\\nu} - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty}) - ``` -""" function get_evolution(param::ArmstrongFrederick, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) param.Hkin * ((2.0/3.0) * ν - βᵢ/param.β∞) end -# Delobelle (Combination of Armstrong-Frederick and Burlet-Cailletaud) +""" + Delobelle(Hkin, β∞, δ) + +Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law +(doi: 10.1016/S0749-6419(95)00001-1) + +```math +g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) = H_\\mathrm{kin} \\left[\\frac{2}{3}\\boldsymbol{\\nu} + - \\delta\\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty} + - \\frac{2}{3\\beta_\\infty}\\left[1 - \\delta\\right]\\left[\\boldsymbol{\\nu}:\\boldsymbol{\\beta}_i\\right]\\boldsymbol{\\nu} + \\right] +``` + +# Arguments +- `Hkin`: Kinematic hardening modulus, ``H_\\mathrm{kin}`` +- `β∞`: Effective back-stress saturation value, ``\\beta_\\infty`` +- `δ`: Amount of Armstrong-Frederick type of kinematic hardening, ``\\delta`` + +""" struct Delobelle{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress @@ -29,54 +50,48 @@ struct Delobelle{T} <: AbstractKinematicHardening end Delobelle(;Hkin, β∞, δ) = Delobelle(Hkin, β∞, δ) # Keyword argument constructor -""" - get_evolution(param::Delobelle, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) - - Kinematic hardening law according to Delobelle, which combines the Armstrong-Frederick law with the Burlet-Cailletaud law - (doi: 10.1016/S0749-6419(95)00001-1) - - ```math - g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin \\left[\\frac{2}{3}\\boldsymbol{\\nu} - - \\delta\\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty} - - \\frac{2}{3\\beta_\\infty}\\left[1 - \\delta\\right]\\left[\\boldsymbol{\\nu}:\\boldsymbol{\\beta}_i\right]\\boldsymbol{\\nu} - \\right] - ``` - -""" function get_evolution(param::Delobelle, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) AF_Term = (param.δ/param.β∞) * βᵢ # Armstrong Frederick term BC_Term = (2.0/3.0) * (1.0-param.δ)*((ν⊡βᵢ)/param.β∞)*ν # Burlet Cailletaud term return param.Hkin * ((2.0/3.0) * ν - AF_Term - BC_Term) # Complete evolution end -# Ohno-Wang + +""" + OhnoWang(Hkin, β∞, m) + +Kinematic hardening law according to Ohno-Wang (doi: 10.1016/0749-6419(93)90042-O) + +```math +g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) = H_\\mathrm{kin} \\left[\\frac{2}{3}\\boldsymbol{\\nu} + - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty} + \\frac{\\langle \\boldsymbol{\\nu}:\\boldsymbol{\\beta}_i \\rangle}{\\beta_\\infty} + \\left[\\frac{\\beta_i^\\mathrm{vM}}{\\beta_\\infty}\\right]^m + \\right] +``` +where ``\\langle x \\rangle`` is 0 if ``x\\leq 0`` and ``x`` if ``x>0``. +``\\beta_i^\\mathrm{vM} = \\sqrt{2\\boldsymbol{\\beta}_i:\\boldsymbol{\\beta}_i/3}``, noting that +``\\boldsymbol{\\beta}_i`` is deviatoric. + +# Arguments +- `Hkin`: Kinematic hardening modulus, ``H_\\mathrm{kin}`` +- `β∞`: Effective back-stress saturation value, ``\\beta_\\infty`` +- `m`: Exponent in the OhnoWang equation, ``m`` + +""" struct OhnoWang{T} <: AbstractKinematicHardening Hkin::T # Initial hardening modulus β∞::T # Saturation stress - mexp::T # Ohno Wang exponent + m::T # Ohno Wang exponent end -OhnoWang(;Hkin, β∞, mexp) = OhnoWang(Hkin, β∞, mexp) # Keyword argument constructor +OhnoWang(;Hkin, β∞, m) = OhnoWang(Hkin, β∞, m) # Keyword argument constructor -""" - get_evolution(param::OhnoWang{Tp}, ν::SecondOrderTensor, βᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} - - Kinematic hardening law according to Ohno-Wang (doi: 10.1016/0749-6419(93)90042-O) - - ```math - g_{\\mathrm{kin},i}(\\nu, \\beta_i) = Hkin \\left[\\frac{2}{3}\\boldsymbol{\\nu} - - \\frac{\\boldsymbol{\\beta}_i}{\\beta_\\infty} - \\frac{\\langle \\boldsymbol{\\nu}:\\boldsymbol{\\beta}_i \\rangle}{\\beta_\\infty} - \\left[\\frac{\\beta_\\mathrm{vM}}{\\beta_\\infty}\right]^\\mathrm{mexp} - \\right] - ``` - -""" function get_evolution(param::OhnoWang{Tp}, ν::SecondOrderTensor, βᵢ::SecondOrderTensor{dim,Tβ}) where{Tp,Tβ,dim} β_vm = vonmises_dev(βᵢ) if β_vm < param.β∞ * eps(promote_type(Tp,Tβ)) return param.Hkin * (2.0/3.0) * ν + 0*βᵢ end mac_term = (macaulay(ν⊡βᵢ) /param.β∞) - exp_term = (β_vm/param.β∞)^param.mexp + exp_term = (β_vm/param.β∞)^param.m return param.Hkin * ((2.0/3.0) * ν - βᵢ * mac_term * exp_term / β_vm ) end \ No newline at end of file diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index 71cbb8f..e7b5789 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -2,7 +2,16 @@ include("IsotropicHardening.jl") include("KinematicHardening.jl") -# Definition of material properties +""" + VonMisesPlasticity(elastic, σ_y0, isotropic, kinematic) +von Mises plasticity with modular elastic law, and multiple modular isotropic and kinematic +hardening contributions. +# Arguments +- `elastic::AbstractMaterial`: Elastic law, see e.g. [`LinearElastic`](@ref) +- `σ_y0::Float64`: Initial yield limit +- `isotropic::NTuple{Niso,AbstractIsotropicHardening}`: Isotropic hardening laws, see e.g. [`Voce`](@ref) +- `kinematic::NTuple{Niso,AbstractKinematicHardening}`: Kinematic hardening laws, see e.g. [`ArmstrongFrederick`](@ref) +""" struct VonMisesPlasticity{T,ElasticType,IsoType,KinType} <:AbstractMaterial elastic::ElasticType # Elastic definition σ_y0::T # Initial yield limit @@ -61,14 +70,14 @@ function get_cache(material::VonMisesPlasticity{T,ElType,IsoType,KinType}) where end """ - material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, state::VonMisesPlasticityState, Δt; ) + material_response(m::VonMisesPlasticity, ϵ::SymmetricTensor{2,3}, old::VonMisesPlasticityState; kwargs...) Return the stress tensor, stress tangent, the new `MaterialState` and a boolean specifying if local iterations converged. The total strain ε and previous material state `state` are given as input, Δt has no influence as the material is rate independent. By specifying different laws in `m.elastic`, `m.isotropic`, and `m.kinematic`, different models can be obtained, within the general equations specified below. Note that `m.isotropic` and `m.kinematic` are of type Tuple, in which each -element contain a hardening law of type ``AbstractIsotropicHardening`` and ``AbstractKinematicHardening``, respectively. +element contain a hardening law of type `AbstractIsotropicHardening` and `AbstractKinematicHardening`, respectively. The stress is calculated from the elastic strains, ``\\boldsymbol{\\epsilon}_\\mathrm{e}``, obtained via the additive decomposition, ``\\boldsymbol{\\epsilon} = \\boldsymbol{\\epsilon}_\\mathrm{e} + \\boldsymbol{\\epsilon}_\\mathrm{p}``. @@ -78,12 +87,12 @@ Von Mises yield function: ```math \\Phi = \\sqrt{\\frac{3}{2}} \\left| \\text{dev} \\left( \\boldsymbol{\\sigma} - \\boldsymbol{\\beta} \\right) \\right| - \\sigma_y - \\kappa ``` -where ``\\boldsymbol{\\beta} = \\sum_{i=1}^{N_\\mathrm{kin} \\boldsymbol{\\beta}_i`` is the total back-stress. +where ``\\boldsymbol{\\beta} = \\sum_{i=1}^{N_\\mathrm{kin}} \\boldsymbol{\\beta}_i`` is the total back-stress. The evolution of ``\\boldsymbol{\\beta}_i`` is given by the kinematic hardening, specified below and using `m.kinematic` Associative plastic flow is used to obtain the plastic strains, ```math -\\dot{\\epsilon}_{\\mathrm{p}} = \\dot{\\lambda} \\frac{\\partial \\Phi}{\\boldsymbol{\\sigma}} +\\dot{\\epsilon}_{\\mathrm{p}} = \\dot{\\lambda} \\frac{\\partial \\Phi}{\\partial \\boldsymbol{\\sigma}} = \\dot{\\lambda} \\boldsymbol{\\nu} ``` @@ -100,7 +109,7 @@ where ``g_{\\mathrm{iso},i}(\\lambda)`` is specified by `m.isotropic[i]` Kinematic hardening is formulated as ```math -\\boldsymbol{\\beta}_i = \\dot{\\lambda} g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) +\\dot{\\boldsymbol{\\beta}}_i = \\dot{\\lambda} g_{\\mathrm{kin},i}(\\nu, \\boldsymbol{\\beta}_i) ``` where ``g_{\\mathrm{kin},i}(\\boldsymbol{\\nu}, \\boldsymbol{\\beta}_i)`` is specified by `m.kinematic[i]` and ``i\\in[1,N_\\mathrm{kin}]``. diff --git a/test/test_vonmises_plasticity.jl b/test/test_vonmises_plasticity.jl index 13c0093..79fb7f3 100644 --- a/test/test_vonmises_plasticity.jl +++ b/test/test_vonmises_plasticity.jl @@ -109,7 +109,7 @@ get_evolution(af, db, ow, ν, β) = MaterialModels.get_evolution.((af, db, ow), Hkin=10.e3; β∞=30.0 af = ArmstrongFrederick(Hkin=Hkin, β∞=β∞) db = Delobelle(Hkin=Hkin, β∞=β∞, δ=0.5) - ow = OhnoWang(Hkin=Hkin, β∞=β∞, mexp=3.0) + ow = OhnoWang(Hkin=Hkin, β∞=β∞, m=3.0) σ_red_dev = dev(rand(SymmetricTensor{2,3})) ν = (3/2)*σ_red_dev/MaterialModels.vonmises(σ_red_dev) β = zero(SymmetricTensor{2,3}) From caa78758327c7a641cb716140fb5e80254978a4e Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 25 Jun 2021 08:49:42 +0200 Subject: [PATCH 55/64] Added docstrings to utility functions and cosmetic fixes --- src/utility_functions.jl | 41 ++++++++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 8 deletions(-) diff --git a/src/utility_functions.jl b/src/utility_functions.jl index 6a5cf78..7b7cb42 100644 --- a/src/utility_functions.jl +++ b/src/utility_functions.jl @@ -1,14 +1,39 @@ -# Should be fixed in ForwardDiff? +# DiffResults have a bug when being used with MArrays, which is fixed by this specialization DiffResults.DiffResult(value::MArray, derivs::Tuple{Vararg{MArray}}) = DiffResults.MutableDiffResult(value, derivs) -# Generic functions, should be defined elsewhere? -function vonmises(𝛔::SymmetricTensor{2,3}) - 𝛔_dev = dev(𝛔) - return vonmises_dev(𝛔_dev) +""" + function vonmises(σ::SymmetricTensor{2,3}) + +Calculate the von Mises effective stress for a symmetric tensor + +# Arguments +- `σ::SymmetricTensor{2,3}`: Symmetric stress tensor + +""" +function vonmises(σ::SymmetricTensor{2,3}) + σ_dev = dev(σ) + return vonmises_dev(σ_dev) end -function vonmises_dev(𝛔_dev::SymmetricTensor{2,3}) - return sqrt((3.0/2.0) * (𝛔_dev ⊡ 𝛔_dev)) +""" + function vonmises_dev(σ_dev::SymmetricTensor{2,3}) + +Calculate the von Mises effective stress for a symmetric and deviatoric tensor + +# Arguments +- `σ_dev::SymmetricTensor{2,3}`: Symmetric and deviatoric stress tensor + +""" +function vonmises_dev(σ_dev::SymmetricTensor{2,3}) + return sqrt((3.0/2.0) * (σ_dev ⊡ σ_dev)) end -macaulay(x::T) where {T} = x > 0.0 ? x : zero(T) \ No newline at end of file +""" + function macaulay(x) + +Calculate the macaulay bracket of x, ``\\langle x \\rangle`` +```math +\\langle x \\rangle = \\left\\lbrace \\begin{matrix} 0 & x\\leq 0 \\\\ x & x>0\\end{matrix} +``` +""" +macaulay(x::T) where {T} = x > zero(T) ? x : zero(T) \ No newline at end of file From 4ae37b5992220e7882156331dcc2f81842ce7a1e Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 25 Jun 2021 08:57:58 +0200 Subject: [PATCH 56/64] Added unit tests for utility functions --- test/runtests.jl | 1 + test/test_utility_functions.jl | 16 ++++++++++++++++ 2 files changed, 17 insertions(+) create mode 100644 test/test_utility_functions.jl diff --git a/test/runtests.jl b/test/runtests.jl index e170067..f159230 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using JLD2 using ForwardDiff include("test_utils.jl") +include("test_utility_functions.jl") include("test_linear_elastic.jl") include("test_plastic.jl") include("test_vonmises_plasticity.jl") diff --git a/test/test_utility_functions.jl b/test/test_utility_functions.jl new file mode 100644 index 0000000..90db9a8 --- /dev/null +++ b/test/test_utility_functions.jl @@ -0,0 +1,16 @@ +@testset "utility functions" begin + s = rand() + 1.0 + σ = SymmetricTensor{2,3}((i,j)-> i==j && i==1 ? s : 0.0) + @test vonmises(σ) ≈ s # Uniaxial stress should give the correct value + + σ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? s : 0.0) + @test vonmises(σ) ≈ √3*s # Shear stress have factor √3 + + σ = rand(SymmetricTensor{2,3}) + @test vonmises(σ) ≈ vonmises_dev(dev(σ)) # Check that the deviatoric version works as intended + + x = 1.0 + rand() + @test macaulay(-x) ≈ zero(typeof(x)) + @test macaulay(x) ≈ x + +end \ No newline at end of file From d45464de4e1719729abb8152f86719bd056fcf61 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Fri, 25 Jun 2021 09:02:51 +0200 Subject: [PATCH 57/64] Fixed test errors due to using internal functions --- test/test_utility_functions.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_utility_functions.jl b/test/test_utility_functions.jl index 90db9a8..4e24584 100644 --- a/test/test_utility_functions.jl +++ b/test/test_utility_functions.jl @@ -1,16 +1,16 @@ @testset "utility functions" begin s = rand() + 1.0 σ = SymmetricTensor{2,3}((i,j)-> i==j && i==1 ? s : 0.0) - @test vonmises(σ) ≈ s # Uniaxial stress should give the correct value + @test MaterialModels.vonmises(σ) ≈ s # Uniaxial stress should give the correct value σ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? s : 0.0) - @test vonmises(σ) ≈ √3*s # Shear stress have factor √3 + @test MaterialModels.vonmises(σ) ≈ √3*s # Shear stress have factor √3 σ = rand(SymmetricTensor{2,3}) - @test vonmises(σ) ≈ vonmises_dev(dev(σ)) # Check that the deviatoric version works as intended + @test MaterialModels.vonmises(σ) ≈ MaterialModels.vonmises_dev(dev(σ)) # Check that the deviatoric version works as intended x = 1.0 + rand() - @test macaulay(-x) ≈ zero(typeof(x)) - @test macaulay(x) ≈ x + @test MaterialModels.macaulay(-x) ≈ zero(typeof(x)) + @test MaterialModels.macaulay(x) ≈ x end \ No newline at end of file From 045209b3f9eccc0eebdd596576c0e2f445685b28 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 26 Jun 2021 18:10:06 +0200 Subject: [PATCH 58/64] Moved calculate_sigma to LinearElastic.jl (where it belongs) --- src/LinearElastic.jl | 6 ++++-- src/VonMisesPlasticity/VonMisesPlasticity.jl | 8 ++------ 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index fdf1846..9de41a2 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -50,6 +50,8 @@ Return the stress tensor and the stress tangent for the given strain ε such tha No `MaterialState` is needed for the stress computation, thus if a state is handed over to `material_response`, the same state is returned. """ function material_response(m::LinearElastic, ε::SymmetricTensor{2,3}, state::LinearElasticState=LinearElasticState(), Δt=nothing; cache=nothing, options=nothing) - σ = m.Eᵉ ⊡ ε + σ = calculate_sigma(m, ε) return σ, m.Eᵉ, state -end \ No newline at end of file +end + +calculate_sigma(m::LinearElastic, ε) = m.Eᵉ ⊡ ε diff --git a/src/VonMisesPlasticity/VonMisesPlasticity.jl b/src/VonMisesPlasticity/VonMisesPlasticity.jl index e7b5789..f4982b5 100644 --- a/src/VonMisesPlasticity/VonMisesPlasticity.jl +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -125,7 +125,7 @@ m = VonMisesPlasticity(elastic=LinearElastic(E=210.e3, ν=0.3), σ_y0=100.0, isotropic=(Voce(Hiso=-100.e3, κ∞=-100.0),Voce(Hiso=10.e3, κ∞=200.0)), kinematic=(ArmstrongFrederick(Hkin=200.e3, β∞=200.0), - OhnoWang(Hkin=1000.e3, β∞=200.0, mexp=3.0)) + OhnoWang(Hkin=1000.e3, β∞=200.0, m=3.0)) ``` """ @@ -193,8 +193,4 @@ end function calculate_plastic_strain(old::VonMisesPlasticityState, ν, Δλ) return old.ϵₚ + Δλ*ν -end - -# This function belongs in LinearElastic.jl, but to avoid conflict with ka/use_full_strain -# it is kept here until that branch is merged with main. -calculate_sigma(m::LinearElastic, ε) = m.Eᵉ ⊡ ε \ No newline at end of file +end \ No newline at end of file From 1834346d9bc159c54a1e83868bfa1833e9389cd7 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 26 Jun 2021 18:17:25 +0200 Subject: [PATCH 59/64] Removed Delta input for Plastic comparison --- test/test_vonmises_plasticity.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/test/test_vonmises_plasticity.jl b/test/test_vonmises_plasticity.jl index 79fb7f3..c9ea59b 100644 --- a/test/test_vonmises_plasticity.jl +++ b/test/test_vonmises_plasticity.jl @@ -1,6 +1,4 @@ -# Functions required as Plastic use the strain increment -get_ramp_value(::Plastic, val_init, val_max, i, num_steps) = (val_max-val_init)/num_steps -get_ramp_value(::VonMisesPlasticity, val_init, val_max, i, num_steps) = val_init + (val_max-val_init)*i/num_steps +get_ramp_value(val_init, val_max, i, num_steps) = val_init + (val_max-val_init)*i/num_steps function uniaxial_loading(m, ϵ11_init, ϵ11_max, num_steps, t_max, options=Dict{Symbol, Any}()) state = initial_material_state(m) @@ -9,7 +7,7 @@ function uniaxial_loading(m, ϵ11_init, ϵ11_max, num_steps, t_max, options=Dict ϵ = SymmetricTensor{2,1}((ϵ11_init,)) σ, dσdϵ, state = material_response(dim, m, ϵ, state, t_max/num_steps; cache=cache, options=options) for i=1:num_steps - ϵ = SymmetricTensor{2,1}((get_ramp_value(m, ϵ11_init, ϵ11_max, i, num_steps),)) + ϵ = SymmetricTensor{2,1}((get_ramp_value(ϵ11_init, ϵ11_max, i, num_steps),)) σ, dσdϵ, state = material_response(dim, m, ϵ, state, t_max/num_steps; cache=cache, options=options) end return σ, dσdϵ, state @@ -21,7 +19,7 @@ function shear_loading(m, ϵ21_init, ϵ21_max, num_steps, t_max, options=Dict{Sy ϵ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? ϵ21_init : zero(typeof(ϵ21_init))) σ, dσdϵ, state = material_response(m, ϵ, state, t_max/num_steps; cache=cache, options=options) for k=1:num_steps - ϵ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? get_ramp_value(m, ϵ21_init, ϵ21_max, k, num_steps) : zero(typeof(ϵ21_max))) + ϵ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? get_ramp_value(ϵ21_init, ϵ21_max, k, num_steps) : zero(typeof(ϵ21_max))) σ, dσdϵ, state = material_response(m, ϵ, state, t_max/num_steps; cache=cache, options=options) end return σ, dσdϵ, state From d19550f768fc9d2b3546c248c4ba8b77f4683091 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 26 Jun 2021 18:18:38 +0200 Subject: [PATCH 60/64] Updated VonMisesPlasticity jld2 due to changes in main branch get_Plastic_loading --- test/jld2_files/VonMisesPlasticity1.jld2 | Bin 15133 -> 13792 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/test/jld2_files/VonMisesPlasticity1.jld2 b/test/jld2_files/VonMisesPlasticity1.jld2 index 59a83a7eb4936326b9627bb64cd103bfc0393881..701b6a604896069e2ba2468c0e2d4bfe668c9c90 100644 GIT binary patch delta 1357 zcmbPR_8@zM43o$MQwA`wHHOg7w#i0rHefPf)I2K$5g|xV?qs%?^b~=pZi3Pc+*hzE;8%-95_1;rrz0w{eRO6!S3#CxIi8z}890TEvX zr8y-b{GVUkHf;Vc9LB`+To591RWeUva-!HVMw!X&;`S_tt~c~1D@scPxzB*yJd;Dk zP_CA+PqS1Z^e#0B&7}dMqckD(DlG{8OB+Ia=|bpPdJy`RK7_V1gwVzU-?xHY z!FYz5xtZC)3}}x6BLlW$RCD zLpO@m2q>YRAF>Wz0;f@5O@TTr{&5GR>u@y(N}SZ${T5wfs)ZgSch%u&h#d?oCO2Bi zPJV5v?gdMuY`|Q>D8R^J00G7_!jP~<6E}f{IhwdBgWTl*R#GgMC!$y;S6PKJWk^oG aXeGkcWdt#XL1bc}_+&O~VU~sb!uJ4P9{=J1 delta 1714 zcmc(gZAep57{|}OUUlNkcQ?1i5~ptFbh<4?%R9qVEE6qX0z=Ep4;B@ry^?*%*Ho}u zMwXD7#Dp>#m`S2&YDo%aU@w?WGy}z8paqI{<~c_}EkAm`{P_RRbDrnV;c)NwoSM8B z653^Z7+79cLJ+4Ds78Z|R#A!xC**8OuEc_MKjwm}t1-!`A#G|E~-=fx`sQ?TNz6Ku`iq`2@KdMgRzGblQq2cLq zET^DUN=zhTtUJqWlO8nUT|JbknPqkxhKOc=_02wt(A8QM=q1hS4iZ8rd>m!8J-re_ zgtfQxT^rn&Ct{iB4sDHCW<%Lipr6LF2TD^(5xFNXF)QlYJJJ;&6PZYZ|2+$V3}wjn zLVW3sXe~!fQ6N6`Mf|-JF}SdHkTMZCBRk(ns;oj1rOQ`Y=k z`0mYXQLi*EKZ^shYC2A8_$l+JnEt~_54k>>WJ4pw8&*tLwY!jtVjpOz9}Mw`@ktI1 zx49TOF*a|yuHvm9>@jn-;VbDO^oNG(iH0>B`EgPdpC_Z>1KvUH9t)?^nwbcYnIy0G z?e9Rl02nJWmJ$@u$BSqh*WHdSNNNf*>B#<{Ae z<>#8lA-k2amo`#>sqLC_d6CE61sRSk`Kt9T(-s@YF3Idbh9S6Z#%W;6miS8P@m+BhjE)ry!q o@A)^mJ=!E;4gRcH!QT$$Kv|_ Date: Sat, 26 Jun 2021 18:34:42 +0200 Subject: [PATCH 61/64] Cosmetic change, zero(x) instead of zero(typeof(x)) --- test/test_utility_functions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_utility_functions.jl b/test/test_utility_functions.jl index 4e24584..4faa995 100644 --- a/test/test_utility_functions.jl +++ b/test/test_utility_functions.jl @@ -10,7 +10,7 @@ @test MaterialModels.vonmises(σ) ≈ MaterialModels.vonmises_dev(dev(σ)) # Check that the deviatoric version works as intended x = 1.0 + rand() - @test MaterialModels.macaulay(-x) ≈ zero(typeof(x)) + @test MaterialModels.macaulay(-x) ≈ zero(x) @test MaterialModels.macaulay(x) ≈ x end \ No newline at end of file From 4cb88a3639b3883a665ff95660d620dc6da9563d Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 26 Jun 2021 18:58:56 +0200 Subject: [PATCH 62/64] Modify Project.toml to minimize changes wrt main --- Project.toml | 5 ++--- docs/Project.toml | 3 +-- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index c2ea052..856e27a 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ version = "0.1.0" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" @@ -14,6 +13,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" [compat] +julia = "1" DiffResults = "1.0" ForwardDiff = "0.10" NLsolve = "4.5" @@ -21,11 +21,10 @@ Reexport = "1.0" Rotations = "1.0" StaticArrays = "1.2" Tensors = "1.4" -julia = "1" [extras] -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" [targets] test = ["Test", "JLD2"] diff --git a/docs/Project.toml b/docs/Project.toml index 1fe7fea..3879473 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -MaterialModels = "d7e62f90-dad9-4fee-a9a8-264218a34fd2" [compat] -Documenter = "0.26" +Documenter = "0.26" \ No newline at end of file From 798e6b00525c0929ac3705378627f1f6edf965b0 Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 26 Jun 2021 19:04:48 +0200 Subject: [PATCH 63/64] Minimize changes to LinearElastic.jl --- src/LinearElastic.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index 9de41a2..20b4f65 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -1,5 +1,7 @@ + """ - LinearElastic(E, ν) +LinearElastic(E, ν) + Isotropic linear elasticity. # Arguments - `E::Float64`: Young's modulus From c8c0d3171efd81db7d9592727ff8f9f5ecaa221b Mon Sep 17 00:00:00 2001 From: Knut Andreas Meyer Date: Sat, 26 Jun 2021 19:07:42 +0200 Subject: [PATCH 64/64] Minimize whitespace changes --- src/LinearElastic.jl | 2 +- src/MaterialModels.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/LinearElastic.jl b/src/LinearElastic.jl index 20b4f65..5964e1b 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -1,6 +1,6 @@ """ -LinearElastic(E, ν) + LinearElastic(E, ν) Isotropic linear elasticity. # Arguments diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 1641398..953afa5 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -9,6 +9,7 @@ import DiffResults using StaticArrays + """ AbstractMaterial