Generalized propulsion functions

docs/src/aero/drag.md

@@ -109,7 +109,7 @@ aerodynamics.airfun(cl, τ, Mach, air::aerodynamics.airfoil)
 
 ## Total drag calculation
 ```@docs
-aerodynamics.cdsum!(parg, para, pare, wing, htail, vtail, icdfun)
+aerodynamics.cdsum!(ac, imission, ip, icdfun)
 ```
 ---
 

docs/src/sizing/sizing.md

@@ -17,11 +17,11 @@ TASOPT.woper(ac, itermax)
 A sized aircraft's mission performance can be obtained (`mission!`), along with operation constraints via a pitch trim calculation (`balance`) and balanced field length calculation (`takeoff!`).
 
 ```@docs
-TASOPT.mission!(pari, parg, parm, para, pare, fuse, wing, htail, vtail, Ldebug)
+TASOPT.mission!(ac, imission, Ldebug)
 
 TASOPT.takeoff!(ac; printTO)
 
-TASOPT.balance(pari, parg, para, fuse, wing, htail, vtail, rfuel, rpay, ξpay, itrim)
+TASOPT.balance(ac, imission, ip, rfuel, rpay, ξpay, itrim)
 
 ```
 

src/IO/read_input.jl

@@ -125,6 +125,7 @@ fuselage = Fuselage()
 wing = Wing()
 htail = Tail()
 vtail = Tail()
+engine = Engine()
 
 # Setup mission variables
 ranges = readmis("range")
@@ -180,13 +181,17 @@ pari[iiopt] = read_input("optimize", options, doptions)
 propsys = read_input("prop_sys_arch", options, doptions)
 if lowercase(propsys) == "tf"
     pari[iiengtype] = 1
+    engine.model_name = "turbofan_md" #TODO: add more engine models
+    engine.weight_model_name = "turbofan"
 elseif lowercase(propsys) == "te"
     pari[iiengtype] = 0
 else
+
     error("Propulsion system \"$propsys\" specified. Choose between
     > TF - turbo-fan
     > TE - turbo-electric" )
 end
+store_engine_model!(engine) #Produce and store the engine functions 
 
 engloc = read_input("engine_location", options, doptions)
 
@@ -1052,9 +1057,8 @@ dHEx = dprop["HeatExchangers"]
     pare[ieTurbCepsilon, :, :] .= read_input("turbine_cooler_effectiveness", HEx, dHEx)
     pare[ieTurbCMp, :, :] .= read_input("turbine_cooler_inlet_mach", HEx, dHEx)
 
-
 return TASOPT.aircraft(name, description,
-    pari, parg, parm, para, pare, [false], fuse_tank, fuselage, wing, htail, vtail)
+    pari, parg, parm, para, pare, [false], fuse_tank, fuselage, wing, htail, vtail, engine)
 
 end
 

src/aero/cdsum.jl

@@ -1,5 +1,5 @@
 """
-    cdsum!(parg,para,pare, wing, htail, vtail, icdfun)
+   cdsum!(ac, imission, ip, icdfun)
 
 Calculates aircraft `CD` components for operating point, ipoint.
 If `icdfun=1`, computes wing `cdf`,`cdp` from airfoil database # `iairf`,
@@ -25,26 +25,25 @@ where:
 
 
 !!! details "🔃 Inputs and Outputs"
-      **Inputs:**
-      - `parg::AbstractArray{Float64}`: Vector of `aircraft` model geometry parameters.
-      - `para::AbstractArray{Float64}`: Vector of `aircraft` model aerodynamic parameters.
-      - `pare::AbstractArray{Float64}`: Vector of `aircraft` model engine parameters.
-      - `Wing::TASOPT.Wing`: Wing Structure.
-      - `Htail::TASOPT.Tail`: Htail Structure.
-      - `Vtail::TASOPT.Tail`: Vtail Structure.
-      - `icdfun::Integer`: Flag if drag should be computed (=1) or if para values should be used (=0).
+**Inputs:**
+      - `ac::aircraft`: aircraft data storage object
+      - `imission::Int64`: mission index
+      - `icdfun::Int64`: Flag if drag should be computed (=1) or if para values should be used (=0).
 
       **Outputs:**
       - No explicit outputs. Computed drag values are saved to `para` of `aircraft` model.
 
 See Section 2.14 of the [TASOPT Technical Desc](@ref dreladocs).
 See also [`trefftz1`](@ref), [`fusebl!`](@ref), [`surfcd2`](@ref), [`surfcd`](@ref), [`cfturb`](@ref), and `cditrp`.
-
-!!! compat "Future Changes"
-      In an upcoming revision, an `aircraft` struct and auxiliary indices will be passed in lieu of pre-sliced `par` arrays.
-
 """
