How to control both the generator torque and motor torque varying with both speed and time? #31855

Wang-Yihu · 2025-11-06T07:36:44Z

Wang-Yihu
Nov 6, 2025

Let me give you an recuperated brayton cycle by me. It is imitated from https://mooseframework.inl.gov/modules/thermal_hydraulics/examples/recuperated_brayton_cycle/recuperated_brayton_cycle.html. The generator torque equals to -0.00115*omega. The motor torque varying with time is from the file
tau_m_data_with_time.csv
The input file is

# 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   - 56789 s: Motor torque decreases from 95.7000348 N-m to 0 according to tau_m_data_with_time.csv
#   * 1000 - 8600 s: Power in main heat source increases from 0 - 1000 kW
#   * 56789 s - ?  : The system will be steady with 0 motor torque.

####################################################################################
# Fuel parameters for snerdi
# Diameter:    16.38 mm    , or 17.4 mm
# Height  :       1.2  m
# Fuel number: 30*36 = 1080, or 30*42 = 1260

# Now I choose diameter = 16.38 mm 
# Radial size: fuel,       gas gap, clad,       monolith
#              16.38e-3/2,          16.72e-3/2, 19.0e-3/2

# Volume: 
# Total volume:                pi/4*(16.38e-3)**2*1.2*1080 =  0.273100615 m3
# Total heat transfer surface: pi*(16.38e-3)*1.2*1260 = 82.651532805 m3
####################################################################################

####################################################################################
# inlet - pipe1 - compressor - pipe2    outlet
#
#                               |         |
#                                       
#                             cold_leg  hot_leg  recuperator
#                               |         |
#   --pipe6--turbine--pipe7---  |   -------
#   |                           |
#                               
# pipe5                        pipe3
#
#   |       reactor             |
# pipe4_3 - pipe4_2- pipe4_1---- 
####################################################################################

I_motor = 1.0

I_generator = 1.0

generator_torque_per_shaft_speed = -0.00115

t_loss_generator = 180000

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.60
D1 = 0.65
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.002
d_i = ${fparse 0.015}

L_reactor = 1.2

num_rods_reactor  = 1080
num_rods_recuperator = 570

k_uo2   = 3.4
k_he = 0.3

fuel_radius =           ${fparse 16.38e-3/2}
clad_radius =           ${fparse 16.72e-3/2}
monolith_radius =       ${fparse 19.0e-3/2}
monolith_radius_outer = ${fparse 20.8e-3/2}
th_1 = ${fparse clad_radius-fuel_radius}
th_2 = ${fparse monolith_radius-clad_radius}
th_3 = ${fparse monolith_radius_outer-monolith_radius}

