Large strain material model with extra output.

docs/make.jl

@@ -8,7 +8,8 @@ makedocs(
         "Interface" => "interface.md",
         "Materials" => [
             "materials/LinearElastic.md",
-            "materials/Plastic.md"
+            "materials/Plastic.md",
+            "materials/LargeDefPlastic.md"
         ]
     ]
 )

docs/src/materials/LargeDefPlastic.md

@@ -0,0 +1,6 @@
+# LargeDefPlastic
+```@docs
+MatHyperElasticPlastic
+material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2, 3, T, 6}, state::MatHyperElasticPlasticState) where T
+material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2, 3, T, 6}, state::MatHyperElasticPlasticState, Δt, ::Symbol; kwargs...) where T
+```

src/FiniteStrain/largedef_plastic.jl

@@ -0,0 +1,229 @@
+export MatHyperElasticPlastic, MatHyperElasticPlasticState
+
+"""
+    MatHyperElasticPlastic(elastic_material, σ_y, H)
+
+Large strain plasticity model with kinematic hardening.
+
+# Arguments
+- `elastic_material::AbstractMaterial`: Hyperelastic material, for exampel NeoHook
+- `σ_y:` yield limit
+- `H:` hardening modulus
+
+# Reference
+J.C. Simo, 1998, Numerical analysis and simulation of plasticity 
+"""
+struct MatHyperElasticPlastic{M <: AbstractMaterial} <: AbstractMaterial
+    elastic_material::M #Elastic material, Yeoh or Neo-hook
+    σ_y::Float64		#Yield stress
+    H::Float64		    #Hardening
+    density::Float64
+end
+
+struct MatHyperElasticPlasticState <: AbstractMaterialState
+    ϵᵖ::Float64					       #Plastic strain
+    Fᵖ::Tensor{2,3,Float64,9}		   #Plastic deformation grad
+    ν::Tensor{2,3,Float64,9}	
+end
+
+struct MatHyperElasticPlasticExtras <: AbstractExtras
+    D  ::Float64 #Dissipation
+    ∂D ::SymmetricTensor{2,3,Float64,6}
+end
+
+# # # # # # #
+# Constructors
+# # # # # # #
+
+function initial_material_state(::MatHyperElasticPlastic)
+    Fᵖ = one(Tensor{2,3,Float64})
+    ν = zero(Tensor{2,3,Float64})
+    return MatHyperElasticPlasticState(0.0, Fᵖ, ν)
+end
+
+function MatHyperElasticPlastic(; 
+    elastic_material::AbstractMaterial,
+    σ_y  ::T,
+    H   ::T,
+    density ::T = NaN, 
+    ) where T
+
+    return MatHyperElasticPlastic(elastic_material, σ_y, H, density)
+end
+
+# # # # # # #
+# Drivers
+# # # # # # #
+
+function _compute_2nd_PK(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,dim,T}, state::MatHyperElasticPlasticState, compute_dissipation::Bool) where {dim,T}
+    emat, σ_y, H = (mp.elastic_material, mp.σ_y, mp.H)
+
+    I = one(C)
+    Iᵈᵉᵛ = otimesu(I,I) - 1/3*otimes(I,I)
+
+    ϵᵖ = state.ϵᵖ
+    Fᵖ = state.Fᵖ
+    ν = state.ν
+
+    dγ = 0.0
+
+    phi, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ = _hyper_yield_function(C, state.Fᵖ, state.ϵᵖ, emat, σ_y, H)
+    dFᵖdC = zero(Tensor{4,dim,T})
+    dΔγdC = zero(Tensor{2,dim,T})
+    dgdC = zero(Tensor{2,dim,T})
+    g = 0.0
+
+    if phi > 0 #Plastic step
+        #Setup newton varibles
+        TOL = 1e-9
+        newton_error, counter = TOL + 1, 0
+        dγ = 0.0
+
+        local ∂ϕ∂Mᵈᵉᵛ, ∂Mᵈᵉᵛ∂M, ∂M∂Cₑ, dϕdΔγ, ∂ϕ∂M, ∂M∂S̃
+        while true; counter +=1 
+
+            Fᵖ = state.Fᵖ - dγ*state.Fᵖ⋅state.ν
+            ϵᵖ = state.ϵᵖ + dγ
+
+            R, Mᵈᵉᵛ, Cₑ, S̃, ∂S̃∂Cₑ = _hyper_yield_function(C, Fᵖ, ϵᵖ, emat, σ_y, H) 
+
+            #Numerical diff
+            #h = 1e-6
+            #J = (yield_function(C, state.Fᵖ - (dγ+h)*state.Fᵖ⋅state.ν, state.ϵᵖ + (dγ+h), μ, λ, σ_y, H)[1] - yield_function(C, state.Fᵖ - dγ*state.Fᵖ⋅state.ν, state.ϵᵖ + dγ, μ, λ, σ_y, H)[1]) / h
+
+            ∂ϕ∂M = sqrt(3/2) * (Mᵈᵉᵛ/norm(Mᵈᵉᵛ)) #⊡ Iᵈᵉᵛ #  state.ν#
+            ∂M∂Cₑ = otimesu(I, S̃)
+            ∂Cₑ∂Fᵖ = otimesl(I, transpose(Fᵖ)⋅C) + otimesu(transpose(Fᵖ)⋅C, I)
+            ∂Fᵖ∂Δγ = -state.Fᵖ⋅state.ν
+            ∂M∂S̃ = otimesu(Cₑ, I)
+
+            dϕdΔγ = ∂ϕ∂M ⊡ ((∂M∂Cₑ + ∂M∂S̃ ⊡ ∂S̃∂Cₑ) ⊡ ∂Cₑ∂Fᵖ ⊡ ∂Fᵖ∂Δγ) - H
+
+            ddγ = dϕdΔγ\-R
+            dγ += ddγ
+
+            newton_error = norm(R)
+
+            if newton_error < TOL
+                break
+            end
+
+            if counter > 6
+                #@warn("Could not find equilibrium in material")
+                break
+            end
+        end
+
+        ∂Cₑ∂C = otimesu(transpose(Fᵖ),transpose(Fᵖ))
+        ∂ϕ∂C = ∂ϕ∂M ⊡ ((∂M∂Cₑ + ∂M∂S̃ ⊡ ∂S̃∂Cₑ) ⊡ ∂Cₑ∂C) 
+        dΔγdC = -inv(dϕdΔγ) * ∂ϕ∂C
+        dFᵖdC = -(state.Fᵖ⋅state.ν) ⊗ dΔγdC
+
+    end
+
+
+    ν = norm(Mᵈᵉᵛ)==0.0 ? zero(Mᵈᵉᵛ) : sqrt(3/2)*Mᵈᵉᵛ/(norm(Mᵈᵉᵛ))
+    S = (Fᵖ ⋅ S̃ ⋅ transpose(Fᵖ))
+
+    ∂S∂Fᵖ = otimesu(I, (Fᵖ ⋅ S̃)) + otimesl((Fᵖ ⋅ S̃), I)
+    ∂S∂S̃ = otimesu(Fᵖ,Fᵖ)
+    ∂Cₑ∂Fᵖ = otimesl(I, transpose(Fᵖ) ⋅ C) + otimesu(transpose(Fᵖ) ⋅ C, I)
+    ∂Cₑ∂C = otimesu(transpose(Fᵖ),transpose(Fᵖ))
+    dCₑdC = ∂Cₑ∂C + ∂Cₑ∂Fᵖ ⊡ dFᵖdC
+    dSdC = ∂S∂Fᵖ ⊡ dFᵖdC + ∂S∂S̃ ⊡ ∂S̃∂Cₑ ⊡ dCₑdC
+
+    g = 0.0 
+    dgdC = zero(SymmetricTensor{2,3,T,6})
+    if compute_dissipation
+        if phi > 0.0
+            M = Cᵉ ⋅ S̃
+            g, dgdC = _compute_dissipation(M, ν, dγ, dCₑdC, ∂S̃∂Cₑ, dΔγdC, ∂M∂Cₑ, ∂M∂S̃, Iᵈᵉᵛ, I)
+        end
+    end
+
+    return symmetric(S), dSdC, ϵᵖ, ν, Fᵖ, g, dgdC
+end
+
+function _hyper_yield_function(C::SymmetricTensor, Fᵖ::Tensor, ϵᵖ::T, emat::AbstractMaterial, σ_y::Float64, H::Float64) where T
+    Cᵉ = symmetric(transpose(Fᵖ)⋅C⋅Fᵖ)
+
+    S̃, ∂S̃∂Cₑ = material_response(emat, Cᵉ)
+
+    #Mandel stress
+    Mbar = Cᵉ⋅S̃
+    Mᵈᵉᵛ = dev(Mbar)
+
+    yield_func = sqrt(3/2)*norm(Mᵈᵉᵛ) - (σ_y + H*ϵᵖ)
+
+    return yield_func, Mᵈᵉᵛ, Cᵉ, S̃, ∂S̃∂Cₑ
+end
+
+function _compute_dissipation(M, ν, Δγ, dCₑdC, ∂S̃∂Cₑ, dΔγdC, ∂M∂Cₑ, ∂M∂S̃, Iᵈᵉᵛ, I)
+
+    Mᵈᵉᵛ = dev(M)
+
+    ∂g∂M = ν * Δγ
+    ∂g∂ν = M * Δγ
+    ∂g∂Δγ = M ⊡ ν
+
+    ∂ν∂Mᵈᵉᵛ = √(3/2) * (1/norm(Mᵈᵉᵛ)) * (Iᵈᵉᵛ - (Mᵈᵉᵛ ⊗ Mᵈᵉᵛ)/(norm(Mᵈᵉᵛ)^2) )
+    ∂Mᵈᵉᵛ∂M = Iᵈᵉᵛ
+    dMdC = (∂M∂Cₑ + ∂M∂S̃ ⊡ ∂S̃∂Cₑ) ⊡ dCₑdC
+
+    g = M ⊡ ν*Δγ
+    dgdC = (∂g∂M + ∂g∂ν ⊡ ∂ν∂Mᵈᵉᵛ ⊡ ∂Mᵈᵉᵛ∂M) ⊡ dMdC  +  ∂g∂Δγ * dΔγdC 
+
+    return g, dgdC
+end
+
+"""
+    material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0;  <keyword arguments>)
+
+The material model assumes a multiplicative split of the deformation gradient
+
+```math
+\\boldsymbol{F} = \\boldsymbol{F}_e \\boldsymbol{F}_p
+```
+The yield function is of von Mises type is formulated in terms of the Mandel stress, \$\\boldsymbol{M}\$
+
+```math
+\\phi = \\sqrt{\\frac{3}{2}} \\left| \\boldsymbol{M}_{\\text{dev}} \\right| - \\sigma_y - H \\varepsilon_p
+```
+
+The evolution equation of associative type for the plastic velocity gradient is
+
+```math
+\\boldsymbol{L}_p = \\dot \\boldsymbol{F}_p \\boldsymbol{F}_p^{-1} = \\dot \\gamma \\frac{\\partial \\phi}{\\partial \\boldsymbol{M}} = \\dot \\gamma \\boldsymbol{\\nu}
+```
+
+As an simplification for the local backward Euler integration in the material routine, it is assumed that
+\$ \\boldsymbol{\\nu} = ^n\\boldsymbol{\\nu} \$. This reduces the system of equations to only one equation 
+(namely the yield function).
+
+"""
+function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt=0.0; 
+    cache=nothing, options=nothing) where T
+
+    S, ∂S∂C, ϵᵖ, ν, Fᵖ, _, _ = _compute_2nd_PK(mp, C, state, false)
+
+    return S, ∂S∂C, MatHyperElasticPlasticState(ϵᵖ, Fᵖ, ν)
+end
+
+"""
+    material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt, extras::Symbol;  <keyword arguments>)
+
+In addition to returning the stress, tangent and state, it also return extra data related to the amount 
+of dissipation energy, according to:
+
+```math
+D = \\dot \\boldsymbol{M} : \\boldsymbol{L}_p (\\geq 0)
+```
+
+"""
+function material_response(mp::MatHyperElasticPlastic, C::SymmetricTensor{2,3,T,6}, state::MatHyperElasticPlasticState, Δt, ::Symbol
+    ; cache=nothing, options=nothing) where T 
+
+    S, ∂S∂C, ϵᵖ, ν, Fᵖ, g, dgdC = _compute_2nd_PK(mp, C, state, true)
+
+    return S, ∂S∂C, MatHyperElasticPlasticState(ϵᵖ, Fᵖ, ν), MatHyperElasticPlasticExtras(g, dgdC)
+end

