Merge pull request #237 from kaipartmann/constitutive_models

Constitutive models

Author: kaipartmann

Project.toml

@@ -1,7 +1,7 @@
 name = "Peridynamics"
 uuid = "4dc47793-80f3-4232-b30e-ca78ca9d621b"
 authors = ["Kai Partmann"]
-version = "0.4.3-DEV"
+version = "0.5.0-DEV"
 
 [deps]
 AbaqusReader = "bc6b9049-e460-56d6-94b4-a597b2c0390d"

docs/src/public_api_reference.md

@@ -30,7 +30,7 @@ ZEMSilling
 ZEMWan
 LinearElastic
 NeoHooke
-MooneyRivlin
+NeoHookePenalty
 SaintVenantKirchhoff
 const_one_kernel
 linear_kernel

src/Peridynamics.jl

@@ -13,7 +13,7 @@ export BBMaterial, DHBBMaterial, OSBMaterial, CMaterial, CRMaterial, BACMaterial
        CKIMaterial, RKCMaterial, RKCRMaterial
 
 # CMaterial related types
-export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling, ZEMWan
+export LinearElastic, NeoHooke, NeoHookePenalty, SaintVenantKirchhoff, ZEMSilling, ZEMWan
 
 # Kernels
 export const_one_kernel, linear_kernel, cubic_b_spline_kernel, cubic_b_spline_kernel_norm

src/physics/ba_correspondence.jl

@@ -8,7 +8,7 @@ correspondence formulation of Chen and Spencer (2019).
 - `kernel::Function`: Kernel function used for weighting the interactions between points.
     (default: `linear_kernel`)
 - `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior.
-    (default: `LinearElastic()`)
+    (default: `SaintVenantKirchhoff()`)
 - `dmgmodel::AbstractDamageModel`: Damage model defining the fracture behavior.
     (default: `CriticalStretch()`)
 - `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. If this value is
@@ -20,7 +20,7 @@ correspondence formulation of Chen and Spencer (2019).
 
 ```julia-repl
 julia> mat = BACMaterial()
-BACMaterial{LinearElastic, typeof(linear_kernel), CriticalStretch}()
+BACMaterial{SaintVenantKirchhoff, typeof(linear_kernel), CriticalStretch}()
 ```
 ---
 
@@ -76,7 +76,7 @@ struct BACMaterial{CM,K,DM} <: AbstractBondAssociatedSystemMaterial
 end
 
 function BACMaterial(; kernel::Function=linear_kernel,
-                     model::AbstractConstitutiveModel=LinearElastic(),
+                     model::AbstractConstitutiveModel=SaintVenantKirchhoff(),
                      dmgmodel::AbstractDamageModel=CriticalStretch(),
                      maxdmg::Real=0.85)
     return BACMaterial(kernel, model, dmgmodel, maxdmg)

src/physics/constitutive_models.jl