L1 = 5.0
L2 = ${L1}
L3 = ${fparse 2 * L1}
L4 = ${fparse 2 * L1}
L5 = ${L1}
L6 = ${L1}
#L7 = ${fparse L1 + recuperator_width}
L7 = ${fparse L1 + 0.1*L1}
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}

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 = 3.6

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 = 1000000
[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
  []
[]

[HeatStructureMaterials]
  [steel]
    type = SolidMaterialProperties
    rho = 8050
    k = 45
    cp = 466
  []
  [uo2]
    type = SolidMaterialProperties
    rho = 10980
    k = ${k_uo2}
    cp = 330
  []
  [he]
    type = SolidMaterialProperties
    rho = 0.0001216
    k = ${k_he}
    cp = 5192
  []
[]

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

[Functions]

  [if_generator_on]
    type = ParsedFunction
    symbol_names = 'omega'
    symbol_values = 'shaft:omega'
    expression = 'if(t>${t_loss_generator}, 0, omega*${generator_torque_per_shaft_speed})'
  []

  ##########################
  # Motor
  ##########################
  
  [tau_m_fn]
    type = PiecewiseLinear
    data_file = 'tau_m_data_with_time.csv'
    format = columns
    extrap = true
    scale_factor = 1.0
  []

  # 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, not active for this example.
  [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
    x = '0.2868480 0.3084640 0.3300760 0.3516920 0.3733040 0.3949200 0.4165200 0.4381600 0.4597600 0.4813600 0.5030000 0.5187200'
    y = '1.5153000 1.5088000 1.5032000 1.4984000 1.4914000 1.4863000 1.4846000 1.4756000 1.4540000 1.4237000 1.3856000 1.3506000'
    extrap = true
  []
  [rp_comp2]
    type = PiecewiseLinear
    x = '0.3674120 0.3890240 0.4106400 0.4322400 0.4523200 0.4754800 0.4970800 0.5177200 0.5442400 0.5658800 0.5874800 0.6090800 0.6307200 0.6484000 0.6621600 0.6759200 0.6886800 0.6994800 0.7093200 0.7181600 0.7253600 0.7338800 0.7407600 0.7466400 0.7517600 0.7595600 0.7633600'
    y = '1.7741000 1.7626000 1.7550000 1.7536000 1.7497000 1.7455000 1.7376000 1.7276000 1.7026000 1.6780000 1.6405000 1.5914000 1.5233000 1.4564000 1.3890000 1.3187000 1.2318000 1.1631000 1.0851000 1.0143000 0.9425300 0.8601900 0.7844300 0.7160400 0.6450700 0.5733400 0.5341100'    
    extrap = true
  []
  [rp_comp3]
    type = PiecewiseLinear
    x = '0.4067200 0.4283200 0.4499200 0.4715600 0.4931600 0.5147600 0.5364000 0.5580000 0.5796400 0.6012400 0.6228400 0.6444800 0.6660800 0.6872000 0.7093200 0.7299200 0.7476400 0.7643200 0.7800400 0.7928000 0.8036400 0.8134400 0.8222800 0.8301600 0.8380000 0.8448800 0.8507600 0.8566800 0.8635600 0.8704400 0.8782800'
    y = '2.1535000 2.1352000 2.1251000 2.1192000 2.1135000 2.1088000 2.1071000 2.1068000 2.1043000 2.0975000 2.0877000 2.0742000 2.0539000 2.0272000 1.9813000 1.9245000 1.8643000 1.7959000 1.7227000 1.6551000 1.5896000 1.5110000 1.4403000 1.3714000 1.2969000 1.2326000 1.1655000 1.0985000 1.0177000 0.9385600 0.8578700' 
    extrap = true
  []
  [rp_comp4]
    type = PiecewiseLinear
    x = '0.4872800 0.5088800 0.5304800 0.5521200 0.5737200 0.5953600 0.6169600 0.6385600 0.6602000 0.6818000 0.7034000 0.7250400 0.7466400 0.7682400 0.7898800 0.8114800 0.8330800 0.8547200 0.8753600 0.8940000 0.9097200 0.9224800 0.9342800 0.9450800 0.9568800'
    y = '2.6634000 2.6485000 2.6366000 2.6248000 2.6184000 2.6124000 2.6110000 2.6102000 2.6110000 2.6102000 2.6110000 2.6099000 2.6060000 2.5989000 2.5874000 2.5642000 2.5309000 2.4824000 2.4202000 2.3500000 2.2832000 2.2195000 2.1437000 2.0706000 1.9947000' 
    extrap = true
  []
  [rp_comp5]
    type = PiecewiseLinear
    x = '0.5658800 0.5874800 0.6090800 0.6307200 0.6523200 0.6739200 0.6955600 0.7171600 0.7387600 0.7604000 0.7820000 0.8036400 0.8252400 0.8468400 0.8684800 0.8900800 0.9116800 0.9333200 0.9549200 0.9765200 0.9981600 1.0197600 1.0413600 1.0600400 1.0757600 1.0904800 1.1032800 1.1150800 1.1248800 1.1327600 1.1406000 1.1484800 1.1563200 1.1632000 1.1690800 1.1750000 1.1808800 1.1867600 1.1916800 1.1976000'
    y = '3.4718000 3.4428000 3.4189000 3.3970000 3.3728000 3.3455000 3.3249000 3.3106000 3.2976000 3.2886000 3.2819000 3.2799000 3.2844000 3.2892000 3.2898000 3.2892000 3.2889000 3.2830000 3.2703000 3.2551000 3.2281000 3.1873000 3.1291000 3.0660000 2.9999000 2.9281000 2.8556000 2.7823000 2.7101000 2.6451000 2.5778000 2.5099000 2.4280000 2.3562000 2.2850000 2.2128000 2.1344000 2.0580000 1.9920000 1.9198000' 
    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
    x = '0.2082520 0.2298640 0.2514800 0.2730960 0.2947080 0.3163240 0.3385240 0.3595520 0.3811640 0.4018000 0.4188400 0.4327600 0.4434800 0.4511200 0.4578000 0.4648800 0.4676400 0.4794000'
    y = '1.1197000 1.1292000 1.1423000 1.1599000 1.1751000 1.1934000 1.2179000 1.2523000 1.2972000 1.3575000 1.4183000 1.4967000 1.5555000 1.6311000 1.7075000 1.7896000 1.8465000 1.9957000'
    extrap = true
  []
  [rp_turb2]
    type = PiecewiseLinear
    x = '0.3674120 0.3890240 0.4106400 0.4322400 0.4538800 0.4705600 0.4922000 0.5079200 0.5196800 0.5275600 0.5354000 0.5422800 0.5491600 0.5590000'
    y = '1.2029000 1.2314000 1.2615000 1.2993000 1.3482000 1.3956000 1.4717000 1.5365000 1.6164000 1.6763000 1.7622000 1.8331000 1.9107000 1.9973000'
    extrap = true
  []
  [rp_turb3]
    type = PiecewiseLinear
    x = '0.4894400 0.5080800 0.5342000 0.5514400 0.5776400 0.5953600 0.6161200 0.6299200 0.6395600 0.6503600 0.6592000 0.6660800 0.6700000 0.6778800 0.6818000 0.6876800 0.6916400 0.6960400 0.6994800 0.7034000 0.7102800 0.7112800 0.7132400 0.7171600'
    y = '1.2856000 1.3172000 1.3532000 1.3853000 1.4625000 1.5312000 1.6149000 1.6750000 1.7413000 1.8149000 1.8875000 1.9559000 2.0013000 2.1282000 2.2221000 2.3390000 2.4564000 2.6010000 2.7379000 2.8729000 3.1802000 3.2705000 3.3605000 3.4468000'
    extrap = true
  []
  [rp_turb4]
    type = PiecewiseLinear
    x = '0.6071200 0.6287600 0.6493600 0.6759200 0.6975200 0.7191200 0.7407600 0.7604000 0.7741600 0.7888800 0.7996800 0.8085200 0.8164000 0.8232800 0.8291600 0.8350800 0.8390000 0.8422800 0.8478400 0.8511600 0.8547200 0.8586400 0.8635600 0.8664800 0.8684800 0.8694400 0.8724000 0.8733600 0.8763200 0.8792800'
    y = '1.3228000 1.3484000 1.3735000 1.4176000 1.4604000 1.5073000 1.5648000 1.6329000 1.6868000 1.7667000 1.8383000 1.9146000 1.9866000 2.0642000 2.1375000 2.2102000 2.2992000 2.3995000 2.5516000 2.6760000 2.8060000 2.9616000 3.1321000 3.2404000 3.3299000 3.4034000 3.4827000 3.5489000 3.6441000 3.7392000'
    extrap = true
  []
  [rp_turb5]
    type = PiecewiseLinear
    x = '0.6975200 0.7191200 0.7407600 0.7623600 0.7839600 0.8056000 0.8304800 0.8488000 0.8704400 0.8920400 0.9136400 0.9352800 0.9549200 0.9716400 0.9853600 0.9971600 1.0079600 1.0168000 1.0244800 1.0296000 1.0364800 1.0404000 1.0443200 1.0459600 1.0480400 1.0492400 1.0531600 1.0576400 1.0610400 1.0639600 1.0669200 1.0705600 1.0736000 1.0755600 1.0787200'
    y = '1.2932000 1.3060000 1.3271000 1.3501000 1.3780000 1.4050000 1.4415000 1.4684000 1.5015000 1.5409000 1.5900000 1.6541000 1.7158000 1.7865000 1.8472000 1.9161000 1.9997000 2.0694000 2.1477000 2.2175000 2.2945000 2.3562000 2.4547000 2.5110000 2.5962000 2.6544000 2.8277000 2.9733000 3.1269000 3.1945000 3.3560000 3.4621000 3.5716000 3.6771000 3.8214000'
    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
    materials = steel
    inner_radius = ${fparse d_i/2}
    num_rods = ${num_rods_recuperator}
  []
  # heat transfer from recuperator to cold leg
  [heat_transfer_cold_leg]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = cold_leg
    hs = recuperator
    hs_side = OUTER
    Hw = 500
  []
  # heat transfer from hot leg to recuperator
  [heat_transfer_hot_leg]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = hot_leg
    hs = recuperator
    hs_side = INNER
    Hw = 500
  []

  [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_1]
    type = JunctionOneToOne1Phase
    connections = 'pipe3:out pipe4_1:in'
  []

  [pipe4_1]
    type = FlowChannel1Phase
    position = '${x4} ${y4} 0'
    orientation = '-1 0 0'
    length = ${fparse (L4-L_reactor)/2}
    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 = '${fparse x4-(L4-L_reactor)/2} ${fparse y4+monolith_radius_outer} 0'
    length = ${L_reactor}
    widths = '${fuel_radius} ${th_1} ${th_2} ${th_3}'
    n_elems = ${n_elems4}
    n_part_elems = '20 2 10 10'
    names =     'fuel  gap  clad    monolith'
    materials = 'uo2   he   steel   steel'
    inner_radius = 0.0
    num_rods = ${fparse num_rods_reactor}
  []
  [total_power]
    type = TotalPower
    power = 0        # controlled
  []
  [heat_generation]
    type = HeatSourceFromTotalPower
    power = total_power
    hs = reactor
    regions = fuel
  []
  [heat_transfer]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = pipe4_2
    hs = reactor
    hs_side = OUTER
    Hw = 1700
  []

  [junction4_1_4_2]
    type = JunctionOneToOne1Phase
    connections = 'pipe4_1:out pipe4_2:in'
  []

  [pipe4_2]
    type = FlowChannel1Phase
    position = '${fparse x4-(L4-L_reactor)/2} ${y4} 0'
    orientation = '-1 0 0'
    length = ${L_reactor}
    n_elems = ${n_elems4}
    A = ${A4}
  []

  [junction4_2_4_3]
    type = JunctionOneToOne1Phase
    connections = 'pipe4_2:out pipe4_3:in'
  []

  [pipe4_3]
    type = FlowChannel1Phase
    position = '${fparse x4-(L4+L_reactor)/2} ${y4} 0'
    orientation = '-1 0 0'
    length = ${fparse (L4-L_reactor)/2}
    n_elems = ${n_elems4}
    A = ${A4}
  []
  
  [junction4_3_5]
    type = JunctionOneToOne1Phase
    connections = 'pipe4_3:out pipe5:in'
  []
  
  # 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
  # motor_control 7
  [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_control 1
  [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
#    torque = 0 # controlled
  []
[]


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

[ControlLogic]
   [motor_torque_logic]
    type = ParsedFunctionControl
    function = 'tau_m_fn'
    symbol_names = 'tau_m_fn'
    symbol_values = 'tau_m_fn'
  []
  [motor_ctrl]
    type = SetComponentRealValueControl
    component = motor
    parameter = torque
    value = motor_torque_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
  []
  
  [generator_torque_logic]
    type = ParsedFunctionControl
    function = 'if_generator_on'
    symbol_names = 'if_generator_on'
    symbol_values = 'if_generator_on'
  []
#  [generator_ctrl]
#    type = SetComponentRealValueControl
#    component = generator
#    parameter = torque
#    value = generator_torque_logic:value
#  []
[]

[Postprocessors]

  [a_if_generator_on]
    type = FunctionValuePostprocessor
    function = if_generator_on
    outputs = csv_2
  []
  
  ##########################
  # 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
  # motor_control 6
  [shaft_speed]
    type = ScalarVariable
    variable = 'shaft:omega'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  # speed in RPM
  # motor_control 5
  [shaft_RPM]
    type = ParsedPostprocessor
    pp_names = 'shaft_speed'
    function = '(shaft_speed * 60) /( 2 * ${fparse pi})'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  ##########################
  # Compressor
  ##########################
  [comp_dissipation_torque]
    type = ScalarVariable
    variable = 'compressor:dissipation_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [comp_isentropic_torque]
    type = ScalarVariable
    variable = 'compressor:isentropic_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [comp_friction_torque]
    type = ScalarVariable
    variable = 'compressor:friction_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [compressor_torque]
    type = ParsedPostprocessor
    pp_names = 'comp_dissipation_torque comp_isentropic_torque comp_friction_torque'
    function = 'comp_dissipation_torque + comp_isentropic_torque + comp_friction_torque'
  []

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

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

  [turb_dissipation_torque]
    type = ScalarVariable
    variable = 'turbine:dissipation_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turb_isentropic_torque]
    type = ScalarVariable
    variable = 'turbine:isentropic_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turb_friction_torque]
    type = ScalarVariable
    variable = 'turbine:friction_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turbine_torque]
    type = ParsedPostprocessor
    pp_names = 'turb_dissipation_torque turb_isentropic_torque turb_friction_torque'
    function = 'turb_dissipation_torque + turb_isentropic_torque + turb_friction_torque'
  []

  [mfr_turb]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = mass
    junction = turbine
  []

  ##########################
  #These parameters are not correct for transient
  #Only for steady!
  # C: compressor, R: recuperator, H: Heat section, T: turbine
  #       6
  #       |
  # 1-C-2-R-3-H-4-T
  #       |       |
  #       ---5----|
  ##########################
  [a_T_1]
    type = SideAverageValue
    variable = T
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_T_2]
    type = SideAverageValue
    variable = T
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_T_3]
    type = SideAverageValue
    variable = T
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_T_4]
    type = SideAverageValue
    variable = T
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_T_5]
    type = SideAverageValue
    variable = T
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_T_6]
    type = SideAverageValue
    variable = T
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_p_1]
    type = SideAverageValue
    variable = p
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_p_2]
    type = SideAverageValue
    variable = p
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_p_3]
    type = SideAverageValue
    variable = p
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_p_4]
    type = SideAverageValue
    variable = p
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_p_5]
    type = SideAverageValue
    variable = p
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_p_6]
    type = SideAverageValue
    variable = p
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_v_1]
    type = SideAverageValue
    variable = v
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_v_2]
    type = SideAverageValue
    variable = v
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_v_3]
    type = SideAverageValue
    variable = v
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_v_4]
    type = SideAverageValue
    variable = v
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_v_5]
    type = SideAverageValue
    variable = v
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_v_6]
    type = SideAverageValue
    variable = v
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_H_1]
    type = SideAverageValue
    variable = H
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_H_2]
    type = SideAverageValue
    variable = H
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_H_3]
    type = SideAverageValue
    variable = H
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_H_4]
    type = SideAverageValue
    variable = H
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_H_5]
    type = SideAverageValue
    variable = H
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_H_6]
    type = SideAverageValue
    variable = H
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_T_recuperator_w_inner]
    type = SideAverageValue
    variable = T_solid
    boundary = recuperator:inner
    outputs = csv_2
  []
  [a_T_recuperator_w_outer]
    type = SideAverageValue
    variable = T_solid
    boundary = recuperator:outer
    outputs = csv_2
  []
  [a_T_reactor_w_inner]
    type = SideAverageValue
    variable = T_solid
    boundary = reactor:inner
    outputs = csv_2
  []
  [a_T_reactor_w_outer]
    type = SideAverageValue
    variable = T_solid
    boundary = reactor:outer
    outputs = csv_2
  []  
  [a_vel_x_1]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_vel_x_2]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_vel_x_3]
    type = SideAverageValue
    variable = vel_x
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_vel_x_4]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_vel_x_5]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_vel_x_6]
    type = SideAverageValue
    variable = vel_x
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_vel_y_1]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_vel_y_2]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_vel_y_3]
    type = SideAverageValue
    variable = vel_y
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_vel_y_4]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_vel_y_5]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_vel_y_6]
    type = SideAverageValue
    variable = vel_y
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_vel_z_1]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_vel_z_2]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_vel_z_3]
    type = SideAverageValue
    variable = vel_z
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_vel_z_4]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_vel_z_5]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_vel_z_6]
    type = SideAverageValue
    variable = vel_z
    boundary = hot_leg:out
    outputs = csv_2     
  []
