diff --git a/Project.toml b/Project.toml index 2485fb2..edf1a15 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Alex Fleming 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" diff --git a/src/ShockwaveProperties.jl b/src/ShockwaveProperties.jl index 25d4f0f..446f452 100644 --- a/src/ShockwaveProperties.jl +++ b/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 diff --git a/src/billig.jl b/src/billig.jl index 5aa2edb..b42f546 100644 --- a/src/billig.jl +++ b/src/billig.jl @@ -1,3 +1,78 @@ -module BilligShock - -end \ No newline at end of file +""" +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 \ No newline at end of file diff --git a/src/cpg.jl b/src/cpg.jl index cd78022..82efb29 100644 --- a/src/cpg.jl +++ b/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 \ No newline at end of file diff --git a/src/normal_shocks.jl b/src/normal_shocks.jl index 797c22c..7d3a2d5 100644 --- a/src/normal_shocks.jl +++ b/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 -""" \ No newline at end of file +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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index c5d3409..255593a 100644 --- a/test/runtests.jl +++ b/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 \ No newline at end of file