@@ -1,18 +1,82 @@
+@doc raw"""
+    SaintVenantKirchhoff
+
+Saint-Venant-Kirchhoff constitutive model that can be specified when using a
+[`CMaterial`](@ref) and [`BACMaterial`](@ref).
+
+The strain energy density ``\Psi`` is given by
+```math
+\Psi = \frac{1}{2} \lambda \, \mathrm{tr}(\boldsymbol{E})^2 + \mu \, \mathrm{tr}(\boldsymbol{E} \cdot \boldsymbol{E}) \; ,
+```
+with the first and second Lamé parameters ``\lambda`` and ``\mu``, and the Green-Lagrange
+strain tensor
+```math
+\boldsymbol{E} = \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I}
+                              \right) \; .
+```
+
+The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+```math
+\begin{aligned}
+\boldsymbol{S} &= \lambda \, \mathrm{tr}(\boldsymbol{E}) \, \boldsymbol{I}
+                + 2 \mu \boldsymbol{E} \; , \\
+\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
+\end{aligned}
+```
+with the deformation gradient ``\boldsymbol{F}`` and the second Piola-Kirchhoff stress
+``\boldsymbol{S}``.
+
+!!! note
+    This model is equivalent to the [`LinearElastic`](@ref) model, both using the same
+    strain energy density function.
+"""
+struct SaintVenantKirchhoff <: AbstractConstitutiveModel end
+
+function first_piola_kirchhoff(::SaintVenantKirchhoff, storage::AbstractStorage,
+                               params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
+    E = 0.5 .* (F' * F - I)
+    S = params.λ * tr(E) * I + 2 * params.μ * E
+    P = F * S
+    return P
+end
+
+function strain_energy_density(::SaintVenantKirchhoff, storage::AbstractStorage,
+                               params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
+    E = 0.5 .* (F' * F - I)
+    Ψ = 0.5 * params.λ * tr(E)^2 + params.μ * tr(E * E)
+    return Ψ
+end
+
 @doc raw"""
     LinearElastic
 
 Linear elastic constitutive model that can be specified when using a [`CMaterial`](@ref) and
 [`BACMaterial`](@ref).
-The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+
+The strain energy density ``\Psi`` is given by
 ```math
-\boldsymbol{P} = \mathbb{C} : \boldsymbol{E} \; ,
+\Psi = \frac{1}{2} \lambda \, \mathrm{tr}(\boldsymbol{E})^2 + \mu \, \mathrm{tr}(\boldsymbol{E} \cdot \boldsymbol{E}) \; ,
 ```
-with the elastic stiffness tensor  ``\mathbb{C}`` and the Green-Lagrange strain tensor
-``\boldsymbol{E}`` with
+with the first and second Lamé parameters ``\lambda`` and ``\mu``, and the Green-Lagrange
+strain tensor ``\boldsymbol{E}``
 ```math
 \boldsymbol{E} = \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I}
                              \right) \; .
 ```
+
+The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+```math
+\begin{aligned}
+\boldsymbol{S} &= \mathbb{C} : \boldsymbol{E} \; , \\
+\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
+\end{aligned}
+```
+with the deformation gradient ``\boldsymbol{F}``, the elastic stiffness tensor ``\mathbb{C}``,
+and the second Piola-Kirchhoff stress ``\boldsymbol{S}``.
+
+!!! note
+    This model is equivalent to the Saint-Venant-Kirchhoff model, but uses a
+    different implementation based on the elastic stiffness tensor.
 """
 struct LinearElastic <: AbstractConstitutiveModel end
 
@@ -21,13 +85,25 @@ function first_piola_kirchhoff(::LinearElastic, storage::AbstractStorage,
     E = 0.5 .* (F' * F - I)
     Evoigt = SVector{6,Float64}(E[1,1], E[2,2], E[3,3], 2 * E[2,3], 2 * E[3,1], 2 * E[1,2])
     Cvoigt = get_hooke_matrix_voigt(params.nu, params.λ, params.μ)
-    Pvoigt = Cvoigt * Evoigt
-    P = SMatrix{3,3,Float64,9}(Pvoigt[1], Pvoigt[6], Pvoigt[5],
-                               Pvoigt[6], Pvoigt[2], Pvoigt[4],
-                               Pvoigt[5], Pvoigt[4], Pvoigt[3])
+    Svoigt = Cvoigt * Evoigt
+    # Convert second Piola-Kirchhoff stress from Voigt to tensor form
+    S = SMatrix{3,3,Float64,9}(Svoigt[1], Svoigt[6], Svoigt[5],
+                               Svoigt[6], Svoigt[2], Svoigt[4],
+                               Svoigt[5], Svoigt[4], Svoigt[3])
+    # First Piola-Kirchhoff stress: P = F * S
+    P = F * S
     return P
 end
 
+function strain_energy_density(::LinearElastic, storage::AbstractStorage,
+                               params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
+    E = 0.5 .* (F' * F - I)
+    # For energy density, use the standard form: Ψ = 0.5 * λ * tr(E)^2 + μ * tr(E*E)
+    # This is equivalent to the Saint-Venant-Kirchhoff model
+    Ψ = 0.5 * params.λ * tr(E)^2 + params.μ * tr(E * E)
+    return Ψ
+end
+
 function get_hooke_matrix_voigt(nu, λ, μ)
     a = (1 - nu) * λ / nu
     Cvoigt = SMatrix{6,6,Float64,36}(
@@ -56,12 +132,21 @@ function get_hooke_matrix(nu, λ, μ)
     return C
 end
 
-
 @doc raw"""
     NeoHooke
 
-Neo-Hookean constitutive model that can be specified when using a [`CMaterial`](@ref) and
-[`BACMaterial`](@ref).
+Compressible Neo-Hookean hyperelastic constitutive model that can be specified when using
+a [`CMaterial`](@ref) and [`BACMaterial`](@ref).
+
+The strain energy density ``\Psi`` is given by
+```math
+\Psi = \frac{1}{2} \mu \left( I_1 - 3 \right) - \mu \log(J) + \frac{1}{2} \lambda \log(J)^2 \; ,
+```
+with the first invariant ``I_1 = \mathrm{tr}(\boldsymbol{C})`` of the right Cauchy-Green
+deformation tensor ``\boldsymbol{C} = \boldsymbol{F}^{\top} \boldsymbol{F}``, the Jacobian
+``J = \mathrm{det}(\boldsymbol{F})``, and the first and second Lamé parameters ``\lambda``
+and ``\mu``.
+
 The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
 ```math
 \begin{aligned}
@@ -71,10 +156,13 @@ The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
 \boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
 \end{aligned}
 ```
-with the deformation gradient ``\boldsymbol{F}``, the right Cauchy-Green deformation tensor
-``\boldsymbol{C}``, the Jacobian ``J = \mathrm{det}(\boldsymbol{F})``, the second
-Piola-Kirchhoff stress ``\boldsymbol{S}``, and the first and second Lamé parameters
-``\lambda`` and ``\mu``.
+with the deformation gradient ``\boldsymbol{F}`` and the second Piola-Kirchhoff stress
+``\boldsymbol{S}``.
+
+# Reference
+Treloar, L. R. G. (1943). "The elasticity of a network of long-chain molecules—II."
+*Transactions of the Faraday Society*, 39, 241–246.
+DOI: [10.1039/TF9433900241](https://doi.org/10.1039/TF9433900241)
 """
 struct NeoHooke <: AbstractConstitutiveModel end
 
@@ -87,33 +175,62 @@ function first_piola_kirchhoff(::NeoHooke, storage::AbstractStorage,
     return P
 end
 
+function strain_energy_density(::NeoHooke, storage::AbstractStorage,
+                               params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
+    J = det(F)
+    C = F' * F
+    I₁ = tr(C)
+    Ψ = 0.5 * params.μ * (I₁ - 3) - params.μ * log(J) + 0.5 * params.λ * log(J)^2
+    return Ψ
+end
+
 @doc raw"""
-    MooneyRivlin
+    NeoHookePenalty
+
+Compressible Neo-Hookean hyperelastic model with a penalty-type volumetric formulation,
+suitable for modeling rubber-like and biological materials. Can be specified when using a
+[`CMaterial`](@ref) and [`BACMaterial`](@ref).
+
+The strain energy density ``\Psi`` is given by
+```math
+\Psi = \frac{1}{2} G \left( \bar{I}_1 - 3 \right) + \frac{K}{8} \left( J^2 + J^{-2} - 2 \right) \; ,
+```
+with the modified first invariant ``\bar{I}_1 = I_1 J^{-2/3}`` where
+``I_1 = \mathrm{tr}(\boldsymbol{C})`` is the first invariant of the right Cauchy-Green
+deformation tensor ``\boldsymbol{C} = \boldsymbol{F}^{\top} \boldsymbol{F}``, the Jacobian
+``J = \mathrm{det}(\boldsymbol{F})``, the shear modulus ``G``, and the bulk modulus ``K``.
 
-Mooney-Rivlin constitutive model that can be specified when using a [`CMaterial`](@ref) and
-[`BACMaterial`](@ref).
 The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
 ```math
 \begin{aligned}
 \boldsymbol{C} &= \boldsymbol{F}^{\top} \boldsymbol{F} \; , \\
 \boldsymbol{S} &= G \left( \boldsymbol{I} - \frac{1}{3} \mathrm{tr}(\boldsymbol{C})
-                           \boldsymbol{C}^{-1} \right) \cdot J^{-\frac{2}{3}}
+                           \boldsymbol{C}^{-1} \right) J^{-\frac{2}{3}}
                 + \frac{K}{4} \left( J^2 - J^{-2} \right) \boldsymbol{C}^{-1} \; , \\
 \boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
 \end{aligned}
 ```
-with the deformation gradient ``\boldsymbol{F}``, the right Cauchy-Green deformation tensor
-``\boldsymbol{C}``, the Jacobian ``J = \mathrm{det}(\boldsymbol{F})``, the second
-Piola-Kirchhoff stress ``\boldsymbol{S}``, the shear modulus ``G``, and the
-bulk modulus ``K``.
+with the deformation gradient ``\boldsymbol{F}`` and the second Piola-Kirchhoff stress
+``\boldsymbol{S}``.
+
+!!! note "Penalty-type volumetric formulation"
+    This model uses a penalty-type volumetric term ``\Psi_{\mathrm{vol}} = \frac{K}{8}(J^2 + J^{-2} - 2)``,
+    which is computationally efficient and widely used in commercial finite element codes
+    for nearly-incompressible materials. The term penalizes volume changes from the
+    reference configuration (``J = 1``).
+
+    This differs from other volumetric formulations such as:
+    - Standard Neo-Hookean (logarithmic): ``\Psi_{\mathrm{vol}} = -\mu \ln J + \frac{\lambda}{2}\ln^2 J``
+    - Simo-Miehe (polyconvex): ``\Psi_{\mathrm{vol}} = \frac{K}{4}(J^2 - 1 - 2\ln J)``
 
 # Error handling
-If the Jacobian ``J`` is smaller than the machine precision `eps()` or a `NaN`, the first
-Piola-Kirchhoff stress tensor is defined as ``\boldsymbol{P} = \boldsymbol{0}``.
+If the Jacobian ``J`` is smaller than the machine precision `eps()` or a `NaN`, the strain
+energy density and first Piola-Kirchhoff stress tensor are defined as zero:
+``\Psi = 0`` and ``\boldsymbol{P} = \boldsymbol{0}``.
 """
-struct MooneyRivlin <: AbstractConstitutiveModel end
+struct NeoHookePenalty <: AbstractConstitutiveModel end
 
-function first_piola_kirchhoff(::MooneyRivlin, storage::AbstractStorage,
+function first_piola_kirchhoff(::NeoHookePenalty, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
     J = det(F)
     J < eps() && return zero(SMatrix{3,3,T,9})
@@ -126,31 +243,13 @@ function first_piola_kirchhoff(::MooneyRivlin, storage::AbstractStorage,
     return P
 end
 
-@doc raw"""
-    SaintVenantKirchhoff
-
-Saint-Venant-Kirchhoff constitutive model that can be specified when using a
-[`CMaterial`](@ref) and [`BACMaterial`](@ref).
-The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
-```math
-\begin{aligned}
-\boldsymbol{E} &= \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I}
-                              \right) \; , \\
-\boldsymbol{S} &= \lambda \, \mathrm{tr}(\boldsymbol{E}) \, \boldsymbol{I}
-                + 2 \mu \boldsymbol{E} \; , \\
-\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
-\end{aligned}
-```
-with the deformation gradient ``\boldsymbol{F}``, the Green-Lagrange strain tensor
-``\boldsymbol{E}``, the second Piola-Kirchhoff stress ``\boldsymbol{S}``, and the first and
-second Lamé parameters ``\lambda`` and ``\mu``.
-"""
-struct SaintVenantKirchhoff <: AbstractConstitutiveModel end
-
-function first_piola_kirchhoff(::SaintVenantKirchhoff, storage::AbstractStorage,
+function strain_energy_density(::NeoHookePenalty, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
-    E = 0.5 .* (F' * F - I)
-    S = params.λ * tr(E) * I + 2 * params.μ * E
-    P = F * S
-    return P
+    J = det(F)
+    J < eps() && return zero(T)
+    isnan(J) && return zero(T)
+    C = F' * F
+    I₁ = tr(C)
+    Ψ = 0.5 * params.G * (I₁ * J^(-2 / 3) - 3) + params.K / 8 * (J^2 + J^(-2) - 2)
+    return Ψ
 end

src/physics/correspondence.jl

@@ -12,12 +12,12 @@ consistent (correspondence) formulation of non-ordinary state-based peridynamics
     - [`linear_kernel`](@ref)
     - [`cubic_b_spline_kernel`](@ref)
 - `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior. \\
-    (default: `LinearElastic()`) \\
+    (default: `SaintVenantKirchhoff()`) \\
     The following models can be used:
+    - [`SaintVenantKirchhoff`](@ref)
     - [`LinearElastic`](@ref)
     - [`NeoHooke`](@ref)
-    - [`MooneyRivlin`](@ref)
-    - [`SaintVenantKirchhoff`](@ref)
+    - [`NeoHookePenalty`](@ref)
 - `zem::AbstractZEMStabilization`: Algorithm of zero-energy mode stabilization. \\
     (default: [`ZEMSilling`](@ref)) \\
     The following algorithms can be used:
@@ -39,7 +39,7 @@ consistent (correspondence) formulation of non-ordinary state-based peridynamics
 
 ```julia-repl
 julia> mat = CMaterial()
-CMaterial{LinearElastic, ZEMSilling, typeof(linear_kernel), CriticalStretch}(maxdmg=0.85)
+CMaterial{SaintVenantKirchhoff, ZEMSilling, typeof(linear_kernel), CriticalStretch}(maxdmg=0.85)
 ```
 
 ---
@@ -119,7 +119,7 @@ struct CMaterial{CM,ZEM,K,DM} <: AbstractCorrespondenceMaterial{CM,ZEM}
 end
 
 function CMaterial(; kernel::Function=linear_kernel,
-                    model::AbstractConstitutiveModel=LinearElastic(),
+                    model::AbstractConstitutiveModel=SaintVenantKirchhoff(),
                     zem::AbstractZEMStabilization=ZEMSilling(),
                     dmgmodel::AbstractDamageModel=CriticalStretch(), maxdmg::Real=0.85)
     return CMaterial(kernel, model, zem, dmgmodel, maxdmg)