############################################################
# The precise way to retrieve parameters
############################################################
# rho*u*A = m_dot
# m_dot*u + p*A = n_dot
# m_dot*cp*T + 0.5*m_dot*u*u = q_not
# p/rho = R/M*T

# So, p = (n_dot-m_dot*u)/A, T = p*M/R/rho, rho = m_dot/u/A
# input them into Eq 3
# (0.5*m_dot - cp*M/R*m_dot)*u^2 + (cp*M/R_u*n_dot)*u - q_dot = 0
# (0.5*m_dot - 1.4/0.4*m_dot)*u^2 + (1.4/0.4*n_dot)*u - q_dot = 0       
  [a_m_dot_1]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = mass
    junction = compressor
    outputs = csv_3
  []
  [a_momentum_dot_1]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = momentum
    junction = compressor
    outputs = csv_3
  []
  [a_energy_dot_1]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = energy
    junction = compressor
    outputs = csv_3
  []
  [a_a_1]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_1'
    function = '0.5*a_m_dot_1 - 1.4/0.4*a_m_dot_1'
    outputs = csv_3
  []
  [a_b_1]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_1'
    function = '(1.4/0.4*a_momentum_dot_1)'
    outputs = csv_3
  []
  [a_c_1]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_1'
    function = '-a_energy_dot_1'
    outputs = csv_3
  []
  [a_delta_1]
    type = ParsedPostprocessor
    pp_names = 'a_a_1 a_b_1 a_c_1'
    function = 'a_b_1^2 - 4*a_a_1*a_c_1'
    outputs = csv_3
  []
  [a_a_u_1]
    type = ParsedPostprocessor
    pp_names = 'a_a_1 a_b_1 a_c_1 a_delta_1'
    function = '(-a_b_1 + sqrt(a_delta_1))/(2*a_a_1)'
    outputs = csv_3
  []
  [a_a_rho_1]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_1 a_a_u_1'
    function = 'a_m_dot_1/a_a_u_1/${A1}'
    outputs = csv_3
  []
  [a_a_p_1]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_1 a_m_dot_1 a_a_u_1'
    function = '(a_momentum_dot_1 - a_m_dot_1*a_a_u_1)/${A1}'
    outputs = csv_3
  []
  [a_a_T_1]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_1 a_a_rho_1'
     function = 'a_a_p_1*0.029/8.3144598/a_a_rho_1'
     outputs = csv_3
  []
  [a_a_H_1]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_1 a_m_dot_1'
    function = 'a_energy_dot_1/a_m_dot_1'
    outputs = csv_3
  []
#########################################################################
    
  [a_m_dot_2]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe2:in
    connection_index = 1
    equation = mass
    junction = compressor
    outputs = csv_3
  []
  [a_momentum_dot_2]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe2:in
    connection_index = 1
    equation = momentum
    junction = compressor
    outputs = csv_3
  []
  [a_energy_dot_2]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe2:in
    connection_index = 1
    equation = energy
    junction = compressor
    outputs = csv_3
  []
  [a_a_2]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_2'
    function = '0.5*a_m_dot_2 - 1.4/0.4*a_m_dot_2'
    outputs = csv_3
  []
  [a_b_2]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_2'
    function = '(1.4/0.4*a_momentum_dot_2)'
    outputs = csv_3
  []
  [a_c_2]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_2'
    function = '-a_energy_dot_2'
    outputs = csv_3
  []
  [a_delta_2]
    type = ParsedPostprocessor
    pp_names = 'a_a_2 a_b_2 a_c_2'
    function = 'a_b_2^2 - 4*a_a_2*a_c_2'
    outputs = csv_3
  []
  [a_a_u_2]
    type = ParsedPostprocessor
    pp_names = 'a_a_2 a_b_2 a_c_2 a_delta_2'
    function = '(-a_b_2 + sqrt(a_delta_2))/(2*a_a_2)'
    outputs = csv_3
  []
  [a_a_rho_2]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_2 a_a_u_2'
    function = 'a_m_dot_2/a_a_u_2/${A2}'
    outputs = csv_3
  []
  [a_a_p_2]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_2 a_m_dot_2 a_a_u_2'
    function = '(a_momentum_dot_2 - a_m_dot_2*a_a_u_2)/${A2}'
    outputs = csv_3
  []
  [a_a_T_2]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_2 a_a_rho_2'
     function = 'a_a_p_2*0.029/8.3144598/a_a_rho_2'
     outputs = csv_3
  []
  [a_a_H_2]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_2 a_m_dot_2'
    function = 'a_energy_dot_2/a_m_dot_2'
    outputs = csv_3
  []  
  
#########################################################################  
  
  [a_m_dot_3]
    type = ADFlowJunctionFlux1Phase
    boundary = cold_leg:out
    connection_index = 1
    equation = mass
    junction = junction_cold_leg_3
    outputs = csv_3
  []
  [a_momentum_dot_3]
    type = ADFlowJunctionFlux1Phase
    boundary = cold_leg:out
    connection_index = 1
    equation = momentum
    junction = junction_cold_leg_3
    outputs = csv_3
  []
  [a_energy_dot_3]
    type = ADFlowJunctionFlux1Phase
    boundary = cold_leg:out
    connection_index = 1
    equation = energy
    junction = junction_cold_leg_3
    outputs = csv_3
  []
  [a_a_3]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_3'
    function = '0.5*a_m_dot_3 - 1.4/0.4*a_m_dot_3'
    outputs = csv_3
  []
  [a_b_3]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_3'
    function = '(1.4/0.4*a_momentum_dot_3)'
    outputs = csv_3
  []
  [a_c_3]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_3'
    function = '-a_energy_dot_3'
    outputs = csv_3
  []
  [a_delta_3]
    type = ParsedPostprocessor
    pp_names = 'a_a_3 a_b_3 a_c_3'
    function = 'a_b_3^2 - 4*a_a_3*a_c_3'
    outputs = csv_3
  []
  [a_a_u_3]
    type = ParsedPostprocessor
    pp_names = 'a_a_3 a_b_3 a_c_3 a_delta_3'
    function = '(-a_b_3 + sqrt(a_delta_3))/(2*a_a_3)'
    outputs = csv_3
  []
  [a_a_rho_3]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_3 a_a_u_3'
    function = 'a_m_dot_3/a_a_u_3/${A3}'
    outputs = csv_3
  []
  [a_a_p_3]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_3 a_m_dot_3 a_a_u_3'
    function = '(a_momentum_dot_3 - a_m_dot_3*a_a_u_3)/${A3}'
    outputs = csv_3
  []
  [a_a_T_3]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_3 a_a_rho_3'
     function = 'a_a_p_3*0.029/8.3144598/a_a_rho_3'
     outputs = csv_3
  []
  [a_a_H_3]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_3 a_m_dot_3'
    function = 'a_energy_dot_3/a_m_dot_3'
    outputs = csv_3
  []

######################################################################### 
  
  [a_m_dot_4]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = mass
    junction = turbine
    outputs = csv_3
  []
  [a_momentum_dot_4]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = momentum
    junction = turbine
    outputs = csv_3
  []
  [a_energy_dot_4]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = energy
    junction = turbine
    outputs = csv_3
  []
  [a_a_4]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_4'
    function = '0.5*a_m_dot_4 - 1.4/0.4*a_m_dot_4'
    outputs = csv_3
  []
  [a_b_4]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_4'
    function = '(1.4/0.4*a_momentum_dot_4)'
    outputs = csv_3
  []
  [a_c_4]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_4'
    function = '-a_energy_dot_4'
    outputs = csv_3
  []
  [a_delta_4]
    type = ParsedPostprocessor
    pp_names = 'a_a_4 a_b_4 a_c_4'
    function = 'a_b_4^2 - 4*a_a_4*a_c_4'
    outputs = csv_3
  []
  [a_a_u_4]
    type = ParsedPostprocessor
    pp_names = 'a_a_4 a_b_4 a_c_4 a_delta_4'
    function = '(-a_b_4 + sqrt(a_delta_4))/(2*a_a_4)'
    outputs = csv_3
  []
  [a_a_rho_4]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_4 a_a_u_4'
    function = 'a_m_dot_4/a_a_u_4/${A6}'
    outputs = csv_3
  []
  [a_a_p_4]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_4 a_m_dot_4 a_a_u_4'
    function = '(a_momentum_dot_4 - a_m_dot_4*a_a_u_4)/${A6}'
    outputs = csv_3
  []
  [a_a_T_4]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_4 a_a_rho_4'
     function = 'a_a_p_4*0.029/8.3144598/a_a_rho_4'
     outputs = csv_3
  []
  [a_a_H_4]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_4 a_m_dot_4'
    function = 'a_energy_dot_4/a_m_dot_4'
    outputs = csv_3
  []
  
