Possible bug in flight condition:static pressure calculation (fc:stat:P)

[flight-test-engineering]

When running the high bypass turbofan example at OFF_DESIGN conditions, I get an anomaly in thrust at 19,000ft. The figure below is for Mach=0.75 and PC=1.0:

Further, I tried to search upstream where this anomaly originates and my best guess is that something is off in the calculation of static pressure at the flight conditions (fc:stat:P). This is the plot I get:

Which is unexpected.

This is the file I used to generate the data. It is the original HBTF example file with the following changes:

finer OFF_DESIGN flight envelope
addition of a logger function to dump the data to a file

import sys
import numpy as np
import openmdao.api as om
import pycycle.api as pyc
from csv import writer
import time

class HBTF(pyc.Cycle):

    def initialize(self):
        # Initialize the model here by setting option variables such as a switch for design vs off-des cases
        self.options.declare('throttle_mode', default='T4', values=['T4', 'percent_thrust'])

        super().initialize()


    def setup(self):
        #Setup the problem by including all the relavant components here - comp, burner, turbine etc
        
        #Create any relavent short hands here:
        design = self.options['design']
        
        USE_TABULAR = False
        if USE_TABULAR: 
            self.options['thermo_method'] = 'TABULAR'
            self.options['thermo_data'] = pyc.AIR_JETA_TAB_SPEC
            FUEL_TYPE = 'FAR'
        else: 
            self.options['thermo_method'] = 'CEA'
            self.options['thermo_data'] = pyc.species_data.janaf
            FUEL_TYPE = 'Jet-A(g)'

        
        #Add subsystems to build the engine deck:
        self.add_subsystem('fc', pyc.FlightConditions())
        self.add_subsystem('inlet', pyc.Inlet())
        
        # Note variable promotion for the fan -- 
        # the LP spool speed and the fan speed are INPUTS that are promoted:
        # Note here that promotion aliases are used. Here Nmech is being aliased to LP_Nmech
        # check out: http://openmdao.org/twodocs/versions/latest/features/core_features/grouping_components/add_subsystem.html?highlight=alias
        self.add_subsystem('fan', pyc.Compressor(map_data=pyc.FanMap,
                                        bleed_names=[], map_extrap=True), promotes_inputs=[('Nmech','LP_Nmech')])
        self.add_subsystem('splitter', pyc.Splitter())
        self.add_subsystem('duct4', pyc.Duct())
        self.add_subsystem('lpc', pyc.Compressor(map_data=pyc.LPCMap,
                                        map_extrap=True),promotes_inputs=[('Nmech','LP_Nmech')])
        self.add_subsystem('duct6', pyc.Duct())
        self.add_subsystem('hpc', pyc.Compressor(map_data=pyc.HPCMap,
                                        bleed_names=['cool1','cool2','cust'], map_extrap=True),promotes_inputs=[('Nmech','HP_Nmech')])
        self.add_subsystem('bld3', pyc.BleedOut(bleed_names=['cool3','cool4']))
        self.add_subsystem('burner', pyc.Combustor(fuel_type=FUEL_TYPE))
        self.add_subsystem('hpt', pyc.Turbine(map_data=pyc.HPTMap,
                                        bleed_names=['cool3','cool4'], map_extrap=True),promotes_inputs=[('Nmech','HP_Nmech')])
        self.add_subsystem('duct11', pyc.Duct())
        self.add_subsystem('lpt', pyc.Turbine(map_data=pyc.LPTMap,
                                        bleed_names=['cool1','cool2'], map_extrap=True),promotes_inputs=[('Nmech','LP_Nmech')])
        self.add_subsystem('duct13', pyc.Duct())
        self.add_subsystem('core_nozz', pyc.Nozzle(nozzType='CV', lossCoef='Cv'))

        self.add_subsystem('byp_bld', pyc.BleedOut(bleed_names=['bypBld']))
        self.add_subsystem('duct15', pyc.Duct())
        self.add_subsystem('byp_nozz', pyc.Nozzle(nozzType='CV', lossCoef='Cv'))
        
        #Create shaft instances. Note that LP shaft has 3 ports! => no gearbox
        self.add_subsystem('lp_shaft', pyc.Shaft(num_ports=3),promotes_inputs=[('Nmech','LP_Nmech')])
        self.add_subsystem('hp_shaft', pyc.Shaft(num_ports=2),promotes_inputs=[('Nmech','HP_Nmech')])
        self.add_subsystem('perf', pyc.Performance(num_nozzles=2, num_burners=1))
    
        # Now use the explicit connect method to make connections -- connect(<from>, <to>)
        
        #Connect the inputs to perf group
        self.connect('inlet.Fl_O:tot:P', 'perf.Pt2')
        self.connect('hpc.Fl_O:tot:P', 'perf.Pt3')
        self.connect('burner.Wfuel', 'perf.Wfuel_0')
        self.connect('inlet.F_ram', 'perf.ram_drag')
        self.connect('core_nozz.Fg', 'perf.Fg_0')
        self.connect('byp_nozz.Fg', 'perf.Fg_1')
        
        #LP-shaft connections
        self.connect('fan.trq', 'lp_shaft.trq_0')
        self.connect('lpc.trq', 'lp_shaft.trq_1')
        self.connect('lpt.trq', 'lp_shaft.trq_2')
        #HP-shaft connections
        self.connect('hpc.trq', 'hp_shaft.trq_0')
        self.connect('hpt.trq', 'hp_shaft.trq_1')
        #Ideally expanding flow by conneting flight condition static pressure to nozzle exhaust pressure
        self.connect('fc.Fl_O:stat:P', 'core_nozz.Ps_exhaust')
        self.connect('fc.Fl_O:stat:P', 'byp_nozz.Ps_exhaust')
        
        #Create a balance component
        # Balances can be a bit confusing, here's some explanation -
        #   State Variables:
        #           (W)        Inlet mass flow rate to implictly balance thrust
        #                      LHS: perf.Fn  == RHS: Thrust requirement (set when TF is instantiated)
        #
        #           (FAR)      Fuel-air ratio to balance Tt4
        #                      LHS: burner.Fl_O:tot:T  == RHS: Tt4 target (set when TF is instantiated)
        #
        #           (lpt_PR)   LPT press ratio to balance shaft power on the low spool
        #           (hpt_PR)   HPT press ratio to balance shaft power on the high spool
        # Ref: look at the XDSM diagrams in the pyCycle paper and this:
        # http://openmdao.org/twodocs/versions/latest/features/building_blocks/components/balance_comp.html

        balance = self.add_subsystem('balance', om.BalanceComp())
        if design:
            balance.add_balance('W', units='lbm/s', eq_units='lbf')
            #Here balance.W is implicit state variable that is the OUTPUT of balance object
            self.connect('balance.W', 'fc.W') #Connect the output of balance to the relevant input
            self.connect('perf.Fn', 'balance.lhs:W')       #This statement makes perf.Fn the LHS of the balance eqn.
            self.promotes('balance', inputs=[('rhs:W', 'Fn_DES')])

            balance.add_balance('FAR', eq_units='degR', lower=1e-4, val=.017)
            self.connect('balance.FAR', 'burner.Fl_I:FAR')
            self.connect('burner.Fl_O:tot:T', 'balance.lhs:FAR')
            self.promotes('balance', inputs=[('rhs:FAR', 'T4_MAX')])
            
            # Note that for the following two balances the mult val is set to -1 so that the NET torque is zero
            balance.add_balance('lpt_PR', val=1.5, lower=1.001, upper=8,
                                eq_units='hp', use_mult=True, mult_val=-1)
            self.connect('balance.lpt_PR', 'lpt.PR')
            self.connect('lp_shaft.pwr_in_real', 'balance.lhs:lpt_PR')
            self.connect('lp_shaft.pwr_out_real', 'balance.rhs:lpt_PR')

            balance.add_balance('hpt_PR', val=1.5, lower=1.001, upper=8,
                                eq_units='hp', use_mult=True, mult_val=-1)
            self.connect('balance.hpt_PR', 'hpt.PR')
            self.connect('hp_shaft.pwr_in_real', 'balance.lhs:hpt_PR')
            self.connect('hp_shaft.pwr_out_real', 'balance.rhs:hpt_PR')

        else:
            
            #In OFF-DESIGN mode we need to redefine the balances:
            #   State Variables:
            #           (W)        Inlet mass flow rate to balance core flow area
            #                      LHS: core_nozz.Throat:stat:area == Area from DESIGN calculation 
            #
            #           (FAR)      Fuel-air ratio to balance Thrust req.
            #                      LHS: perf.Fn  == RHS: Thrust requirement (set when TF is instantiated)
            #
            #           (BPR)      Bypass ratio to balance byp. noz. area
            #                      LHS: byp_nozz.Throat:stat:area == Area from DESIGN calculation
            #
            #           (lp_Nmech)   LP spool speed to balance shaft power on the low spool
            #           (hp_Nmech)   HP spool speed to balance shaft power on the high spool

            if self.options['throttle_mode'] == 'T4': 
                balance.add_balance('FAR', val=0.017, lower=1e-4, eq_units='degR')
                self.connect('balance.FAR', 'burner.Fl_I:FAR')
                self.connect('burner.Fl_O:tot:T', 'balance.lhs:FAR')
                self.promotes('balance', inputs=[('rhs:FAR', 'T4_MAX')])

            elif self.options['throttle_mode'] == 'percent_thrust': 
                balance.add_balance('FAR', val=0.017, lower=1e-4, eq_units='lbf', use_mult=True)
                self.connect('balance.FAR', 'burner.Fl_I:FAR')
                self.connect('perf.Fn', 'balance.rhs:FAR')
                self.promotes('balance', inputs=[('mult:FAR', 'PC'), ('lhs:FAR', 'Fn_max')])


            balance.add_balance('W', units='lbm/s', lower=10., upper=1000., eq_units='inch**2')
            self.connect('balance.W', 'fc.W')
            self.connect('core_nozz.Throat:stat:area', 'balance.lhs:W')

            balance.add_balance('BPR', lower=2., upper=10., eq_units='inch**2')
            self.connect('balance.BPR', 'splitter.BPR')
            self.connect('byp_nozz.Throat:stat:area', 'balance.lhs:BPR')

            # Again for the following two balances the mult val is set to -1 so that the NET torque is zero
            balance.add_balance('lp_Nmech', val=1.5, units='rpm', lower=500., eq_units='hp', use_mult=True, mult_val=-1)
            self.connect('balance.lp_Nmech', 'LP_Nmech')
            self.connect('lp_shaft.pwr_in_real', 'balance.lhs:lp_Nmech')
            self.connect('lp_shaft.pwr_out_real', 'balance.rhs:lp_Nmech')

            balance.add_balance('hp_Nmech', val=1.5, units='rpm', lower=500., eq_units='hp', use_mult=True, mult_val=-1)
            self.connect('balance.hp_Nmech', 'HP_Nmech')
            self.connect('hp_shaft.pwr_in_real', 'balance.lhs:hp_Nmech')
            self.connect('hp_shaft.pwr_out_real', 'balance.rhs:hp_Nmech')
            
            # Specify the order in which the subsystems are executed:
            
            # self.set_order(['balance', 'fc', 'inlet', 'fan', 'splitter', 'duct4', 'lpc', 'duct6', 'hpc', 'bld3', 'burner', 'hpt', 'duct11',
            #                 'lpt', 'duct13', 'core_nozz', 'byp_bld', 'duct15', 'byp_nozz', 'lp_shaft', 'hp_shaft', 'perf'])
        
        # Set up all the flow connections:
        self.pyc_connect_flow('fc.Fl_O', 'inlet.Fl_I')
        self.pyc_connect_flow('inlet.Fl_O', 'fan.Fl_I')
        self.pyc_connect_flow('fan.Fl_O', 'splitter.Fl_I')
        self.pyc_connect_flow('splitter.Fl_O1', 'duct4.Fl_I')
        self.pyc_connect_flow('duct4.Fl_O', 'lpc.Fl_I')
        self.pyc_connect_flow('lpc.Fl_O', 'duct6.Fl_I')
        self.pyc_connect_flow('duct6.Fl_O', 'hpc.Fl_I')
        self.pyc_connect_flow('hpc.Fl_O', 'bld3.Fl_I')
        self.pyc_connect_flow('bld3.Fl_O', 'burner.Fl_I')
        self.pyc_connect_flow('burner.Fl_O', 'hpt.Fl_I')
        self.pyc_connect_flow('hpt.Fl_O', 'duct11.Fl_I')
        self.pyc_connect_flow('duct11.Fl_O', 'lpt.Fl_I')
        self.pyc_connect_flow('lpt.Fl_O', 'duct13.Fl_I')
        self.pyc_connect_flow('duct13.Fl_O','core_nozz.Fl_I')
        self.pyc_connect_flow('splitter.Fl_O2', 'byp_bld.Fl_I')
        self.pyc_connect_flow('byp_bld.Fl_O', 'duct15.Fl_I')
        self.pyc_connect_flow('duct15.Fl_O', 'byp_nozz.Fl_I')

        #Bleed flows:
        self.pyc_connect_flow('hpc.cool1', 'lpt.cool1', connect_stat=False)
        self.pyc_connect_flow('hpc.cool2', 'lpt.cool2', connect_stat=False)
        self.pyc_connect_flow('bld3.cool3', 'hpt.cool3', connect_stat=False)
        self.pyc_connect_flow('bld3.cool4', 'hpt.cool4', connect_stat=False)
        
        #Specify solver settings:
        newton = self.nonlinear_solver = om.NewtonSolver()
        newton.options['atol'] = 1e-8

        # set this very small, so it never activates and we rely on atol
        newton.options['rtol'] = 1e-99 
        newton.options['iprint'] = 2
        newton.options['maxiter'] = 50
        newton.options['solve_subsystems'] = True
        newton.options['max_sub_solves'] = 1000
        newton.options['reraise_child_analysiserror'] = False
        # ls = newton.linesearch = BoundsEnforceLS()
        ls = newton.linesearch = om.ArmijoGoldsteinLS()
        ls.options['maxiter'] = 6 # original was 3
        ls.options['rho'] = 0.75
        # ls.options['print_bound_enforce'] = True

        self.linear_solver = om.DirectSolver()

        super().setup()