-function cdsum!(parg,para,pare, wing, htail, vtail, icdfun)
+function cdsum!(ac, imission, ip, icdfun)
+      #Unpack data storage
+      parg = ac.parg
+      para = view(ac.para, :, ip, imission)
+      pare = view(ac.pare, :, ip, imission)
+      wing = ac.wing
+      htail = ac.htail
+      vtail = ac.vtail
 
       Ldebug = false
 #      Ldebug = true

src/balance/balance.jl

@@ -1,39 +1,30 @@
 """
-balance(pari, parg, para, fuse, wing, htail, vtail, rfuel, rpay, ξpay, itrim)
+      balance(ac, imission, ip, rfuel, rpay, ξpay, itrim)
 
 Makes one of three (or none) changes to achieve pitch trim
 calculates resulting CG, CP, NP locations.
 
 Inputs:
-- `pari[.]`  integer flag array
-- `parg[.]`  geometry parameter array
-- `para[.]`  aero parameter array
-- `fuse`      TASOPT.Fuselage
-- `wing`      TASOPT.Wing
-- `htail`     TASOPT.Tail
-- `vtail`     TASOPT.Tail
-- `rfuel`    fuel fraction   Wfuel_actual/Wfuel_MTOW
+- `ac::aircraft`: structure with aircraft parameters
+- `imission::Int64`: mission index (1 is design mission)
+- `ip::Int64`: mission point index
+- `rfuel`    fuel fraction
 - `rpay`     payload fraction Wpay_actual/Wpay_MTOW
 - `ξpay`    partial-payload packing location
     * = 0.0   all the way in front  of cabin
     * = 0.5   all the way in middle of cabin
     * = 1.0   all the way in back   of cabin
-- `iengloc`  engine location index
 - `itrim`      = 0  no changes
     * = 1  adjust CLh   (horizontal tail cl)
     * = 2  adjust Sh    (horizontal tail area)
     * = 3  adjust xwbox (wing box location)
 
 Outputs: 
-
-- `para[iaxCG]`  center of gravity
-- `para[iaxCP]`  center of pressure ( = xCG if itrim=1,2,3 )
-- `para[iaxNP]`  neutral point location
-
-!!! compat "Future Changes"
-      In an upcoming revision, an `aircraft` struct and auxiliary indices will be passed in lieu of pre-sliced `par` arrays.
+No direct outputs. Fields in `ac` are modified.
 """
-function balance(pari, parg, para, fuse, wing, htail, vtail, rfuel, rpay, ξpay, itrim)
+function balance(ac, imission, ip, rfuel, rpay, ξpay, itrim)
+      #Unpack aircraft
+      pari, parg, parm, para, pare, fuse, fuse_tank, wing, htail, vtail = unpack_ac(ac, imission, ip = ip)
 
       iengloc = pari[iiengloc]
 