#########################################################################
   
  [a_m_dot_5]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe7:in
    connection_index = 1
    equation = mass
    junction = turbine
    outputs = csv_3
  []
  [a_momentum_dot_5]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe7:in
    connection_index = 1
    equation = momentum
    junction = turbine
    outputs = csv_3
  []
  [a_energy_dot_5]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe7:in
    connection_index = 1
    equation = energy
    junction = turbine
    outputs = csv_3
  []
  [a_a_5]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_5'
    function = '0.5*a_m_dot_5 - 1.4/0.4*a_m_dot_5'
    outputs = csv_3
  []
  [a_b_5]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_5'
    function = '(1.4/0.4*a_momentum_dot_5)'
    outputs = csv_3
  []
  [a_c_5]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_5'
    function = '-a_energy_dot_5'
    outputs = csv_3
  []
  [a_delta_5]
    type = ParsedPostprocessor
    pp_names = 'a_a_5 a_b_5 a_c_5'
    function = 'a_b_5^2 - 4*a_a_5*a_c_5'
    outputs = csv_3
  []
  [a_a_u_5]
    type = ParsedPostprocessor
    pp_names = 'a_a_5 a_b_5 a_c_5 a_delta_5'
    function = '(-a_b_5 + sqrt(a_delta_5))/(2*a_a_5)'
    outputs = csv_3
  []
  [a_a_rho_5]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_5 a_a_u_5'
    function = 'a_m_dot_5/a_a_u_5/${A7}'
    outputs = csv_3
  []
  [a_a_p_5]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_5 a_m_dot_5 a_a_u_5'
    function = '(a_momentum_dot_5 - a_m_dot_5*a_a_u_5)/${A7}'
    outputs = csv_3
  []
  [a_a_T_5]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_5 a_a_rho_5'
     function = 'a_a_p_5*0.029/8.3144598/a_a_rho_5'
     outputs = csv_3
  []
  [a_a_H_5]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_5 a_m_dot_5'
    function = 'a_energy_dot_5/a_m_dot_5'
    outputs = csv_3
  []

#########################################################################
  
  [a_m_dot_6]
    type = ADFlowBoundaryFlux1Phase
    boundary = outlet
    equation = mass
    outputs = csv_3
  []
  [a_momentum_dot_6]
    type = ADFlowBoundaryFlux1Phase
    boundary = outlet
    equation = momentum
    outputs = csv_3
  []
  [a_energy_dot_6]
    type = ADFlowBoundaryFlux1Phase
    boundary = outlet
    equation = energy
    outputs = csv_3
  []
  [a_a_6]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_6'
    function = '0.5*a_m_dot_6 - 1.4/0.4*a_m_dot_6'
    outputs = csv_3
  []
  [a_b_6]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_6'
    function = '(1.4/0.4*a_momentum_dot_6)'
    outputs = csv_3
  []
  [a_c_6]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_6'
    function = '-a_energy_dot_6'
    outputs = csv_3
  []
  [a_delta_6]
    type = ParsedPostprocessor
    pp_names = 'a_a_6 a_b_6 a_c_6'
    function = 'a_b_6^2 - 4*a_a_6*a_c_6'
    outputs = csv_3
  []
  [a_a_u_6]
    type = ParsedPostprocessor
    pp_names = 'a_a_6 a_b_6 a_c_6 a_delta_6'
    function = '(-a_b_6 + sqrt(a_delta_6))/(2*a_a_6)'
    outputs = csv_3
  []
  [a_a_rho_6]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_6 a_a_u_6'
    function = 'a_m_dot_6/a_a_u_6/${A8}'
    outputs = csv_3
  []
  [a_a_p_6]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_6 a_m_dot_6 a_a_u_6'
    function = '(a_momentum_dot_6 - a_m_dot_6*a_a_u_6)/${A8}'
    outputs = csv_3
  []
  [a_a_T_6]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_6 a_a_rho_6'
     function = 'a_a_p_6*0.029/8.3144598/a_a_rho_6'
     outputs = csv_3
  []
  [a_a_H_6]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_6 a_m_dot_6'
    function = 'a_energy_dot_6/a_m_dot_6'
    outputs = csv_3
  []
  
#########################################################################
# We can prove the accuracy
#########################################################################

  [b_energy_dot_4-a_energy_dot_3]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_4 a_energy_dot_3'
    function = 'a_energy_dot_4 - a_energy_dot_3'
    outputs = csv_3
  []  
  
  # compressor_torque
  [b_(-m_dot_div_omega)*(h_out-h_in)_c]                   
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_1 shaft_speed a_a_H_2 a_a_H_1'
    function = '-a_m_dot_1/shaft_speed*(a_a_H_2 - a_a_H_1)'
    outputs = csv_3
  []
  
  # turbine_torque
  [b_(-m_dot_div_omega)*(h_out-h_in)_t]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_4 shaft_speed a_a_H_5 a_a_H_4'
    function = '-a_m_dot_4/shaft_speed*(a_a_H_5 - a_a_H_4)'
    outputs = csv_3
  []       
[]

