Add a function to calculate a thermal diffusivity

Author: MasanoriKanamaru

src/heat_conduction.jl

@@ -85,7 +85,8 @@ function forward_euler!(stpm::SingleTPM, Δt)
             Cₚ = stpm.thermo_params.heat_capacity[i_face]
             Δz = stpm.thermo_params.Δz
 
-            λ = k * Δt / (ρ * Cₚ * Δz^2)
+            α = thermal_diffusivity(k, ρ, Cₚ)  # Thermal diffusivity [m²/s]
+            λ = α * Δt / Δz^2
             λ ≥ 0.5 && error("The forward Euler method is unstable because λ = $λ. This should be less than 0.5.")
 
             for i_depth in 2:(n_depth-1)

src/thermo_params.jl

@@ -33,6 +33,20 @@ thermal_skin_depth(P, k, ρ, Cₚ) = @. √(4π * P * k / (ρ * Cₚ))
 thermal_inertia(k, ρ, Cₚ) = @. √(k * ρ * Cₚ)
 
 
+"""
+    thermal_diffusivity(k, ρ, Cp) -> α
+
+# Arguments
+- `k`  : Thermal conductivity [W/m/K]
+- `ρ`  : Material density [kg/m³]
+- `Cₚ` : Heat capacity [J/kg/K]
+
+# Return
+- `α` : Thermal diffusivity [m²/s]
+"""
+thermal_diffusivity(k, ρ, Cₚ) = @. k / (ρ * Cₚ)
+
+
 # ****************************************************************
 #               Struct for thermophysical parameters
 # ****************************************************************