@@ -69,7 +60,7 @@ function balance(pari, parg, para, fuse, wing, htail, vtail, rfuel, rpay, ξpay,
 
       xwbox = wing.layout.box_x
 
-      rfuelF, rfuelB, rpayF, rpayB, xcgF, xcgB = cglpay(pari, parg, fuse, wing, htail, vtail)
+      rfuelF, rfuelB, rpayF, rpayB, xcgF, xcgB = cglpay(ac)
 
       #---- wing centroid offset from wingbox, assumed fixed in CG calculations
       dxwing = wing.layout.x - wing.layout.box_x
@@ -90,7 +81,6 @@ function balance(pari, parg, para, fuse, wing, htail, vtail, rfuel, rpay, ξpay,
       Sh1 = Sh
 
       #---- total weight, weight moment, and derivatives
-
       W = rpay * Wpay +
           rfuel * Wfuel +
           Wfuse +
@@ -273,20 +263,18 @@ Calculates resulting CG, CP, NP locations
 
 Inputs:  
 
-- `pari[.]`  integer fla array
+- `ac`       aircraft object
 - `parg[.]`  geometry parameter array
 - `paraF[.]` aero parameter array for fwdCG case
 - `paraB[.]` aero parameter array for aft CG case
 - `paraC[.]` aero parameter array for cruise tail CL case
 
 Outputs: 
-
-- `parg[igSh]`    HT area
-- `parg[igxwbox]` wingbox location
-- `parg[igxwing]` wing centroid location
-
+No direct outputs. Fields in `ac` are modified.
 """
-function htsize(pari, parg, paraF, paraB, paraC,fuse,wing, htail, vtail)
+function htsize(ac, paraF, paraB, paraC)
+      #TODO find a way to remove the para inputs and use ac instead
+      pari, parg, fuse, fuse_tank, wing, htail, vtail = unpack_ac_components(ac)
 
       itmax = 10
       toler = 1.0e-7
@@ -298,7 +286,7 @@ function htsize(pari, parg, paraF, paraB, paraC,fuse,wing, htail, vtail)
       cosL = cosd(sweep)
 
       #---- set CG limits with worst-case payload arrangements
-      rfuelF, rfuelB, rpayF, rpayB, xcgF, xcgB = cglpay(pari, parg,fuse, wing, htail, vtail)
+      rfuelF, rfuelB, rpayF, rpayB, xcgF, xcgB = cglpay(ac)
 
       rpayC = 1.0
 
@@ -490,7 +478,7 @@ function htsize(pari, parg, paraF, paraB, paraC,fuse,wing, htail, vtail)
                   lhtail = htail.layout.x - (xwbox + dxwing)
                   lhtail_xw = -1.0
 
-                  Vh = htail.volume      # Tail volume is specified here!
+                  Vh = htail.volume     # Tail volume is specified here!
                   Sh = Vh * S * cma / lhtail # Corresponding tail surf area
 
                   r[1] = Sh * lhtail - Vh * S * cma  # This is 0 (residual is 0 ∵ Vh = Sh*lh/(S*cma))
@@ -665,7 +653,8 @@ which gives an explicit solution for `rpayF`,`rpayB`.
 The alternative 2D search for `rfuel`,`rpay` is kinda ugly, 
 and unwarranted in practice.
 """
-function cglpay(pari, parg, fuse, wing, htail, vtail)
+function cglpay(ac)
+      pari, parg, fuse, fuse_tank, wing, htail, vtail = unpack_ac_components(ac)
 
       Wpay = parg[igWpay]
       Wfuel = parg[igWfuel]

src/engine/engine.jl

@@ -8,8 +8,10 @@ using NLopt
 using Roots
 using LinearAlgebra
 
-export tfcalc!, mcool, Tmcalc, gas_tset, gaschem
-export tfweight, ddct, ddat, gct, gat, tfsize!, Ncmap, ecmap, Ncmap1, ecmap1, etmap, Pimap, tfoper!
+export Engine, store_engine_model!, extract_engine_model
+
+export tfwrap!, tfcalc!, mcool, Tmcalc, gas_tset, gaschem
+export tfweightwrap!, tfweight, ddct, ddat, gct, gat, tfsize!, Ncmap, ecmap, Ncmap1, ecmap1, etmap, Pimap, tfoper!
 
 export gassum, gassumd, gas_prat, gas_delh, gas_delhd, gas_burn, gas_burnd, gas_mach, gas_machd, gas_mass, gasfuel, fuelLHV, gasPr
 export hxdesign!, radiator_design!, hxweight, resetHXs
@@ -18,17 +20,20 @@ import ..TASOPT: __TASOPTindices__, __TASOPTroot__, StructuralAlloy
 
 include(__TASOPTindices__)
 include(joinpath(__TASOPTroot__,"misc/constants.jl"))
+include("../misc/engine.jl")
 include("gasfun.jl")
 include("gascalc.jl")
 # include("tfan.jl")
 include("tfmap.jl")
-include("tfcool.jl")
-include("tfsize.jl")
+include("turbofan/tfcool.jl")
+include("turbofan/tfsize.jl")
 include("gaussn.jl")
 include("compare.jl")
-include("tfoper.jl")
-include("tfcalc.jl")
-include("tfweight.jl")
+include("turbofan/tfoper.jl")
+include("turbofan/tfcalc.jl")
+include("turbofan/tfweight.jl")
+include("turbofan/tfwrap.jl")
+include("turbofan/tfweightwrap.jl")
 include("hxfun.jl")
 include("PEMfuelcell.jl")
 

src/engine/tfcalc.jl → src/engine/turbofan/tfcalc.jl

@@ -326,6 +326,9 @@ function tfcalc!(pari, parg, para, pare, wing, ip,
                 pare[iepihtD] = pihtD
                 pare[iepiltD] = piltD
 
+                #Fuel mass flow rate
+                pare[iemfuel] = ff * mcore * neng
+
                 #--------------------------------------------------------------------------
         else
                 #----- off-design operation case
@@ -701,6 +704,8 @@ function tfcalc!(pari, parg, para, pare, wing, ip,
         pare[ieetaht] = etaht
         pare[ieetalt] = etalt
 
+        #Fuel mass flow rate
+        pare[iemfuel] = ff * mcore * neng
 
         if (M5 <= 0.999999)
                 ichoke5 = 0

src/engine/tfweight.jl → src/engine/turbofan/tfweight.jl

@@ -1,32 +1,37 @@
 """
-      tfweight(iengwgt, Gearf, OPR, BPR, mdotc, dfan, rSnace,
-      dlcomp, neng, feadd, fpylon, HXs, HX_add_mass_frac)
+    tfweight(ac, HXs)
 
 Engine weight estimation function using Giulia Pantalone, Drela, or Fitzgerald model.
 
 !!! details "🔃 Inputs and Outputs"
     **Input:**
-    - `iengwgt`: Engine model index, Drela=0, Fitzgerald=1, and Pantalone>=3, 
-    - `OPR`: Overall pressure ratio.
-    - `BPR`: By-pass ratio.
-    - `mdotc`: Engine core mass flow rate.
-    - `dfan`: Fan diameter.
-    - `rSnace`: 
-    - `dlcomp`:
-    - `neng`: Number of engines.
-    - `feadd`: Fuel system weight ratio.
-    - `fpylon`: Pylon weight fraction.
+    - `ac::aircraft`: aircraft object
     - `HXs`: vector with heat exchanger performance data
-    - `HX_add_mass_frac`: added mass fraction to heat exchangers
 
     **Output:**
     - `Weng`: Total engine weight.
     - `Wnac`: Nacelle weight. 
     - `Webare`: Bare engine weight.
     - `Snace1`: Nacelle area.
 """
-function tfweight(iengwgt, Gearf, OPR, BPR, mdotc, dfan, rSnace,
-    dlcomp, neng, feadd, fpylon, HXs, HX_add_mass_frac)
+function tfweight(ac, HXs)
+
+    TSL = Tref
+    pSL = pref
+    ip = ipcruise1 #OPR and BPR designed for start-of-cruise
+
+    iengwgt = ac.pari[iiengwgt]
+    Gearf = ac.parg[igGearf]
+    mdotc = ac.pared[iemblcD, ip] * sqrt(Tref / TSL) * (pSL / pref)
+    BPR = ac.pared[ieBPR, ip]
+    OPR = ac.pared[iepilc, ip] * ac.pared[iepihc, ip]
+    dfan = ac.parg[igdfan]
+    dlcomp = ac.parg[igdlcomp]
+    rSnace = ac.parg[igrSnace]
+    neng = ac.parg[igneng]
+    feadd = ac.parg[igfeadd]
+    fpylon = ac.parg[igfpylon]
+    HX_add_mass_frac = ac.parg[igHXaddmassfrac]
 
     # include("constants.inc")
 

src/engine/turbofan/tfweightwrap.jl

@@ -0,0 +1,42 @@
+"""
+    tfweightwrap!(ac, HXs)
+
+General function to estimate and store the weight of a turbofan engine. 
+This function is basically a wrapper on tfweight, going from the
+basic aircraft inputs to those required by the function and storing the outputs.
+
+!!! details "🔃 Inputs and Outputs"
+    **Input:**
+    - `ac::aircraft`: aircraft object
+    - `HXs`: vector with heat exchanger performance data
+
+    **Output:**
+    No direct outputs. The `ac` object gets modified with the engine weights.
+"""
+function tfweightwrap!(ac, HXs)
+    pari = ac.pari
+    parg = ac.parg
+    wing = ac.wing
+    neng = parg[igneng]
+
+    Weng, Wnace, Webare, W_HXs, Snace1 = tfweight(ac, HXs)
+
+    parg[igWeng] = Weng
+    parg[igWebare] = Webare
+    parg[igWnace] = Wnace
+    parg[igWHXs] = W_HXs
+
+    # set new nacelle area / reference area  fraction fSnace
+    S = wing.layout.S
+
+    Snace = Snace1 * neng
+    fSnace = Snace / S
+    parg[igfSnace] = fSnace
+
+    # set new nacelle area / reference area  fraction fSnace
+    Snace = Snace1 * neng
+    fSnace = Snace / S
+    parg[igfSnace] = fSnace
+    lnace = parg[igdfan] * parg[igrSnace] * 0.15
+    parg[iglnace] = lnace
+end