diff --git a/docs/src/index.md b/docs/src/index.md index d4d5c77..f3b22d1 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,7 +1,24 @@ # Materials +## Elastic materials ```@docs LinearElastic -material_response(m::LinearElastic, Δε::SymmetricTensor{2,3}, state::LinearElasticState) +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}, 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/LinearElastic.jl b/src/LinearElastic.jl index 96c73bc..5964e1b 100644 --- a/src/LinearElastic.jl +++ b/src/LinearElastic.jl @@ -52,6 +52,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/MaterialModels.jl b/src/MaterialModels.jl index 10d0b7e..953afa5 100644 --- a/src/MaterialModels.jl +++ b/src/MaterialModels.jl @@ -71,7 +71,9 @@ include("Plastic.jl") include("CrystalViscoPlastic/slipsystems.jl") include("CrystalViscoPlastic/CrystalViscoPlastic.jl") include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl") +include("VonMisesPlasticity/VonMisesPlasticity.jl") +include("utility_functions.jl") include("nonlinear_solver.jl") include("wrappers.jl") @@ -83,4 +85,9 @@ export LinearElastic, Plastic export LinearElasticState, PlasticState export OneD, UniaxialStrain, UniaxialStress, PlaneStrain, PlaneStress +export VonMisesPlasticity +export LinearIsotropicElasticity +export Voce, Swift +export ArmstrongFrederick, Delobelle, OhnoWang + end diff --git a/src/VonMisesPlasticity/IsotropicHardening.jl b/src/VonMisesPlasticity/IsotropicHardening.jl new file mode 100644 index 0000000..6f7777b --- /dev/null +++ b/src/VonMisesPlasticity/IsotropicHardening.jl @@ -0,0 +1,58 @@ +# Isotropic hardening +abstract type AbstractIsotropicHardening end + +""" + 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 + +function get_hardening(param::Voce, λ::Number) + param.κ∞ * (1.0 - exp(-param.Hiso * λ / param.κ∞)) +end + +""" + 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 + n::T +end +Swift(;K, λ0, n) = Swift(K, λ0, n) # Keyword argument constructor + + +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 new file mode 100644 index 0000000..16439c2 --- /dev/null +++ b/src/VonMisesPlasticity/KinematicHardening.jl @@ -0,0 +1,97 @@ +# Kinematic hardening +abstract type AbstractKinematicHardening end + +""" + 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 + +function get_evolution(param::ArmstrongFrederick, ν::SecondOrderTensor, βᵢ::SecondOrderTensor) + param.Hkin * ((2.0/3.0) * ν - βᵢ/param.β∞) +end + +""" + 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 + δ::T # Amount of Armstrong-Frederick hardening +end +Delobelle(;Hkin, β∞, δ) = Delobelle(Hkin, β∞, δ) # Keyword argument constructor + +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 + + +""" + 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 + m::T # Ohno Wang exponent +end +OhnoWang(;Hkin, β∞, m) = OhnoWang(Hkin, β∞, m) # Keyword argument constructor + +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.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 new file mode 100644 index 0000000..f4982b5 --- /dev/null +++ b/src/VonMisesPlasticity/VonMisesPlasticity.jl @@ -0,0 +1,196 @@ +include("IsotropicHardening.jl") +include("KinematicHardening.jl") + + +""" + 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 + isotropic::IsoType # Tuple of isotropic hardening definitions + kinematic::KinType # Tuple of kinematic hardening definitions +end +VonMisesPlasticity(;elastic, σ_y0, isotropic, kinematic) = VonMisesPlasticity(elastic, σ_y0, isotropic, kinematic) + +# Definition of material state +struct VonMisesPlasticityState{Nkin,T,N} <:AbstractMaterialState + ϵₚ::SymmetricTensor{2,3,T,N} + λ::T + β::NTuple{Nkin, SymmetricTensor{2,3,T,N}} +end + +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{<:VonMisesPlasticityResidual{NKin}}) where{NKin} = VonMisesPlasticityResidual{NKin} # needed for frommandel + +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 + tomandel!(v, r.β[i], offset=1+N_tens*i) + end + return v +end + +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 VonMisesPlasticityResidual(σ,λ,β) +end + +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 + +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 + state_tmp = initial_material_state(material) + ϵ = 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, ϵ) + 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: + return OnceDifferentiable(rf!, X0, X0; autodiff = :forward) +end + +""" + 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. + +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 is specified by `m.elastic` and is evaluated by giving it the elastic strain. + +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. +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}{\\partial \\boldsymbol{\\sigma}} += \\dot{\\lambda} \\boldsymbol{\\nu} +``` + +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 `m.isotropic[i]` + +Kinematic hardening is formulated as +```math +\\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}]``. + +# 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. + +# 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, m=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} + + σ_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 + 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 + update_cache!(cache, rf!) + + 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, cache.x_f; nlsolve_options...) + + if result.f_converged + x = frommandel(VonMisesPlasticityResidual{Nkin}, result.zero::MVector{7 + 6*Nkin, T}) + dRdx = cache.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ϵ, VonMisesPlasticityState(ϵₚ, x.λ, x.β) + else + error("Material model not converged. Could not find material state.") + end + end + +end + +# General residual function +function residual(x::VonMisesPlasticityResidual{NKin}, m::VonMisesPlasticity, old::VonMisesPlasticityState, ϵ) 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.λ)) + + 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) + + return VonMisesPlasticityResidual(Rσ, Rλ, Rβ) +end + +function initial_guess(m::VonMisesPlasticity, old::VonMisesPlasticityState{Nkin}, ϵ) where {Nkin} + σ_trial = calculate_sigma(m.elastic, ϵ-old.ϵₚ) + λ = old.λ + β = ntuple(i->old.β[i], Nkin) + return VonMisesPlasticityResidual(σ_trial,λ,β) +end + +function calculate_elastic_strain(old::VonMisesPlasticityState, ϵ, ν, Δλ) + return ϵ - calculate_plastic_strain(old, ν, Δλ) +end + +function calculate_plastic_strain(old::VonMisesPlasticityState, ν, Δλ) + return old.ϵₚ + Δλ*ν +end \ No newline at end of file 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) diff --git a/src/utility_functions.jl b/src/utility_functions.jl new file mode 100644 index 0000000..7b7cb42 --- /dev/null +++ b/src/utility_functions.jl @@ -0,0 +1,39 @@ +# 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) + +""" + 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}) + +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 + +""" + 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 diff --git a/test/jld2_files/VonMisesPlasticity1.jld2 b/test/jld2_files/VonMisesPlasticity1.jld2 new file mode 100644 index 0000000..701b6a6 Binary files /dev/null and b/test/jld2_files/VonMisesPlasticity1.jld2 differ diff --git a/test/runtests.jl b/test/runtests.jl index ec7e4d5..f159230 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,10 +3,13 @@ using Test using Rotations using StaticArrays 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") include("test_crystal_visco_plastic.jl") include("test_crystal_visco_plastic_red.jl") include("test_wrappers.jl") diff --git a/test/test_utility_functions.jl b/test/test_utility_functions.jl new file mode 100644 index 0000000..4faa995 --- /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 MaterialModels.vonmises(σ) ≈ s # Uniaxial stress should give the correct value + + σ = SymmetricTensor{2,3}((i,j)-> i==2 && j==1 ? s : 0.0) + @test MaterialModels.vonmises(σ) ≈ √3*s # Shear stress have factor √3 + + σ = rand(SymmetricTensor{2,3}) + @test MaterialModels.vonmises(σ) ≈ MaterialModels.vonmises_dev(dev(σ)) # Check that the deviatoric version works as intended + + x = 1.0 + rand() + @test MaterialModels.macaulay(-x) ≈ zero(x) + @test MaterialModels.macaulay(x) ≈ x + +end \ No newline at end of file diff --git a/test/test_vonmises_plasticity.jl b/test/test_vonmises_plasticity.jl new file mode 100644 index 0000000..c9ea59b --- /dev/null +++ b/test/test_vonmises_plasticity.jl @@ -0,0 +1,146 @@ +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) + cache = get_cache(m) + dim = UniaxialStress() + ϵ = 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(ϵ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 +end + +function shear_loading(m, ϵ21_init, ϵ21_max, num_steps, t_max, options=Dict{Symbol, Any}()) + state = initial_material_state(m) + cache = get_cache(m) + ϵ = 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(ϵ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 + 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 + ) + 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.λ)) + + # 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) + +end + + +@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) + + @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 + + # 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, β∞=β∞, m=3.0) + σ_red_dev = dev(rand(SymmetricTensor{2,3})) + ν = (3/2)*σ_red_dev/MaterialModels.vonmises(σ_red_dev) + β = zero(SymmetricTensor{2,3}) + + # 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 + + # 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κdλ = ForwardDiff.derivative(λarg->MaterialModels.get_hardening(voce, λarg), λ) + @test dκdλ ≈ Hiso*(1 - κ/κ∞) + + # 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