src/MaterialModels.jl

@@ -27,6 +27,14 @@ Store state variables here. For now, this should **not** be mutable, a new objec
 abstract type AbstractMaterialState end
 abstract type AbstractResiduals end
 
+"""
+    AbstractExtras
+
+Outputs from the material routine in addition to what is usual or strictly necessary (stress, tangent and state)
+"""
+abstract type AbstractExtras end
+struct EmptyExtras <: AbstractExtras end
+
 """
     material_response(m::AbstractMaterial, Δε::SymmetricTensor{2,3}, state::AbstractMaterialState, Δt; cache, options)
 
@@ -40,6 +48,16 @@ This function signature must be the same for all material models, even if they d
 """
 function material_response end
 
+"""
+    material_response(m::AbstractMaterial, Δε::SymmetricTensor{2,3}, state::AbstractMaterialState, Δt, ::Symbol; cache, options)
+
+Fallback for materials that has no extra outputs
+"""
+function material_response(m::AbstractMaterial, strain, state::AbstractMaterialState, Δt, ::Symbol; cache, options)
+    stress, tangent, state = material_response(m, strain, state, Δt; cache=cache, options=options)
+    return stress, tangent, state, EmptyExtras()
+end
+
 """
     initial_material_state(::AbstractMaterial)
 
@@ -68,6 +86,7 @@ function update_cache! end
 
 include("LinearElastic.jl")
 include("Plastic.jl")
+include("FiniteStrain/largedef_plastic.jl")
 include("CrystalViscoPlastic/slipsystems.jl")
 include("CrystalViscoPlastic/CrystalViscoPlastic.jl")
 include("CrystalViscoPlastic/CrystalViscoPlasticRed.jl")
@@ -85,6 +104,7 @@ export material_response
 export AbstractMaterial
 export LinearElastic, Plastic
 export LinearElasticState, PlasticState
+export MatHyperElasticPlastic, MatHyperElasticPlasticState
 export OneD, UniaxialStrain, UniaxialStress, PlaneStrain, PlaneStress
 
 export NeoHook, Yeoh, StVenant

test/runtests.jl

@@ -12,4 +12,5 @@ include("test_crystal_visco_plastic_red.jl")
 include("test_wrappers.jl")
 include("test_neohook.jl")
 include("test_yeoh.jl")
-include("test_stvenant.jl")
+include("test_stvenant.jl")
+include("test_largedef_plastic.jl")

test/test_largedef_plastic.jl

@@ -0,0 +1,73 @@
+@testset "MatHyperElasticPlastic" begin
+    # constructor
+    E=200e3
+    ν=0.3
+    neohook = NeoHook(
+        μ = E / (2(1 + ν)),
+        λ = (E * ν) / ((1 + ν) * (1 - 2ν))
+    )
+
+    m = MatHyperElasticPlastic(
+        elastic_material = neohook,
+        σ_y = 200.0, 
+        H = 50.0,
+    )
+
+    # initial state
+    state = initial_material_state(m)
+
+    # elastic branch
+    ε11_yield = m.σ_y / E 
+    Δt = 0.0
+
+    F = Tensor{2,3}((0.5ε11_yield + 1, 0.0, 0.0, 0.0, (-0.5ε11_yield*ν + 1.0), 0.0, 0.0, 0.0, (-0.5ε11_yield*ν + 1.0)))
+    C = tdot(F)
+    S_neo, ∂S_neo, _ = material_response(neohook, C)
+    S, ∂S, temp_state = material_response(m, C, state)
+    _,_,_,extras = material_response(m, C, state, Δt , :extras)
+
+    @test S_neo ≈ S
+    @test ∂S == ∂S_neo
+    @test temp_state.ϵᵖ == 0.0
+    @test temp_state.Fᵖ == one(Tensor{2,3,Float64})
+    @test extras.D == 0.0 # No dissipation
+
+    # Plastic branch
+    F = Tensor{2,3}((1.1ε11_yield + 1, 0.0, 0.0, 0.0, (-1.1ε11_yield*ν + 1.0), 0.0, 0.0, 0.0, (-1.1ε11_yield*ν + 1.0)))
+    C = tdot(F)
+
+    S, ∂S, temp_state, extras = material_response(m, C, state, Δt, :extras)
+
+    @test temp_state.ϵᵖ > 0.0
+    @test extras.D > 0.0
+
+end
+
+
+function get_MatHyperElasticPlastic_loading()
+    strain1 = range(0.0,  0.002, length=5)
+    strain2 = range(0.002, 0.001, length=5)
+    strain3 = range(0.001, 0.007, length=5)
+
+    _x = [strain1[1:end]..., strain2[1:end]..., strain3[1:end]...]
+    _F = [Tensor{2,3}((x + 1.0, x/10, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) for x in _x]
+    C = [tdot(F) for F in _F]
+
+    return C
+end  
+
+@testset "MatHyperElasticPlastic jld2" begin
+
+    E=200e3
+    ν=0.3
+    m = MatHyperElasticPlastic(
+        elastic_material = MaterialModels.NeoHook(
+            μ = E / (2(1 + ν)),
+            λ = (E * ν) / ((1 + ν) * (1 - 2ν))
+        ),
+        σ_y = 200.0, 
+        H = 50.0,
+    )
+    loading = get_MatHyperElasticPlastic_loading()
+    check_jld2(m, loading, "MatHyperElasticPlastic1")#, debug_print=true, OVERWRITE_JLD2=true)
+end