Reopen #23: Small strain plasticity

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
+```

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
+end
+
+calculate_sigma(m::LinearElastic, ε) =  m.Eᵉ ⊡ ε

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

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
+

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