Constitutive model cleanup

benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl

@@ -5,7 +5,7 @@ using BenchmarkTools
 
 function benchmark_viscous_model()
   elasto = NeoHookean3D(λ=1e6, μ=1e3)
-  visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=1e3), τ=10.)
+  visco = ViscousIncompressible(IsochoricNeoHookean3D(μ=1e3), τ=10.)
   model = GeneralizedMaxwell(elasto, visco)
   update_time_step!(model, 1e-2)
   Ψ, ∂Ψu, ∂Ψuu = model()

src/Exports.jl

@@ -33,7 +33,6 @@ end
 @publish PhysicalModels IncompressibleNeoHookean3D
 @publish PhysicalModels IncompressibleNeoHookean2D
 @publish PhysicalModels IncompressibleNeoHookean2D_CV
-@publish PhysicalModels IncompressibleNeoHookean3D_2dP
 @publish PhysicalModels VolumetricEnergy
 @publish PhysicalModels MooneyRivlin3D
 @publish PhysicalModels MooneyRivlin2D

src/PhysicalModels/MechanicalModels.jl

@@ -4,6 +4,13 @@
 # Coercive volumetric Mechanical models
 # ============================================
 
+"""
+Coercive volumetric energy term of the form:
+
+```math
+\\Psi = \\frac{1}{\\kappa} (J-1)^2
+```
+"""
 struct VolumetricEnergy <: Volumetric
   λ::Float64
   function VolumetricEnergy(; λ::Float64)
@@ -190,6 +197,13 @@ end
 # Mechanical models
 # ===================
 
+"""
+Yeoh constitutive model.
+
+```math
+\\Psi = \\sum_{i=1}^3 C_i (I_1 - 3)^i
+```
+"""
 struct Yeoh3D <: IsoElastic
   λ::Float64
   C10::Float64
@@ -625,6 +639,14 @@ I1iso(F) = det(F)^(-2 / 3) * F ⊙ F
 ∂I1iso_∂F∂Ftotal(F) = ∂I1iso_∂F∂F(F) + ∂I1iso_∂F∂J(F) ⊗ cof(F) + cof(F) ⊗ ∂I1iso_∂F∂J(F) + ∂I1iso_∂J∂J(F)*cof(F) ⊗ cof(F) + ∂I1iso_∂J(F)*×ᵢ⁴(F)
 
 
+"""
+Simplified eight-chain model by Arruda and Boyce.
+the implementation uses the first five terms of the inverse Langevin function.
+
+```math
+\\Psi = C_1 \\sum_{i=1}^{3} \\alpha_i \\beta^{i-1} (I_1^i - 3^i)
+```
+"""
 struct EightChain <: IsoElastic
   μ::Float64
   N::Float64
@@ -846,22 +868,6 @@ struct IncompressibleNeoHookean3D <: IsoElastic
 
 end
 