class MPhbtf(pyc.MPCycle):

    def setup(self):

        self.pyc_add_pnt('DESIGN', HBTF(thermo_method='CEA')) # Create an instace of the High Bypass ratio Turbofan

        self.set_input_defaults('DESIGN.inlet.MN', 0.751)
        self.set_input_defaults('DESIGN.fan.MN', 0.4578)
        self.set_input_defaults('DESIGN.splitter.BPR', 5.105) ### 5.105     #### BPR !!!! ####
        self.set_input_defaults('DESIGN.splitter.MN1', 0.3104)
        self.set_input_defaults('DESIGN.splitter.MN2', 0.4518)
        self.set_input_defaults('DESIGN.duct4.MN', 0.3121)
        self.set_input_defaults('DESIGN.lpc.MN', 0.3059)
        self.set_input_defaults('DESIGN.duct6.MN', 0.3563)
        self.set_input_defaults('DESIGN.hpc.MN', 0.2442)
        self.set_input_defaults('DESIGN.bld3.MN', 0.3000)
        self.set_input_defaults('DESIGN.burner.MN', 0.1025)
        self.set_input_defaults('DESIGN.hpt.MN', 0.3650)
        self.set_input_defaults('DESIGN.duct11.MN', 0.3063)
        self.set_input_defaults('DESIGN.lpt.MN', 0.4127)
        self.set_input_defaults('DESIGN.duct13.MN', 0.4463)
        self.set_input_defaults('DESIGN.byp_bld.MN', 0.4489)
        self.set_input_defaults('DESIGN.duct15.MN', 0.4589)
        self.set_input_defaults('DESIGN.LP_Nmech', 4666.1, units='rpm') ### 4666.1 -> CF34 5452
        self.set_input_defaults('DESIGN.HP_Nmech', 14705.7, units='rpm') ### 14705.7 -> CF34 18616

        # --- Set up bleed values -----
        
        self.pyc_add_cycle_param('inlet.ram_recovery', 0.9990)
        self.pyc_add_cycle_param('duct4.dPqP', 0.0048)
        self.pyc_add_cycle_param('duct6.dPqP', 0.0101)
        self.pyc_add_cycle_param('burner.dPqP', 0.0540)
        self.pyc_add_cycle_param('duct11.dPqP', 0.0051)
        self.pyc_add_cycle_param('duct13.dPqP', 0.0107)
        self.pyc_add_cycle_param('duct15.dPqP', 0.0149)
        self.pyc_add_cycle_param('core_nozz.Cv', 0.9933)
        self.pyc_add_cycle_param('byp_bld.bypBld:frac_W', 0.005)
        self.pyc_add_cycle_param('byp_nozz.Cv', 0.9939) ### 0.9939
        self.pyc_add_cycle_param('hpc.cool1:frac_W', 0.050708) ### 0.050708
        self.pyc_add_cycle_param('hpc.cool1:frac_P', 0.5)
        self.pyc_add_cycle_param('hpc.cool1:frac_work', 0.5)
        self.pyc_add_cycle_param('hpc.cool2:frac_W', 0.020274) ### 0.020274
        self.pyc_add_cycle_param('hpc.cool2:frac_P', 0.55)
        self.pyc_add_cycle_param('hpc.cool2:frac_work', 0.5)
        self.pyc_add_cycle_param('bld3.cool3:frac_W', 0.067214) ### 0.067214
        self.pyc_add_cycle_param('bld3.cool4:frac_W', 0.101256)
        self.pyc_add_cycle_param('hpc.cust:frac_P', 0.5)
        self.pyc_add_cycle_param('hpc.cust:frac_work', 0.5)
        self.pyc_add_cycle_param('hpc.cust:frac_W', 0.0445)
        self.pyc_add_cycle_param('hpt.cool3:frac_P', 1.0)
        self.pyc_add_cycle_param('hpt.cool4:frac_P', 0.0)
        self.pyc_add_cycle_param('lpt.cool1:frac_P', 1.0)
        self.pyc_add_cycle_param('lpt.cool2:frac_P', 0.0)
        self.pyc_add_cycle_param('hp_shaft.HPX', 250.0, units='hp') ### 250.0

        self.od_pts = ['OD_full_pwr', 'OD_part_pwr'] 

        self.od_MNs = [0.8, 0.8] ###[0.8, 0.8]
        self.od_alts = [35000.0, 35000.0]
        self.od_Fn_target = [5500.0, 5300.0] ###[5500.0, 5300.0]
        self.od_dTs = [0.0, 0.0]

        self.pyc_add_pnt('OD_full_pwr', HBTF(design=False, thermo_method='CEA', throttle_mode='T4'))

        self.set_input_defaults('OD_full_pwr.fc.MN', 0.8) ###0.8
        self.set_input_defaults('OD_full_pwr.fc.alt', 35000, units='ft')
        self.set_input_defaults('OD_full_pwr.fc.dTs', 0., units='degR')

        self.pyc_add_pnt('OD_part_pwr', HBTF(design=False, thermo_method='CEA', throttle_mode='percent_thrust'))

        self.set_input_defaults('OD_part_pwr.fc.MN', 0.8) ###0.8
        self.set_input_defaults('OD_part_pwr.fc.alt', 35000, units='ft')
        self.set_input_defaults('OD_part_pwr.fc.dTs', 0., units='degR')

        self.connect('OD_full_pwr.perf.Fn', 'OD_part_pwr.Fn_max')

        self.pyc_use_default_des_od_conns()

        #Set up the RHS of the balances!
        self.pyc_connect_des_od('core_nozz.Throat:stat:area','balance.rhs:W')
        self.pyc_connect_des_od('byp_nozz.Throat:stat:area','balance.rhs:BPR')

        super().setup()
        
        