[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'
  []
  [csv_2]
    type = CSV
    file_base = 'canshu'
  []
  [csv_3]
    type = CSV
    file_base = 'xincanshu'
  []
  [console]
    type = Console
    show = 'shaft_speed compressor:pressure_ratio turbine:pressure_ratio'
  []
[]

The input file above will give you some promsing result.
What if I want to make generator torque equal to 0 when t > t_loss_generator to simulate load rejection?
I add this into [Functions]

  [if_generator_on]
    type = ParsedFunction
    symbol_names = 'omega'
    symbol_values = 'shaft:omega'
    expression = 'if(t>${t_loss_generator}, 0, omega*${generator_torque_per_shaft_speed})'
  []

And these into [ControlLogic]

  [generator_torque_logic]
    type = ParsedFunctionControl
    function = 'if_generator_on'
    symbol_names = 'if_generator_on'
    symbol_values = 'if_generator_on'
  []
  [generator_ctrl]
    type = SetComponentRealValueControl
    component = generator
    parameter = torque
    value = generator_torque_logic:value
  []

And this into [Components]

  [generator]
    type = ShaftConnectedMotor
    inertia = ${I_generator}
#     torque = generator_torque_fn
    torque = 0 # controlled
  []

But the result shows that the generator torque and motor torque will be both 0. The control for generator torque fails! And it involves motor torque failing! I want to know the reason!
Also, I want to know how to control both the generator torque and motor torque varying with both speed and time?

Answered by Wang-Yihu

Nov 9, 2025

You are right sir. We can set different "placeholder" values for both torque in motor and torque in generator and give them controllogic about shaft:omega, but there will be mistakes for torque value: the 1st step: "placeholder" values; the 2nd step: function about shaft:omega in time step 1 with time; the 3rd step: function about shaft:omega in time step 2 with time... - This is not what I want.

What I want is : torque value = the 1st step: function about shaft:omega in time step 1 with time, the 2nd step: function about shaft:omega in time step 2 with time; the 3rd step: function about shaft:omega in time step 3 with time...

Specifically, I want generator torque = f*omega, t<= t_loss_ge…

View full answer

Wang-Yihu · 2025-11-06T12:22:21Z

Wang-Yihu
Nov 6, 2025
Author

And I also try this:

  [generator_ctrl]
    type = TimeFunctionComponentControl
    component = generator
    parameter = torque
    function = if_generator_on
  []

But the result is the same. The control for generator and motor torque fail again! So how to control both the generator torque and motor torque varying with both speed and time?

0 replies

lindsayad · 2025-11-07T00:06:58Z

lindsayad
Nov 7, 2025
Maintainer

@joshuahansel FYI

0 replies

joshuahansel · 2025-11-07T14:18:48Z

joshuahansel
Nov 7, 2025
Collaborator

Can you provide tau_m_data_with_time.csv?

3 replies

joshuahansel Nov 7, 2025
Collaborator

Also, how old is your executable? This input file seems to be using old code.

Wang-Yihu Nov 7, 2025
Author

Haha, it is quite old. Maybe 2023. tau_m_data_with_time.csv is here

joshuahansel Nov 7, 2025
Collaborator

Oh, sorry I'm blind 👵

joshuahansel · 2025-11-07T16:40:36Z

joshuahansel
Nov 7, 2025
Collaborator

BUG! 🐛 I'm still identifying its nature, but it happens when your "placeholder" torque value in your motor and generator components (which one usually sets to zero) are the same. What is happening is that your motor torque control actually ends up controlling both of the torques. Under the hood, what is happening is that there is a Function object created under the hood when a user supplies a constant for a Function parameter, but apparently these are the same object when they have the same value!

Workaround: use different placeholder values.

I will work on a fix to this and post here afterward, so if you want the proper solution (maybe a week or two from now), maybe you should update your application/input file (it's not 2023 anymore 😉 ). Thanks for telling us the bad behavior!

4 replies

joshuahansel Nov 7, 2025
Collaborator

Let me know if you do update your input file and need help. Here's one tip to start:

change these

  [turb_dissipation_torque]
    type = ScalarVariable
    variable = 'turbine:dissipation_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []

to

  [turb_dissipation_torque]
    type = ElementAverageValue
    variable = dissipation_torque
    block = 'turbine'
    execute_on = 'INITIAL TIMESTEP_END'
  []

joshuahansel Nov 7, 2025
Collaborator

See #31869 for the issue report.

Wang-Yihu Nov 9, 2025
Author

You are right sir. We can set different "placeholder" values for both torque in motor and torque in generator and give them controllogic about shaft:omega, but there will be mistakes for torque value: the 1st step: "placeholder" values; the 2nd step: function about shaft:omega in time step 1 with time; the 3rd step: function about shaft:omega in time step 2 with time... - This is not what I want.

What I want is : torque value = the 1st step: function about shaft:omega in time step 1 with time, the 2nd step: function about shaft:omega in time step 2 with time; the 3rd step: function about shaft:omega in time step 3 with time...

Specifically, I want generator torque = f*omega, t<= t_loss_generator; 0, t > t_loss_generator.

The trouble is I want both time and shaft:omega in that function, however, as https://mooseframework.inl.gov/source/components/ShaftConnectedMotor.html, the time in function will be converted to shaft:omega forcibly. But, I have the solution: I can make something equal to time but pretend not to be time.

So I can add this in [Functions]

  [if_generator_on]
    type = ParsedFunction
    symbol_names = ''
    symbol_values = ''
    expression = 'if(t<=${t_loss_generator}, 1, 0)'
  []

  [generator_torque_fn]
    type = ParsedFunction
    expression = 'slope * if_generator_on * t'
    symbol_names = 'slope if_generator_on'
    symbol_values = '${generator_torque_per_shaft_speed} a_if_generator_on'
  []

add this in [Postprocessors]

  [a_if_generator_on]
    type = FunctionValuePostprocessor
    function = if_generator_on
    outputs = csv_2
  []

Here is the input file, you can get some promising result...

# 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   - 56789 s: Motor torque decreases from 95.7000348 N-m to 0 according to tau_m_data_with_time.csv
#   * 1000 - 8600 s: Power in main heat source increases from 0 - 1000 kW
#   * 56789 s - ?  : The system will be steady with 0 motor torque.

####################################################################################
# Fuel parameters for snerdi
# Diameter:    16.38 mm    , or 17.4 mm
# Height  :       1.2  m
# Fuel number: 30*36 = 1080, or 30*42 = 1260

# Now I choose diameter = 16.38 mm 
# Radial size: fuel,       gas gap, clad,       monolith
#              16.38e-3/2,          16.72e-3/2, 19.0e-3/2

# Volume: 
# Total volume:                pi/4*(16.38e-3)**2*1.2*1080 =  0.273100615 m3
# Total heat transfer surface: pi*(16.38e-3)*1.2*1260 = 82.651532805 m3
####################################################################################

####################################################################################
# inlet - pipe1 - compressor - pipe2    outlet
#
#                               |         |
#                                       
#                             cold_leg  hot_leg  recuperator
#                               |         |
#   --pipe6--turbine--pipe7---  |   -------
#   |                           |
#                               
# pipe5                        pipe3
#
#   |       reactor             |
# pipe4_3 - pipe4_2- pipe4_1---- 
####################################################################################

I_motor = 1.0

I_generator = 1.0

generator_torque_per_shaft_speed = -0.00115

t_loss_generator = 180000

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.60
D1 = 0.65
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.002
d_i = ${fparse 0.015}

L_reactor = 1.2

num_rods_reactor  = 1080
num_rods_recuperator = 570

k_uo2   = 3.4
k_he = 0.3

fuel_radius =           ${fparse 16.38e-3/2}
clad_radius =           ${fparse 16.72e-3/2}
monolith_radius =       ${fparse 19.0e-3/2}
monolith_radius_outer = ${fparse 20.8e-3/2}
th_1 = ${fparse clad_radius-fuel_radius}
th_2 = ${fparse monolith_radius-clad_radius}
th_3 = ${fparse monolith_radius_outer-monolith_radius}

L1 = 5.0
L2 = ${L1}
L3 = ${fparse 2 * L1}
L4 = ${fparse 2 * L1}
L5 = ${L1}
L6 = ${L1}
#L7 = ${fparse L1 + recuperator_width}
L7 = ${fparse L1 + 0.1*L1}
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}

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 = 3.6

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 = 1000000
[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
  []
[]

[HeatStructureMaterials]
  [steel]
    type = SolidMaterialProperties
    rho = 8050
    k = 45
    cp = 466
  []
  [uo2]
    type = SolidMaterialProperties
    rho = 10980
    k = ${k_uo2}
    cp = 330
  []
  [he]
    type = SolidMaterialProperties
    rho = 0.0001216
    k = ${k_he}
    cp = 5192
  []
[]

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

[Functions]

  [time_fn]
    type = ParsedFunction
    expression = t
  []
  
  [if_generator_on]
    type = ParsedFunction
    symbol_names = ''
    symbol_values = ''
    expression = 'if(t<=${t_loss_generator}, 1, 0)'
  []

  ##########################
  # Motor
  ##########################
  
  [tau_m_fn]
    type = PiecewiseLinear
    data_file = 'tau_m_data_with_time.csv'
    format = columns
    extrap = true
    scale_factor = 1.0
  []

  # 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, not active for this example.
  [generator_torque_fn]
    type = ParsedFunction
    expression = 'slope * if_generator_on * t'
    symbol_names = 'slope if_generator_on'
    symbol_values = '${generator_torque_per_shaft_speed} a_if_generator_on'
  []
  # 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
    x = '0.2868480 0.3084640 0.3300760 0.3516920 0.3733040 0.3949200 0.4165200 0.4381600 0.4597600 0.4813600 0.5030000 0.5187200'
    y = '1.5153000 1.5088000 1.5032000 1.4984000 1.4914000 1.4863000 1.4846000 1.4756000 1.4540000 1.4237000 1.3856000 1.3506000'
    extrap = true
  []
  [rp_comp2]
    type = PiecewiseLinear
    x = '0.3674120 0.3890240 0.4106400 0.4322400 0.4523200 0.4754800 0.4970800 0.5177200 0.5442400 0.5658800 0.5874800 0.6090800 0.6307200 0.6484000 0.6621600 0.6759200 0.6886800 0.6994800 0.7093200 0.7181600 0.7253600 0.7338800 0.7407600 0.7466400 0.7517600 0.7595600 0.7633600'
    y = '1.7741000 1.7626000 1.7550000 1.7536000 1.7497000 1.7455000 1.7376000 1.7276000 1.7026000 1.6780000 1.6405000 1.5914000 1.5233000 1.4564000 1.3890000 1.3187000 1.2318000 1.1631000 1.0851000 1.0143000 0.9425300 0.8601900 0.7844300 0.7160400 0.6450700 0.5733400 0.5341100'    
    extrap = true
  []
  [rp_comp3]
    type = PiecewiseLinear
    x = '0.4067200 0.4283200 0.4499200 0.4715600 0.4931600 0.5147600 0.5364000 0.5580000 0.5796400 0.6012400 0.6228400 0.6444800 0.6660800 0.6872000 0.7093200 0.7299200 0.7476400 0.7643200 0.7800400 0.7928000 0.8036400 0.8134400 0.8222800 0.8301600 0.8380000 0.8448800 0.8507600 0.8566800 0.8635600 0.8704400 0.8782800'
    y = '2.1535000 2.1352000 2.1251000 2.1192000 2.1135000 2.1088000 2.1071000 2.1068000 2.1043000 2.0975000 2.0877000 2.0742000 2.0539000 2.0272000 1.9813000 1.9245000 1.8643000 1.7959000 1.7227000 1.6551000 1.5896000 1.5110000 1.4403000 1.3714000 1.2969000 1.2326000 1.1655000 1.0985000 1.0177000 0.9385600 0.8578700' 
    extrap = true
  []
  [rp_comp4]
    type = PiecewiseLinear
    x = '0.4872800 0.5088800 0.5304800 0.5521200 0.5737200 0.5953600 0.6169600 0.6385600 0.6602000 0.6818000 0.7034000 0.7250400 0.7466400 0.7682400 0.7898800 0.8114800 0.8330800 0.8547200 0.8753600 0.8940000 0.9097200 0.9224800 0.9342800 0.9450800 0.9568800'
    y = '2.6634000 2.6485000 2.6366000 2.6248000 2.6184000 2.6124000 2.6110000 2.6102000 2.6110000 2.6102000 2.6110000 2.6099000 2.6060000 2.5989000 2.5874000 2.5642000 2.5309000 2.4824000 2.4202000 2.3500000 2.2832000 2.2195000 2.1437000 2.0706000 1.9947000' 
    extrap = true
  []
  [rp_comp5]
    type = PiecewiseLinear
    x = '0.5658800 0.5874800 0.6090800 0.6307200 0.6523200 0.6739200 0.6955600 0.7171600 0.7387600 0.7604000 0.7820000 0.8036400 0.8252400 0.8468400 0.8684800 0.8900800 0.9116800 0.9333200 0.9549200 0.9765200 0.9981600 1.0197600 1.0413600 1.0600400 1.0757600 1.0904800 1.1032800 1.1150800 1.1248800 1.1327600 1.1406000 1.1484800 1.1563200 1.1632000 1.1690800 1.1750000 1.1808800 1.1867600 1.1916800 1.1976000'
    y = '3.4718000 3.4428000 3.4189000 3.3970000 3.3728000 3.3455000 3.3249000 3.3106000 3.2976000 3.2886000 3.2819000 3.2799000 3.2844000 3.2892000 3.2898000 3.2892000 3.2889000 3.2830000 3.2703000 3.2551000 3.2281000 3.1873000 3.1291000 3.0660000 2.9999000 2.9281000 2.8556000 2.7823000 2.7101000 2.6451000 2.5778000 2.5099000 2.4280000 2.3562000 2.2850000 2.2128000 2.1344000 2.0580000 1.9920000 1.9198000' 
    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
    x = '0.2082520 0.2298640 0.2514800 0.2730960 0.2947080 0.3163240 0.3385240 0.3595520 0.3811640 0.4018000 0.4188400 0.4327600 0.4434800 0.4511200 0.4578000 0.4648800 0.4676400 0.4794000'
    y = '1.1197000 1.1292000 1.1423000 1.1599000 1.1751000 1.1934000 1.2179000 1.2523000 1.2972000 1.3575000 1.4183000 1.4967000 1.5555000 1.6311000 1.7075000 1.7896000 1.8465000 1.9957000'
    extrap = true
  []
  [rp_turb2]
    type = PiecewiseLinear
    x = '0.3674120 0.3890240 0.4106400 0.4322400 0.4538800 0.4705600 0.4922000 0.5079200 0.5196800 0.5275600 0.5354000 0.5422800 0.5491600 0.5590000'
    y = '1.2029000 1.2314000 1.2615000 1.2993000 1.3482000 1.3956000 1.4717000 1.5365000 1.6164000 1.6763000 1.7622000 1.8331000 1.9107000 1.9973000'
    extrap = true
  []
  [rp_turb3]
    type = PiecewiseLinear
    x = '0.4894400 0.5080800 0.5342000 0.5514400 0.5776400 0.5953600 0.6161200 0.6299200 0.6395600 0.6503600 0.6592000 0.6660800 0.6700000 0.6778800 0.6818000 0.6876800 0.6916400 0.6960400 0.6994800 0.7034000 0.7102800 0.7112800 0.7132400 0.7171600'
    y = '1.2856000 1.3172000 1.3532000 1.3853000 1.4625000 1.5312000 1.6149000 1.6750000 1.7413000 1.8149000 1.8875000 1.9559000 2.0013000 2.1282000 2.2221000 2.3390000 2.4564000 2.6010000 2.7379000 2.8729000 3.1802000 3.2705000 3.3605000 3.4468000'
    extrap = true
  []
  [rp_turb4]
    type = PiecewiseLinear
    x = '0.6071200 0.6287600 0.6493600 0.6759200 0.6975200 0.7191200 0.7407600 0.7604000 0.7741600 0.7888800 0.7996800 0.8085200 0.8164000 0.8232800 0.8291600 0.8350800 0.8390000 0.8422800 0.8478400 0.8511600 0.8547200 0.8586400 0.8635600 0.8664800 0.8684800 0.8694400 0.8724000 0.8733600 0.8763200 0.8792800'
    y = '1.3228000 1.3484000 1.3735000 1.4176000 1.4604000 1.5073000 1.5648000 1.6329000 1.6868000 1.7667000 1.8383000 1.9146000 1.9866000 2.0642000 2.1375000 2.2102000 2.2992000 2.3995000 2.5516000 2.6760000 2.8060000 2.9616000 3.1321000 3.2404000 3.3299000 3.4034000 3.4827000 3.5489000 3.6441000 3.7392000'
    extrap = true
  []
  [rp_turb5]
    type = PiecewiseLinear
    x = '0.6975200 0.7191200 0.7407600 0.7623600 0.7839600 0.8056000 0.8304800 0.8488000 0.8704400 0.8920400 0.9136400 0.9352800 0.9549200 0.9716400 0.9853600 0.9971600 1.0079600 1.0168000 1.0244800 1.0296000 1.0364800 1.0404000 1.0443200 1.0459600 1.0480400 1.0492400 1.0531600 1.0576400 1.0610400 1.0639600 1.0669200 1.0705600 1.0736000 1.0755600 1.0787200'
    y = '1.2932000 1.3060000 1.3271000 1.3501000 1.3780000 1.4050000 1.4415000 1.4684000 1.5015000 1.5409000 1.5900000 1.6541000 1.7158000 1.7865000 1.8472000 1.9161000 1.9997000 2.0694000 2.1477000 2.2175000 2.2945000 2.3562000 2.4547000 2.5110000 2.5962000 2.6544000 2.8277000 2.9733000 3.1269000 3.1945000 3.3560000 3.4621000 3.5716000 3.6771000 3.8214000'
    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
    materials = steel
    inner_radius = ${fparse d_i/2}
    num_rods = ${num_rods_recuperator}
  []
  # heat transfer from recuperator to cold leg
  [heat_transfer_cold_leg]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = cold_leg
    hs = recuperator
    hs_side = OUTER
    Hw = 500
  []
  # heat transfer from hot leg to recuperator
  [heat_transfer_hot_leg]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = hot_leg
    hs = recuperator
    hs_side = INNER
    Hw = 500
  []

  [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_1]
    type = JunctionOneToOne1Phase
    connections = 'pipe3:out pipe4_1:in'
  []

  [pipe4_1]
    type = FlowChannel1Phase
    position = '${x4} ${y4} 0'
    orientation = '-1 0 0'
    length = ${fparse (L4-L_reactor)/2}
    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 = '${fparse x4-(L4-L_reactor)/2} ${fparse y4+monolith_radius_outer} 0'
    length = ${L_reactor}
    widths = '${fuel_radius} ${th_1} ${th_2} ${th_3}'
    n_elems = ${n_elems4}
    n_part_elems = '20 2 10 10'
    names =     'fuel  gap  clad    monolith'
    materials = 'uo2   he   steel   steel'
    inner_radius = 0.0
    num_rods = ${fparse num_rods_reactor}
  []
  [total_power]
    type = TotalPower
    power = 0        # controlled
  []
  [heat_generation]
    type = HeatSourceFromTotalPower
    power = total_power
    hs = reactor
    regions = fuel
  []
  [heat_transfer]
    type = HeatTransferFromHeatStructure1Phase
    flow_channel = pipe4_2
    hs = reactor
    hs_side = OUTER
    Hw = 1700
  []

  [junction4_1_4_2]
    type = JunctionOneToOne1Phase
    connections = 'pipe4_1:out pipe4_2:in'
  []

  [pipe4_2]
    type = FlowChannel1Phase
    position = '${fparse x4-(L4-L_reactor)/2} ${y4} 0'
    orientation = '-1 0 0'
    length = ${L_reactor}
    n_elems = ${n_elems4}
    A = ${A4}
  []

  [junction4_2_4_3]
    type = JunctionOneToOne1Phase
    connections = 'pipe4_2:out pipe4_3:in'
  []

  [pipe4_3]
    type = FlowChannel1Phase
    position = '${fparse x4-(L4+L_reactor)/2} ${y4} 0'
    orientation = '-1 0 0'
    length = ${fparse (L4-L_reactor)/2}
    n_elems = ${n_elems4}
    A = ${A4}
  []
  
  [junction4_3_5]
    type = JunctionOneToOne1Phase
    connections = 'pipe4_3:out pipe5:in'
  []
  
  # 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
  # motor_control 7
  [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_control 1
  [motor]
    type = ShaftConnectedMotor
    inertia = ${I_motor}
    torque = 0.0000002015151314 # 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]
   [motor_torque_logic]
    type = ParsedFunctionControl
    function = 'tau_m_fn'
    symbol_names = 'tau_m_fn'
    symbol_values = 'tau_m_fn'
  []
  [motor_ctrl]
    type = SetComponentRealValueControl
    component = motor
    parameter = torque
    value = motor_torque_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
  []
  
[]

[Postprocessors]

  [a_if_generator_on]
    type = FunctionValuePostprocessor
    function = if_generator_on
    outputs = csv_2
  []
  
  ##########################
  # 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
  # motor_control 6
  [shaft_speed]
    type = ScalarVariable
    variable = 'shaft:omega'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [shaft_speed_another]
    type = ScalarVariable
    variable = 'shaft:omega'
    execute_on = 'TIMESTEP_BEGIN'    
  []
  # speed in RPM
  # motor_control 5
  [shaft_RPM]
    type = ParsedPostprocessor
    pp_names = 'shaft_speed'
    function = '(shaft_speed * 60) /( 2 * ${fparse pi})'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  ##########################
  # Compressor
  ##########################
  [comp_dissipation_torque]
    type = ScalarVariable
    variable = 'compressor:dissipation_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [comp_isentropic_torque]
    type = ScalarVariable
    variable = 'compressor:isentropic_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [comp_friction_torque]
    type = ScalarVariable
    variable = 'compressor:friction_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [compressor_torque]
    type = ParsedPostprocessor
    pp_names = 'comp_dissipation_torque comp_isentropic_torque comp_friction_torque'
    function = 'comp_dissipation_torque + comp_isentropic_torque + comp_friction_torque'
  []

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

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

  [turb_dissipation_torque]
    type = ScalarVariable
    variable = 'turbine:dissipation_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turb_isentropic_torque]
    type = ScalarVariable
    variable = 'turbine:isentropic_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turb_friction_torque]
    type = ScalarVariable
    variable = 'turbine:friction_torque'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [turbine_torque]
    type = ParsedPostprocessor
    pp_names = 'turb_dissipation_torque turb_isentropic_torque turb_friction_torque'
    function = 'turb_dissipation_torque + turb_isentropic_torque + turb_friction_torque'
  []

  [mfr_turb]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = mass
    junction = turbine
  []

  ##########################
  #These parameters are not correct for transient
  #Only for steady!
  # C: compressor, R: recuperator, H: Heat section, T: turbine
  #       6
  #       |
  # 1-C-2-R-3-H-4-T
  #       |       |
  #       ---5----|
  ##########################
  [a_T_1]
    type = SideAverageValue
    variable = T
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_T_2]
    type = SideAverageValue
    variable = T
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_T_3]
    type = SideAverageValue
    variable = T
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_T_4]
    type = SideAverageValue
    variable = T
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_T_5]
    type = SideAverageValue
    variable = T
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_T_6]
    type = SideAverageValue
    variable = T
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_p_1]
    type = SideAverageValue
    variable = p
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_p_2]
    type = SideAverageValue
    variable = p
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_p_3]
    type = SideAverageValue
    variable = p
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_p_4]
    type = SideAverageValue
    variable = p
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_p_5]
    type = SideAverageValue
    variable = p
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_p_6]
    type = SideAverageValue
    variable = p
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_v_1]
    type = SideAverageValue
    variable = v
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_v_2]
    type = SideAverageValue
    variable = v
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_v_3]
    type = SideAverageValue
    variable = v
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_v_4]
    type = SideAverageValue
    variable = v
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_v_5]
    type = SideAverageValue
    variable = v
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_v_6]
    type = SideAverageValue
    variable = v
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_H_1]
    type = SideAverageValue
    variable = H
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_H_2]
    type = SideAverageValue
    variable = H
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_H_3]
    type = SideAverageValue
    variable = H
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_H_4]
    type = SideAverageValue
    variable = H
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_H_5]
    type = SideAverageValue
    variable = H
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_H_6]
    type = SideAverageValue
    variable = H
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_T_recuperator_w_inner]
    type = SideAverageValue
    variable = T_solid
    boundary = recuperator:inner
    outputs = csv_2
  []
  [a_T_recuperator_w_outer]
    type = SideAverageValue
    variable = T_solid
    boundary = recuperator:outer
    outputs = csv_2
  []
  [a_T_reactor_w_inner]
    type = SideAverageValue
    variable = T_solid
    boundary = reactor:inner
    outputs = csv_2
  []
  [a_T_reactor_w_outer]
    type = SideAverageValue
    variable = T_solid
    boundary = reactor:outer
    outputs = csv_2
  []  
  [a_vel_x_1]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_vel_x_2]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_vel_x_3]
    type = SideAverageValue
    variable = vel_x
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_vel_x_4]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_vel_x_5]
    type = SideAverageValue
    variable = vel_x
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_vel_x_6]
    type = SideAverageValue
    variable = vel_x
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_vel_y_1]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_vel_y_2]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_vel_y_3]
    type = SideAverageValue
    variable = vel_y
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_vel_y_4]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_vel_y_5]
    type = SideAverageValue
    variable = vel_y
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_vel_y_6]
    type = SideAverageValue
    variable = vel_y
    boundary = hot_leg:out
    outputs = csv_2     
  []
  [a_vel_z_1]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe1:in
    outputs = csv_2
  []
  [a_vel_z_2]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe2:in
    outputs = csv_2     
  []
  [a_vel_z_3]
    type = SideAverageValue
    variable = vel_z
    boundary = cold_leg:out
    outputs = csv_2     
  []
  [a_vel_z_4]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe4_2:out
    outputs = csv_2     
  []
  [a_vel_z_5]
    type = SideAverageValue
    variable = vel_z
    boundary = pipe7:in
    outputs = csv_2     
  []
  [a_vel_z_6]
    type = SideAverageValue
    variable = vel_z
    boundary = hot_leg:out
    outputs = csv_2     
  []
