Updated the stres function definitions

Project.toml

@@ -1,7 +1,7 @@
 name = "Hyperelastics"
 uuid = "c20c605d-adbe-46e4-9db6-e2b04c3ba023"
 authors = ["Carson Farmer"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -29,7 +29,6 @@ DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [extensions]
@@ -47,6 +46,12 @@ Integrals = "4"
 InverseLangevinApproximations = "0.2"
 LabelledArrays = "1"
 LossFunctions = "0.11"
+<<<<<<< HEAD
+=======
+MakieCore = "0.8"
+NaNMath = "1"
+Optimization = "3"
+>>>>>>> 602d4c70ec14c8c4d41dfcc0d828b3ebcda2db1a
 PackageExtensionCompat = "1"
 QuadGK = "2"
 RecursiveArrayTools = "3"
@@ -56,7 +61,6 @@ SpecialFunctions = "2"
 Statistics = "1"
 StructArrays = "0.6"
 Unitful = "1"
-julia = "1.9"
 
 [extras]
 CPUSummary = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"

docs/Project.toml

@@ -6,18 +6,16 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Hyperelastics = "c20c605d-adbe-46e4-9db6-e2b04c3ba023"
-MakiePublication = "dde8697e-0d61-460d-88dd-856f66710dd1"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 
 [compat]
-CairoMakie = "0.10, 0.11"
-ComponentArrays = "0.15"
-DataInterpolations = "4"
-Documenter = "1"
-DocumenterCitations = "1"
-ForwardDiff = "0.10"
+DataInterpolations = "5"
+OptimizationOptimJL = "0.3"
 Hyperelastics = "0.1"
-MakiePublication = "0.3"
+ForwardDiff = "0.10"
+DocumenterCitations = "1"
+CairoMakie = "0.12"
+Documenter = "1"
 Optimization = "3"
-OptimizationOptimJL = "0.2"
+ComponentArrays = "0.15"

docs/src/example.md

@@ -4,9 +4,15 @@
 using Hyperelastics
 using Optimization, OptimizationOptimJL
 using ComponentArrays: ComponentVector
+<<<<<<< HEAD
 using DifferentiationInterface
 import ForwardDiff
 using CairoMakie, MakiePublication
+=======
+import ForwardDiff
+using DifferentiationInterface
+using CairoMakie
+>>>>>>> 602d4c70ec14c8c4d41dfcc0d828b3ebcda2db1a
 set_theme!(theme_latexfonts())
 ```
 

src/isotropic_compressible_models.jl

@@ -383,32 +383,44 @@ end
 
 
 struct CHIPFoam <: AbstractCompressibleModel
-
+    Wg::Function
     """
     $(SIGNATURES)
 
     CHIPFoam Model
+
     # Model:
-    ```math
-    \\psi =
+    Refer to: Lewis M. A robust, compressible, hyperelastic constitutive model for the mechanical response of foamed rubber. Technische Mechanik-European Journal of Engineering Mechanics. 2016;36(1-2):88-101.
 
-    ```
     # Arguments:
-
+    - isothermal
+        - Boolean to determine if the model is isothermal or adiabatic
 
     # Parameters:
-    - μ
-    - α
-    - β
+    - Ĝ
+    - K̂
+    - Jb
+    - pg
+    - C10
+    - φ₀,
+    - K
+    - p₀
+    - γ
 
     """
-    function CHIPFoam()
-        new()
+    function CHIPFoam(isothermal::Bool=true)
+        if isothermal
+            Wg(Jg, p₀, φ₀, γ) = p₀ * φ₀ * (Jg - log(Jg) - 1)
+        else
+            Wg(Jg, p₀, φ₀, γ) = p₀ * φ₀ * (Jg - 1 / (γ - 1) * (γ - Jg^(1 - γ)))
+        end
+        new(Wg)
     end
 end
+
 function ContinuumMechanicsBase.StrainEnergyDensity(
-    ::CHIPFoam,
-    I,
+    ψ::CHIPFoam,
+    I::Vector,
     (; Ĝ, K̂, Jb, pg, C10, φ₀, K, p₀, γ)
 )
     I1 = I[1]
@@ -418,33 +430,33 @@ function ContinuumMechanicsBase.StrainEnergyDensity(
     J̄ = J / Jm
     Jg = Jm * (J̄ - 1 + φ₀) / (φ₀)
     f = (2J - 1) / (cbrt(J)) + (2 - 2J + φ₀) * cbrt((φ₀) / (J - (1 - φ₀)))
+
     cbrt_φ₀ = cbrt(φ₀)
     cbrt_φ₀2 = cbrt_φ₀^(2)
-    dJ̄dJ_1 = (-C10 - 4C10 * (cbrt_φ₀2) + K * (cbrt_φ₀2)) / (K * exp((5.0C10 - pg - 5.0C10 * (cbrt_φ₀)) / K) * (cbrt_φ₀2))q
+    dJ̄dJ_1 = (-C10 - 4C10 * (cbrt_φ₀2) + K * (cbrt_φ₀2)) / (K * exp((5.0C10 - pg - 5.0C10 * (cbrt_φ₀)) / K) * (cbrt_φ₀2))
+
     W_LB = Ĝ / 2 * (I1 - 3) + K̂ * (
         (Jb - 1) * (J - (Jb + 1) / 2) + ((J - Jb) >= 0) * ((J - 1)^2 / 2 - (Jb - 1) * (J - (Jb + 1) / 2))
     )
 
     W_D = C10 * (Jm * (I1 * f - 3 * (1 - φ₀)) - J * (I1 - 3) * (1 - φ₀) * dJ̄dJ_1)
 
     W_M = (1 - φ₀) * K * (Jm * log(Jm) - Jm + 1)
-    W_g = p₀ * φ₀ * (Jg - log(Jg) - 1)
-    # Or
-    W_g = p₀ * φ₀ * (Jb - 1 / (γ - 1) * (γ - Jg^(1 - γ)))
+    W_g = ψ.Wg(Jg, p₀, φ₀, γ)
     return W_LB + +W_D + W_M + W_g
 end
 
 function ContinuumMechanicsBase.CauchyStressTensor(
     ψ::CHIPFoam,
-    F,
+    F::Matrix,
     p;
     ad_type=nothing,
     kwargs...
 )
     J = det(F)
     F̄ = F ./ cbrt(J)
-    C = F̄' * F̄
-    Ī₁ = tr(F̄)
+    C̄ = F̄' * F̄
+    Ī₁ = tr(C̄)
     Ī₂ = zero(Ī₁)
     I₃ = J^2
     ∂W∂Ī₁, _, ∂W∂I₃ = ∂ψ(ψ, [Ī₁, Ī₂, I₃], p, ad_type)
@@ -459,7 +471,7 @@ end
 
 function ContinuumMechanicsBase.SecondPiolaKirchoffStressTensor(
     ψ::CHIPFoam,
-    F,
+    F::Matrix,
     p;
     ad_type=nothing,
     kwargs...

src/stress_functions.jl

@@ -115,9 +115,10 @@ function ContinuumMechanicsBase.SecondPiolaKirchoffStressTensor(
     ad_type = nothing,
     kwargs...,
 ) where {T<:InvariantForm,R}
-    I1 = I₁(F)
-    I2 = I₂(F)
-    I3 = I₃(F)
+    C = F'*F
+    I1 = I₁(C)
+    I2 = I₂(C)
+    I3 = I₃(C)
     ∂ψ∂I = ∂ψ(ψ, [I1, I2, I3], p, ad_type; kwargs...)
     S = 2∂ψ∂I[1] * F' + 2∂ψ∂I[2] * (I1 * F' + F' * F * F') + 2I3 * ∂ψ∂I[3] * inv(F)
     return S
@@ -190,9 +191,10 @@ function ContinuumMechanicsBase.CauchyStressTensor(
     ad_type,
     kwargs...,
 ) where {T<:InvariantForm,S}
-    I1 = I₁(F)
-    I2 = I₂(F)
-    I3 = I₃(F)
+    C = F'*F
+    I1 = I₁(C)
+    I2 = I₂(C)
+    I3 = I₃(C)
     J = sqrt(I3)
     B = F * F'
     ∂ψ∂I = ∂ψ(ψ, [I1, I2, I3], p, ad_type; kwargs...)
@@ -201,3 +203,33 @@ function ContinuumMechanicsBase.CauchyStressTensor(
         2 * J * ∂ψ∂I[3] * I
     return σ
 end
+
+function ContinuumMechanicsBase.MaterialElasticStiffnessTensor(
+ ψ::Hyperelastics.AbstractHyperelasticModel{T},
+    F::Matrix{S},
+    p;
+    ad_type,
+    kwargs...,
+) where {T<:InvariantForm,S}
+    I1 = I₁(F)
+    I2 = I₂(F)
+    I3 = I₃(F)
+    J = sqrt(I3)
+    B = F * F'
+    ∂²ψ∂I² = ∂²ψ(ψ, [I1, I2, I3], p, ad_type; kwargs...)
+    C = zeros(3, 3, 3, 3)
+    for i in 1:3
+        for j in 1:3
+            for k in 1:3
+                for l in 1:3
+                    C[i, j, k, l] =
+                        2 * inv(J) * ∂²ψ∂I²[1][i, j, k, l] +
+                        2 * I1 * ∂²ψ∂I²[2][i, j, k, l] +
+                        2 * I2 * ∂²ψ∂I²[3][i, j, k, l] +
+                        2 * J * ∂²ψ∂I²[4][i, j, k, l]
+                end
+            end
+        end
+    end
+    return C
+end