def viewer(prob, pt, file=sys.stdout):
    """
    print a report of all the relevant cycle properties
    """

    if pt == 'DESIGN':
        MN = prob['DESIGN.fc.Fl_O:stat:MN']
        LPT_PR = prob['DESIGN.balance.lpt_PR']
        HPT_PR = prob['DESIGN.balance.hpt_PR']
        FAR = prob['DESIGN.balance.FAR']
    else:
        MN = prob[pt+'.fc.Fl_O:stat:MN']
        LPT_PR = prob[pt+'.lpt.PR']
        HPT_PR = prob[pt+'.hpt.PR']
        FAR = prob[pt+'.balance.FAR']

    summary_data = (MN[0], prob[pt+'.fc.alt'][0], prob[pt+'.inlet.Fl_O:stat:W'][0], prob[pt+'.perf.Fn'][0],
                        prob[pt+'.perf.Fg'][0], prob[pt+'.inlet.F_ram'][0], prob[pt+'.perf.OPR'][0],
                        prob[pt+'.perf.TSFC'][0], prob[pt+'.splitter.BPR'][0])
  
    print(file=file, flush=True)
    print(file=file, flush=True)
    print(file=file, flush=True)
    print("----------------------------------------------------------------------------", file=file, flush=True)
    print("                              POINT:", pt, file=file, flush=True)
    print("----------------------------------------------------------------------------", file=file, flush=True)
    print("                       PERFORMANCE CHARACTERISTICS", file=file, flush=True)
    print("    Mach      Alt       W      Fn      Fg    Fram     OPR     TSFC      BPR ", file=file, flush=True)
    print(" %7.5f  %7.1f %7.3f %7.1f %7.1f %7.1f %7.3f  %7.5f  %7.3f" %summary_data, file=file, flush=True)


    fs_names = ['fc.Fl_O', 'inlet.Fl_O', 'fan.Fl_O', 'splitter.Fl_O1', 'splitter.Fl_O2',
                'duct4.Fl_O', 'lpc.Fl_O', 'duct6.Fl_O', 'hpc.Fl_O', 'bld3.Fl_O', 'burner.Fl_O',
                'hpt.Fl_O', 'duct11.Fl_O', 'lpt.Fl_O', 'duct13.Fl_O', 'core_nozz.Fl_O', 'byp_bld.Fl_O',
                'duct15.Fl_O', 'byp_nozz.Fl_O']
    fs_full_names = [f'{pt}.{fs}' for fs in fs_names]
    pyc.print_flow_station(prob, fs_full_names, file=file)

    comp_names = ['fan', 'lpc', 'hpc']
    comp_full_names = [f'{pt}.{c}' for c in comp_names]
    pyc.print_compressor(prob, comp_full_names, file=file)

    pyc.print_burner(prob, [f'{pt}.burner'], file=file)

    turb_names = ['hpt', 'lpt']
    turb_full_names = [f'{pt}.{t}' for t in turb_names]
    pyc.print_turbine(prob, turb_full_names, file=file)

    noz_names = ['core_nozz', 'byp_nozz']
    noz_full_names = [f'{pt}.{n}' for n in noz_names]
    pyc.print_nozzle(prob, noz_full_names, file=file)

    shaft_names = ['hp_shaft', 'lp_shaft']
    shaft_full_names = [f'{pt}.{s}' for s in shaft_names]
    pyc.print_shaft(prob, shaft_full_names, file=file)

    bleed_names = ['hpc', 'bld3', 'byp_bld']
    bleed_full_names = [f'{pt}.{b}' for b in bleed_names]
    pyc.print_bleed(prob, bleed_full_names, file=file)


