Recuperated Brayton Cycle

Introduction

As shown in the Brayton Cycle modelingguide, and the given input files closed_brayton_cycle.i and open_brayton_cycle.i, a Brayton Cycle Power Conversion Unit (PCU) consists of a motor, compressor, turbine, and generator all coupled by a single shaft.

Detailed descriptions of the compressor and turbine components used in this example can be found in the Brayton Cycle modeling guide. In the aforementioned example, a simplified startup transient with a simplified heat source was conducted which demonstrated the Thermal Hydraulics module’s capability to produce torque, power, mass flow rate, and pressure ratios for all PCU components based on the shaft speed.

In this example a different heat source and recuperator were added to the same open Brayton Cycle PCU to demonstrate a more nuanced model along with a more realistic piping structure. Also, a Proportional-Integral-Derivative (PID) controller was added to the motor to conduct a more realistic startup transient of the system.

A heat rate of 105.75 kW is applied as a uniform volumetric source to a 10-m-long steel tube with an outer radius of 15 cm. Gas which has been compressed in the compressor, and preheated in the recuperator, passes over the exterior of the heat structure and removes heat via forced convection.

The recuperator is represented by an annular cylindrical heat structure of the same material and width as the main heat source but only 5 m in length and without internal heat generation. Instead, hot exhaust gas leaving the turbine transfers heat to the inside of the recuperator, and cooler gas from the compressor outlet removes heat from the outer surface. This preheats the gas entering the main heat source and improves thermal efficiency of the cycle.

Transient Description

Initially the shaft system is at rest, and the working fluid and heat structures are the ambient temperature and pressure. At t = 0 s the motor PID activates and ramps the PCU shaft to approximately 87,000 RPM over the course of 1600 s. When t = 1000 s, heat generation in the main heat structure is activated and linearly increases from 0 kW to a maximum power of 105.75 kW by t = 8600 s.

As the working fluid begins to heat, torque supplied to the shaft from the turbine increases. The motor PID control logic slowly decreases the torque supplied by the motor to 0 N·m once the turbine torque is greater than the torque supplied by the motor. This allows for a graceful transition of PCU control from the motor PID to the working fluid and heat source, and prevents drastically over-speeding the system above the rated speed of 96,000 RPM.

Input File description

The recuperated Brayton Cycle example is executed with the input file recuperated_brayton_cycle.i. Shown below in Figure 1 is a diagram of the modeled system.

Figure 1: Diagram of the recuperated Brayton Cycle.

All 90° pipe bends are modeled using VolumeJunction1Phase as shown in the example below:

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [junction2_cold_leg]
    type = VolumeJunction1Phase<<<{"description": "Junction between 1-phase flow channels that has a non-zero volume", "href": "../../../../source/components/VolumeJunction1Phase.html"}>>>
    connections<<<{"description": "Junction connections"}>>> = 'pipe2:out cold_leg:in'
    position<<<{"description": "Spatial position of the center of the junction [m]"}>>> = '${x3} ${y3} 0'
    volume<<<{"description": "Volume of the junction [m^3]"}>>> = '${fparse A2*0.1}'
  []
[]

The main heat source is modeled as a heat structure with internal heat generation and forced convection heat transfer to Pipe 4.

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [reactor]
    type = HeatStructureCylindrical<<<{"description": "Cylindrical heat structure", "href": "../../../../source/components/HeatStructureCylindrical.html"}>>>
    orientation<<<{"description": "Direction of axis from start position to end position (no need to normalize)"}>>> = '-1 0 0'
    position<<<{"description": "Start position of axis in 3-D space [m]"}>>> = '${x4} ${y4} 0'
    length<<<{"description": "Length of each axial section [m]"}>>> = ${L4}
    widths<<<{"description": "Width of each radial region [m]"}>>> = 0.15
    n_elems<<<{"description": "Number of elements in each axial section"}>>> = ${n_elems4}
    n_part_elems<<<{"description": "Number of elements of each radial region"}>>> = 2
    names<<<{"description": "Name of each radial region"}>>> = core
    solid_properties<<<{"description": "Solid properties object name for each radial region"}>>> = steel
    solid_properties_T_ref<<<{"description": "Density reference temperatures for each radial region. This is required if 'solid_properties' is provided. The density in each region will be a constant value computed by evaluating the density function at the reference temperature."}>>> = '300'
  []
[]

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [total_power]
    type = TotalPower<<<{"description": "Prescribes total power via a user supplied value", "href": "../../../../source/components/TotalPower.html"}>>>
    power<<<{"description": "Total power [W]"}>>> = 0
  []
[]

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [heat_generation]
    type = HeatSourceFromTotalPower<<<{"description": "Heat generation from total power", "href": "../../../../source/components/HeatSourceFromTotalPower.html"}>>>
    power<<<{"description": "Component that provides total power"}>>> = total_power
    hs<<<{"description": "Heat structure in which to apply heat source"}>>> = reactor
    regions<<<{"description": "Heat structure regions where heat generation is to be applied"}>>> = core
  []
[]

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [heat_transfer]
    type = HeatTransferFromHeatStructure1Phase<<<{"description": "Connects a 1-phase flow channel and a heat structure", "href": "../../../../source/components/HeatTransferFromHeatStructure1Phase.html"}>>>
    flow_channel<<<{"description": "Name of flow channel component to connect to"}>>> = pipe4
    hs<<<{"description": "Heat structure name"}>>> = reactor
    hs_side<<<{"description": "Heat structure side"}>>> = OUTER
    Hw<<<{"description": "Convective heat transfer coefficient [W/(m^2-K)]"}>>> = 10000
  []
[]

The Recuperator is modeled as a heat structure without internal heat generation and heat transfer to both cold_leg and hot_leg.

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [recuperator]
    type = HeatStructureCylindrical<<<{"description": "Cylindrical heat structure", "href": "../../../../source/components/HeatStructureCylindrical.html"}>>>
    orientation<<<{"description": "Direction of axis from start position to end position (no need to normalize)"}>>> = '0 -1 0'
    position<<<{"description": "Start position of axis in 3-D space [m]"}>>> = '${x3} ${y3} 0'
    length<<<{"description": "Length of each axial section [m]"}>>> = '${fparse L3/2}'
    widths<<<{"description": "Width of each radial region [m]"}>>> = ${recuperator_width}
    n_elems<<<{"description": "Number of elements in each axial section"}>>> = '${fparse n_elems3/2}'
    n_part_elems<<<{"description": "Number of elements of each radial region"}>>> = 2
    names<<<{"description": "Name of each radial region"}>>> = recuperator
    solid_properties<<<{"description": "Solid properties object name for each radial region"}>>> = steel
    solid_properties_T_ref<<<{"description": "Density reference temperatures for each radial region. This is required if 'solid_properties' is provided. The density in each region will be a constant value computed by evaluating the density function at the reference temperature."}>>> = '300'
    inner_radius<<<{"description": "Inner radius of the heat structure [m]"}>>> = ${D1}
  []
[]

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [heat_transfer_cold_leg]
    type = HeatTransferFromHeatStructure1Phase<<<{"description": "Connects a 1-phase flow channel and a heat structure", "href": "../../../../source/components/HeatTransferFromHeatStructure1Phase.html"}>>>
    flow_channel<<<{"description": "Name of flow channel component to connect to"}>>> = cold_leg
    hs<<<{"description": "Heat structure name"}>>> = recuperator
    hs_side<<<{"description": "Heat structure side"}>>> = OUTER
    Hw<<<{"description": "Convective heat transfer coefficient [W/(m^2-K)]"}>>> = 10000
  []
[]

[Components<<<{"href": "../../../../syntax/Components/index.html"}>>>]
  [heat_transfer_hot_leg]
    type = HeatTransferFromHeatStructure1Phase<<<{"description": "Connects a 1-phase flow channel and a heat structure", "href": "../../../../source/components/HeatTransferFromHeatStructure1Phase.html"}>>>
    flow_channel<<<{"description": "Name of flow channel component to connect to"}>>> = hot_leg
    hs<<<{"description": "Heat structure name"}>>> = recuperator
    hs_side<<<{"description": "Heat structure side"}>>> = INNER
    Hw<<<{"description": "Convective heat transfer coefficient [W/(m^2-K)]"}>>> = 10000
  []
[]

Control Logic

One of the main changes from the Brayton Cycle modeling guide is the addition of a PID controller to the motor for the start up transient. This PID control operates a 3-Phase electric motor which applies torque to the shaft to reach a desired speed set by the user, in this case 87,000 RPM. A speed lower than the rated turbine speed of 96,000 RPM was chosen to prevent significant over-speeding of the turbine during the start up transient.

As the shaft increases in speed, the compressor begins to compress the inlet atmospheric air and produces a coolant flow through the system. Once the heat source begins heating the coolant, the gas begins conducting work on the system via the turbine. However, until a fully developed and fully heated flow is established, the working fluid still requires an assist from the motor to maintain shaft speed and adequate coolant mass flow rate. At this point, a direct shutoff of the motor would cause the system to stall. Conversely, if the motor were to supply a constant torque the shaft would over speed due to the increasing torque supplied by the turbine from the heating coolant.