-function SecondPiola(obj::IncompressibleNeoHookean3D, Λ::Float64=1.0)
-  Ψ(C) = obj.μ / 2 * tr(C) * det(C)^(-1 / 3) - 3 * obj.μ / 2
-  S(C) = begin
-    detC = det(C)
-    invC = inv(C)
-    obj.μ * detC^(-1 / 3) * I3 - obj.μ / 3 * tr(C) * detC^(-1 / 3) * invC
-  end
-  ∂S∂C(C) = begin
-    detC = det(C)
-    trC = tr(C)
-    invC = inv(C)
-    IinvC = I3 ⊗ invC
-    1 / 3 * obj.μ * detC^(-1 / 3) * (4 / 3 * trC * invC ⊗ invC - (IinvC + IinvC') - trC / detC * ×ᵢ⁴(C))
-  end
-  return (Ψ, S, ∂S∂C)
-end
 
 struct IncompressibleNeoHookean2D <: IsoElastic
   λ::Float64
@@ -1004,7 +1010,6 @@ struct ARAP2D <: IsoElastic
 end
 
 
-
 struct NonlinearARAP2D <: IsoElastic
   μ::Float64
   p::Float64
@@ -1043,9 +1048,13 @@ struct NonlinearARAP2D <: IsoElastic
 end
 
 
+"""
+Neo-Hooke hyperelastic model
 
-
-
+```math
+W = \\frac{1}{2}\\mu (I_1 - 3)
+```
+"""
 struct IsochoricNeoHookean3D <: IsoElastic
   μ::Float64
 end
@@ -1055,33 +1064,21 @@ function IsochoricNeoHookean3D(; μ::Real)
 end
 
 function (obj::IsochoricNeoHookean3D)(::Float64=1.0)
-  Ψ(F) = obj.μ / 2 * (F⊙F * det(F)^(-2/3) - 3)
-  ∂Ψ∂F(F) = begin
-    μ = obj.μ
-    J = det(F)
-    Ic = F⊙F
-    obj.μ * J^(-2/3) * (F - 1/3*Ic*inv(F)')
-  end
-  ∂Ψ∂FF(F) = begin
-    μ = obj.μ
-    J = det(F)
-    Ic = F⊙F
-    invF = inv(F)
-    H = cof(F)
-    TensorValue(ForwardDiff.jacobian(∂Ψ∂F, get_array(F)))
-  end
-  return (Ψ, ∂Ψ∂F, ∂Ψ∂FF)
+  W(I) = obj.μ / 2 * (I - 3)
+  ∂W∂I(I) = obj.μ / 2
+  Ψ(F) = W(I1iso(F))
+  ∂Ψ∂F(F) = ∂W∂I(I1iso(F)) * ∂I1iso_∂Ftotal(F)
+  ∂∂Ψ∂FF(F) = ∂W∂I(I1iso(F)) * ∂I1iso_∂F∂Ftotal(F)
+  return Ψ, ∂Ψ∂F, ∂∂Ψ∂FF
 end
 
 function SecondPiola(obj::IsochoricNeoHookean3D)
   μ = obj.μ
   H(F) = cof(F)
-  Ψ(C) = μ / 2 * tr(C) * (det(C))^(-1 / 3)
+  Ψ(C) = μ / 2 * tr(C) * (det(C))^(-1 / 3) -3*μ/2
   ∂Ψ∂C(C) = μ / 2 * I3 * (det(C))^(-1 / 3)
   ∂Ψ∂dC(C) = -μ / 6 * tr(C) * (det(C))^(-4 / 3)
-  S(C) = let HC = H(C)
-    2 * (∂Ψ∂C(C) + ∂Ψ∂dC(C) * HC)
-  end
+  S(C) = 2 * (∂Ψ∂C(C) + ∂Ψ∂dC(C) * H(C))
   ∂2Ψ∂CdC(C) = -μ / 6 * I3 * (det(C))^(-4 / 3)
   ∂2Ψ∂2dC(C) = 2 * μ / 9 * tr(C) * (det(C))^(-7 / 3)
   ∂S∂C(C) = let HC = H(C)
@@ -1090,32 +1087,8 @@ function SecondPiola(obj::IsochoricNeoHookean3D)
   return (Ψ, S, ∂S∂C)
 end
 
-
-struct IncompressibleNeoHookean3D_2dP <: Mechano
-  μ::Float64
-  τ::Float64
-  Δt::Float64
-  ρ::Float64
-
-  function IncompressibleNeoHookean3D_2dP(; μ::Float64, τ::Float64, Δt::Float64, ρ::Float64=0.0)
-    new(μ, τ, Δt, ρ)
-  end
-
-  function (obj::IncompressibleNeoHookean3D_2dP)(Λ::Float64=1.0; Threshold=0.01)
-    μ = obj.μ
-    H(F) = det(F) * inv(F)'
-    Ψ(Ce) = μ / 2 * tr(Ce) * (det(Ce))^(-1 / 3)
-    ∂Ψ∂Ce(Ce) = μ / 2 * I3 * (det(Ce))^(-1 / 3)
-    ∂Ψ∂dCe(Ce) = -μ / 6 * tr(Ce) * (det(Ce))^(-4 / 3)
-    Se(Ce) = let HCe = H(Ce)
-      2 * (∂Ψ∂Ce(Ce) + ∂Ψ∂dCe(Ce) * HCe)
-    end
-    ∂2Ψ∂CedCe(Ce) = -μ / 6 * I3 * (det(Ce))^(-4 / 3)
-    ∂2Ψ∂2dCe(Ce) = 2 * μ / 9 * tr(Ce) * (det(Ce))^(-7 / 3)
-    ∂Se∂Ce(Ce) = let HCe = H(Ce)
-      2 * (∂2Ψ∂2dCe(Ce) * (HCe ⊗ HCe) + ∂2Ψ∂CedCe(Ce) ⊗ HCe + HCe ⊗ ∂2Ψ∂CedCe(Ce) + ∂Ψ∂dCe(Ce) * ×ᵢ⁴(Ce))
-    end
-
-    return (Ψ, Se, ∂Se∂Ce)
-  end
+function SecondPiola(obj::IncompressibleNeoHookean3D)
+  @warn "[DEPRECATED] SecondPiola(::IncompressibleNeoHookean3D) is deprecated. Use SecondPiola(::IsochoricNeoHookean3D) instead."
+  obj2 = IsochoricNeoHookean3D(obj.μ)
+  SecondPiola(obj2)
 end

src/PhysicalModels/PhysicalModels.jl

@@ -21,7 +21,6 @@ export IsochoricNeoHookean3D
 export IncompressibleNeoHookean3D
 export IncompressibleNeoHookean2D
 export IncompressibleNeoHookean2D_CV
-export IncompressibleNeoHookean3D_2dP
 export ARAP2D
 export ARAP2D_regularized
 export NonlinearARAP2D

test/TestConstitutiveModels/ElectroMechanicalTests.jl

@@ -119,7 +119,7 @@ end
 #     157 μs        Histogram: log(frequency) by time        391 μs <
 #  Memory estimate: 240.03 KiB, allocs estimate: 3069.
   hyper_elastic = NeoHookean3D(λ=1000., μ=10.)
-  short_term = IncompressibleNeoHookean3D(μ=5., λ=0.)
+  short_term = IsochoricNeoHookean3D(μ=5.0)
   viscous_branch1 = ViscousIncompressible(short_term, τ=6.)
   visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1)
   dielectric = IdealDielectric(ε=1.0)
@@ -141,7 +141,7 @@ end
 
 @testset "ViscoElectricModel 2-branch" begin
   hyper_elastic = NeoHookean3D(λ=1000., μ=10.)
-  short_term = IncompressibleNeoHookean3D(μ=5., λ=0.)
+  short_term = IsochoricNeoHookean3D(μ=5.)
   viscous_branch1 = ViscousIncompressible(short_term, τ=6.)
   viscous_branch2 = ViscousIncompressible(short_term, τ=60.)
   visco_elastic = GeneralizedMaxwell(hyper_elastic, viscous_branch1, viscous_branch2)

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -227,25 +227,6 @@ end
 end
 
 
-@testset "IncompressibleNeoHookean3D_2dP" begin
-  #  Memory estimate: 0 bytes, allocs estimate: 0.
-  Ce = TensorValue(0.01 + 1.0, 0.02, 0.03, 0.04, 0.05 + 1.0, 0.06, 0.07, 0.08, 0.09 + 1.0)
-  model = IncompressibleNeoHookean3D_2dP(μ=1.0, τ=1.0, Δt=1.0)
-  Ψ, Se, ∂Se = model()
-
-
-  # Se_(Ce) =2*TensorValue(ForwardDiff.gradient(x -> Ψ(x), get_array(Ce)))
-  # ∂Se_(Ce) =2*TensorValue(ForwardDiff.hessian(x -> Ψ(x), get_array(Ce)))
-
-  #  norm(Se_(Ce)) -norm(Se(Ce))
-  #  norm(∂Se_(Ce)) -norm(∂Se(Ce))
-  @test (Ψ(Ce)) == 1.5040930711508358
-  @test norm(Se(Ce)) == 0.12632997589595116
-  @test norm(∂Se(Ce)) == 2.616897862779383
-
-end
-
-
 @testset "LinearElasticity2D" begin
   #  Memory estimate: 0 bytes, allocs estimate: 0.
   model = LinearElasticity2D(λ=3.0, μ=1.0)

test/TestConstitutiveModels/ViscousModelsTests.jl

@@ -124,12 +124,12 @@ end;
 end;
 
 @testset "ViscousIncompressible" begin
-  visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ=τ1)
+  visco = ViscousIncompressible(IsochoricNeoHookean3D(μ=μ1), τ=τ1)
   test_viscous_derivatives_numerical(visco, rtolP=1e-3, rtolH=1e-3)
 end
 
 @testset "ViscousIncompressible2" begin
-  visco = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=1.0), τ=10.0)
+  visco = ViscousIncompressible(IsochoricNeoHookean3D(μ=1.0), τ=10.0)
   update_time_step!(visco, 0.1)
   Ψ, ∂Ψu, ∂Ψuu = visco()
   F    =  1e-2*TensorValue(1,2,3,4,5,6,7,8,9) + I3