def logger(prob, file=sys.stdout):
    """
    log data of relevant cycle properties
    """
    
    pt = 'OD_part_pwr'

    MN = prob[pt+'.fc.Fl_O:stat:MN'][0]
    alt = prob[pt+'.fc.alt'][0]
    Fn = prob[pt+'.perf.Fn'][0]
    Fg = prob[pt+'.perf.Fg'][0]
    F_ram = prob[pt+'.inlet.F_ram'][0]
    BPR = prob[pt+'.splitter.BPR'][0]
    OPR = prob[pt+'.perf.OPR'][0]
    
    LPT_PR = prob[pt+'.lpt.PR'][0]
    HPT_PR = prob[pt+'.hpt.PR'][0]
    
    FAR = prob[pt+'.balance.FAR'][0]
    PC = prob[pt+'.PC'][0]
    TSFC = prob[pt+'.perf.TSFC'][0]

    Wf = prob[pt+'.burner.mix_fuel.thermo_add.mix:W'][0]
    
    line_props = [alt, MN, LPT_PR, HPT_PR,
                 FAR, PC, Fn, Fg, F_ram,
                 OPR, TSFC, BPR, Wf]

    ### flight (static) conditions
    st_props = ['stat:P', 'stat:T', 'stat:h', 'stat:S', 'stat:rho']

    st = 'fc.conv.fs.Fl_O'

    line_props += [prob[pt+'.'+st+':'+st_p][0] for st_p in st_props]

    ###
    
    # stations which we want to log at
    st_list = ['inlet.real_flow.Fl_O', 'fan.real_flow.Fl_O', 'lpc.real_flow.Fl_O', 
               'hpc.real_flow.Fl_O', 'burner.Fl_O',
                'hpt.real_flow.Fl_O', 'lpt.real_flow.Fl_O']
    
    st_names = [st_name.split('.')[0] for st_name in st_list]

    # properties we want to log
    st_props = ['tot:P', 'tot:T', 'tot:h', 'tot:S', 'tot:rho']
    
    for st_idx, st in enumerate(st_list):
        props = [prob[pt+'.'+st+':'+st_p][0] for st_p in st_props]
        props.append(prob[pt+'.'+st_names[st_idx]+'.out_stat.flow_static.Fl_O:stat:W'][0])
        [line_props.append(p) for p in props]

    return line_props
   