To prevent stalling or significant over-speeding of the shaft while the coolant is reaching a fully developed and fully heated flow, a shutdown function is applied to the motor. A linearly decreasing torque is applied to the motor starting from the time the turbine first provides more torque and reaches 0 N·m after 35,000 s. This function uses ControlLogic to determine when the torque supplied by the turbine is greater than the motor, and initiates a graceful ramp down of the motor. The proceeding sections discuss in detail how this control logic is implemented in the input file.

PID Control

Shown below is the ControlLogic block containing the PID block utilizing PIDControl. Here, the $var element = document.getElementById("moose-equation-9ba4963b-1633-4d66-b305-6efe12bc41c5");katex.render("K_p", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , $var element = document.getElementById("moose-equation-89ff1977-9528-4593-b794-8205df6a98de");katex.render("K_i", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , and $var element = document.getElementById("moose-equation-7a665359-76cd-4739-bef0-93905d7479e6");katex.render("K_d", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ variables are the proportional, integral, and derivative gains respectively. As with any PID control, gains must be determined for a specific system and are not universal. In this specific PID, the output signal is used to govern the torque supplied by the motor to the system shaft based upon the current shaft speed.

[ControlLogic<<<{"href": "../../../../syntax/ControlLogic/index.html"}>>>]
  [initial_motor_PID]
    type = PIDControl<<<{"description": "Declares a control data named 'output' and uses Proportional Integral Derivative logic on the 'value' control data to set it.", "href": "../../../../source/controllogic/PIDControl.html"}>>>
    set_point<<<{"description": "The name of the control data with the set point."}>>> = set_point:value
    input<<<{"description": "The name of the control data that we read in."}>>> = shaft_RPM
    initial_value<<<{"description": "The initial value for the integral part."}>>> = 0
    K_p<<<{"description": "The coefficient for the proportional term."}>>> = 0.0011
    K_i<<<{"description": "The coefficient for the integral term."}>>> = 0.00000004
    K_d<<<{"description": "The coefficient for the derivative term."}>>> = 0
  []
[]

The PID attempts to reach a desired shaft speed, in RPM, which is determined by the set_point block using GetFunctionValueControl. The rated speed of the turbine is set at the beginning of the input file as a global parameter speed_rated_rpm. To prevent over-speeding of the turbine and compressor, the set point of the PID is adjusted to be lower than the rated turbine speed.

[ControlLogic<<<{"href": "../../../../syntax/ControlLogic/index.html"}>>>]
  [set_point]
    type = GetFunctionValueControl<<<{"description": "Sets a ControlData named 'value' with the value of a function", "href": "../../../../source/controllogic/GetFunctionValueControl.html"}>>>
    function<<<{"description": "The name of the function prescribing a value."}>>> = '${fparse speed_rated_rpm - 9000}'
  []
[]

A ParsedFunctionControl is applied to determine whether to send the PID output to the motor or shutdown the motor because the turbine is supplying enough torque.

[ControlLogic<<<{"href": "../../../../syntax/ControlLogic/index.html"}>>>]
  [logic]
    type = ParsedFunctionControl<<<{"description": "Control that evaluates a parsed function", "href": "../../../../source/controllogic/ParsedFunctionControl.html"}>>>
    function<<<{"description": "The function to be evaluated by this control."}>>> = 'if(motor+0.5 > turb, PID, shutdown_fn)'
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'motor turb PID shutdown_fn'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'motor_torque turbine_torque initial_motor_PID:output motor_torque_fn_shutdown'
  []
[]

If logic determines to send the PID signal to the motor, then SetComponentRealValueControl motor_PID is used to apply the signal directly to the motor component torque.

[ControlLogic<<<{"href": "../../../../syntax/ControlLogic/index.html"}>>>]
  [motor_PID]
    type = SetComponentRealValueControl<<<{"description": "Control to set a floating point (Real) value of a component parameter with control data boolean", "href": "../../../../source/controllogic/SetComponentRealValueControl.html"}>>>
    component<<<{"description": "The name of the component to be controlled."}>>> = motor
    parameter<<<{"description": "The name of the parameter in the component to be controlled."}>>> = torque
    value<<<{"description": "The name of control data to be set in the component."}>>> = logic:value
  []
[]

Motor Shutdown

Should the turbine provide more torque than the motor, the logic activates a ParsedFunction, motor_torque_fn_shutdown, to ramp down the motor allowing for a graceful transition of the system from the motor to the turbine.

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [motor_torque_fn_shutdown]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'PID_trip_status time_trip'
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'PID_trip_status time_trip'
    expression<<<{"description": "The user defined function."}>>> = 'if(PID_trip_status = 1, max(2.4 - (2.4 * ((t - time_trip) / 35000)),0.0), 1)'
  []
[]

The shutdown function utilizes AuxVariables values PID_trip_status and time_trip. These values are dependent upon a number of factors including the motor PID, startup rate of the heat source, heat source total power, etc. The user has no way of knowing these values while writing the input file, so they must be determined during input file execution. The following steps outline one method to determine, store, and apply initially unknown variables within an input file.

AuxVariables

First, AuxVariables are used to create two variables, time_trip and PID_trip_status. time_trip is used to store the time at which turbine torque is greater than the motor torque. PID_trip_status is used to change a trip status from 0 to 1.

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [time_trip]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = SCALAR
  []
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [PID_trip_status]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = SCALAR
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
[]

Functions

Next, four Functions are created which are used to output the desired values for the previously mentioned AuxVariables. These functions are PID_tripped_status_fn, PID_tripped_constant_value, time_fn, and is_tripped_fn.

time_fn records the current simulation time.

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [time_fn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = t
  []
[]

is_tripped_fn determines if the turbine torque is greater than the motor torque using the Postprocessors values turbine_torque and motor_torque.

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [is_tripped_fn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'motor_torque turbine_torque'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'motor_torque turbine_torque'
    expression<<<{"description": "The user defined function."}>>> = 'turbine_torque > motor_torque'
  []
[]

PID_tripped_constant_value creates a constant function with a value of 1. This is later used to change the PID trip status from 0 to 1. 0 indicates the motor is supplying more torque than the turbine, 1 indicates the turbine is supplying more torque than the motor.

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [PID_tripped_constant_value]
    type = ConstantFunction<<<{"description": "A function that returns a constant value as defined by an input parameter.", "href": "../../../../source/functions/ConstantFunction.html"}>>>
    value<<<{"description": "The constant value"}>>> = 1
  []
[]

PID_tripped_status_fn stores the value from the Control PID_trip_status (which will be discussed next) into a separate variable of the same name.

[Functions<<<{"href": "../../../../syntax/Functions/index.html"}>>>]
  [PID_tripped_status_fn]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../../source/functions/MooseParsedFunction.html"}>>>
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'PID_trip_status'
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'PID_trip_status'
    expression<<<{"description": "The user defined function."}>>> = 'PID_trip_status'
  []
[]

Controls

Multiple variables and functions have been discussed which are now used by Controls and AuxScalarKernels to actually record the desired values. Two Controls use the previously mentioned functions to activate two AuxScalarKernels which then apply the desired values to the previously mentioned AuxVariables.

PID_trip_status is a ConditionalFunctionEnableControl within the Controls block. When the function is_tripped_fn becomes true, PID_trip_status activates the AuxScalarKernel set_PID_tripped.

[Controls<<<{"href": "../../../../syntax/Controls/index.html"}>>>]
  [PID_trip_status]
    type = ConditionalFunctionEnableControl<<<{"description": "Control for enabling/disabling objects when a function value is true", "href": "../../../../source/controls/ConditionalFunctionEnableControl.html"}>>>
    conditional_function<<<{"description": "The function to give a true or false value"}>>> = is_tripped_fn
    enable_objects<<<{"description": "A list of object tags to enable."}>>> = 'AuxScalarKernels::PID_trip_status_aux'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_END'
  []
[]

time_PID is also a ConditionalFunctionEnableControl which activates when PID_trip_status activates. This control then disables the AuxScalarKernel set_time_PID which was recording the time up until the point of the PID trip.

[Controls<<<{"href": "../../../../syntax/Controls/index.html"}>>>]
  [time_PID]
    type = ConditionalFunctionEnableControl<<<{"description": "Control for enabling/disabling objects when a function value is true", "href": "../../../../source/controls/ConditionalFunctionEnableControl.html"}>>>
    conditional_function<<<{"description": "The function to give a true or false value"}>>> = PID_tripped_status_fn
    disable_objects<<<{"description": "A list of object tags to disable."}>>> = 'AuxScalarKernels::time_trip_aux'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_END'
  []
[]

AuxScalarKernels

Finally, the AuxScalarKernels apply all of the values and conditions discussed up to this point.

time_trip_aux takes the function time_fn and applies it to the AuxVariable time_trip. When the time_PID Control activates, it disables this AuxScalarVariable which results in time_trip remaining constant at the time the PID trip occurred. This value can now be applied to motor_torque_fn_shutdown.

[AuxScalarKernels<<<{"href": "../../../../syntax/AuxScalarKernels/index.html"}>>>]
  [time_trip_aux]
    type = FunctionScalarAux<<<{"description": "Sets a value of a scalar variable based on a function.", "href": "../../../../source/auxscalarkernels/FunctionScalarAux.html"}>>>
    function<<<{"description": "The functions to set the scalar variable components."}>>> = time_fn
    variable<<<{"description": "The name of the variable that this kernel operates on"}>>> = time_trip
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_END'
  []
[]

PID_trip_status_aux, once activated, takes the function PID_tripped_constant_value and applies it to the AuxVariable PID_trip_status. This causes a simple switch from 0 to 1 meaning the turbine went from supplying less torque than the motor to supplying more torque than the motor.

[AuxScalarKernels<<<{"href": "../../../../syntax/AuxScalarKernels/index.html"}>>>]
  [PID_trip_status_aux]
    type = FunctionScalarAux<<<{"description": "Sets a value of a scalar variable based on a function.", "href": "../../../../source/auxscalarkernels/FunctionScalarAux.html"}>>>
    function<<<{"description": "The functions to set the scalar variable components."}>>> = PID_tripped_constant_value
    variable<<<{"description": "The name of the variable that this kernel operates on"}>>> = PID_trip_status
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'TIMESTEP_END'
    enable<<<{"description": "Set the enabled status of the MooseObject."}>>> = false
  []
[]

Results

Figure 2 shows the first 2000 seconds of the transient where the PID controlled motor increased shaft speed from 0 to approximately 85,000 RPM. This is followed by Figure 3 which displays the shaft speed over the course of the entire transient.

Figure 2: PID start-up of shaft system.

Figure 3: Shaft speed during transient.

As shown in Figure 3, shaft speed is quickly ramped up by the PID and then linearly increases as the working fluid is heated and begins working on the turbine to produce shaft torque. At approximately 20,000 seconds the motor and turbine provide the same amount of torque which initiates the motor shutdown function. A comparison of the motor and turbine torques is shown below in Figure 4.

Figure 4: Comparison of the motor and turbine torques.

Figure 5 visualizes the pressure ratios of both the compressor and turbine during the transient. A slight difference is seen between the compressor and turbine steady state pressure ratios due to frictional losses across the piping system and 90° pipe bends, modeled by VolumeJunction1Phase.

Figure 5: Compressor and turbine pressure ratios.

Finally, coolant temperatures across key points of the system are displayed in Figure 6. The reactor inlet and outlet temperatures were omitted as there was no difference between the cold_leg outlet and turbine inlet respectively.

Also, the compressor inlet was not plotted as it remained constant at ambient temperature of 300 K. The compressor outlet temperature sees an immediate increase corresponding to the startup of the PID system while the rest of the temperature responses are dependent upon the heat source power output.

Figure 6: Coolant temperatures across key components.

(modules/thermal_hydraulics/test/tests/problems/brayton_cycle/closed_brayton_cycle.i)

# This input file is used to demonstrate a simple closed, air Brayton cycle using
# a compressor, turbine, shaft, motor, and generator.
# The flow length is divided into 6 segments as illustrated below, where
#   - "(C)" denotes the compressor
#   - "(T)" denotes the turbine
#   - "*" denotes a fictitious junction
#
#                Heated section               Cooled section
# *-----(C)-----*--------------*-----(T)-----*--------------*
#    1       2         3          4       5         6
#
# Initially the fluid is at rest at ambient conditions, the shaft speed is zero,
# and no heat transfer occurs with the system.
# The transient is controlled as follows:
#   * 0   - 100 s: motor ramps up torque linearly from zero
#   * 100 - 200 s: motor ramps down torque linearly to zero, HTC ramps up linearly from zero.
#   * 200 - 300 s: (no changes; should approach steady condition)

I_motor = 1.0
motor_torque_max = 400.0

I_generator = 1.0
generator_torque_per_shaft_speed = -0.00025

motor_ramp_up_duration = 100.0
motor_ramp_down_duration = 100.0
post_motor_time = 100.0
t1 = ${motor_ramp_up_duration}
t2 = ${fparse t1 + motor_ramp_down_duration}
t3 = ${fparse t2 + post_motor_time}

D1 = 0.15
D2 = ${D1}
D3 = ${D1}
D4 = ${D1}
D5 = ${D1}
D6 = ${D1}

A1 = ${fparse 0.25 * pi * D1^2}
A2 = ${fparse 0.25 * pi * D2^2}
A3 = ${fparse 0.25 * pi * D3^2}
A4 = ${fparse 0.25 * pi * D4^2}
A5 = ${fparse 0.25 * pi * D5^2}
A6 = ${fparse 0.25 * pi * D6^2}

L1 = 10.0
L2 = ${L1}
L3 = ${L1}
L4 = ${L1}
L5 = ${L1}
L6 = ${L1}

x1 = 0.0
x2 = ${fparse x1 + L1}
x3 = ${fparse x2 + L2}
x4 = ${fparse x3 + L3}
x5 = ${fparse x4 + L4}
x6 = ${fparse x5 + L5}

x2_minus = ${fparse x2 - 0.001}
x2_plus = ${fparse x2 + 0.001}
x5_minus = ${fparse x5 - 0.001}
x5_plus = ${fparse x5 + 0.001}

n_elems1 = 10
n_elems2 = ${n_elems1}
n_elems3 = ${n_elems1}
n_elems4 = ${n_elems1}
n_elems5 = ${n_elems1}
n_elems6 = ${n_elems1}

A_ref_comp = ${fparse 0.5 * (A1 + A2)}
V_comp = ${fparse A_ref_comp * 1.0}
I_comp = 1.0

A_ref_turb = ${fparse 0.5 * (A4 + A5)}
V_turb = ${fparse A_ref_turb * 1.0}
I_turb = 1.0

c0_rated_comp = 351.6925137
rho0_rated_comp = 1.146881112

rated_mfr = 0.25

speed_rated_rpm = 96000
speed_rated = ${fparse speed_rated_rpm * 2 * pi / 60.0}

speed_initial = 0

eff_comp = 0.79
eff_turb = 0.843

T_hot = 1000
T_cold = 300
T_ambient = 300
p_ambient = 1e5

[GlobalParams]
  orientation = '1 0 0'
  gravity_vector = '0 0 0'

  initial_p = ${p_ambient}
  initial_T = ${T_ambient}
  initial_vel = 0
  initial_vel_x = 0
  initial_vel_y = 0
  initial_vel_z = 0

  fp = fp_air
  closures = closures
  f = 0

  scaling_factor_1phase = '1 1 1e-5'
  scaling_factor_rhoV = 1
  scaling_factor_rhouV = 1
  scaling_factor_rhovV = 1
  scaling_factor_rhowV = 1
  scaling_factor_rhoEV = 1e-5

  rdg_slope_reconstruction = none
[]

[Functions]
  [motor_torque_fn]
    type = PiecewiseLinear
    x = '0 ${t1} ${t2}'
    y = '0 ${motor_torque_max} 0'
  []
  [motor_power_fn]
    type = ParsedFunction
    expression = 'torque * speed'
    symbol_names = 'torque speed'
    symbol_values = 'motor_torque shaft:omega'
  []
  [generator_torque_fn]
    type = ParsedFunction
    expression = 'slope * t'
    symbol_names = 'slope'
    symbol_values = '${generator_torque_per_shaft_speed}'
  []
  [generator_power_fn]
    type = ParsedFunction
    expression = 'torque * speed'
    symbol_names = 'torque speed'
    symbol_values = 'generator_torque shaft:omega'
  []
  [htc_wall_fn]
    type = PiecewiseLinear
    x = '0 ${t1} ${t2}'
    y = '0 0 1e3'
  []
[]

[FluidProperties]
  [fp_air]
    type = IdealGasFluidProperties
    emit_on_nan = none
  []
[]

[Closures]
  [closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [shaft]
    type = Shaft
    connected_components = 'motor compressor turbine generator'
    initial_speed = ${speed_initial}
  []
  [motor]
    type = ShaftConnectedMotor
    inertia = ${I_motor}
    torque = 0 # controlled
  []
  [generator]
    type = ShaftConnectedMotor
    inertia = ${I_generator}
    torque = generator_torque_fn
  []

  [pipe1]
    type = FlowChannel1Phase
    position = '${x1} 0 0'
    length = ${L1}
    n_elems = ${n_elems1}
    A = ${A1}
  []
  [compressor]
    type = ShaftConnectedCompressor1Phase
    position = '${x2} 0 0'
    inlet = 'pipe1:out'
    outlet = 'pipe2:in'
    A_ref = ${A_ref_comp}
    volume = ${V_comp}

    omega_rated = ${speed_rated}
    mdot_rated = ${rated_mfr}
    c0_rated = ${c0_rated_comp}
    rho0_rated = ${rho0_rated_comp}

    speeds = '0.5208 0.6250 0.7292 0.8333 0.9375'
    Rp_functions = 'rp_comp1 rp_comp2 rp_comp3 rp_comp4 rp_comp5'
    eff_functions = 'eff_comp1 eff_comp2 eff_comp3 eff_comp4 eff_comp5'

    min_pressure_ratio = 1.0

    speed_cr_I = 0
    inertia_const = ${I_comp}
    inertia_coeff = '${I_comp} 0 0 0'

    # assume no shaft friction
    speed_cr_fr = 0
    tau_fr_const = 0
    tau_fr_coeff = '0 0 0 0'
  []
  [pipe2]
    type = FlowChannel1Phase
    position = '${x2} 0 0'
    length = ${L2}
    n_elems = ${n_elems2}
    A = ${A2}
  []
  [junction2_3]
    type = JunctionOneToOne1Phase
    connections = 'pipe2:out pipe3:in'
  []
  [pipe3]
    type = FlowChannel1Phase
    position = '${x3} 0 0'
    length = ${L3}
    n_elems = ${n_elems3}
    A = ${A3}
  []
  [junction3_4]
    type = JunctionOneToOne1Phase
    connections = 'pipe3:out pipe4:in'
  []
  [pipe4]
    type = FlowChannel1Phase
    position = '${x4} 0 0'
    length = ${L4}
    n_elems = ${n_elems4}
    A = ${A4}
  []
  [turbine]
    type = ShaftConnectedCompressor1Phase
    position = '${x5} 0 0'
    inlet = 'pipe4:out'
    outlet = 'pipe5:in'
    A_ref = ${A_ref_turb}
    volume = ${V_turb}

    treat_as_turbine = true

    omega_rated = ${speed_rated}
    mdot_rated = ${rated_mfr}
    c0_rated = ${c0_rated_comp}
    rho0_rated = ${rho0_rated_comp}

    speeds = '0 0.5208 0.6250 0.7292 0.8333 0.9375'
    Rp_functions = 'rp_turb0 rp_turb1 rp_turb2 rp_turb3 rp_turb4 rp_turb5'
    eff_functions = 'eff_turb1 eff_turb1 eff_turb2 eff_turb3 eff_turb4 eff_turb5'

    min_pressure_ratio = 1.0

    speed_cr_I = 0
    inertia_const = ${I_turb}
    inertia_coeff = '${I_turb} 0 0 0'

    # assume no shaft friction
    speed_cr_fr = 0
    tau_fr_const = 0
    tau_fr_coeff = '0 0 0 0'
  []
  [pipe5]
    type = FlowChannel1Phase
    position = '${x5} 0 0'
    length = ${L5}
    n_elems = ${n_elems5}
    A = ${A5}
  []
  [junction5_6]
    type = JunctionOneToOne1Phase
    connections = 'pipe5:out pipe6:in'
  []
  [pipe6]
    type = FlowChannel1Phase
    position = '${x6} 0 0'
    length = ${L6}
    n_elems = ${n_elems6}
    A = ${A6}
  []
  [junction6_1]
    type = JunctionOneToOne1Phase
    connections = 'pipe6:out pipe1:in'
  []

  [heating]
    type = HeatTransferFromSpecifiedTemperature1Phase
    flow_channel = pipe3
    T_wall = ${T_hot}
    Hw = htc_wall_fn
  []
  [cooling]
    type = HeatTransferFromSpecifiedTemperature1Phase
    flow_channel = pipe6
    T_wall = ${T_cold}
    Hw = htc_wall_fn
  []
[]

[ControlLogic]
  [motor_ctrl]
    type = TimeFunctionComponentControl
    component = motor
    parameter = torque
    function = motor_torque_fn
  []
[]

[Postprocessors]
  [heating_rate]
    type = ADHeatRateConvection1Phase
    block = 'pipe3'
    T = T
    T_wall = T_wall
    Hw = Hw
    P_hf = P_hf
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [cooling_rate]
    type = ADHeatRateConvection1Phase
    block = 'pipe6'
    T = T
    T_wall = T_wall
    Hw = Hw
    P_hf = P_hf
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [motor_torque]
    type = RealComponentParameterValuePostprocessor
    component = motor
    parameter = torque
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [motor_power]
    type = FunctionValuePostprocessor
    function = motor_power_fn
    execute_on = 'INITIAL TIMESTEP_END'
    indirect_dependencies = 'motor_torque shaft:omega'
  []

  [generator_torque]
    type = ShaftConnectedComponentPostprocessor
    quantity = torque
    shaft_connected_component_uo = generator:shaftconnected_uo
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [generator_power]
    type = FunctionValuePostprocessor
    function = generator_power_fn
    execute_on = 'INITIAL TIMESTEP_END'
    indirect_dependencies = 'generator_torque shaft:omega'
  []

  [shaft_speed]
    type = ScalarVariable
    variable = 'shaft:omega'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [p_in_comp]
    type = PointValue
    variable = p
    point = '${x2_minus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_out_comp]
    type = PointValue
    variable = p
    point = '${x2_plus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_ratio_comp]
    type = ParsedPostprocessor
    pp_names = 'p_in_comp p_out_comp'
    expression = 'p_out_comp / p_in_comp'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [p_in_turb]
    type = PointValue
    variable = p
    point = '${x5_minus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_out_turb]
    type = PointValue
    variable = p
    point = '${x5_plus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_ratio_turb]
    type = ParsedPostprocessor
    pp_names = 'p_in_turb p_out_turb'
    expression = 'p_in_turb / p_out_turb'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [mfr_comp]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = mass
    junction = compressor
  []
  [mfr_turb]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe4:out
    connection_index = 0
    equation = mass
    junction = turbine
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  end_time = ${t3}
  dt = 0.1
  abort_on_solve_fail = true

  solve_type = NEWTON
  nl_rel_tol = 1e-50
  nl_abs_tol = 1e-11
  nl_max_its = 15

  l_tol = 1e-4
  l_max_its = 10
[]

[Outputs]
  [csv]
    type = CSV
    file_base = 'closed_brayton_cycle'
    execute_vector_postprocessors_on = 'INITIAL'
  []
  [console]
    type = Console
    show = 'shaft_speed p_ratio_comp p_ratio_turb compressor:pressure_ratio turbine:pressure_ratio'
  []
[]

[Functions]
  # compressor pressure ratio
  [rp_comp1]
    type = PiecewiseLinear
    data_file = 'rp_comp1.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp2]
    type = PiecewiseLinear
    data_file = 'rp_comp2.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp3]
    type = PiecewiseLinear
    data_file = 'rp_comp3.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp4]
    type = PiecewiseLinear
    data_file = 'rp_comp4.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp5]
    type = PiecewiseLinear
    data_file = 'rp_comp5.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []

  # compressor efficiency
  [eff_comp1]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp2]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp3]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp4]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp5]
    type = ConstantFunction
    value = ${eff_comp}
  []

  # turbine pressure ratio
  [rp_turb0]
    type = ConstantFunction
    value = 1
  []
  [rp_turb1]
    type = PiecewiseLinear
    data_file = 'rp_turb1.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb2]
    type = PiecewiseLinear
    data_file = 'rp_turb2.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb3]
    type = PiecewiseLinear
    data_file = 'rp_turb3.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb4]
    type = PiecewiseLinear
    data_file = 'rp_turb4.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb5]
    type = PiecewiseLinear
    data_file = 'rp_turb5.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []

  # turbine efficiency
  [eff_turb1]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb2]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb3]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb4]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb5]
    type = ConstantFunction
    value = ${eff_turb}
  []
