Add dimensionless number calculations and bump version

Project.toml

@@ -1,7 +1,7 @@
 name = "ShockwaveProperties"
 uuid = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50"
 authors = ["Alex Fleming <[email protected]> and contributors"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/ShockwaveProperties.jl

@@ -19,6 +19,10 @@ export temperature, pressure, speed_of_sound
 export specific_static_enthalpy, specific_static_internal_energy
 export static_internal_energy_density
 export total_internal_energy_density, total_enthalpy_density
+export kinematic_viscosity, thermal_diffusivity
+
+# dimensionless numbers
+export reynolds_number, froude_number, fourier_number, prandtl_number
 
 # representations
 export ConservedProps, PrimitiveProps
@@ -36,7 +40,7 @@ export primitive_state_vector, primitive_state_behind
 include("cpg.jl")
 
 # At 300K
-const DRY_AIR = CaloricallyPerfectGas(1004.9u"J/kg/K", 717.8u"J/kg/K", 0.0289647u"kg/mol")
+const DRY_AIR = CaloricallyPerfectGas(1004.9u"J/kg/K", 717.8u"J/kg/K", 0.0289647u"kg/mol", 0.02624u"W/m/K", 1.846e-5u"kg/m/s")
 
 include("state_variables.jl")
 include("billig.jl")

src/cpg.jl

@@ -1,5 +1,7 @@
+using Unitful: Length, Time, MolarMass, DynamicViscosity
 @derived_dimension HeatCapacity 𝐋^2 * 𝐓^-2 * 𝚯^-1 true
-@derived_dimension MolarMass 𝐌 * 𝐍^-1 true
+@derived_dimension ThermalConductivity 𝐌 * 𝐋 * 𝐓^-3 * 𝚯^-1 true
+# @derived_dimension MolarMass 𝐌 * 𝐍^-1 true
 @derived_dimension MomentumDensity 𝐌 * 𝐋^-2 * 𝐓^-1 true
 @derived_dimension SpecificEnergy 𝐋^2 * 𝐓^-2 true
 @derived_dimension EnergyDensity 𝐌 * 𝐋^-1 * 𝐓^-2 true
@@ -11,6 +13,8 @@ const _units_ρ = u"kg/m^3"
 const _units_v = u"m/s"
 const _units_T = u"K"
 const _units_P = u"Pa"
+const _units_k = u"W/m/K"
+const _units_μ = u"Pa*s"
 
 const _units_ρv = _units_ρ * _units_v
 const _units_int_e = _units_cvcp * _units_T
@@ -25,33 +29,61 @@ Provides the properties of a calorically perfect gas (or mixture of gases).
  - ``γ``: Heat capacity ratio
  - ``R``: Specific gas constant
 """
-struct CaloricallyPerfectGas{U1<:HeatCapacity,U2<:MolarMass}
+struct CaloricallyPerfectGas{
+    U1<:HeatCapacity,
+    U2<:MolarMass,
+    U3<:ThermalConductivity,
+    U4<:DynamicViscosity,
+}
     c_p::U1
     c_v::U1
     ℳ::U2
 
     γ::Float64
     R::U1
+
+    k::U3
+    μ::U4
 end
 
-function CaloricallyPerfectGas(c_p, c_v, ℳ)
+function CaloricallyPerfectGas(c_p, c_v, ℳ, k, μ)
     q_cp = Quantity(c_p, _units_cvcp)
     q_cv = Quantity(c_v, _units_cvcp)
     q_ℳ = Quantity(ℳ, _units_ℳ)
     q_R = q_cp - q_cv
-    return CaloricallyPerfectGas{typeof(q_cp),typeof(q_ℳ)}(q_cp, q_cv, q_ℳ, c_p / c_v, q_R)
+    q_k = Quantity(k, _units_k)
+    q_μ = Quantity(k, _units_μ)
+    return CaloricallyPerfectGas{typeof(q_cp),typeof(q_ℳ),typeof(q_k),typeof(q_μ)}(
+        q_cp,
+        q_cv,
+        q_ℳ,
+        c_p / c_v,
+        q_R,
+        q_k,
+        q_μ,
+    )
 end
 
 function CaloricallyPerfectGas(
     c_p::T,
     c_v::T,
     ℳ::U,
-) where {T<:HeatCapacity,U<:Unitful.MolarMass}
+    k::V,
+    μ::W,
+) where {T<:HeatCapacity,U<:MolarMass,V<:ThermalConductivity,W<:DynamicViscosity}
     R = c_p - c_v
     γ = c_p / c_v
-    return CaloricallyPerfectGas{T,U}(c_p, c_v, ℳ, γ, R)
+    return CaloricallyPerfectGas{T,U,V,W}(c_p, c_v, ℳ, γ, R, k, μ)
 end
 
+## DIMENSIONLESS NUMBERS THAT DEPEND ON GAS PROPERTIES
+
+function prandtl_number(gas::CaloricallyPerfectGas)
+    return uconvert(Unitful.NoUnits, gas.c_p * gas.μ / gas.k)
+end
+
+## STATES
+
 """
     PrimitiveProps{N, DTYPE, Q1<:Density, Q2<:Temperature}
 

src/state_variables.jl

@@ -1,3 +1,5 @@
+using Unitful: Length, Acceleration
+
 ## DENSITY
 
 """
@@ -270,4 +272,45 @@ end
 
 function speed_of_sound(u::ConservedProps, gas::CaloricallyPerfectGas)
     return speed_of_sound(density(u), pressure(u, gas), gas)
-end
+end
+
+## KINEMATIC VISCOSITY
+
+function kinematic_viscosity(s, gas::CaloricallyPerfectGas)
+    return gas.μ / density(s)
+end
+
+## THERMAL DIFFUSIVITY
+
+function thermal_diffusivity(s, gas::CaloricallyPerfectGas)
+    return gas.k / (gas.c_p * density(s))
+end
+
+## DIMENSIONLESS NUMBERS THAT DEPEND ON FLUID STATES
+
+function reynolds_number(s, gas::CaloricallyPerfectGas, D::T) where {T<:Length}
+    return uconvert(Unitful.NoUnits, density(s) * norm(velocity(s, gas)) * D / gas.μ)
+end
+
+function froude_number(
+    s,
+    gas,
+    external_force_field::U1,
+    L::U2,
+) where {U1<:Acceleration,U2<:Length}
+    return uconvert(
+        Unitful.NoUnits,
+        norm(velocity(s, gas)) / sqrt(norm(external_force_field) * L),
+    )
+end
+
+function fourier_number(
+    s,
+    gas::CaloricallyPerfectGas,
+    t::U1,
+    L::U2,
+) where {U1<:Time,U2<:Length}
+    return uconvert(Unitful.NoUnits, thermal_diffusivity(s, gas) * t / L^2)
+end
+
+

test/runtests.jl

@@ -33,21 +33,40 @@ using Unitful
             (v1, v2) = test_items[i]
             @test pressure(v1, DRY_AIR) ≈ pressure(v2, DRY_AIR)
             @test temperature(v1, DRY_AIR) ≈ temperature(v2, DRY_AIR)
-            @test all(
-                momentum_density(v1, DRY_AIR) ≈ momentum_density(v2, DRY_AIR),
-            )
+            @test all(momentum_density(v1, DRY_AIR) ≈ momentum_density(v2, DRY_AIR))
             @test all(velocity(v1, DRY_AIR) ≈ velocity(v2, DRY_AIR))
             @test (
                 specific_static_internal_energy(v1, DRY_AIR) ≈
                 specific_static_internal_energy(v2, DRY_AIR)
             )
             @test speed_of_sound(v1, DRY_AIR) ≈ speed_of_sound(v2, DRY_AIR)
+            @test thermal_diffusivity(v1, DRY_AIR) ≈ thermal_diffusivity(v2, DRY_AIR)
+            @test kinematic_viscosity(v1, DRY_AIR) ≈ kinematic_viscosity(v2, DRY_AIR)
+        end
+
+        length_scale = 1.0u"cm"
+        time_scale = 1.0u"ms"
+        g = 9.8086u"m/s^2"
+
+        @testset "Dimensionless Parameters Preserved: $i" for i ∈ eachindex(test_items)
+            (v1, v2) = test_items[i]
+            @test reynolds_number(v1, DRY_AIR, length_scale) ≈ reynolds_number(v2, DRY_AIR, length_scale)
+            @test froude_number(v1, DRY_AIR, g, length_scale) ≈ froude_number(v2, DRY_AIR, g, length_scale)
+            @test fourier_number(v1, DRY_AIR, time_scale, length_scale) ≈ fourier_number(v2, DRY_AIR, time_scale, length_scale)
+
         end
     end
 
     @testset "Convert Between Units" begin
         v1 = fill(ConservedProps(1.0, (2.0, 3.0), 4.0), 2)
-        v2 = fill(ConservedProps(1.0u"lb/ft^3", Quantity.((1.0,1.0), u"lb/ft^2/s"), 4.0u"J/mm^3"), 2)
+        v2 = fill(
+            ConservedProps(
+                1.0u"lb/ft^3",
+                Quantity.((1.0, 1.0), u"lb/ft^2/s"),
+                4.0u"J/mm^3",
+            ),
+            2,
+        )
         v1 .= v2
         @test all(v1 .≈ v2)
 
@@ -56,7 +75,7 @@ using Unitful
 
         @test all(v1 .≈ v3)
     end
-    
+
     @testset "Convert Primitve ↔ Conserved" begin
         s1 = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
         u = ConservedProps(s1, DRY_AIR)
@@ -65,7 +84,7 @@ using Unitful
         @test all(s1.M ≈ s2.M)
         @test s1.T ≈ s2.T
     end
-    
+
     @testset "Equivalency Between Unitful and Unitless" begin
         n = @SVector [-1.0, 0]
         t = @SVector [0.0, 1.0]
@@ -77,17 +96,15 @@ using Unitful
             @test sR_nounits[1] ≈ ustrip(sR.ρ)
             @test all(sR_nounits[2:end-1] .≈ sR.M)
             @test sR_nounits[end] ≈ ustrip(sR.T)
-            @test pressure(sL_nounits[1], sL_nounits[end], DRY_AIR) ≈
-                  pressure(sL, DRY_AIR)
-            @test pressure(sR_nounits[1], sR_nounits[end], DRY_AIR) ≈
-                  pressure(sR, DRY_AIR)
+            @test pressure(sL_nounits[1], sL_nounits[end], DRY_AIR) ≈ pressure(sL, DRY_AIR)
+            @test pressure(sR_nounits[1], sR_nounits[end], DRY_AIR) ≈ pressure(sR, DRY_AIR)
         end
-    
+
         uL = ConservedProps(sL, DRY_AIR)
         uL_nounits = state_to_vector(uL)
         uR = state_behind(uL, n, t, DRY_AIR)
         uR_nounits = conserved_state_behind(uL_nounits, n, t, DRY_AIR)
-    
+
         @testset "Conserved State Quantities" begin
             @test uR_nounits[1] ≈ ustrip(uR.ρ)
             @test all(uR_nounits[2:end-1] .≈ ustrip.(uR.ρv))
@@ -100,7 +117,7 @@ using Unitful
             @test ρe_nounits ≈ ustrip(static_internal_energy_density(uR))
             @test pressure(ρe_nounits, DRY_AIR) ≈ pressure(uR, DRY_AIR)
         end
-    
+
         @testset "Conversion" begin
             uR_nounits2 = conserved_state_vector(sR_nounits, DRY_AIR)
             sR_nounits2 = primitive_state_vector(uR_nounits, DRY_AIR)
@@ -109,14 +126,14 @@ using Unitful
             @test all(sR_nounits .≈ sR_nounits2)
         end
     end
-    
+
     # flux for the euler equations
-    function F(u::ConservedProps{N, T, U1, U2, U3}) where {N, T, U1, U2, U3}
+    function F(u::ConservedProps{N,T,U1,U2,U3}) where {N,T,U1,U2,U3}
         v = velocity(u)
         P = pressure(u, DRY_AIR)
         return vcat(u.ρv', u.ρv * v' + I * P, (v * (u.ρE + P))')
     end
-    
+
     @testset "Rankine-Hugoniot Condition" begin
         free_stream = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
         u_L = ConservedProps(free_stream, DRY_AIR)
@@ -131,21 +148,31 @@ using Unitful
             @test all(F(u_L) * n .≈ F(u_R) * n)
         end
     end
-    
+
     @testset "Speed of Sound" begin
         gas = DRY_AIR
         s = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
         u = ConservedProps(s, gas)
-    
+
         a_s = speed_of_sound(s, gas)
         a_u = speed_of_sound(u, gas)
-    
+
         a_s_pressure = speed_of_sound(density(s), pressure(s, gas), gas)
         a_u_pressure = speed_of_sound(density(u), pressure(u, gas), gas)
-
-        a_s_ie = speed_of_sound(density(s), momentum_density(s, gas), total_internal_energy_density(s, gas), gas)
-        a_u_ie = speed_of_sound(density(u), momentum_density(u), total_internal_energy_density(u), gas)
-
+
+        a_s_ie = speed_of_sound(
+            density(s),
+            momentum_density(s, gas),
+            total_internal_energy_density(s, gas),
+            gas,
+        )
+        a_u_ie = speed_of_sound(
+            density(u),
+            momentum_density(u),
+            total_internal_energy_density(u),
+            gas,
+        )
+
         @test a_s ≈ a_u
         @test a_s_pressure ≈ a_u_pressure
         @test a_s_ie ≈ a_u_ie