#flight envelope definition
alternate_altitudes = np.arange(35000, 1000, -250)
flight_env = [(0.8, 35000), (0.75, 35000)]
flight_env = flight_env+ [(0.75, a) for a in alternate_altitudes]
PCs = np.array([1.0])

prob = om.Problem()

prob.model = mp_hbtf = MPhbtf()

prob.setup()

prob.set_val('DESIGN.fan.PR', 1.685)
prob.set_val('DESIGN.fan.eff', 0.8948)

prob.set_val('DESIGN.lpc.PR', 1.935)  #1.935 is the original; CF34 is 2.0062
prob.set_val('DESIGN.lpc.eff', 0.9243)

prob.set_val('DESIGN.hpc.PR', 9.369)  #9.369 is the original; CF34 is 8.0755
prob.set_val('DESIGN.hpc.eff', 0.8707)

prob.set_val('DESIGN.hpt.eff', 0.888) # 0.888
prob.set_val('DESIGN.lpt.eff', 0.8996) # 0.8996

prob.set_val('DESIGN.fc.alt', 35000., units='ft')
prob.set_val('DESIGN.fc.MN', 0.78)

prob.set_val('DESIGN.T4_MAX', 2857, units='degR') # note, this is not EGT
prob.set_val('DESIGN.Fn_DES', 5900.0, units='lbf') # tweak this until EASA CO2 is obtained at 794ft ## original was 5900lbf