[]

(modules/thermal_hydraulics/test/tests/problems/brayton_cycle/open_brayton_cycle.i)

# This input file is used to demonstrate a simple open-air Brayton cycle using
# a compressor, turbine, shaft, motor, and generator.
# The flow length is divided into 5 segments as illustrated below, where
#   - "(I)" denotes the inlet
#   - "(C)" denotes the compressor
#   - "(T)" denotes the turbine
#   - "(O)" denotes the outlet
#   - "*" denotes a fictitious junction
#
#                  Heated section
# (I)-----(C)-----*--------------*-----(T)-----(O)
#      1       2         3          4       5
#
# Initially the fluid is at rest at ambient conditions, the shaft speed is zero,
# and no heat transfer occurs with the system.
# The transient is controlled as follows:
#   * 0   - 100 s: motor ramps up torque linearly from zero
#   * 100 - 200 s: motor ramps down torque linearly to zero, HTC ramps up linearly from zero.
#   * 200 - 300 s: (no changes; should approach steady condition)

I_motor = 1.0
motor_torque_max = 400.0

I_generator = 1.0
generator_torque_per_shaft_speed = -0.00025

motor_ramp_up_duration = 100.0
motor_ramp_down_duration = 100.0
post_motor_time = 100.0
t1 = ${motor_ramp_up_duration}
t2 = ${fparse t1 + motor_ramp_down_duration}
t3 = ${fparse t2 + post_motor_time}

