diff --git a/Project.toml b/Project.toml index 45530be..e8430d8 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/example.md b/docs/src/example.md index 31c39d1..1a96b83 100644 --- a/docs/src/example.md +++ b/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()) ``` diff --git a/examples/CHIPFoamModel.jl b/examples/CHIPFoamModel.jl new file mode 100644 index 0000000..0863d9f --- /dev/null +++ b/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, )) diff --git a/examples/Project.toml b/examples/Project.toml index daa1781..1c7beb4 100644 --- a/examples/Project.toml +++ b/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" diff --git a/examples/TreloarModelFitting.jl b/examples/TreloarModelFitting.jl new file mode 100644 index 0000000..011c427 --- /dev/null +++ b/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 diff --git a/examples/ferrite_example.jl b/examples/ferrite_example.jl index df5215a..9569a5f 100644 --- a/examples/ferrite_example.jl +++ b/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 diff --git a/examples/gridap_example.jl b/examples/gridap_example.jl index bbdf089..f0f0d21 100644 --- a/examples/gridap_example.jl +++ b/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) diff --git a/examples/treloar.jl b/examples/treloar.jl index 52dc516..fd9a282 100644 --- a/examples/treloar.jl +++ b/examples/treloar.jl @@ -1,16 +1,17 @@ # # 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 @@ -18,7 +19,7 @@ treloar_data = HyperelasticUniaxialTest(λ₁, s₁, name = "test") # ## 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), diff --git a/ferrite_picture.png b/ferrite_picture.png new file mode 100644 index 0000000..82ffce3 Binary files /dev/null and b/ferrite_picture.png differ diff --git a/gripad_picture.png b/gripad_picture.png index b3f3209..1a70b3d 100644 Binary files a/gripad_picture.png and b/gripad_picture.png differ diff --git a/hyperelasticity.vtu b/hyperelasticity.vtu new file mode 100644 index 0000000..cf0b374 Binary files /dev/null and b/hyperelasticity.vtu differ diff --git a/results_step.vtu b/results_step.vtu index b7c64d9..7d08312 100644 Binary files a/results_step.vtu and b/results_step.vtu differ diff --git a/src/isotropic_compressible_models.jl b/src/isotropic_compressible_models.jl index ae77bc0..db22cca 100644 --- a/src/isotropic_compressible_models.jl +++ b/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 diff --git a/src/stress_functions.jl b/src/stress_functions.jl index 3c6a623..12c19c4 100644 --- a/src/stress_functions.jl +++ b/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