prob.set_val('OD_full_pwr.T4_MAX', 2857, units='degR') # note, this is not EGT
prob.set_val('OD_part_pwr.PC', 0.8)


# Set initial guesses for balances
prob['DESIGN.balance.FAR'] = 0.025
prob['DESIGN.balance.W'] = 100.
prob['DESIGN.balance.lpt_PR'] = 4.0
prob['DESIGN.balance.hpt_PR'] = 3.0
prob['DESIGN.fc.balance.Pt'] = 5.2
prob['DESIGN.fc.balance.Tt'] = 440.0


for pt in ['OD_full_pwr', 'OD_part_pwr']:

    # initial guesses
    prob[pt+'.balance.FAR'] = 0.02467
    prob[pt+'.balance.W'] = 300
    prob[pt+'.balance.BPR'] = 5.105
    prob[pt+'.balance.lp_Nmech'] = 5000 
    prob[pt+'.balance.hp_Nmech'] = 15000 
    prob[pt+'.hpt.PR'] = 3.
    prob[pt+'.lpt.PR'] = 4.
    prob[pt+'.fan.map.RlineMap'] = 2.0
    prob[pt+'.lpc.map.RlineMap'] = 2.0
    prob[pt+'.hpc.map.RlineMap'] = 2.0



prob.set_solver_print(level=-1)
prob.set_solver_print(level=2, depth=1)


##################

# create dataframe labels/header
# collect variables we want to save for LUT
header = ['alt', 'MN', 'LPT_PR', 'HPT_PR',
                 'FAR', 'PC', 'Fn', 'Fg', 'F_ram',
                 'OPR', 'TSFC', 'BPR', 'Wf']

st_list = ['fc']
st = st_list[0]

# static inlet conditions
st_props = ['stat:P', 'stat:T', 'stat:h', 'stat:S', 'stat:rho']

header = header + [[st+':'+st_p][0] for st_p in st_props]

# engine stations flow properties
st_list = ['inlet.real_flow.Fl_O', 'fan.real_flow.Fl_O', 'lpc.real_flow.Fl_O', 
           'hpc.real_flow.Fl_O', 'burner.Fl_O',
           'hpt.real_flow.Fl_O', 'lpt.real_flow.Fl_O']
           
st_props = ['tot:P', 'tot:T', 'tot:h', 'tot:S', 'tot:rho']

st_names = [st_name.split('.')[0] for st_name in st_list]

for st_idx, st in enumerate(st_list):
    props = [st_names[st_idx]+':'+st_p for st_p in st_props]
    props.append(st_names[st_idx]+':stat:W')
    header = header + props

##################




tot_time = 0 # time counter

viewer_file = open('03_hbtf_view_19k_orig_pcycle.out', 'w')
data_log_file = open('pycycle_HBTF_deck_19k.csv', 'w')
writer_object = writer(data_log_file)
writer_object.writerow(header)
first_pass = True

#first_pass = True
for MN, alt in flight_env: 

    # NOTE: You never change the MN,alt for the 
    # design point because that is a fixed reference condition.
    
    start_time = time.time()

    print('***'*10)
    print(f'* MN: {MN}, alt: {alt}')
    print('***'*10)
    prob['OD_full_pwr.fc.MN'] = MN
    prob['OD_full_pwr.fc.alt'] = alt

    prob['OD_part_pwr.fc.MN'] = MN
    prob['OD_part_pwr.fc.alt'] = alt

    for PC in PCs:
        print(f'* MN: {MN}, alt: {alt}, PC: {PC}, Run time= {(time.time() - start_time):.5f}', end="\r")
        prob['OD_part_pwr.PC'] = PC
        
        prob.run_model()
        
        run_data = logger(prob)

        writer_object.writerow(run_data)
        

        if first_pass: 
            first_pass = False
            viewer(prob, 'OD_full_pwr')
            viewer(prob, 'OD_part_pwr')
            
            viewer(prob, 'DESIGN', file=viewer_file)

            viewer(prob, 'OD_full_pwr', file=viewer_file)
            viewer(prob, 'OD_part_pwr', file=viewer_file)
        else:
            #viewer(prob, 'DESIGN', file=viewer_file)
            viewer(prob, 'OD_full_pwr')
            viewer(prob, 'OD_full_pwr', file=viewer_file)
            viewer(prob, 'OD_part_pwr', file=viewer_file)


    runtime = time.time() - start_time
    tot_time += runtime

    print()