@@ -156,13 +156,13 @@ end
 
 @testset "GeneralizedMaxwell NeoHookean 1-branch" begin
   hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
-  branch1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ=τ1)
+  branch1 = ViscousIncompressible(IsochoricNeoHookean3D(μ=μ1), τ=τ1)
   cons_model = GeneralizedMaxwell(hyper_elastic_model, branch1)
   test_viscous_derivatives_numerical(cons_model, rtolP=1e-3, rtolH=1e-2)
 end
 
 @testset "Dissipation ViscousIncompressible" begin
-  branch1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ1), τ=τ1)
+  branch1 = ViscousIncompressible(IsochoricNeoHookean3D(μ=μ1), τ=τ1)
   D = Dissipation(branch1, 0.1)
   F = I3
   Fn = I3
@@ -174,7 +174,7 @@ end
   Δt = 0.1
   μ = 1.0
   τ = 1.234
-  short_term = IncompressibleNeoHookean3D(λ=0., μ=μ)
+  short_term = IsochoricNeoHookean3D(μ=μ)
   branch1 = ViscousIncompressible(short_term, τ=τ)
   branch1.Δt[] = Δt
 
@@ -210,7 +210,7 @@ end
 
 @testset "GeneralizedMaxell energy at rest" begin
   hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
-  short_term = IncompressibleNeoHookean3D(λ=0., μ=μ1)
+  short_term = IsochoricNeoHookean3D(μ=μ1)
   visco_branch = ViscousIncompressible(short_term, τ=τ1)
   cons_model = GeneralizedMaxwell(hyper_elastic_model, visco_branch)
   update_time_step!(cons_model, 0.1)

