implement Billig formulae.

Project.toml

@@ -4,6 +4,7 @@ authors = ["Alex Fleming <[email protected]> and contributors"]
 version = "1.0.0-DEV"
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 

src/ShockwaveProperties.jl

@@ -1,6 +1,14 @@
 module ShockwaveProperties
 
+export DRY_AIR
+export BilligShockParametrization
+export shock_density_ratio, shock_pressure_ratio, shock_temperature_ratio
+export conserved_state_behind
+
 include("cpg.jl")
+include("billig.jl")
 include("normal_shocks.jl")
 
+const DRY_AIR = CaloricallyPerfectGas(1.0049, 0.7178, 0.0289647)
+
 end

src/billig.jl

@@ -1,3 +1,78 @@
-module BilligShock
-
-end
+"""
+Implementation of the shockwave parametrization about cylindrical blunt bodies developed by Billig in
+*Shock-wave shapes around spherical-and cylindrical-nosed bodies.*
+"""
+module BilligShockParametrization
+using LinearAlgebra
+
+export shock_front, shock_normal, shock_tangent
+
+const BILLIG_CONSTANTS = (A=0.386, B=4.67, C=1.386, D=1.8)
+
+# compute standoff distance relative to a disk of radius 1 in the plane
+@inline function standoff_distance(M_inf)
+    return BILLIG_CONSTANTS.A * exp(BILLIG_CONSTANTS.B / M_inf^2)
+end
+
+# compute critical radius relative to a disk of radius 1 in the plane
+@inline function critical_radius(M_inf)
+    return BILLIG_CONSTANTS.C * exp(BILLIG_CONSTANTS.D / (M_inf - 1)^0.75)
+end
+
+# compute sine and cosine of the critical angle 
+# formed between the x-axis and the tangent to the shock front at x=0
+@inline function critical_shock_angle(M_inf)
+    sinθ = 1 / M_inf
+    cosθ = sqrt(1 - sinθ^2)
+    return (sinθ, cosθ)
+end
+
+"""
+    shock_front(α, M_inf, R_b)
+Maps a parameter ``α`` to the shock wave generated about a body with radius ``R_b`` 
+by free-stream flow to the right with mach number ``M_infty``.
+
+The shock wave parametrization developed by *Billig* parametrizes the shock wave on the y-coordinate.
+"""
+function shock_front(α, M_inf, R_b)
+    y̅ = α / R_b
+    Δ = standoff_distance(M_inf)
+    R̅c = critical_radius(M_inf)
+    sinθ, cosθ = critical_shock_angle(M_inf)
+    cotθ = (cosθ / sinθ)
+    tanθ = (sinθ / cosθ)
+    return [R_b * (-1 - Δ + R̅c * cotθ^2 * (sqrt(1 + (y̅ * tanθ / R̅c)^2) - 1)), α]
+end
+
+# we'll manually take derivatives here... Zygote doesn't like 2nd derivatives.
+# we can also optimize out cot^2 * tan^2 ;)
+
+"""
+    shock_normal(α, M_inf, R_b)
+Maps a parameter ``α`` to the outward-facing normal to a shock wave generated about a body with radius ``R_b``
+by a free-stream flow to the right with nach number ``M_infty``
+
+The shock wave parametrization developed by *Billig* parametrizes the shock wave on the y-coordinate.
+"""
+function shock_normal(α, M_inf, R_b)
+    R_c = R_b * critical_radius(M_inf)
+    sinθ, cosθ = critical_shock_angle(M_inf)
+    tanθ = (sinθ / cosθ)
+    dxdα = α / (R_c * sqrt(1 + (α * tanθ / R_c)^2))
+    n = [-1, dxdα]
+    return n ./ norm(n)
+end
+
+"""
+    shock_normal(α, M_inf, R_b)
+Maps a parameter ``α`` to the tangent to a shock wave generated about a body with radius ``R_b``
+by a free-stream flow to the right with nach number ``M_infty``
+
+The shock wave parametrization developed by *Billig* parametrizes the shock wave on the y-coordinate.
+"""
+function shock_tangent(α, M_inf, R_b)
+    n = shock_normal(α, M_inf, R_b)
+    return [n[2], 1]
+end
+
+end # module

src/cpg.jl

@@ -1,24 +1,36 @@
 using Unitful
+using LinearAlgebra
 
 const _units_cvcp = u"J/kg/K"
 const _units_ℳ = u"kg/mol"
 const _units_ρ = u"kg/m^3"