print(f'total run time: {tot_time}')

viewer_file.close()
data_log_file.close()

[flight-test-engineering]

Further investigating this issue, I get the following results when running the code below derived directly from pycycle/elements/flight_conditions.py (I will post a separate issue regarding an AttributeError thrown if you try to run it straight from the repo):

import openmdao.api as om

from pycycle.thermo.cea import species_data
from pycycle.thermo.cea import thermo_data
from pycycle.constants import THERMO_DEFAULT_COMPOSITIONS
from pycycle.elements.ambient import Ambient
from pycycle.elements.flow_start import FlowStart
from pycycle.thermo.thermo import ThermoAdd, Thermo
from pycycle.element_base import Element


class FlightConditions(Element):
    """Determines total and static flow properties given an altitude and Mach number using the input atmosphere model"""

    def initialize(self):
        self.options.declare('composition', default=None,
                              desc='composition of the flow. If None, default for thermo package is used')
        self.options.declare('reactant', default=False, types=(bool, str), 
                              desc='If False, flow matches base composition. If a string, then that reactant '
                                   'is mixed into the flow at at the ratio set by the `mix_ratio` input')
        self.options.declare('mix_ratio_name', default='mix:ratio', 
                             desc='The name of the input that governs the mix ratio of the reactant to the primary flow')

        super().initialize()
        
    def pyc_setup_output_ports(self): 
        composition = self.options['composition']
        thermo_method = self.options['thermo_method']
        thermo_data = self.options['thermo_data']
        reactant = self.options['reactant']

        if reactant is not False: 
            # just make a local ThermoAdd, since the FS will make the real one for us later
            thermo_add = ThermoAdd(method=thermo_method, mix_mode='reactant', 
                                   thermo_kwargs={'spec':thermo_data, 
                                                  'inflow_composition':composition, 
                                                  'mix_composition':reactant, })
          
            self.init_output_flow('Fl_O', thermo_add)

        else: 
            if composition is None: 
                composition = THERMO_DEFAULT_COMPOSITIONS[thermo_method]
            self.init_output_flow('Fl_O', composition)
        

    def setup(self):
        thermo_method = self.options['thermo_method']
        thermo_data = self.options['thermo_data']
        reactant = self.options['reactant']
        mix_ratio_name = self.options['mix_ratio_name']

        # composition = self.Fl_O_data['Fl_O']
        composition = self.options['composition']

        self.add_subsystem('ambient', Ambient(), promotes=('alt', 'dTs'))  # inputs

        conv = self.add_subsystem('conv', om.Group(), promotes=['*'])
        if reactant is not False:
            proms = ['Fl_O:*', 'MN', 'W', mix_ratio_name]
        else:
            proms = ['Fl_O:*', 'MN', 'W']
        fs_start = conv.add_subsystem('fs', FlowStart(thermo_method=thermo_method,
                                                      thermo_data=thermo_data, 
                                                      composition=composition, 
                                                      reactant=reactant, 
                                                      mix_ratio_name=mix_ratio_name), 
                                      promotes=proms)

        # need to manually call this in this setup, because we have an element within an element
        fs_start.pyc_setup_output_ports() 

        balance = conv.add_subsystem('balance', om.BalanceComp())
        balance.add_balance('Tt', val=500.0, lower=1e-4, units='degR', desc='Total temperature', eq_units='degR')
        balance.add_balance('Pt', val=14.696, lower=1e-4, units='psi', desc='Total pressure', eq_units='psi')
        # sub.set_order(['fs','balance'])

        newton = conv.nonlinear_solver = om.NewtonSolver()
        newton.options['atol'] = 1e-10
        newton.options['rtol'] = 1e-10
        newton.options['maxiter'] = 10
        newton.options['iprint'] = -1
        newton.options['solve_subsystems'] = True
        newton.options['reraise_child_analysiserror'] = False
        newton.linesearch = om.BoundsEnforceLS()
        newton.linesearch.options['bound_enforcement'] = 'scalar'

        newton.linesearch.options['iprint'] = -1
        # newton.linesearch.options['solve_subsystems'] = True

        conv.linear_solver = om.DirectSolver(assemble_jac=True)

        self.connect('ambient.Ps', 'balance.rhs:Pt')
        self.connect('ambient.Ts', 'balance.rhs:Tt')

        self.connect('balance.Pt', 'fs.P')
        self.connect('balance.Tt', 'fs.T')

        self.connect('Fl_O:stat:P', 'balance.lhs:Pt')
        self.connect('Fl_O:stat:T', 'balance.lhs:Tt')

        # self.set_order(['ambient', 'subgroup'])

        super().setup()