D1 = 0.15
D2 = ${D1}
D3 = ${D1}
D4 = ${D1}
D5 = ${D1}

A1 = ${fparse 0.25 * pi * D1^2}
A2 = ${fparse 0.25 * pi * D2^2}
A3 = ${fparse 0.25 * pi * D3^2}
A4 = ${fparse 0.25 * pi * D4^2}
A5 = ${fparse 0.25 * pi * D5^2}

L1 = 10.0
L2 = ${L1}
L3 = ${L1}
L4 = ${L1}
L5 = ${L1}

x1 = 0.0
x2 = ${fparse x1 + L1}
x3 = ${fparse x2 + L2}
x4 = ${fparse x3 + L3}
x5 = ${fparse x4 + L4}

x2_minus = ${fparse x2 - 0.001}
x2_plus = ${fparse x2 + 0.001}
x5_minus = ${fparse x5 - 0.001}
x5_plus = ${fparse x5 + 0.001}

n_elems1 = 10
n_elems2 = ${n_elems1}
n_elems3 = ${n_elems1}
n_elems4 = ${n_elems1}
n_elems5 = ${n_elems1}

A_ref_comp = ${fparse 0.5 * (A1 + A2)}
V_comp = ${fparse A_ref_comp * 1.0}
I_comp = 1.0