############################################################
# The precise way to retrieve parameters
############################################################
# rho*u*A = m_dot
# m_dot*u + p*A = n_dot
# m_dot*cp*T + 0.5*m_dot*u*u = q_not
# p/rho = R/M*T

# So, p = (n_dot-m_dot*u)/A, T = p*M/R/rho, rho = m_dot/u/A
# input them into Eq 3
# (0.5*m_dot - cp*M/R*m_dot)*u^2 + (cp*M/R_u*n_dot)*u - q_dot = 0
# (0.5*m_dot - 1.4/0.4*m_dot)*u^2 + (1.4/0.4*n_dot)*u - q_dot = 0       
  [a_m_dot_1]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = mass
    junction = compressor
    outputs = csv_3
  []
  [a_momentum_dot_1]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = momentum
    junction = compressor
    outputs = csv_3
  []
  [a_energy_dot_1]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe1:out
    connection_index = 0
    equation = energy
    junction = compressor
    outputs = csv_3
  []
  [a_a_1]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_1'
    function = '0.5*a_m_dot_1 - 1.4/0.4*a_m_dot_1'
    outputs = csv_3
  []
  [a_b_1]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_1'
    function = '(1.4/0.4*a_momentum_dot_1)'
    outputs = csv_3
  []
  [a_c_1]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_1'
    function = '-a_energy_dot_1'
    outputs = csv_3
  []
  [a_delta_1]
    type = ParsedPostprocessor
    pp_names = 'a_a_1 a_b_1 a_c_1'
    function = 'a_b_1^2 - 4*a_a_1*a_c_1'
    outputs = csv_3
  []
  [a_a_u_1]
    type = ParsedPostprocessor
    pp_names = 'a_a_1 a_b_1 a_c_1 a_delta_1'
    function = '(-a_b_1 + sqrt(a_delta_1))/(2*a_a_1)'
    outputs = csv_3
  []
  [a_a_rho_1]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_1 a_a_u_1'
    function = 'a_m_dot_1/a_a_u_1/${A1}'
    outputs = csv_3
  []
  [a_a_p_1]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_1 a_m_dot_1 a_a_u_1'
    function = '(a_momentum_dot_1 - a_m_dot_1*a_a_u_1)/${A1}'
    outputs = csv_3
  []
  [a_a_T_1]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_1 a_a_rho_1'
     function = 'a_a_p_1*0.029/8.3144598/a_a_rho_1'
     outputs = csv_3
  []
  [a_a_H_1]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_1 a_m_dot_1'
    function = 'a_energy_dot_1/a_m_dot_1'
    outputs = csv_3
  []