if __name__ == "__main__":

  p1 = om.Problem()
  p1.model = om.Group()
  
  des_vars = p1.model.add_subsystem('des_vars', om.IndepVarComp())
  des_vars.add_output('W', 0.0, units='lbm/s')
  des_vars.add_output('alt', 1., units='ft')
  des_vars.add_output('MN', 0.5)
  des_vars.add_output('dTs', 0.0, units='degR')
  
  
  fc = p1.model.add_subsystem("fc", FlightConditions(thermo_data=thermo_data.wet_air))
  
  p1.model.connect('des_vars.W', 'fc.W')
  p1.model.connect('des_vars.alt', 'fc.alt')
  p1.model.connect('des_vars.MN', 'fc.MN')
  p1.model.connect('des_vars.dTs', 'fc.dTs')
  1.set_solver_print(level=-1)
  
  p1.setup()
  
  # p1.root.list_connections()
  
  altitudes = np.arange(17000, 21000, 50)
  
  ps = []
  
  for alt in altitudes:
          
      
      p1['des_vars.alt'] = alt
      p1['des_vars.MN'] = 0.5
      p1['des_vars.dTs'] = 0.0
      p1['des_vars.W'] = 1.0
      
      p1.run_model()
      
      print(f'alt: {alt}, Ps_set: {p1['fc.Fl_O:stat:P']}')
      ps.append(p1['fc.Fl_O:stat:P'][0])
  
  import matplotlib.pyplot as plt
  plt. plot(altitudes, ps)
    

RESULTS confirming there is something going on in the 19000 ft band:

starting setup....
alt: 17000, Ps_set: [7.64664999]
alt: 17050, Ps_set: [7.6307941]
alt: 17100, Ps_set: [7.61500891]
alt: 17150, Ps_set: [7.59929084]
alt: 17200, Ps_set: [7.58363633]
alt: 17250, Ps_set: [7.56804178]
alt: 17300, Ps_set: [7.55250364]
alt: 17350, Ps_set: [7.53701833]
alt: 17400, Ps_set: [7.52158227]
alt: 17450, Ps_set: [7.50619188]
alt: 17500, Ps_set: [7.4908436]
alt: 17550, Ps_set: [7.47553385]
alt: 17600, Ps_set: [7.46025905]
alt: 17650, Ps_set: [7.44501562]
alt: 17700, Ps_set: [7.4298]
alt: 17750, Ps_set: [7.41460861]
alt: 17800, Ps_set: [7.39943788]
alt: 17850, Ps_set: [7.38428422]
alt: 17900, Ps_set: [7.36914407]
alt: 17950, Ps_set: [7.35401386]
alt: 18000, Ps_set: [7.33888999]
alt: 18050, Ps_set: [7.32647053]
alt: 18100, Ps_set: [7.31908477]
alt: 18150, Ps_set: [7.31617604]
alt: 18200, Ps_set: [7.31718768]
alt: 18250, Ps_set: [7.32156301]
alt: 18300, Ps_set: [7.32874536]
alt: 18350, Ps_set: [7.33817807]
alt: 18400, Ps_set: [7.34930446]
alt: 18450, Ps_set: [7.36156785]
alt: 18500, Ps_set: [7.37441159]
alt: 18550, Ps_set: [7.38727899]
alt: 18600, Ps_set: [7.39961339]
alt: 18650, Ps_set: [7.41085812]
alt: 18700, Ps_set: [7.4204565]
alt: 18750, Ps_set: [7.42785187]
alt: 18800, Ps_set: [7.43248755]
alt: 18850, Ps_set: [7.43380687]
alt: 18900, Ps_set: [7.43125317]
alt: 18950, Ps_set: [7.42426977]
alt: 19000, Ps_set: [7.41229999]
alt: 19050, Ps_set: [7.39497502]
alt: 19100, Ps_set: [7.37266803]
alt: 19150, Ps_set: [7.34593534]
alt: 19200, Ps_set: [7.31533329]
alt: 19250, Ps_set: [7.2814182]
alt: 19300, Ps_set: [7.24474639]
alt: 19350, Ps_set: [7.2058742]
alt: 19400, Ps_set: [7.16535796]
alt: 19450, Ps_set: [7.12375398]
alt: 19500, Ps_set: [7.08161861]
alt: 19550, Ps_set: [7.03950817]
alt: 19600, Ps_set: [6.99797898]
alt: 19650, Ps_set: [6.95758737]
alt: 19700, Ps_set: [6.91888967]
alt: 19750, Ps_set: [6.88244222]
alt: 19800, Ps_set: [6.84880133]
alt: 19850, Ps_set: [6.81852333]
alt: 19900, Ps_set: [6.79216456]
alt: 19950, Ps_set: [6.77028134]
alt: 20000, Ps_set: [6.75342999]
alt: 20050, Ps_set: [6.73928174]
alt: 20100, Ps_set: [6.72513907]
alt: 20150, Ps_set: [6.71100524]
alt: 20200, Ps_set: [6.69688348]
alt: 20250, Ps_set: [6.68277706]
alt: 20300, Ps_set: [6.66868922]
alt: 20350, Ps_set: [6.65462321]
alt: 20400, Ps_set: [6.64058228]
alt: 20450, Ps_set: [6.62656968]
alt: 20500, Ps_set: [6.61258865]
alt: 20550, Ps_set: [6.59864246]
alt: 20600, Ps_set: [6.58473433]
alt: 20650, Ps_set: [6.57086754]
alt: 20700, Ps_set: [6.55704531]
alt: 20750, Ps_set: [6.54327092]
alt: 20800, Ps_set: [6.52954759]
alt: 20850, Ps_set: [6.51587858]
alt: 20900, Ps_set: [6.50226715]
alt: 20950, Ps_set: [6.48871654]

I also ran with set_solver_print(level=2) and I could not find any non-convergence messages.

[flight-test-engineering]

I found the bug. It is in pycycle/elements/US1976.py

The tabulated pressure for 19000ft is currently 7.4123, whereas it should be 7.04123.