A_ref_turb = ${fparse 0.5 * (A4 + A5)}
V_turb = ${fparse A_ref_turb * 1.0}
I_turb = 1.0

c0_rated_comp = 351.6925137
rho0_rated_comp = 1.146881112

rated_mfr = 0.25

speed_rated_rpm = 96000
speed_rated = ${fparse speed_rated_rpm * 2 * pi / 60.0}

speed_initial = 0

eff_comp = 0.79
eff_turb = 0.843

T_hot = 1000
T_ambient = 300
p_ambient = 1e5

[GlobalParams]
  orientation = '1 0 0'
  gravity_vector = '0 0 0'

  initial_p = ${p_ambient}
  initial_T = ${T_ambient}
  initial_vel = 0
  initial_vel_x = 0
  initial_vel_y = 0
  initial_vel_z = 0

  fp = fp_air
  closures = closures
  f = 0

  scaling_factor_1phase = '1 1 1e-5'
  scaling_factor_rhoV = 1
  scaling_factor_rhouV = 1
  scaling_factor_rhovV = 1
  scaling_factor_rhowV = 1
  scaling_factor_rhoEV = 1e-5

  rdg_slope_reconstruction = none
[]

[Functions]
  [motor_torque_fn]
    type = PiecewiseLinear
    x = '0 ${t1} ${t2}'
    y = '0 ${motor_torque_max} 0'
  []
  [motor_power_fn]
    type = ParsedFunction
    expression = 'torque * speed'
    symbol_names = 'torque speed'
    symbol_values = 'motor_torque shaft:omega'
  []
  [generator_torque_fn]
    type = ParsedFunction
    expression = 'slope * t'
    symbol_names = 'slope'
    symbol_values = '${generator_torque_per_shaft_speed}'
  []
  [generator_power_fn]
    type = ParsedFunction
    expression = 'torque * speed'
    symbol_names = 'torque speed'
    symbol_values = 'generator_torque shaft:omega'
  []
  [htc_wall_fn]
    type = PiecewiseLinear
    x = '0 ${t1} ${t2}'
    y = '0 0 1e3'
  []
[]

[FluidProperties]
  [fp_air]
    type = IdealGasFluidProperties
    emit_on_nan = none
  []
[]

[Closures]
  [closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [shaft]
    type = Shaft
    connected_components = 'motor compressor turbine generator'
    initial_speed = ${speed_initial}
  []
  [motor]
    type = ShaftConnectedMotor
    inertia = ${I_motor}
    torque = 0 # controlled
  []
  [generator]
    type = ShaftConnectedMotor
    inertia = ${I_generator}
    torque = generator_torque_fn
  []

  [inlet]
    type = InletStagnationPressureTemperature1Phase
    input = 'pipe1:in'
    p0 = ${p_ambient}
    T0 = ${T_ambient}
  []
  [pipe1]
    type = FlowChannel1Phase
    position = '${x1} 0 0'
    length = ${L1}
    n_elems = ${n_elems1}
    A = ${A1}
  []
  [compressor]
    type = ShaftConnectedCompressor1Phase
    position = '${x2} 0 0'
    inlet = 'pipe1:out'
    outlet = 'pipe2:in'
    A_ref = ${A_ref_comp}
    volume = ${V_comp}

    omega_rated = ${speed_rated}
    mdot_rated = ${rated_mfr}
    c0_rated = ${c0_rated_comp}
    rho0_rated = ${rho0_rated_comp}

    speeds = '0.5208 0.6250 0.7292 0.8333 0.9375'
    Rp_functions = 'rp_comp1 rp_comp2 rp_comp3 rp_comp4 rp_comp5'
    eff_functions = 'eff_comp1 eff_comp2 eff_comp3 eff_comp4 eff_comp5'

    min_pressure_ratio = 1.0

    speed_cr_I = 0
    inertia_const = ${I_comp}
    inertia_coeff = '${I_comp} 0 0 0'

    # assume no shaft friction
    speed_cr_fr = 0
    tau_fr_const = 0
    tau_fr_coeff = '0 0 0 0'
  []
  [pipe2]
    type = FlowChannel1Phase
    position = '${x2} 0 0'
    length = ${L2}
    n_elems = ${n_elems2}
    A = ${A2}
  []
  [junction2_3]
    type = JunctionOneToOne1Phase
    connections = 'pipe2:out pipe3:in'
  []
  [pipe3]
    type = FlowChannel1Phase
    position = '${x3} 0 0'
    length = ${L3}
    n_elems = ${n_elems3}
    A = ${A3}
  []
  [junction3_4]
    type = JunctionOneToOne1Phase
    connections = 'pipe3:out pipe4:in'
  []
  [pipe4]
    type = FlowChannel1Phase
    position = '${x4} 0 0'
    length = ${L4}
    n_elems = ${n_elems4}
    A = ${A4}
  []
  [turbine]
    type = ShaftConnectedCompressor1Phase
    position = '${x5} 0 0'
    inlet = 'pipe4:out'
    outlet = 'pipe5:in'
    A_ref = ${A_ref_turb}
    volume = ${V_turb}

    treat_as_turbine = true

    omega_rated = ${speed_rated}
    mdot_rated = ${rated_mfr}
    c0_rated = ${c0_rated_comp}
    rho0_rated = ${rho0_rated_comp}

    speeds = '0 0.5208 0.6250 0.7292 0.8333 0.9375'
    Rp_functions = 'rp_turb0 rp_turb1 rp_turb2 rp_turb3 rp_turb4 rp_turb5'
    eff_functions = 'eff_turb1 eff_turb1 eff_turb2 eff_turb3 eff_turb4 eff_turb5'

    min_pressure_ratio = 1.0

    speed_cr_I = 0
    inertia_const = ${I_turb}
    inertia_coeff = '${I_turb} 0 0 0'

    # assume no shaft friction
    speed_cr_fr = 0
    tau_fr_const = 0
    tau_fr_coeff = '0 0 0 0'
  []
  [pipe5]
    type = FlowChannel1Phase
    position = '${x5} 0 0'
    length = ${L5}
    n_elems = ${n_elems5}
    A = ${A5}
  []
  [outlet]
    type = Outlet1Phase
    input = 'pipe5:out'
    p = ${p_ambient}
  []

  [heating]
    type = HeatTransferFromSpecifiedTemperature1Phase
    flow_channel = pipe3
    T_wall = ${T_hot}
    Hw = htc_wall_fn
  []
[]

[ControlLogic]
  [motor_ctrl]
    type = TimeFunctionComponentControl
    component = motor
    parameter = torque
    function = motor_torque_fn
  []
[]

[Postprocessors]
  [heating_rate]
    type = ADHeatRateConvection1Phase
    block = 'pipe3'
    T = T
    T_wall = T_wall
    Hw = Hw
    P_hf = P_hf
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [motor_torque]
    type = RealComponentParameterValuePostprocessor
    component = motor
    parameter = torque
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [motor_power]
    type = FunctionValuePostprocessor
    function = motor_power_fn
    execute_on = 'INITIAL TIMESTEP_END'
    indirect_dependencies = 'motor_torque shaft:omega'
  []

  [generator_torque]
    type = ShaftConnectedComponentPostprocessor
    quantity = torque
    shaft_connected_component_uo = generator:shaftconnected_uo
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [generator_power]
    type = FunctionValuePostprocessor
    function = generator_power_fn
    execute_on = 'INITIAL TIMESTEP_END'
    indirect_dependencies = 'generator_torque shaft:omega'
  []

  [shaft_speed]
    type = ScalarVariable
    variable = 'shaft:omega'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [p_in_comp]
    type = PointValue
    variable = p
    point = '${x2_minus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_out_comp]
    type = PointValue
    variable = p
    point = '${x2_plus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_ratio_comp]
    type = ParsedPostprocessor
    pp_names = 'p_in_comp p_out_comp'
    expression = 'p_out_comp / p_in_comp'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [p_in_turb]
    type = PointValue
    variable = p
    point = '${x5_minus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_out_turb]
    type = PointValue
    variable = p
    point = '${x5_plus} 0 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_ratio_turb]
    type = ParsedPostprocessor
    pp_names = 'p_in_turb p_out_turb'
    expression = 'p_in_turb / p_out_turb'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [mfr_comp]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = mass
    junction = compressor
  []
  [mfr_turb]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe4:out
    connection_index = 0
    equation = mass
    junction = turbine
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  end_time = ${t3}
  dt = 0.1
  abort_on_solve_fail = true

  solve_type = NEWTON
  nl_rel_tol = 1e-50
  nl_abs_tol = 1e-11
  nl_max_its = 15

  l_tol = 1e-4
  l_max_its = 10
[]

[Outputs]
  [csv]
    type = CSV
    file_base = 'open_brayton_cycle'
    execute_vector_postprocessors_on = 'INITIAL'
  []
  [console]
    type = Console
    show = 'shaft_speed p_ratio_comp p_ratio_turb compressor:pressure_ratio turbine:pressure_ratio'
  []
[]

[Functions]
  # compressor pressure ratio
  [rp_comp1]
    type = PiecewiseLinear
    data_file = 'rp_comp1.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp2]
    type = PiecewiseLinear
    data_file = 'rp_comp2.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp3]
    type = PiecewiseLinear
    data_file = 'rp_comp3.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp4]
    type = PiecewiseLinear
    data_file = 'rp_comp4.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp5]
    type = PiecewiseLinear
    data_file = 'rp_comp5.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []

  # compressor efficiency
  [eff_comp1]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp2]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp3]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp4]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp5]
    type = ConstantFunction
    value = ${eff_comp}
  []

  # turbine pressure ratio
  [rp_turb0]
    type = ConstantFunction
    value = 1
  []
  [rp_turb1]
    type = PiecewiseLinear
    data_file = 'rp_turb1.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb2]
    type = PiecewiseLinear
    data_file = 'rp_turb2.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb3]
    type = PiecewiseLinear
    data_file = 'rp_turb3.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb4]
    type = PiecewiseLinear
    data_file = 'rp_turb4.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb5]
    type = PiecewiseLinear
    data_file = 'rp_turb5.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []

  # turbine efficiency
  [eff_turb1]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb2]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb3]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb4]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb5]
    type = ConstantFunction
    value = ${eff_turb}
  []
