Added simplified Aruda Boyce model

src/Exports.jl

@@ -45,6 +45,7 @@ end
 @publish PhysicalModels NonlinearMooneyRivlin_CV
 @publish PhysicalModels NonlinearIncompressibleMooneyRivlin2D_CV
 @publish PhysicalModels EightChain
+@publish PhysicalModels EightChain5Terms
 @publish PhysicalModels TransverseIsotropy3D
 @publish PhysicalModels TransverseIsotropy2D
 @publish PhysicalModels ThermalModel

src/PhysicalModels/MechanicalModels.jl

@@ -678,6 +678,37 @@ struct EightChain <: IsoElastic
 end
 
 
+struct EightChain5Terms <: IsoElastic
+  μ::Float64
+  N::Float64
+  EightChain5Terms(; μ::Float64, N::Float64) = new(μ, N)
+end
+
+function (obj::EightChain5Terms)()
+  (; μ, N) = obj
+  α = (1/2, 1/20, 11/1050, 19/7000, 519/673750)
+  β = 1 / N
+  ∂W∂I(I) = sum(i*αi*β^(i-1)*I^(i-1) for (i, αi) in enumerate(α))
+  ∂∂W∂II(I) = sum(i*(i-1)*α[i]*β^(i-1)*I^(i-2) for i in 2:length(α))
+  C1 = μ / 2 / ∂W∂I(3)
+  Ψ(F) = begin
+    I1F = F⊙F
+    C1 * sum(αi*β^(i-1)*(I1F^i - 3^i) for (i, αi) in enumerate(α))
+  end
+  ∂Ψ∂I(I) = C1 * ∂W∂I(I)
+  ∂Ψ∂F(F) = begin
+    I1F = F⊙F
+    ∂Ψ∂I(I1F) * 2F
+  end
+  ∂∂Ψ∂II(I) = C1 * ∂∂W∂II(I)
+  ∂∂Ψ∂FF(F) = begin
+    I1F = F⊙F
+    4*∂∂Ψ∂II(I1F)*(F⊗F) + 2*∂Ψ∂I(I1F)*δᵢₖδⱼₗ3D
+  end
+  return (Ψ, ∂Ψ∂F, ∂∂Ψ∂FF)
+end
+
+
 
 struct TransverseIsotropy3D <: AnisoElastic
   μ::Float64

src/PhysicalModels/PhysicalModels.jl

@@ -35,6 +35,7 @@ export NonlinearNeoHookean_CV
 export NonlinearMooneyRivlin_CV
 export NonlinearIncompressibleMooneyRivlin2D_CV
 export EightChain
+export EightChain5Terms
 export TransverseIsotropy3D
 export TransverseIsotropy2D
 export LinearElasticity3D

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -355,6 +355,14 @@ end
 end
 
 
+@testset "EightChain5Terms" begin
+  #  Memory estimate: 0 bytes, allocs estimate: 0.
+  model = EightChain5Terms(μ=μParams[1], N=μParams[2])
+  test_derivatives_3D_(model, Kinematics(Mechano, Solid), rtol=1e-13)
+  test_equilibrium_at_rest_3D(model)
+end
+
+
 @testset "TransverseIsotropy2D" begin
   #  Memory estimate: 0 bytes, allocs estimate: 0.
   ∇u = TensorValue(1.0, 2.0, 3.0, 4.0) * 1e-3