Fix viscoelasticity tutorial material parameterization

KnutAM · Aug 8, 2024 · cf7e8bf · cf7e8bf
1 parent 9c8651d
commit cf7e8bf
Showing 1 changed file with 6 additions and 6 deletions.
diff --git a/docs/src/literate_tutorials/viscoelasticity.jl b/docs/src/literate_tutorials/viscoelasticity.jl
@@ -11,7 +11,7 @@
 # ```math
 # \begin{aligned}
 # \boldsymbol{\sigma} &= \boldsymbol{\sigma}^\mathrm{vol} + \boldsymbol{\sigma}^\mathrm{dev} \\
-# \boldsymbol{\sigma}^\mathrm{vol} &= K\boldsymbol{\epsilon}^\mathrm{vol}\\
+# \boldsymbol{\sigma}^\mathrm{vol} &= 3K\boldsymbol{\epsilon}^\mathrm{vol}\\
 # \boldsymbol{\sigma}^\mathrm{dev} &= 2G_1 \boldsymbol{\epsilon}^\mathrm{dev} + 2G_2 \boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{e} \\
 # 2G_2 \boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{e} &= \eta \dot{\boldsymbol{\epsilon}}_\mathrm{v}^\mathrm{dev} \\
 # \boldsymbol{\epsilon}^\mathrm{dev} &= \boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{e} + \boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{v}
@@ -25,7 +25,7 @@
 # \begin{aligned}
 # \boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{v} &= 
 # \frac{2\Delta t*G_2*\boldsymbol{\epsilon}^\mathrm{dev} + \eta {}^\mathrm{n}\boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{v}}{\eta + 2\Delta t G_2} \\
-# \boldsymbol{\sigma} &= K \boldsymbol{\epsilon}^\mathrm{vol} + 2G_1 \boldsymbol{\epsilon}^\mathrm{dev} + 2G_2 [\boldsymbol{\epsilon}^\mathrm{dev}-\boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{v}]
+# \boldsymbol{\sigma} &= 3 K \boldsymbol{\epsilon}^\mathrm{vol} + 2G_1 \boldsymbol{\epsilon}^\mathrm{dev} + 2G_2 [\boldsymbol{\epsilon}^\mathrm{dev}-\boldsymbol{\epsilon}^\mathrm{dev}_\mathrm{v}]
 # \end{aligned}
 # ```
 # 
@@ -44,7 +44,7 @@ import CairoMakie as CM
 # finite element code. To start the material modeling, we define a material struct with all 
 # parameters. 
 Base.@kwdef struct ZenerMaterial{T}
-    K::T =5.0   # Bulk modulus 
+    K::T =5.0/3 # Bulk modulus 
     G1::T=1.0   # Shear modulus, parallel
     G2::T=50.   # Shear modulus, series 
     η::T =5.0   # Damping modulus 
@@ -68,11 +68,11 @@ function FerriteAssembly.element_residual!(re, state, ae, m::ZenerMaterial, cv::
         dΩ = getdetJdV(cv, q_point)
         ϵ = function_symmetric_gradient(cv, q_point, ae)
         ϵdev = dev(ϵ)
-        ϵv = (Δt*2*m.G2*ϵdev + m.η*old_ϵv)/(m.η + Δt*2*m.G2)
-        σ = (m.G1+m.G2)*2*ϵdev - 2*m.G2*ϵv + m.K*vol(ϵ)
+        ϵv = (Δt * 2 * m.G2 * ϵdev + m.η * old_ϵv)/(m.η + Δt * 2 * m.G2)
+        σ = (m.G1 + m.G2) * 2 * ϵdev - 2 * m.G2 * ϵv + 3 * m.K * vol(ϵ)
         for i in 1:getnbasefunctions(cv)
             δ∇N = shape_symmetric_gradient(cv, q_point, i)
-            re[i] += (δ∇N⊡σ)*dΩ
+            re[i] += (δ∇N ⊡ σ) * dΩ
         end
         ## Note that to save the state by mutation, we need to extract the value from the dual 
         ## number. Consequently, we do this before assigning to the state vector. Note that