[]

(modules/thermal_hydraulics/test/tests/problems/brayton_cycle/recuperated_brayton_cycle.i)

# This input file models an open, recuperated Brayton cycle with a PID
# controlled start up using a coupled motor.
#
# Heat is supplied to the system by a volumetric heat source, and a second heat
# source is used to model a recuperator. The recuperator transfers heat from the
# turbine exhaust gas to the compressor outlet gas.
#
# Initially the fluid and heat structures are at rest at ambient conditions,
# and the shaft speed is zero.
# The transient is controlled as follows:
#   * 0   - 2000 s: Motor increases shaft speed to approx. 85,000 RPM by PID control
#   * 1000 - 8600 s: Power in main heat source increases from 0 - 104 kW
#   * 2000 - 200000 s: Torque supplied by turbine increases to steady state level
#                      as working fluid temperature increases. Torque supplied by
#                      the motor is ramped down to 0 N-m transitioning shaft control
#                      to the turbine at its rated speed of 96,000 RPM.


I_motor = 1.0

I_generator = 1.0
generator_torque_per_shaft_speed = -0.00025

motor_ramp_up_duration = 3605
motor_ramp_down_duration = 1800
post_motor_time = 2160000
t1 = ${motor_ramp_up_duration}
t2 = ${fparse t1 + motor_ramp_down_duration}
t3 = ${fparse t2 + post_motor_time}

D1 = 0.15
D2 = ${D1}
D3 = ${D1}
D4 = ${D1}
D5 = ${D1}
D6 = ${D1}
D7 = ${D1}
D8 = ${D1}

A1 = ${fparse 0.25 * pi * D1^2}
A2 = ${fparse 0.25 * pi * D2^2}
A3 = ${fparse 0.25 * pi * D3^2}
A4 = ${fparse 0.25 * pi * D4^2}
A5 = ${fparse 0.25 * pi * D5^2}
A6 = ${fparse 0.25 * pi * D6^2}
A7 = ${fparse 0.25 * pi * D7^2}
A8 = ${fparse 0.25 * pi * D8^2}

recuperator_width = 0.15

L1 = 5.0
L2 = ${L1}
L3 = ${fparse 2 * L1}
L4 = ${fparse 2 * L1}
L5 = ${L1}
L6 = ${L1}
L7 = ${fparse L1 + recuperator_width}
L8 = ${L1}

x1 = 0.0
x2 = ${fparse x1 + L1}
x3 = ${fparse x2 + L2}
x4 = ${x3}
x5 = ${fparse x4 - L4}
x6 = ${x5}
x7 = ${fparse x6 + L6}
x8 = ${fparse x7 + L7}

y1 = 0
y2 = ${y1}
y3 = ${y2}
y4 = ${fparse y3 - L3}
y5 = ${y4}
y6 = ${fparse y5 + L5}
y7 = ${y6}
y8 = ${y7}

x1_out = ${fparse x1 + L1 - 0.001}
x2_in = ${fparse x2 + 0.001}
y5_in = ${fparse y5 + 0.001}
x6_out = ${fparse x6 + L6 - 0.001}
x7_in = ${fparse x7 + 0.001}
y8_in = ${fparse y8 + 0.001}
y8_out = ${fparse y8 + L8 - 0.001}
hot_leg_in = ${y8_in}
hot_leg_out = ${y8_out}
cold_leg_in = ${fparse y3 - 0.001}
cold_leg_out = ${fparse y3 - (L3/2) - 0.001}

n_elems1 = 5
n_elems2 = ${n_elems1}
n_elems3 = ${fparse 2 * n_elems1}
n_elems4 = ${fparse 2 * n_elems1}
n_elems5 = ${n_elems1}
n_elems6 = ${n_elems1}
n_elems7 = ${n_elems1}
n_elems8 = ${n_elems1}

A_ref_comp = ${fparse 0.5 * (A1 + A2)}
V_comp = ${fparse A_ref_comp * 1.0}
I_comp = 1.0

A_ref_turb = ${fparse 0.5 * (A4 + A5)}
V_turb = ${fparse A_ref_turb * 1.0}
I_turb = 1.0

c0_rated_comp = 351.6925137
rho0_rated_comp = 1.146881112

rated_mfr = 0.25

speed_rated_rpm = 96000
speed_rated = ${fparse speed_rated_rpm * 2 * pi / 60.0}

speed_initial = 0

eff_comp = 0.79
eff_turb = 0.843

T_ambient = 300
p_ambient = 1e5

hs_power = 105750
[GlobalParams]
  gravity_vector = '0 0 0'

  initial_p = ${p_ambient}
  initial_T = ${T_ambient}
  initial_vel = 0
  initial_vel_x = 0
  initial_vel_y = 0
  initial_vel_z = 0

  fp = fp_air
  closures = closures
  f = 0

  scaling_factor_1phase = '1 1 1e-5'
  scaling_factor_rhoV = 1
  scaling_factor_rhouV = 1e-2
  scaling_factor_rhovV = 1e-2
  scaling_factor_rhowV = 1e-2
  scaling_factor_rhoEV = 1e-5
  scaling_factor_temperature = 1e-2
  rdg_slope_reconstruction = none
[]

[FluidProperties]
  [fp_air]
    type = IdealGasFluidProperties
    emit_on_nan = none
  []
[]

[SolidProperties]
  [steel]
    type = ThermalFunctionSolidProperties
    rho = 8050
    k = 45
    cp = 466
  []
[]

[Closures]
  [closures]
    type = Closures1PhaseSimple
  []
[]

[Functions]

  ##########################
  # Motor
  ##########################


  # Functions for control logic that determines when to shut off the PID system
  [is_tripped_fn]
    type = ParsedFunction
    symbol_names = 'motor_torque turbine_torque'
    symbol_values = 'motor_torque turbine_torque'
    expression = 'turbine_torque > motor_torque'
  []
  [PID_tripped_constant_value]
    type = ConstantFunction
    value = 1
  []
  [PID_tripped_status_fn]
    type = ParsedFunction
    symbol_values = 'PID_trip_status'
    symbol_names = 'PID_trip_status'
    expression = 'PID_trip_status'
  []
  [time_fn]
    type = ParsedFunction
    expression = t
  []

  # Shutdown function which ramps down the motor once told by the control logic
  [motor_torque_fn_shutdown]
    type = ParsedFunction
    symbol_values = 'PID_trip_status time_trip'
    symbol_names = 'PID_trip_status time_trip'
    expression = 'if(PID_trip_status = 1, max(2.4 - (2.4 * ((t - time_trip) / 35000)),0.0), 1)'
  []

  # Generates motor power curve
  [motor_power_fn]
    type = ParsedFunction
    expression = 'torque * speed'
    symbol_names = 'torque speed'
    symbol_values = 'motor_torque shaft:omega'
  []

  ##########################
  # Generator
  ##########################

  # Generates generator torque curve
  [generator_torque_fn]
    type = ParsedFunction
    expression = 'slope * t'
    symbol_names = 'slope'
    symbol_values = '${generator_torque_per_shaft_speed}'
  []
  # Generates generator power curve
  [generator_power_fn]
    type = ParsedFunction
    expression = 'torque * speed'
    symbol_names = 'torque speed'
    symbol_values = 'generator_torque shaft:omega'
  []

  ##########################
  # Reactor
  ##########################

  # Ramps up reactor power when activated by control logic
  [power_fn]
    type = PiecewiseLinear
    x = '0 1000 8600'
    y = '0 0 ${hs_power}'
  []

  ##########################
  # Compressor
  ##########################

  # compressor pressure ratios
  [rp_comp1]
    type = PiecewiseLinear
    data_file = 'rp_comp1.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp2]
    type = PiecewiseLinear
    data_file = 'rp_comp2.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp3]
    type = PiecewiseLinear
    data_file = 'rp_comp3.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp4]
    type = PiecewiseLinear
    data_file = 'rp_comp4.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_comp5]
    type = PiecewiseLinear
    data_file = 'rp_comp5.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []

  # compressor efficiencies
  [eff_comp1]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp2]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp3]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp4]
    type = ConstantFunction
    value = ${eff_comp}
  []
  [eff_comp5]
    type = ConstantFunction
    value = ${eff_comp}
  []

  ##########################
  # Turbine
  ##########################

  # turbine pressure ratios
  [rp_turb0]
    type = ConstantFunction
    value = 1
  []
  [rp_turb1]
    type = PiecewiseLinear
    data_file = 'rp_turb1.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb2]
    type = PiecewiseLinear
    data_file = 'rp_turb2.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb3]
    type = PiecewiseLinear
    data_file = 'rp_turb3.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb4]
    type = PiecewiseLinear
    data_file = 'rp_turb4.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []
  [rp_turb5]
    type = PiecewiseLinear
    data_file = 'rp_turb5.csv'
    x_index_in_file = 0
    y_index_in_file = 1
    format = columns
    extrap = true
  []

  # turbine efficiency
  [eff_turb1]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb2]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb3]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb4]
    type = ConstantFunction
    value = ${eff_turb}
  []
  [eff_turb5]
    type = ConstantFunction
    value = ${eff_turb}
  []