#########################################################################
    
  [a_m_dot_2]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe2:in
    connection_index = 1
    equation = mass
    junction = compressor
    outputs = csv_3
  []
  [a_momentum_dot_2]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe2:in
    connection_index = 1
    equation = momentum
    junction = compressor
    outputs = csv_3
  []
  [a_energy_dot_2]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe2:in
    connection_index = 1
    equation = energy
    junction = compressor
    outputs = csv_3
  []
  [a_a_2]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_2'
    function = '0.5*a_m_dot_2 - 1.4/0.4*a_m_dot_2'
    outputs = csv_3
  []
  [a_b_2]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_2'
    function = '(1.4/0.4*a_momentum_dot_2)'
    outputs = csv_3
  []
  [a_c_2]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_2'
    function = '-a_energy_dot_2'
    outputs = csv_3
  []
  [a_delta_2]
    type = ParsedPostprocessor
    pp_names = 'a_a_2 a_b_2 a_c_2'
    function = 'a_b_2^2 - 4*a_a_2*a_c_2'
    outputs = csv_3
  []
  [a_a_u_2]
    type = ParsedPostprocessor
    pp_names = 'a_a_2 a_b_2 a_c_2 a_delta_2'
    function = '(-a_b_2 + sqrt(a_delta_2))/(2*a_a_2)'
    outputs = csv_3
  []
  [a_a_rho_2]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_2 a_a_u_2'
    function = 'a_m_dot_2/a_a_u_2/${A2}'
    outputs = csv_3
  []
  [a_a_p_2]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_2 a_m_dot_2 a_a_u_2'
    function = '(a_momentum_dot_2 - a_m_dot_2*a_a_u_2)/${A2}'
    outputs = csv_3
  []
  [a_a_T_2]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_2 a_a_rho_2'
     function = 'a_a_p_2*0.029/8.3144598/a_a_rho_2'
     outputs = csv_3
  []
  [a_a_H_2]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_2 a_m_dot_2'
    function = 'a_energy_dot_2/a_m_dot_2'
    outputs = csv_3
  []  
  
#########################################################################  
  
  [a_m_dot_3]
    type = ADFlowJunctionFlux1Phase
    boundary = cold_leg:out
    connection_index = 1
    equation = mass
    junction = junction_cold_leg_3
    outputs = csv_3
  []
  [a_momentum_dot_3]
    type = ADFlowJunctionFlux1Phase
    boundary = cold_leg:out
    connection_index = 1
    equation = momentum
    junction = junction_cold_leg_3
    outputs = csv_3
  []
  [a_energy_dot_3]
    type = ADFlowJunctionFlux1Phase
    boundary = cold_leg:out
    connection_index = 1
    equation = energy
    junction = junction_cold_leg_3
    outputs = csv_3
  []
  [a_a_3]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_3'
    function = '0.5*a_m_dot_3 - 1.4/0.4*a_m_dot_3'
    outputs = csv_3
  []
  [a_b_3]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_3'
    function = '(1.4/0.4*a_momentum_dot_3)'
    outputs = csv_3
  []
  [a_c_3]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_3'
    function = '-a_energy_dot_3'
    outputs = csv_3
  []
  [a_delta_3]
    type = ParsedPostprocessor
    pp_names = 'a_a_3 a_b_3 a_c_3'
    function = 'a_b_3^2 - 4*a_a_3*a_c_3'
    outputs = csv_3
  []
  [a_a_u_3]
    type = ParsedPostprocessor
    pp_names = 'a_a_3 a_b_3 a_c_3 a_delta_3'
    function = '(-a_b_3 + sqrt(a_delta_3))/(2*a_a_3)'
    outputs = csv_3
  []
  [a_a_rho_3]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_3 a_a_u_3'
    function = 'a_m_dot_3/a_a_u_3/${A3}'
    outputs = csv_3
  []
  [a_a_p_3]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_3 a_m_dot_3 a_a_u_3'
    function = '(a_momentum_dot_3 - a_m_dot_3*a_a_u_3)/${A3}'
    outputs = csv_3
  []
  [a_a_T_3]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_3 a_a_rho_3'
     function = 'a_a_p_3*0.029/8.3144598/a_a_rho_3'
     outputs = csv_3
  []
  [a_a_H_3]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_3 a_m_dot_3'
    function = 'a_energy_dot_3/a_m_dot_3'
    outputs = csv_3
  []

