Merge pull request #185 from TRACER-LULab/CHIPFoam

Chip foam

Project.toml

@@ -46,9 +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"

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())
 ```
 

examples/CHIPFoamModel.jl

@@ -0,0 +1,23 @@
+using Hyperelastics
+using DifferentiationInterface
+import ForwardDiff
+# ---
+function J̄(J, (;pg, C10, φ₀, K))
+    p̃ = pg + C10 * (φ₀^(1 / 3) * (4J - 4 + 5φ₀) / (J - 1 + φ₀) - (4J - 1) * (4J + 1) / (3J^(4 / 3)))
+    Jm = exp(-p̃ / K)
+    return J / Jm
+end
+
+function dJ̄dJ1((;pg, C10, φ₀, K))
+    cbrt_φ₀ = cbrt(φ₀)
+    cbrt_φ₀2 = cbrt_φ₀^(2)
+    (-C10 - 4C10 * (cbrt_φ₀2) + K * (cbrt_φ₀2)) / (K * exp((5.0C10 - pg - 5.0C10 * (cbrt_φ₀)) / K) * (cbrt_φ₀2))
+end
+#
+using Symbolics
+@variables J pg C10 φ₀ K
+p̃ = pg + C10 * (φ₀^(1 // 3) * (4J - 4 + 5φ₀) / (J - 1 + φ₀) - (4J - 1) * (4J + 1) / (3J^(4 // 3)))
+Jm = exp(-p̃ / K)
+J̄ = J / Jm
+dJ̄dJ = Symbolics.derivative(J̄, J)
+dJ̄dJ_1 = substitute(dJ̄dJ, Dict(J => 1, ))

examples/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"

examples/TreloarModelFitting.jl

@@ -0,0 +1,41 @@
+using Hyperelastics
+using Optimization, OptimizationOptimJL
+using ComponentArrays: ComponentVector
+using DifferentiationInterface
+import ForwardDiff
+using CairoMakie, MakiePublication
+set_theme!(theme_latexfonts())
+f = Figure(size=(800, 800))
+ax = Makie.Axis(f[1, 1], xlabel="Stretch [-]", ylabel="Stress [kg/cm²]")
+treloar_data = Treloar1944Uniaxial()
+scatter!(ax,
+    getindex.(treloar_data.data.λ, 1),
+    getindex.(treloar_data.data.s, 1),
+    label="Treloar 1944 Experimental",
+    color=:black
+)
+axislegend(position=:lt)
+save("treloar_data.png", f)
+models = Dict(
+    Gent => ComponentVector(μ=240e-3, J_m=80.0),
+    EdwardVilgis => ComponentVector(Ns=0.10, Nc=0.20, α=0.001, η=0.001),
+    NeoHookean => ComponentVector(μ=200e-3),
+    Beda => ComponentVector(C1=0.1237, C2=0.0424, C3=7.84e-5, K1=0.0168, α=0.9, β=0.68, ζ=3.015)
+)
+
+sol = Dict{Any,SciMLSolution}()
+for (ψ, p_0) in models
+    HEProblem = HyperelasticProblem(ψ(), treloar_data, p_0, ad_type=AutoForwardDiff())
+    sol[ψ] = solve(HEProblem, NelderMead())
+end
+return sol # hide
+f = Figure(size=(400, 400))
+ax = Makie.Axis(f[1, 1], xlabel="Stretch [-]", ylabel="Stress [kg/cm²]")
+for (ψ, p) in sol
+    pred = predict(ψ(), treloar_data, p.u, ad_type=AutoForwardDiff())
+    lines!(ax, getindex.(pred.data.λ, 1), getindex.(pred.data.s, 1), label=string(ψ))
+end
+scatter!(ax, getindex.(treloar_data.data.λ, 1), getindex.(treloar_data.data.s, 1), label="Treloar 1944 Experimental", color=:black)
+axislegend(position=:lt)
+f
+save("treloar_data_fits.png", f) # hide

examples/ferrite_example.jl

@@ -129,8 +129,8 @@ function solve()
         x, y, z = X
         return t * Vec{3}((
             0.0,
-            L / 2 - y + (y - L / 2) * cos(θ) - (z - L / 2) * sin(θ),
-            L / 2 - z + (y - L / 2) * sin(θ) + (z - L / 2) * cos(θ)
+            L / 2 - y + (y - L / 2) * cos(-θ) - (z - L / 2) * sin(-θ),
+            L / 2 - z + (y - L / 2) * sin(-θ) + (z - L / 2) * cos(-θ)
         ))
     end
 

examples/gridap_example.jl

@@ -23,7 +23,7 @@ end
 S(F) = ∇(ψ)(F)
 
 domain = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
-partition = (5, 5, 5)
+partition = (10, 10, 10)
 model = CartesianDiscreteModel(domain, partition)
 
 labels = get_face_labeling(model)

examples/treloar.jl

@@ -1,24 +1,25 @@
 # # Package Imports
 using Hyperelastics
-using Unitful
-using ForwardDiff, FiniteDiff
+# using Unitful
+using DifferentiationInterface
+import ForwardDiff, FiniteDiff
 using Optimization
 using OptimizationOptimJL
 using ComponentArrays
-
+using CairoMakie
 # # Load the Treloar 1994 Uniaxial Dataset
 # treloar_data = [Treloar1944Uniaxial(),Treloar1944Uniaxial()]
 treloar_data = Treloar1944Uniaxial()
-λ₁ = getindex.(treloar_data.data.λ, 1)u"m/m"
-s₁ = getindex.(treloar_data.data.s, 1)u"MPa"
+λ₁ = getindex.(treloar_data.data.λ, 1)
+s₁ = getindex.(treloar_data.data.s, 1)
 treloar_data = HyperelasticUniaxialTest(λ₁, s₁, name = "test")
 
 # ## Fit the Gent Model
 # $W(\vec{\lambda}) = -\frac{\mu J_m}{2}\log{\bigg(1-\frac{I_1-3}{J_m}\bigg)}$
 #
 ## Initial guess for the parameters
 models = Dict(
-    Gent => ComponentVector((μ = 240e-3, Jₘ = 79.0)),
+    Gent => ComponentVector((μ = 240e-3, J_m = 79.0)),
     # EdwardVilgis => ComponentVector(Ns=0.10, Nc=0.20, α=0.001, η=0.001),
     # ModifiedFloryErman => ComponentVector(μ=0.24, N=50.0, κ=10.0),
     # NeoHookean => ComponentVector(μ=200e-3),

src/isotropic_compressible_models.jl

@@ -380,3 +380,102 @@ function ContinuumMechanicsBase.CauchyStressTensor(
     σ_vol = 2 * J * ∂ψ∂I3 * I
     return σ_dev + σ_vol
 end
+
+
+struct CHIPFoam <: AbstractCompressibleModel
+    Wg::Function
+    """
+    $(SIGNATURES)
+
+    CHIPFoam Model
+
+    # Model:
+    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(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::Vector,
+    (; Ĝ, K̂, Jb, pg, C10, φ₀, K, p₀, γ)
+)
+    I1 = I[1]
+    J = sqrt(I[3])
+    p̃ = pg + C10 * (φ₀^(1 / 3) * (4J - 4 + 5φ₀) / (J - 1 + φ₀) - (4J - 1) * (4J + 1) / (3J^(4 / 3)))
+    Jm = exp(-p̃ / K)
+    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))
+
+    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 = ψ.Wg(Jg, p₀, φ₀, γ)
+    return W_LB + +W_D + W_M + W_g
+end
+
+function ContinuumMechanicsBase.CauchyStressTensor(
+    ψ::CHIPFoam,
+    F::Matrix,
+    p;
+    ad_type=nothing,
+    kwargs...
+)
+    J = det(F)
+    F̄ = F ./ cbrt(J)
+    C̄ = F̄' * F̄
+    Ī₁ = tr(C̄)
+    Ī₂ = zero(Ī₁)
+    I₃ = J^2
+    ∂W∂Ī₁, _, ∂W∂I₃ = ∂ψ(ψ, [Ī₁, Ī₂, I₃], p, ad_type)
+    ∂W∂Ī₁
+    ∂W∂J = ∂W∂I₃ * 2 * J
+
+    B̄ = F̄ * F̄'
+
+    σ = 2 / J * ∂W∂Ī₁ * (B̄ - (LinearAlgebra.I * (Ī₁ / 3))) + ∂W∂J * LinearAlgebra.I
+    return σ
+end
+
+function ContinuumMechanicsBase.SecondPiolaKirchoffStressTensor(
+    ψ::CHIPFoam,
+    F::Matrix,
+    p;
+    ad_type=nothing,
+    kwargs...
+)
+    σ = CauchyStressTensor(ψ, F, p)
+    S = det(F) * inv(F) * σ' * inv(F)'
+end

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