diff --git a/benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl b/benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl index 1569f20..0e84d33 100644 --- a/benchmark/ConstitutiveModelsBenchmark/ViscousModelsBenchmarks.jl +++ b/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() diff --git a/src/Exports.jl b/src/Exports.jl index 78170f3..a2ca1df 100644 --- a/src/Exports.jl +++ b/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 diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl index 850f27f..4176fac 100644 --- a/src/PhysicalModels/MechanicalModels.jl +++ b/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 diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl index 31eb4f7..455a776 100644 --- a/src/PhysicalModels/PhysicalModels.jl +++ b/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 diff --git a/test/TestConstitutiveModels/ElectroMechanicalTests.jl b/test/TestConstitutiveModels/ElectroMechanicalTests.jl index 815e467..e19071c 100644 --- a/test/TestConstitutiveModels/ElectroMechanicalTests.jl +++ b/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) diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl index 0f428b8..36c87e1 100644 --- a/test/TestConstitutiveModels/PhysicalModelTests.jl +++ b/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) diff --git a/test/TestConstitutiveModels/ViscousModelsTests.jl b/test/TestConstitutiveModels/ViscousModelsTests.jl index cd50353..fe1c6c0 100644 --- a/test/TestConstitutiveModels/ViscousModelsTests.jl +++ b/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) diff --git a/test/data/StaggeredViscoElectricSimulation.jl b/test/data/StaggeredViscoElectricSimulation.jl index 74ea0cd..94f228f 100644 --- a/test/data/StaggeredViscoElectricSimulation.jl +++ b/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) diff --git a/test/data/ViscoElasticSimulation.jl b/test/data/ViscoElasticSimulation.jl index d14ea16..0aa2a4c 100644 --- a/test/data/ViscoElasticSimulation.jl +++ b/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)