+const _units_v = u"m/s"
+const _units_T = u"K"
+
+const _units_ρv = _units_ρ * _units_v
+const _units_int_e = _units_cvcp * _units_T
+const _units_ρE = _units_ρ * _units_int_e
+
 const _dimension_ρ = Unitful.dimension(_units_ρ)
+const _dimension_v = Unitful.dimension(_units_v)
+const _dimension_int_e = Unitful.dimension(_units_int_e)
+const _dimension_ρE = Unitful.dimension(_units_ρE)
+const _dimension_ρv = Unitful.dimension(_units_ρv)
 
 """
 Provides the properties of a calorically perfect gas (or mixture of gases). 
 """
-struct CaloricallyPerfectGas{T}
-    c_p::Quantity{T,Unitful.dimension(_units_cvcp),typeof(_units_cvcp)}
-    c_v::Quantity{T,Unitful.dimension(_units_cvcp),typeof(_units_cvcp)}
-    ℳ::Quantity{T,Unitful.dimension(_units_ℳ),typeof(_units_ℳ)}
+struct CaloricallyPerfectGas
+    c_p::Quantity{Float64,Unitful.dimension(_units_cvcp),typeof(_units_cvcp)}
+    c_v::Quantity{Float64,Unitful.dimension(_units_cvcp),typeof(_units_cvcp)}
+    ℳ::Quantity{Float64,Unitful.dimension(_units_ℳ),typeof(_units_ℳ)}
 
-    γ::T
-    R::Quantity{T,Unitful.dimension(_units_cvcp),typeof(_units_cvcp)}
+    γ::Float64
+    R::Quantity{Float64,Unitful.dimension(_units_cvcp),typeof(_units_cvcp)}
 end
 
