Merge pull request #2 from ngomezve/main

new combustion models + changes in Mach number + standard atmosphere

src/Gas.jl

@@ -27,6 +27,7 @@ Constructs `Gas` with given composition `Y`
 """
 function Gas(Y::AbstractVector)
    gas = Gas(); gas.Y = convert(Vector{Float64}, Y)
+   set_TP!(gas, Tstd, Pstd) #setting temperature and pressure to recalculate thermodynamic properties
    return gas
 end
 
@@ -71,7 +72,7 @@ function Gas()
 
 end
 
-# Overload Base.getproperty for convinence
+# Overload Base.getproperty for convenience
 function Base.getproperty(gas::Gas, sym::Symbol)
    if sym === :h_T # dh/dT
       return getfield(gas, :cp)
@@ -110,6 +111,17 @@ function Base.getproperty(gas::Gas, sym::Symbol)
       return cp/(cp - R)
    elseif sym === :gamma
       return getproperty(gas, :γ)
+   elseif sym === :ρ
+      R = getproperty(gas, :R)
+      T = getfield(gas, :T)
+      P = getfield(gas, :P)
+      return P / (R * T)
+   elseif sym === :rho
+      return getproperty(gas, :ρ)
+   elseif sym === :ν     
+      return 1 / getproperty(gas, :ρ)
+   elseif sym === :nu     
+      return getproperty(gas, :ν)
    elseif sym === :X # Get mole fractions
       Y = getfield(gas, :Y)
       MW = spdict.MW
@@ -140,7 +152,7 @@ function Base.setproperty!(gas::Gas, sym::Symbol, val::Float64)
    if sym === :T
       setfield!(gas, :T, val) # first set T
       setfield!(gas, :Tarray, Tarray!(val, getfield(gas, :Tarray))) # update Tarray
-      TT = view(getfield(gas, :Tarray), :) # Just convinence
+      TT = view(getfield(gas, :Tarray), :) # Just convenience
       # Next set the cp, h and s of the gas
       ## Get the right coefficients 
       ## (assumes Tmid is always 1000.0. Check performed in readThermo.jl.):
@@ -176,7 +188,7 @@ function Base.setproperty!(gas::Gas, sym::Symbol, val::Float64)
    ## Setting Pressure
    elseif sym === :P
       setfield!(gas, :P, val)
-      TT = view(getfield(gas, :Tarray), :) # Just convinence
+      TT = view(getfield(gas, :Tarray), :) # Just convenience
       # Next set s of the gas
       ## Get the right coefficients 
       ## (assumes Tmid is always 1000.0. Check performed in readThermo.jl.):

src/Gas1D.jl

@@ -114,6 +114,17 @@ function Base.getproperty(gas::Gas1D, sym::Symbol)
       return cp / (cp - R)
    elseif sym === :gamma
       return getproperty(gas, :γ)
+   elseif sym === :ρ
+      R = getproperty(gas, :R)
+      T = getfield(gas, :T)
+      P = getfield(gas, :P)
+      return P / (R * T)
+   elseif sym === :rho
+      return getproperty(gas, :ρ)
+   elseif sym === :ν     
+      return 1 / getproperty(gas, :ρ)
+   elseif sym === :nu     
+      return getproperty(gas, :ν)
    else
       return getfield(gas, sym)
    end

src/IdealGases.jl

@@ -27,6 +27,7 @@ export print_thermo_table
 include("utils.jl")
 export X2Y, Y2X
 include("idealgasthermo.jl")
+include("atmosphere.jl")
 
 gas = Gas()
 gas.X = Xair

src/atmosphere.jl

@@ -0,0 +1,58 @@
+"""
+    standard_atmosphere(z, alt_type = "geometric")
+Calculate the gas state at a given desired altitude in the atmosphere using the 1976
+US Standard Atmosphere. Model valid between -5 and 86 km in geometric altitude. Accepts 
+geometric altitude ("alt_type == "geometric") or geopotential altitude ("alt_type == "geopotential"). 
+Model from: `https://ntrs.nasa.gov/citations/19770009539`
+"""
+function standard_atmosphere(z, alt_type = "geometric")
+
+    g = 9.80665 #Acceleration of gravity on Earth surface (z = 0)
+    R_e = 6.356766e6; #Radius of Earth
+
+    #Initialize gas (composition does not change)
+    gas = Gas1D()
+    R = gas.R
+
+    #Calculate geopotential altitude
+    if alt_type == "geometric"
+        H = (R_e * z) / (R_e + z) #Geopotential height
+
+    elseif alt_type == "geopotential"
+        H = z
+    end
+
+    #Database with 1976 US Standard Atmosphere temperatures and pressures
+    H0s = [0, 11e3, 20e3, 32e3, 47e3, 51e3, 71e3] #atm. regions
+    λs = [-6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2e-3] #lapse rates
+    T0s = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65] #temperatures at interfaces
+    P0s = [101325, 22632.1, 5474.89, 868.019, 110.906, 66.9389, 3.95642] #pressures at interfaces
+
+    #Table lookup to find atmospheric region
+    ind = 1
+    for i in 1:length(H0s)
+        if H > H0s[i]
+            ind = i
+        end
+    end
+    #Store base height, temperature, pressure and lapse rate
+    H0 = H0s[ind]
+    λ = λs[ind]
+    T0 = T0s[ind]
+    P0 = P0s[ind]
+
+    #Find temperature and pressure using the lapse rate
+    #Use closed-form solutions
+    if λ != 0
+        P = P0 * ((T0 + λ * (H - H0)) / T0)^(-g / (R * λ))
+        T = T0 + λ * (H - H0)
+
+    elseif λ == 0
+        T = T0
+        P = P0 * exp(-g * (H - H0) / (R * T))
+    end
+
+    #set the gas to the correct T and P
+    set_TP!(gas, T, P)
+    return gas
+end

src/combustion.jl

@@ -278,6 +278,10 @@ function vitiated_mixture(fuel::AbstractSpecies, oxidizer::AbstractSpecies,
 
     Xvitiated::Dict{String, Float64} = mergewith(+, Xin, Xdict)
 
+    for key in keys(Xvitiated) #Normalize such that sum(Xi) = 1
+        Xvitiated[key] = Xvitiated[key] / (1 + ηburn * molFAR)
+    end
+
     return Xvitiated
 end  # function vitiated_gas
 
@@ -288,7 +292,7 @@ vitiated_mixture(fuel, oxidizer, stoich_FOR(fuel, oxidizer))
     vitiated_mixture(fuel::AbstractString, oxidizer::AbstractString, 
     FAR::Float64, ηburn::Float64=1.0)
 
-Convinence function that finds fuel and oxidizer from thermo database
+Convenience function that finds fuel and oxidizer from thermo database
 """
 function vitiated_mixture(fuel::AbstractString, oxidizer::AbstractString, 
     FAR::Float64, ηburn::Float64=1.0)