[]

[Components]

  # system inlet pulling air from the open atmosphere
  [inlet]
    type = InletStagnationPressureTemperature1Phase
    input = 'pipe1:in'
    p0 = ${p_ambient}
    T0 = ${T_ambient}
  []

  # Inlet pipe
  [pipe1]
    type = FlowChannel1Phase
    position = '${x1} ${y1} 0'
    orientation = '1 0 0'
    length = ${L1}
    n_elems = ${n_elems1}
    A = ${A1}
  []

  # Compressor as defined in MAGNET PCU document (Guillen 2020)
  [compressor]
    type = ShaftConnectedCompressor1Phase
    position = '${x2} ${y2} 0'
    inlet = 'pipe1:out'
    outlet = 'pipe2:in'
    A_ref = ${A_ref_comp}
    volume = ${V_comp}

    omega_rated = ${speed_rated}
    mdot_rated = ${rated_mfr}
    c0_rated = ${c0_rated_comp}
    rho0_rated = ${rho0_rated_comp}

    # Determines which compression ratio curve and efficiency curve to use depending on ratio of speed/rated_speed
    speeds = '0.5208 0.6250 0.7292 0.8333 0.9375'
    Rp_functions = 'rp_comp1 rp_comp2 rp_comp3 rp_comp4 rp_comp5'
    eff_functions = 'eff_comp1 eff_comp2 eff_comp3 eff_comp4 eff_comp5'

    min_pressure_ratio = 1.0

    speed_cr_I = 0
    inertia_const = ${I_comp}
    inertia_coeff = '${I_comp} 0 0 0'

    # assume no shaft friction
    speed_cr_fr = 0
    tau_fr_const = 0
    tau_fr_coeff = '0 0 0 0'
  []

  # Outlet pipe from the compressor
  [pipe2]
    type = FlowChannel1Phase
    position = '${x2} ${y2} 0'
    orientation = '1 0 0'
    length = ${L2}
    n_elems = ${n_elems2}
    A = ${A2}
  []

  # 90 degree connection between pipe 2 and 3
  [junction2_cold_leg]
    type = VolumeJunction1Phase
    connections = 'pipe2:out cold_leg:in'
    position = '${x3} ${y3} 0'
    volume = ${fparse A2*0.1}
  []

  # Cold leg of the recuperator
  [cold_leg]
    type = FlowChannel1Phase
    position = '${x3} ${y3} 0'
    orientation = '0 -1 0'
    length = ${fparse L3/2}
    n_elems = ${fparse n_elems3/2}
    A = ${A3}
  []

  # Recuperator which transfers heat from exhaust gas to reactor inlet gas to improve thermal efficency
  [recuperator]
    type = HeatStructureCylindrical
    orientation = '0 -1 0'
    position = '${x3} ${y3} 0'
    length = ${fparse L3/2}
    widths = ${recuperator_width}
    n_elems = ${fparse n_elems3/2}
    n_part_elems = 2
    names = recuperator
    solid_properties = steel
    solid_properties_T_ref = '300'
    inner_radius = ${D1}
  []
  # heat transfer from recuperator to cold leg
  [heat_transfer_cold_leg]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = cold_leg
    hs = recuperator
    hs_side = OUTER
    Hw = 10000
  []
  # heat transfer from hot leg to recuperator
  [heat_transfer_hot_leg]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = hot_leg
    hs = recuperator
    hs_side = INNER
    Hw = 10000
  []

  [junction_cold_leg_3]
    type = JunctionOneToOne1Phase
    connections = 'cold_leg:out pipe3:in'
  []
  [pipe3]
    type = FlowChannel1Phase
    position = '${x3} ${fparse y3 - (L3/2)} 0'
    orientation = '0 -1 0'
    length = ${fparse L3/2}
    n_elems = ${fparse n_elems3/2}
    A = ${A3}
  []
  # 90 degree connection between pipe 3 and 4
  [junction3_4]
    type = VolumeJunction1Phase
    connections = 'pipe3:out pipe4:in'
    position = '${x4} ${y4} 0'
    volume = ${fparse A3*0.1}
  []

  # Pipe through the "reactor core"
  [pipe4]
    type = FlowChannel1Phase
    position = '${x4} ${y4} 0'
    orientation = '-1 0 0'
    length = ${L4}
    n_elems = ${n_elems4}
    A = ${A4}
  []

  # "Reactor Core" and it's associated heat transfer to pipe 4
  [reactor]
    type = HeatStructureCylindrical
    orientation = '-1 0 0'
    position = '${x4} ${y4} 0'
    length = ${L4}
    widths = 0.15
    n_elems = ${n_elems4}
    n_part_elems = 2
    names = core
    solid_properties = steel
    solid_properties_T_ref = '300'
  []
  [total_power]
    type = TotalPower
    power = 0
  []
  [heat_generation]
    type = HeatSourceFromTotalPower
    power = total_power
    hs = reactor
    regions = core
  []
  [heat_transfer]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = pipe4
    hs = reactor
    hs_side = OUTER
    Hw = 10000
  []

  # 90 degree connection between pipe 4 and 5
  [junction4_5]
    type = VolumeJunction1Phase
    connections = 'pipe4:out pipe5:in'
    position = '${x5} ${y5} 0'
    volume = ${fparse A4*0.1}
  []

  # Pipe carrying hot gas back to the PCU
  [pipe5]
    type = FlowChannel1Phase
    position = '${x5} ${y5} 0'
    orientation = '0 1 0'
    length = ${L5}
    n_elems = ${n_elems5}
    A = ${A5}
  []

  # 90 degree connection between pipe 5 and 6
  [junction5_6]
    type = VolumeJunction1Phase
    connections = 'pipe5:out pipe6:in'
    position = '${x6} ${y6} 0'
    volume = ${fparse A5*0.1}
  []

  # Inlet pipe to the turbine
  [pipe6]
    type = FlowChannel1Phase
    position = '${x6} ${y6} 0'
    orientation = '1 0 0'
    length = ${L6}
    n_elems = ${n_elems6}
    A = ${A6}
  []

  # Turbine as defined in MAGNET PCU document (Guillen 2020) and (Wright 2006)
  [turbine]
    type = ShaftConnectedCompressor1Phase
    position = '${x7} ${y7} 0'
    inlet = 'pipe6:out'
    outlet = 'pipe7:in'
    A_ref = ${A_ref_turb}
    volume = ${V_turb}

    # A turbine is treated as an "inverse" compressor, this value determines if component is to be treated as turbine or compressor
    # If treat_as_turbine is omitted, code automatically assumes it is a compressor
    treat_as_turbine = true

    omega_rated = ${speed_rated}
    mdot_rated = ${rated_mfr}
    c0_rated = ${c0_rated_comp}
    rho0_rated = ${rho0_rated_comp}

    # Determines which compression ratio curve and efficiency curve to use depending on ratio of speed/rated_speed
    speeds = '0 0.5208 0.6250 0.7292 0.8333 0.9375'
    Rp_functions = 'rp_turb0 rp_turb1 rp_turb2 rp_turb3 rp_turb4 rp_turb5'
    eff_functions = 'eff_turb1 eff_turb1 eff_turb2 eff_turb3 eff_turb4 eff_turb5'

    min_pressure_ratio = 1.0

    speed_cr_I = 0
    inertia_const = ${I_turb}
    inertia_coeff = '${I_turb} 0 0 0'

    # assume no shaft friction
    speed_cr_fr = 0
    tau_fr_const = 0
    tau_fr_coeff = '0 0 0 0'
  []

  # Outlet pipe from turbine
  [pipe7]
    type = FlowChannel1Phase
    position = '${x7} ${y7} 0'
    orientation = '1 0 0'
    length = ${L7}
    n_elems = ${n_elems7}
    A = ${A7}
  []

  # 90 degree connection between pipe 7 and 8
  [junction7_hot_leg]
    type = VolumeJunction1Phase
    connections = 'pipe7:out hot_leg:in'
    position = '${x8} ${y8} 0'
    volume = ${fparse A7*0.1}
  []

  # Hot leg of the recuperator
  [hot_leg]
    type = FlowChannel1Phase
    position = '${x8} ${y8} 0'
    orientation = '0 1 0'
    length = ${L8}
    n_elems = ${n_elems8}
    A = ${A8}
  []

  # System outlet dumping exhaust gas to the atmosphere
  [outlet]
    type = Outlet1Phase
    input = 'hot_leg:out'
    p = ${p_ambient}
  []

  # Roatating shaft connecting motor, compressor, turbine, and generator
  [shaft]
    type = Shaft
    connected_components = 'motor compressor turbine generator'
    initial_speed = ${speed_initial}
  []

  # 3-Phase electircal motor used for system start-up, controlled by PID
  [motor]
    type = ShaftConnectedMotor
    inertia = ${I_motor}
    torque = 0 # controlled
  []

  # Electric generator supplying power to the grid
  [generator]
    type = ShaftConnectedMotor
    inertia = ${I_generator}
    torque = generator_torque_fn
  []