-function CaloricallyPerfectGas(c_p::T, c_v::T, ℳ::T) where {T}
-    return CaloricallyPerfectGas{T}(
+function CaloricallyPerfectGas(c_p::Float64, c_v::Float64, ℳ::Float64)
+    return CaloricallyPerfectGas(
         Quantity(c_p, _units_cvcp),
         Quantity(c_v, _units_cvcp),
         Quantity(ℳ, _units_ℳ),
@@ -27,30 +39,91 @@ function CaloricallyPerfectGas(c_p::T, c_v::T, ℳ::T) where {T}
     )
 end
 
-const DRY_AIR = CaloricallyPerfectGas(1.0049, 0.7178, 0.0289647)
+enthalpy(gas::CaloricallyPerfectGas, T::Float64) = gas.c_p * Quantity(T, u"K")
+enthalpy(gas::CaloricallyPerfectGas, T::Quantity{Float64,Unitful.𝚯,Units}) where {Units} = gas.c_p * T
 
-enthalpy(gas::CaloricallyPerfectGas, T::Real) = gas.c_p * Quantity(T, u"K")
-enthalpy(gas::CaloricallyPerfectGas, T::Quantity{T1,Unitful.𝚯,Units}) where {T1,Units} = gas.c_p * T
+internal_energy(gas::CaloricallyPerfectGas, T::Float64) = gas.c_v * Quantity(T, u"K")
+internal_energy(gas::CaloricallyPerfectGas, T::Quantity{Float64,Unitful.𝚯,Units}) where {Units} = gas.c_v * T
 
-int_energy(gas::CaloricallyPerfectGas, T::Real) = gas.c_v * Quantity(T, u"K")
-int_energy(gas::CaloricallyPerfectGas, T::Quantity{T1,Unitful.𝚯,Units}) where {T1,Units} = gas.c_v * T
 
-function speed_of_sound(gas::CaloricallyPerfectGas, T::Real)
+function speed_of_sound(gas::CaloricallyPerfectGas, T::Float64)
     return sqrt(gas.γ * gas.R * Quantity(T, u"K"))
 end
 
-function speed_of_sound(gas::CaloricallyPerfectGas, T::Quantity{T1,Unitful.𝚯,Units}) where {T1,Units}
+"""
+Compute the speed of sound in an ideal gas at a temperature ``T``. We assume 
+that the gas is a non-dispersive medium.
+"""
+function speed_of_sound(gas::CaloricallyPerfectGas, T::Quantity{Float64,Unitful.𝚯,Units}) where {Units}
     return sqrt(gas.γ * gas.R * T)
 end
 
-function pressure(gas::CaloricallyPerfectGas, ρ::Real, T::Real)
+"""
+Compute the pressure in a calorically perfect gas from its density and temperature.
+"""
+function pressure(gas::CaloricallyPerfectGas, ρ::Float64, T::Float64)
     # calculate ρe then 
     # calculate p from calorically perfect gas relations
-    return (gas.γ - 1) * Quantity(ρ, _units_ρ) * int_energy(gas, T)
+    return (gas.γ - 1) * Quantity(ρ, _units_ρ) * internal_energy(gas, T)
+end
+
+function pressure(
+    gas::CaloricallyPerfectGas,
+    ρ::Quantity{Float64,_dimension_ρ,Units1},
+    T::Quantity{Float64,Unitful.𝚯,Units2}) where {Units1,Units2}
+    return (gas.γ - 1) * ρ * internal_energy(gas, T)
 end
 
-function pressure(gas::CaloricallyPerfectGas,
-    ρ::Quanity{U,_dimension_ρ,Units1},
-    T::Quantity{U,Unitful.𝚯,Units2}) where {U,Units1,Units2}
-    return (gas.γ-1) * ρ * T
+"""
+Compute the pressure in a calorically perfect gas from its internal energy density.
+"""
+pressure(gas::CaloricallyPerfectGas, ρe::Float64) = (gas.γ - 1) * Quantity(ρe, _units_ρE)
+pressure(gas::CaloricallyPerfectGas, ρe::Quantity{Float64,_dimension_ρE,Units}) where {Units} = (gas.γ - 1) * ρe
+
+"""
+Named tuple type of the conserved quantities in the Euler equations.
+"""
+ConservedState = @NamedTuple begin
+    ρ::Quantity{Float64,_dimension_ρ,_units_ρ}
+    ρv::Vector{Quantity{Float64,_dimension_ρv,_units_ρv}}
+    ρE::Quantity{Float64,_dimension_ρE,_units_ρE}
+end
+
+"""
+Compute the internal energy volume density (ρe) from conserved state quantities.
+"""
+internal_energy_density(state::ConservedState) = state.ρE - (state.ρv ⋅ state.ρv) / (2 * state.ρ)
+
+"""
+Named tuple type of the primitive (not conserved) quantities 
+that completely determine the state of a calorically perfect gas.
+"""
+PrimitiveState = @NamedTuple begin
+    ρ::Quantity{Float64,_dimension_ρ,_units_ρ}
+    M::Vector{Float64}
+    T::Quantity{Float64,Unitful.𝚯,_units_T}
 end
+
+"""
+Compute the pressure at a given state in a gas.
+"""
+pressure(gas::CaloricallyPerfectGas, state::ConservedState) = pressure(gas, internal_energy_density(state))
+pressure(gas::CaloricallyPerfectGas, state::PrimitiveState) = pressure(gas, state.ρ, state.T)
+
+"""
+Compute the temperature at a given state in a gas.
+"""
+temperature(gas::CaloricallyPerfectGas, state::ConservedState) = internal_energy_density(state) / gas.c_v
+temperature(::CaloricallyPerfectGas, state::PrimitiveState) = state.T
+
+"""
+Compute the density at a given state in a gas.
+"""
+density(::CaloricallyPerfectGas, state::Union{ConservedState,PrimitiveState}) = state.ρ
+
+"""
+Compute the speed of sound in a gas at a given state. We assume that the gas is a non-dispersive medium.
+"""
+function speed_of_sound(gas::CaloricallyPerfectGas, state::Union{ConservedState,PrimitiveState})
+    return speed_of_sound(gas, temperature(gas, state))
+end

src/normal_shocks.jl

@@ -1,54 +1,54 @@
+using LinearAlgebra
+
 """
 **Equation 4.8** from Anderson&Anderson
 
-Computes the density change across a shock wave where free stream mach number is ``M``
-and the outward (away from body) normal is ``n̂``
+Computes the density change across a shock wave. 
+The incoming flow has Mach number(s) ``M=[M_x, M_y]`` and the outward (away from body) normal is ``n̂``.
 """
-function shock_density_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
-    Mn2 = (M * n̂[1])^2
+@inline function shock_density_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
+    Mn2 = (M ⋅ n̂)^2
     return ((gas.γ + 1) * Mn2) / ((gas.γ - 1) * Mn2 + 2)
 end
 
 """
 **Equation 4.9** from Anderson&Anderson
 
-Computes the pressure ratio across a shock wave where the free stream mach number is ``M``
-and the outward (away from body) normal is ``n̂``.
+Computes the pressure ratio across a shock wave.
+The incoming flow has Mach number(s) ``M=[M_x, M_y]`` and the outward (away from body) normal is ``n̂``.
 """
-function shock_pressure_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
-    Mn2 = (M * n̂[1])^2
+@inline function shock_pressure_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
+    Mn2 = (M ⋅ n̂)^2
     return 1 + (2 * gas.γ) / (gas.γ + 1) * (Mn2 - 1)
 end
 
 """
 **Equation 4.10** from Anderson&Anderson
 
-Computes the normal Mach number ratio across a shock wave where the free stream mach number is ``M``
-and the outward (away from body) normal is ``n̂``.
+Computes the normal Mach number ratio across a shock wave.
+The incoming flow has Mach number(s) ``M=[M_x, M_y]`` and the outward (away from body) normal is ``n̂``.
 """
 function shock_normal_mach_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
-    MnL2 = (M * n̂[1])^2
+    MnL2 = (M ⋅ n̂)^2
     Mn2sqr = (MnL2 + (2 / (gas.γ - 1))) / (2 * gas.γ / (gas.γ - 1) * MnL2 - 1)
     Mn2 = sqrt(Mn2sqr)
     return Mn2
 end
 
-
 """
-Computes the normal velocity ratio across a shock wave where the free stream mach number is ``M``
-and the outward (away from body) normal is ``n̂``.
+Computes the normal velocity ratio across a shock wave.
+The incoming flow has Mach number(s) ``M=[M_x, M_y]`` and the outward (away from body) normal is ``n̂``.
 
 Derived from speed of sound proportional to the square root of temperature.
-
 Useful for computing momentum ratio.
 """
-function shock_normal_velocity_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
+@inline function shock_normal_velocity_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
     return sqrt(1 / shock_temperature_ratio(M, n̂; gas=gas)) * shock_normal_mach_ratio(M, n̂; gas=gas)
 end
 
 """
-Computes the tangential mach number ratio across a shock wave, where the free stream mach number is ``M``
-and the outward (away from body) normal is ``n̂``.
+Computes the tangential mach number ratio across a shock wave.
+The incoming flow has Mach number(s) ``M=[M_x, M_y]`` and the outward (away from body) normal is ``n̂``.
 
 Derived from speed of sound proportional to the square root of temperature.
 """
@@ -59,13 +59,20 @@ end
 """
 **Equation 4.11** from Anderson&Anderson
 
-Computes the temperature across a shock wave where the free stream mach number is ``M``
-and the outward (away from body) normal is ``n̂``.
+Computes the temperature ratio across a shock wave.
+The incoming flow has Mach number(s) ``M=[M_x, M_y]`` and the outward (away from body) normal is ``n̂``.
 """
 function shock_temperature_ratio(M, n̂; gas::CaloricallyPerfectGas=DRY_AIR)
     return shock_pressure_ratio(M, n̂; gas=gas) / shock_density_ratio(M, n̂; gas=gas)
 end
 
 """
-Computes the 
-"""
+Computes the conservation properties behind a shockwave.
+The outward (away from body) normal is ``n̂``.
+"""
+function conserved_state_behind(state_L::ConservedState, n̂, t̂; gas::CaloricallyPerfectGas=DRY_AIR)
+    @assert t̂ ⋅ n̂ == 0. "tangent and normal vectors should be normal."
+    M_L = state_L.ρv / (state_L.ρ * speed_of_sound(gas, state_L))
+    mach_t_ratio = shock_tangent_mach_ratio(M_L, n̂; gas=gas)
+    M_R_n = (M_L⋅n̂) * shock_normal_mach_ratio(M_L, n̂; gas=gas)
+end

test/runtests.jl

@@ -1,6 +1,18 @@
+using LinearAlgebra
 using ShockwaveProperties
 using Test
 
 @testset "ShockwaveProperties.jl" begin
-    # Write your tests here.
-end
+    @testset "BilligShockParametrization" begin
+        using .BilligShockParametrization
+        @testset "Shockwave Normals" begin
+            y = -10:0.1:10
+            for R_b ∈ 1.0:0.25:4.0, M ∈ 1.5:0.5:5.0
+                vals = mapreduce(hcat, y) do α
+                    shock_normal(α, M, R_b) ⋅ shock_tangent(α, M, R_b)
+                end
+                @test all(≈(0), vals)
+            end
+        end
+    end
+end