######################################################################### 
  
  [a_m_dot_4]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = mass
    junction = turbine
    outputs = csv_3
  []
  [a_momentum_dot_4]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = momentum
    junction = turbine
    outputs = csv_3
  []
  [a_energy_dot_4]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe6:out
    connection_index = 0
    equation = energy
    junction = turbine
    outputs = csv_3
  []
  [a_a_4]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_4'
    function = '0.5*a_m_dot_4 - 1.4/0.4*a_m_dot_4'
    outputs = csv_3
  []
  [a_b_4]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_4'
    function = '(1.4/0.4*a_momentum_dot_4)'
    outputs = csv_3
  []
  [a_c_4]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_4'
    function = '-a_energy_dot_4'
    outputs = csv_3
  []
  [a_delta_4]
    type = ParsedPostprocessor
    pp_names = 'a_a_4 a_b_4 a_c_4'
    function = 'a_b_4^2 - 4*a_a_4*a_c_4'
    outputs = csv_3
  []
  [a_a_u_4]
    type = ParsedPostprocessor
    pp_names = 'a_a_4 a_b_4 a_c_4 a_delta_4'
    function = '(-a_b_4 + sqrt(a_delta_4))/(2*a_a_4)'
    outputs = csv_3
  []
  [a_a_rho_4]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_4 a_a_u_4'
    function = 'a_m_dot_4/a_a_u_4/${A6}'
    outputs = csv_3
  []
  [a_a_p_4]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_4 a_m_dot_4 a_a_u_4'
    function = '(a_momentum_dot_4 - a_m_dot_4*a_a_u_4)/${A6}'
    outputs = csv_3
  []
  [a_a_T_4]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_4 a_a_rho_4'
     function = 'a_a_p_4*0.029/8.3144598/a_a_rho_4'
     outputs = csv_3
  []
  [a_a_H_4]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_4 a_m_dot_4'
    function = 'a_energy_dot_4/a_m_dot_4'
    outputs = csv_3
  []
  
#########################################################################
   
  [a_m_dot_5]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe7:in
    connection_index = 1
    equation = mass
    junction = turbine
    outputs = csv_3
  []
  [a_momentum_dot_5]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe7:in
    connection_index = 1
    equation = momentum
    junction = turbine
    outputs = csv_3
  []
  [a_energy_dot_5]
    type = ADFlowJunctionFlux1Phase
    boundary = pipe7:in
    connection_index = 1
    equation = energy
    junction = turbine
    outputs = csv_3
  []
  [a_a_5]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_5'
    function = '0.5*a_m_dot_5 - 1.4/0.4*a_m_dot_5'
    outputs = csv_3
  []
  [a_b_5]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_5'
    function = '(1.4/0.4*a_momentum_dot_5)'
    outputs = csv_3
  []
  [a_c_5]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_5'
    function = '-a_energy_dot_5'
    outputs = csv_3
  []
  [a_delta_5]
    type = ParsedPostprocessor
    pp_names = 'a_a_5 a_b_5 a_c_5'
    function = 'a_b_5^2 - 4*a_a_5*a_c_5'
    outputs = csv_3
  []
  [a_a_u_5]
    type = ParsedPostprocessor
    pp_names = 'a_a_5 a_b_5 a_c_5 a_delta_5'
    function = '(-a_b_5 + sqrt(a_delta_5))/(2*a_a_5)'
    outputs = csv_3
  []
  [a_a_rho_5]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_5 a_a_u_5'
    function = 'a_m_dot_5/a_a_u_5/${A7}'
    outputs = csv_3
  []
  [a_a_p_5]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_5 a_m_dot_5 a_a_u_5'
    function = '(a_momentum_dot_5 - a_m_dot_5*a_a_u_5)/${A7}'
    outputs = csv_3
  []
  [a_a_T_5]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_5 a_a_rho_5'
     function = 'a_a_p_5*0.029/8.3144598/a_a_rho_5'
     outputs = csv_3
  []
  [a_a_H_5]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_5 a_m_dot_5'
    function = 'a_energy_dot_5/a_m_dot_5'
    outputs = csv_3
  []

#########################################################################
  
  [a_m_dot_6]
    type = ADFlowBoundaryFlux1Phase
    boundary = outlet
    equation = mass
    outputs = csv_3
  []
  [a_momentum_dot_6]
    type = ADFlowBoundaryFlux1Phase
    boundary = outlet
    equation = momentum
    outputs = csv_3
  []
  [a_energy_dot_6]
    type = ADFlowBoundaryFlux1Phase
    boundary = outlet
    equation = energy
    outputs = csv_3
  []
  [a_a_6]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_6'
    function = '0.5*a_m_dot_6 - 1.4/0.4*a_m_dot_6'
    outputs = csv_3
  []
  [a_b_6]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_6'
    function = '(1.4/0.4*a_momentum_dot_6)'
    outputs = csv_3
  []
  [a_c_6]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_6'
    function = '-a_energy_dot_6'
    outputs = csv_3
  []
  [a_delta_6]
    type = ParsedPostprocessor
    pp_names = 'a_a_6 a_b_6 a_c_6'
    function = 'a_b_6^2 - 4*a_a_6*a_c_6'
    outputs = csv_3
  []
  [a_a_u_6]
    type = ParsedPostprocessor
    pp_names = 'a_a_6 a_b_6 a_c_6 a_delta_6'
    function = '(-a_b_6 + sqrt(a_delta_6))/(2*a_a_6)'
    outputs = csv_3
  []
  [a_a_rho_6]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_6 a_a_u_6'
    function = 'a_m_dot_6/a_a_u_6/${A8}'
    outputs = csv_3
  []
  [a_a_p_6]
    type = ParsedPostprocessor
    pp_names = 'a_momentum_dot_6 a_m_dot_6 a_a_u_6'
    function = '(a_momentum_dot_6 - a_m_dot_6*a_a_u_6)/${A8}'
    outputs = csv_3
  []
  [a_a_T_6]
     type = ParsedPostprocessor
     pp_names = 'a_a_p_6 a_a_rho_6'
     function = 'a_a_p_6*0.029/8.3144598/a_a_rho_6'
     outputs = csv_3
  []
  [a_a_H_6]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_6 a_m_dot_6'
    function = 'a_energy_dot_6/a_m_dot_6'
    outputs = csv_3
  []
  
#########################################################################
# We can prove the accuracy
#########################################################################

  [b_energy_dot_4-a_energy_dot_3]
    type = ParsedPostprocessor
    pp_names = 'a_energy_dot_4 a_energy_dot_3'
    function = 'a_energy_dot_4 - a_energy_dot_3'
    outputs = csv_3
  []  
  
  # compressor_torque
  [b_(-m_dot_div_omega)*(h_out-h_in)_c]                   
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_1 shaft_speed a_a_H_2 a_a_H_1'
    function = '-a_m_dot_1/shaft_speed*(a_a_H_2 - a_a_H_1)'
    outputs = csv_3
  []
  
  # turbine_torque
  [b_(-m_dot_div_omega)*(h_out-h_in)_t]
    type = ParsedPostprocessor
    pp_names = 'a_m_dot_4 shaft_speed a_a_H_5 a_a_H_4'
    function = '-a_m_dot_4/shaft_speed*(a_a_H_5 - a_a_H_4)'
    outputs = csv_3
  []       
[]

[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'
  []
  [csv_2]
    type = CSV
    file_base = 'canshu'
  []
  [csv_3]
    type = CSV
    file_base = 'xincanshu'
  []
  [console]
    type = Console
    show = 'shaft_speed compressor:pressure_ratio turbine:pressure_ratio'
  []
[]

Answer selected by Wang-Yihu

joshuahansel Nov 10, 2025
Collaborator

I think that should work, but note there is also ParsedFunctionControl, which may be a more elegant solution.

How to control both the generator torque and motor torque varying with both speed and time? #31855

Uh oh!

Uh oh!

Wang-Yihu Nov 6, 2025

Replies: 4 comments · 7 replies

Uh oh!

Wang-Yihu Nov 6, 2025 Author

Uh oh!

lindsayad Nov 7, 2025 Maintainer

Uh oh!

joshuahansel Nov 7, 2025 Collaborator

Uh oh!

joshuahansel Nov 7, 2025 Collaborator

Uh oh!

Wang-Yihu Nov 7, 2025 Author

Uh oh!

joshuahansel Nov 7, 2025 Collaborator

Uh oh!

joshuahansel Nov 7, 2025 Collaborator

Uh oh!

joshuahansel Nov 7, 2025 Collaborator

Uh oh!

joshuahansel Nov 7, 2025 Collaborator

Uh oh!

Uh oh!

Wang-Yihu Nov 9, 2025 Author

Uh oh!

joshuahansel Nov 10, 2025 Collaborator

Wang-Yihu
Nov 6, 2025

Replies: 4 comments 7 replies

Wang-Yihu
Nov 6, 2025
Author

lindsayad
Nov 7, 2025
Maintainer

joshuahansel
Nov 7, 2025
Collaborator

joshuahansel Nov 7, 2025
Collaborator

Wang-Yihu Nov 7, 2025
Author

joshuahansel Nov 7, 2025
Collaborator

joshuahansel
Nov 7, 2025
Collaborator

joshuahansel Nov 7, 2025
Collaborator

joshuahansel Nov 7, 2025
Collaborator

Wang-Yihu Nov 9, 2025
Author

joshuahansel Nov 10, 2025
Collaborator