[]


# Control logics which govern startup of the motor, startup of the "reactor core", and shutdown of the motor
[ControlLogic]

  # Sets desired shaft speed to be reached by motor NOTE: SHOULD BE SET LOWER THAN RATED TURBINE RPM
  [set_point]
    type = GetFunctionValueControl
    function = ${fparse speed_rated_rpm - 9000}
  []

  # PID with gains determined by iterative process NOTE: Gain values are system specific
  [initial_motor_PID]
    type = PIDControl
    set_point = set_point:value
    input = shaft_RPM
    initial_value = 0
    K_p = 0.0011
    K_i = 0.00000004
    K_d = 0
  []

  # Determines when the PID system should be running and when it should begin the shutdown cycle. If needed: PID output, else: shutdown function
  [logic]
    type = ParsedFunctionControl
    function = 'if(motor+0.5 > turb, PID, shutdown_fn)'
    symbol_names = 'motor turb PID shutdown_fn'
    symbol_values = 'motor_torque turbine_torque initial_motor_PID:output motor_torque_fn_shutdown'
  []

  # Takes the output generated in [logic] and applies it to the motor torque
  [motor_PID]
    type = SetComponentRealValueControl
    component = motor
    parameter = torque
    value = logic:value
  []
  # Determines when to turn on heat source
  [power_logic]
    type = ParsedFunctionControl
    function = 'power_fn'
    symbol_names = 'power_fn'
    symbol_values = 'power_fn'
  []
  # Applies heat source to the total_power block
  [power_applied]
    type = SetComponentRealValueControl
    component = total_power
    parameter = power
    value = power_logic:value
  []
[]

[Controls]

  # Enables set_PID_tripped
  [PID_trip_status]
    type = ConditionalFunctionEnableControl
    conditional_function = is_tripped_fn
    enable_objects = 'AuxScalarKernels::PID_trip_status_aux'
    execute_on = 'TIMESTEP_END'
  []

  # Enables set_time_PID
  [time_PID]
    type = ConditionalFunctionEnableControl
    conditional_function = PID_tripped_status_fn
    disable_objects = 'AuxScalarKernels::time_trip_aux'
    execute_on = 'TIMESTEP_END'
  []
[]

[AuxVariables]

  # Creates a variable that will later be set to the time when tau_turbine > tau_motor
  [time_trip]
    order = FIRST
    family = SCALAR
  []

  # Creates variable which indicates if tau_turbine > tau_motor....... If tau_motor > tau_turbine, 0, else 1
  [PID_trip_status]
    order = FIRST
    family = SCALAR
    initial_condition = 0
  []
[]

[AuxScalarKernels]

  # Creates variable from time_fn which indicates when tau_turbine > tau_motor
  [time_trip_aux]
    type = FunctionScalarAux
    function = time_fn
    variable = time_trip
    execute_on = 'TIMESTEP_END'
  []

  # Overwrites variable PID_trip_status to the value from PID_tripped_constant_value (changes 0 to 1)
  [PID_trip_status_aux]
    type = FunctionScalarAux
    function = PID_tripped_constant_value
    variable = PID_trip_status
    execute_on = 'TIMESTEP_END'
    enable = false
  []
[]


[Postprocessors]

  # Indicates when tau_turbine > tau_motor
  [trip_time]
    type = ScalarVariable
    variable = time_trip
    execute_on = 'TIMESTEP_END'
  []

  ##########################
  # Motor
  ##########################

  [motor_torque]
    type = RealComponentParameterValuePostprocessor
    component = motor
    parameter = torque
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [motor_power]
    type = FunctionValuePostprocessor
    function = motor_power_fn
    execute_on = 'INITIAL TIMESTEP_END'
  []

  ##########################
  # generator
  ##########################

  [generator_torque]
    type = ShaftConnectedComponentPostprocessor
    quantity = torque
    shaft_connected_component_uo = generator:shaftconnected_uo
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [generator_power]
    type = FunctionValuePostprocessor
    function = generator_power_fn
    execute_on = 'INITIAL TIMESTEP_END'
  []

  ##########################
  # Shaft
  ##########################

  # Speed in rad/s
  [shaft_speed]
    type = ScalarVariable
    variable = 'shaft:omega'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  # speed in RPM
  [shaft_RPM]
    type = ParsedPostprocessor
    pp_names = 'shaft_speed'
    expression = '(shaft_speed * 60) /( 2 * ${fparse pi})'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  ##########################
  # Compressor
  ##########################
  [comp_dissipation_torque]
    type = ElementAverageValue
    variable = dissipation_torque
    block = 'compressor'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [comp_isentropic_torque]
    type = ElementAverageValue
    variable = isentropic_torque
    block = 'compressor'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [comp_friction_torque]
    type = ElementAverageValue
    variable = friction_torque
    block = 'compressor'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [compressor_torque]
    type = ParsedPostprocessor
    pp_names = 'comp_dissipation_torque comp_isentropic_torque comp_friction_torque'
    expression = 'comp_dissipation_torque + comp_isentropic_torque + comp_friction_torque'
  []
  [p_in_comp]
    type = PointValue
    variable = p
    point = '${x1_out} ${y1} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_out_comp]
    type = PointValue
    variable = p
    point = '${x2_in} ${y2} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_ratio_comp]
    type = ParsedPostprocessor
    pp_names = 'p_in_comp p_out_comp'
    expression = 'p_out_comp / p_in_comp'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_in_comp]
    type = PointValue
    variable = T
    point = '${x1_out} ${y1} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_out_comp]
    type = PointValue
    variable = T
    point = '${x2_in} ${y2} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_ratio_comp]
    type = ParsedPostprocessor
    pp_names = 'T_in_comp T_out_comp'
    expression = '(T_out_comp - T_in_comp) / T_out_comp'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [mfr_comp]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = mass
    junction = compressor
  []

  ##########################
  # turbine
  ##########################

  [turb_dissipation_torque]
    type = ElementAverageValue
    variable = dissipation_torque
    block = 'turbine'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turb_isentropic_torque]
    type = ElementAverageValue
    variable = isentropic_torque
    block = 'turbine'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turb_friction_torque]
    type = ElementAverageValue
    variable = friction_torque
    block = 'turbine'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turbine_torque]
    type = ParsedPostprocessor
    pp_names = 'turb_dissipation_torque turb_isentropic_torque turb_friction_torque'
    expression = 'turb_dissipation_torque + turb_isentropic_torque + turb_friction_torque'
  []
  [p_in_turb]
    type = PointValue
    variable = p
    point = '${x6_out} ${y6} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_out_turb]
    type = PointValue
    variable = p
    point = '${x7_in} ${y7} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [p_ratio_turb]
    type = ParsedPostprocessor
    pp_names = 'p_in_turb p_out_turb'
    expression = 'p_in_turb / p_out_turb'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_in_turb]
    type = PointValue
    variable = T
    point = '${x6_out} ${y6} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_out_turb]
    type = PointValue
    variable = T
    point = '${x7_in} ${y7} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [mfr_turb]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = mass
    junction = turbine
  []

  ##########################
  # Recuperator
  ##########################

  [cold_leg_in]
    type = PointValue
    variable = T
    point = '${x3} ${cold_leg_in} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [cold_leg_out]
    type = PointValue
    variable = T
    point = '${x3} ${cold_leg_out} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [hot_leg_in]
    type = PointValue
    variable = T
    point = '${x8} ${hot_leg_in} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [hot_leg_out]
    type = PointValue
    variable = T
    point = '${x8} ${hot_leg_out} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  ##########################
  # Reactor
  ##########################

  [reactor_inlet]
    type = PointValue
    variable = T
    point = '${x4} ${y4} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [reactor_outlet]
    type = PointValue
    variable = T
    point = '${x5} ${y5_in} 0'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  end_time = ${t3}
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.01
    growth_factor = 1.1
    cutback_factor = 0.9
  []
  dtmin = 1e-5
  dtmax = 1000
  steady_state_detection = true
  steady_state_start_time = 200000

  solve_type = NEWTON
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8
  nl_max_its = 15

  l_tol = 1e-4
  l_max_its = 10

  petsc_options_iname  = '-pc_type'
  petsc_options_value  = ' lu     '
[]

[Outputs]
  [e]
    type = Exodus
    file_base = 'recuperated_brayton_cycle_out'
  []
  [csv]
    type = CSV
    file_base = 'recuperated_brayton_cycle'
    execute_vector_postprocessors_on = 'INITIAL'
  []
  [console]
    type = Console
    show = 'shaft_speed p_ratio_comp p_ratio_turb pressure_ratio pressure_ratio'
  []
[]

Introduction
Transient Description
Input File description
Control Logic
Motor Shutdown
Results