test/data/StaggeredViscoElectricSimulation.jl

@@ -29,8 +29,8 @@ function staggered_visco_electric_simulation(; t_end=2, writevtk=true, verbose=t
 
   # Constitutive model
   hyper_elastic_model = NeoHookean3D(λ=1e6, μ=1.4e4)
-  viscous_branch_1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0.0, μ=5.6e4); τ=0.82)
-  viscous_branch_2 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0.0, μ=3.4e4); τ=10.7)
+  viscous_branch_1 = ViscousIncompressible(IsochoricNeoHookean3D(μ=5.6e4); τ=0.82)
+  viscous_branch_2 = ViscousIncompressible(IsochoricNeoHookean3D(μ=3.4e4); τ=10.7)
   visco_elastic_model = GeneralizedMaxwell(hyper_elastic_model, viscous_branch_1, viscous_branch_2)
   elec_model = IdealDielectric(ε=1.0)
   cons_model = ElectroMechModel(mechano=visco_elastic_model, electro=elec_model)

test/data/ViscoElasticSimulation.jl

@@ -31,9 +31,9 @@ function visco_elastic_simulation(;t_end=15, writevtk=true, verbose=true)
   μ₃ = 1.98e4  # Pa
   τ₃ = 500.0   # s
   hyper_elastic_model = NeoHookean3D(λ=λ, μ=μ)
-  viscous_branch_1 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ₁), τ=τ₁)
-  viscous_branch_2 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ₂), τ=τ₂)
-  viscous_branch_3 = ViscousIncompressible(IncompressibleNeoHookean3D(λ=0., μ=μ₃), τ=τ₃)
+  viscous_branch_1 = ViscousIncompressible(IsochoricNeoHookean3D(μ=μ₁), τ=τ₁)
+  viscous_branch_2 = ViscousIncompressible(IsochoricNeoHookean3D(μ=μ₂), τ=τ₂)
+  viscous_branch_3 = ViscousIncompressible(IsochoricNeoHookean3D(μ=μ₃), τ=τ₃)
   cons_model = GeneralizedMaxwell(hyper_elastic_model, viscous_branch_1, viscous_branch_2, viscous_branch_3)
   k=Kinematics(Mechano,Solid)