@@ -506,3 +510,130 @@ function fixed_fuel_vitiated_species(fuel, oxidizer, ηburn::Float64=1.0)
 
     return burntgas
 end  # function fixed_fuel_vitiated_species
+
+"""
+    fuel_combustion(gas_ox::AbstractGas, fuel::String, Tf::Float64, FAR::Float64, 
+    ηburn::Float64 = 1.0, hvap::Float64 = 0.0)
+
+This function returns an `AbstractGas` with the combustion products at the combustor exit temperature
+for a given fuel type, oxidizer gas at a given enthalpy, and FAR. It includes the combustion 
+efficiency and the fuel enthalpy of vaporization as optional inputs. 
+"""
+function fuel_combustion(gas_ox::AbstractGas, fuel::String, Tf::Float64, FAR::Float64, 
+    ηburn::Float64 = 1.0, hvap::Float64 = 0.0)
+
+    #Create variables corresponding to the oxidizer and fuel species and mixtures
+    fuel_sps = species_in_spdict(fuel)
+
+    #Find oxidizer composition
+    if typeof(gas_ox) == Gas1D
+        gas_sps = gas_ox.comp_sp
+    else
+        if "Air" in keys(gas_ox.Xdict)
+            Xin = Xair
+        else
+            Xin = gas_ox.Xdict
+        end
+        gas_sps = generate_composite_species(IdealGases.Xidict2Array(Xin))
+    end
+
+    #Find the vectors with the fuel mole and mass fractions
+    Xfuel = Xidict2Array(Dict([(fuel, 1.0)])) #Mole fraction
+    Yfuel = X2Y(Xfuel) #Mass fraction
+    gas_fuel = Gas(Yfuel) #Create a fuel gas to calculate enthalpy
+
+    set_TP!(gas_fuel, Tf, gas_ox.P) #Set the fuel gas at the correct conditions
+
+    #Find product composition and mole and mass fractions
+    Xprod_dict =  vitiated_mixture(fuel_sps, gas_sps, FAR, ηburn)
+    Xprod = Xidict2Array(Xprod_dict)
+    Yprod = X2Y(Xprod)
+
+    #Initialize output 
+    gas_prod = Gas(Yprod)
+
+    # Combustion enthalpy calculations
+    #enthalpy before reaction = enthalpy after reaction
+    h0 = (gas_ox.h + FAR * (gas_fuel.h - abs(hvap))) / (1 + FAR)
+
+    #Set the products at the correct enthalpy and pressure
+    set_hP!(gas_prod, h0, gas_ox.P)
+    return gas_prod
+end
+
+"""
+    gas_burn(gas_ox::AbstractGas, fuel::String, Tf::Float64, Tburn::Float64, 
+    ηburn::Float64 = 1.0, hvap::Float64 = 0.0)
+
+This function returns a tuple containing the FAR and an `AbstractGas` with the combustion 
+products for a given fuel type, oxidizer gas at a given enthalpy, and desired combustor exit temperature. 
+It includes the combustion efficiency and the fuel enthalpy of vaporization as optional inputs. 
+"""
+function gas_burn(gas_ox::AbstractGas, fuel::String, Tf::Float64, Tburn::Float64, 
+    ηburn::Float64 = 1.0, hvap::Float64 = 0.0)
+
+    #Create variables corresponding to the oxidizer and fuel species and mixtures
+    fuel_sps = species_in_spdict(fuel)
+
+    #Extract composite species with oxidizer gas composition
+    if typeof(gas_ox) == Gas1D
+        gas_sps = gas_ox.comp_sp
+    else
+        if "Air" in keys(gas_ox.Xdict)
+            Xin = Xair
+        else
+            Xin = gas_ox.Xdict
+        end
+        gas_sps = generate_composite_species(IdealGases.Xidict2Array(Xin))
+    end
+
+    #Find the vectors with the fuel mole and mass fractions
+    Xfuel = Xidict2Array(Dict([(fuel, 1.0)])) #Mole fraction
+    Yfuel = X2Y(Xfuel) #Mass fraction
+    gas_fuel = Gas(Yfuel) #Create a fuel gas to calculate enthalpy
+
+    set_TP!(gas_fuel, Tf, gas_ox.P) #Set the fuel gas at the correct conditions
+
+    #Store enthalpies of oxidizer and fuel at original temperatures
+    ho = gas_ox.h 
+    hf = gas_fuel.h 
+
+    #Find change in gas composition for a FAR of 1
+    nCO2, nN2, nH2O, nO2 = reaction_change_molar_fraction(fuel_sps.name)
+
+    names = ["CO2", "H2O", "N2", "O2"]
+    ΔX    = [nCO2,  nH2O,  nN2,  nO2]
+    Xdict = Dict(zip(names, ΔX))
+
+    Xc = Xidict2Array(Xdict)
+    Yc = X2Y(Xc) #mass fraction change in combustion for FAR = 1
+
+    gas_c = Gas(Yc) #Create a "virtual" gas with the changes in combustion, for enthalpy
+                    #calculations
+    set_TP!(gas_c, Tburn, gas_ox.P)
+
+    hc = gas_c.h #Enthalpy change for FAR = 1
+
+    gas_ox_burnt = deepcopy(gas_ox) #make a copy to avoid modifying the input
+    set_TP!(gas_ox_burnt, Tburn, gas_ox.P)
+    ha = gas_ox_burnt.h #Enthalpy of original oxidizer at final temperature
+
+    set_TP!(gas_fuel, Tburn, gas_ox.P)
+    hff = gas_fuel.h #Enthalpy of fuel at final temperature
+
+    #Find FAR corresponding to Tburn
+    FAR = (ha - ho) / (hf - ηburn * hc - (1 - ηburn) * hff - abs(hvap)) #solve for FAR 
+
+    #Find product composition and mole and mass fractions
+    Xprod_dict = vitiated_mixture(fuel_sps, gas_sps, FAR, ηburn)
+    Xprod = Xidict2Array(Xprod_dict)
+    Yprod = X2Y(Xprod)
+
+    #Initialize output 
+    gas_prod = Gas(Yprod)
+
+    #Set the products at the correct enthalpy and pressure
+    set_TP!(gas_prod, Tburn, gas_ox.P)
+
+    return FAR, gas_prod
+end