Three-Dimensional Transient Simulation of Supercritical RP-3 Pyrolysis and Flow Maldistribution in Parallel Regenerative Cooling Channels

1 Introduction

Under extreme thermal environments and high heat-flux conditions, the hot wall of a regenerative cooling thermal-protection system is subjected to a severe thermal load; the total temperature of the gas stream on the hot side can exceed 2800 K [1,2], so conventional passive thermal protection is insufficient for long-duration service. Constrained by the external coolant supply, system volume, and allowable pressure drop, developing an efficient and reliable thermal protection and management scheme remains a key bottleneck. Regenerative cooling introduces hydrocarbon fuel to the vicinity of a hot wall, utilizing both its physical heat sink and the chemical heat sink provided by high-temperature endothermic reactions [3]. This approach simultaneously achieves structural cooling and fuel preheating and is considered an effective active-cooling route for hot-end components in extreme heat environments.

In typical regenerative cooling structures, hydrocarbon fuel flows through multiple parallel cooling channels under supercritical conditions. Owing to manifold effects, channel coupling, and nonuniform wall heat flux, flow redistribution readily develops among channels [4,5,6,7]. Such maldistribution may drive some channels toward a low-flow/high-heat-flux state, resulting in local overtemperature and elevated thermal failure risk, whereas other channels remain overcooled, reducing overall heat sink utilization. Therefore, clarifying the flow distribution characteristics of parallel channels under nonuniform heat flux is of substantial engineering significance. Moreover, hydrocarbon fuels undergo pronounced thermal cracking under supercritical high-temperature and high-pressure conditions, and their thermophysical properties and heat-transfer characteristics vary sharply with temperature, pressure, and cracking conversion [8,9,10]. Although pyrolysis enhances the chemical heat sink, the associated compositional evolution also alters density, viscosity, and pressure-drop characteristics, thereby making flow redistribution in parallel systems more complicated. Recent studies have examined this issue from the viewpoints of the channel aspect ratio and resistance characteristics [11,12], topology optimization under nonuniform thermal environments [13], interconnection-based adaptive regulation [14], and hybrid flow-configuration design [15]. These works have demonstrated that structural regulation can improve thermal-state uniformity and heat-sink utilization under steady or quasi-steady conditions. Nevertheless, under transient heating conditions relevant to hypersonic flight, the coupled action of pseudocritical property variation and endothermic pyrolysis may establish a positive feedback loop among the flow rate, temperature, and reaction intensity, causing rapid amplification of maldistribution. The transient spatiotemporal evolution of this process, as well as the mechanism by which an interconnection structure interrupts this feedback loop, remains poorly understood.

Therefore, this work aims to move beyond steady-state analysis by establishing a three-dimensional transient coupled heat-transfer–flow–pyrolysis model to capture the dynamic response of supercritical RP-3 under nonuniform heat flux. The coupling mechanism between the pseudocritical property transition and pyrolysis along the time axis is analyzed, and the manner in which a midstream interconnection structure passively suppresses maldistribution through pressure-driven lateral leakage from the high-heat-flux channel to the low-heat-flux channel during transient evolution is investigated.

2 Numerical Model and Method Validation

The commercial CFD software ANSYS Fluent 2020 R2 was used with a pressure-based transient solver to simulate three-dimensional conjugate heat transfer coupled with chemical reactions.

2.1 Computational Model of Parallel Channels

As shown in Fig. 1a, the computational domain comprises a regenerative cooling structure consisting of three parallel rectangular channels. The three channels are arranged horizontally in parallel with a total length of 1100 mm, comprising a 200 mm adiabatic inlet section, a 200 mm adiabatic outlet section to ensure fully developed flow, and a 700 mm heated section in the middle. A heat-flux boundary is imposed on the bottom wall to represent heating from the combustor wall, while the other walls are treated as adiabatic. A 304 stainless-steel solid wall is wrapped around the channels to account for conjugate conduction. The three channels are connected through inlet and outlet manifolds. At the inlet manifold, a constant total mass flow rate and a fixed inlet temperature are prescribed; at the outlet manifold, a given static pressure is imposed. The working fluid is supercritical RP-3 with a total inlet mass flow rate of 6 g/s and an inlet temperature of 300 K, and the system pressure is maintained at approximately 3 MPa.

As shown in Fig. 1b, under nonuniform heat flux, the heated surfaces of the three channels are subjected to different constant heat-flux densities, which increase from the low-heat-flux channel to the high-heat-flux channel. To quantitatively describe the degree of heat-flux nonuniformity, the heat-flux nonuniformity factor θ is defined as follows: θ=qh−qmqm=qm−qlqm(1) where qh, qm, and ql are the heat-flux densities of the high-, medium-, and low-heat-flux channels, respectively, and qm is the weighted-average heat flux of the three channels. The simulation cases without an interconnection structure are listed in Table 1.

The applied average heat-flux levels were set to 1.5 MW/m2 to represent the typical combustor-wall thermal load reported in scramjet supersonic combustor heat-flux measurements and calculations [16]. A lower level of 1.0 MW/m2 was further considered as a baseline case to examine the mass-flow-rate distribution under nonuniform heating with negligible pyrolysis, thereby enabling a clearer comparison between nonreacting and reacting regimes. The nonuniformity factor θ was chosen to range from 0.05–0.20 to cover mild to strong nonuniform heating [17].

Figure 1: Parallel rectangular channels without interconnection. (a) Three-dimensional geometry of the parallel-channel configuration; (b) channel cross-section and boundary conditions: nonuniform heat flux and gravity direction.

In the cases without an interconnection structure, the three channels communicate only through the inlet and outlet manifolds, and there is no direct lateral flow path between channels. In the cases with an interconnection structure, a rectangular interconnection is added to the heated section (see Fig. 2a,b); the interconnection is a rectangular slot with the same cross-sectional width as the channel, allowing lateral exchange of fluid between adjacent channels and enabling adaptive flow redistribution. The interconnection is placed in the heated section; x denotes the distance from the interconnection center to the inlet of the heated section. Five layouts are considered: x = 100 mm, 200 mm, 350 mm, 500 mm, and 600 mm. The interconnection geometry is kept identical so that the differences are solely due to axial location.

All interconnection-structure cases were conducted under the Q1.5–0.2 condition, which ensures sufficiently strong maldistribution and significant pyrolysis, making it easier to identify the differences in suppression performance among interconnection locations.

Figure 2: Parallel rectangular channels with an interconnection structure. (a) Three-dimensional geometry of the parallel channels with an interconnection structure; (b) streamwise layout and definition of the interconnection location parameters.

2.2 Governing Equations

In the fluid domain, the species transport equations, continuity equation, Reynolds-averaged Navier‒Stokes (RANS) equations, and energy equation are solved simultaneously.

The species transport equation describes the spatial distribution of mass fractions of pyrolysis products as follows: ∂ρYm∂t+∂ρuiYm∂xi=∂∂xiρ(Dm+μtρSct)∂Ym∂xi+Rm(2) where ρ is the fluid density, ui is the velocity component, t is time, Dm is the diffusion coefficient of species m, Ym is the mass fraction of species m, and Rm is the mass source term of species m.

The continuity equation ensures mass conservation:

The momentum equation includes the gravity term and the Reynolds stress term: ∂∂tρui+∂∂xjρuiuj=−∂p∂xi+∂∂xjη∂ui∂xj−ρui′uj′¯+ρgi(4) where p is the pressure, gi is the i-th component of gravitational acceleration, and τ i j = − ρ u i ′ u j ′ ¯ is the Reynolds stress tensor, in which the overbar denotes time averaging and u′i and u′j are the fluctuating velocity components in the i- and j-directions, respectively. The energy equation is written in enthalpy form, explicitly accounting for heat conduction, turbulent thermal diffusion, and the endothermic term of chemical reactions: ∂∂tρh+∂∂xiρuih=∂∂xiλcp∂h∂xi+μtPrt∂h∂xi+S(5) where h is the specific enthalpy, λ is the thermal conductivity, Cp is the specific heat at constant pressure, μt is the turbulent viscosity, Prt is the turbulent Prandtl number, and S is the energy source term representing the energy change within the control volume induced by pyrolysis reactions.

The k-ω SST model was adopted because it has shown good performance in predicting near-wall heat transfer, adverse-pressure-gradient effects, and secondary-flow structures in regenerative-cooling channels with supercritical hydrocarbon fuel [18,19,20]. The model reduces to the standard k-ω formulation in the near-wall region and smoothly blends to the k-ε formulation in the core flow, thereby providing accurate near-wall velocity and temperature gradients with good numerical robustness [21,22]. Turbulent kinetic energy equation:

Specific dissipation rate equation: ∂ρω∂t+∂ρω∂xi+∂ρuiω∂xi=∂∂xiμ+μtσω∂ω∂xi+Gω−Yω+Sωb+Gω(7) where k is the turbulent kinetic energy and ω is the specific dissipation rate; Gk and Gω are the production terms of k and ω; Yk and Yω are the dissipation terms of k and ω; and Sk and Sω are the user-defined source terms in the k and ω equations, respectively. σk and σω are the diffusion coefficients in the transport equations of turbulent kinetic energy and the specific dissipation rate, respectively. For the standard SST k-ω model used in this work, the model constants are σk1 = 0.85, σk2 = 1.00, σω1 = 0.50, σω2 = 0.856, β1 = 0.075, β2 = 0.0828, β* = 0.09, γ1 = 5/9, γ2 = 0.44, and α1 = 0.31. Subscripts 1 and 2 denote the near-wall and outer-flow constants, respectively.

In the solid domain, the energy equation reduces to the heat conduction equation: ρscs∂T∂t=∂∂xiks∂T∂xi(8) where ks is the thermal conductivity of the solid and T is the temperature.

The governing equations were solved using a pressure-based transient solver. Pressure–velocity coupling was achieved using the PISO algorithm. Convective terms were discretized with a second-order upwind scheme, and time integration was performed using a second-order implicit scheme. A constant total mass flow rate and a fixed inlet temperature were specified at the inlet, and a prescribed static pressure was applied at the outlet. The transient time step was set to Δt = 5 × 10−4 s, with a maximum of 20 inner iterations per time step. In addition to residual criteria, key quantities, including the outlet mass flow rate of each channel, total pressure drop, and wall temperature were monitored. A time step was considered to have converged when the residuals of continuity, momentum, energy, and species equations were all below 10−4 and when the relative change in monitored quantities between two consecutive inner iterations was less than 0.5%. For postprocessing, field data were sampled and exported every 0.5 s. To exclude the influence of time-discretization error on the transient-maldistribution evolution, a time-step independence study was conducted by comparing the baseline Δt = 5 × 10−4 s with a smaller time step Δt2 = 1 × 10−4 s. At representative times, the maximum flow-deviation indicator βi and the total pressure drop ΔP were compared. The differences in these key quantities were less than 2%, and the transient trends were consistent. Considering accuracy and computational cost, Δt = 5 × 10−4 s was adopted for all subsequent simulations.

2.3 Chemical Kinetic Model

To reduce computational cost while maintaining accuracy, this study adopted a 6R–17S RP-3 fuel pyrolytic reaction mechanism, which was simplified by Tian et al. [23] from the 24-step RP-3 cracking kinetic model proposed by Jiang et al. [24]. The mechanism consists of one primary cracking reaction and several secondary cracking reactions, and it can accurately predict the temperature evolution and the overall chemical heat sink of RP-3 in heated channels over a wide conversion range of up to approximately 90%. Although the 6R–17S mechanism is reduced and cannot fully reproduce detailed product distributions, its impact on the present study is expected to be limited because the density evolution predicted by different pyrolytic mechanisms is almost identical along the channel, and the relative differences in thermophysical properties are negligible (except for some cases at very high conversion). Tian et al. [23] further reported that the simplified model and baseline model are in close agreement for temperature, flow field, and thermophysical properties when the conversion rate is less than 100% and that the primary reaction dominates the macroscopic effect of pyrolysis (70.8% acceleration and 85.5% chemical heat sink at full cracking). Hence, kinetic-model reduction is unlikely to change the main pressure-drop trend in an engineering-oriented maldistribution analysis focusing on bulk properties and integral resistance. All reactions are expressed in the Arrhenius form, and the corresponding preexponential factors and activation energies are listed in Table 2.

Because the thermophysical properties of RP-3 aviation kerosene and its pyrolysis products vary significantly with temperature at supercritical pressure, a componentwise approach was employed to determine the thermophysical properties. For the base fuel RP-3, the correlations for specific heat at constant pressure, density, thermal conductivity, and dynamic viscosity at approximately 3 MPa were taken directly from the fitted expressions reported by the Dalian University of Technology group [25]. For conventional pyrolysis products, the thermophysical properties were obtained from database-based correlations, while an equivalent-substitute treatment was adopted for the heavy lumped species in the reduced reaction mechanism. During the transient calculation, the local mixture properties were dynamically updated through a UDF according to the local temperature and species composition. The RP-3 property correlations from Beihang University [26,27,28,29] and the Dalian University of Technology group are compared in Fig. 3.

Figure 3: Comparison of RP-3 thermophysical properties at 3 MPa [25,26,27,28,29].

2.4 Grid Independence Verification

Unstructured hexahedral meshes were employed in the simulations, balancing geometric adaptability and numerical accuracy; their adequacy has been confirmed in previous studies. For the parallel-channel model, near-wall mesh refinement is performed with a first-layer thickness of 0.001 mm to satisfy y+ < 1 for the first layer, and boundary-layer meshes are arranged in the wall-normal direction with a geometrically increasing thickness, with a growth rate of 1.2.

Taking a parallel-channel case with an interconnection structure and uniform heat flux as an example, the maximum outlet temperature in the heated section and the inlet-outlet pressure drop are compared for six meshes. As shown in Fig. 4, when the mesh size reaches 1,886,048 cells, the variations in the maximum outlet temperature and the pressure drop become negligible, whereas the curve shapes remain consistent, indicating that mesh independence is essentially achieved. Considering both accuracy and computational cost, this mesh is adopted for subsequent simulations.

Figure 4: Grid independence results.

2.5 Validation of the Numerical Method

To validate the applicability of the coupled turbulence flow–heat transfer–pyrolysis kinetics model with temperature-dependent thermophysical properties developed in this study for predicting the heat transfer of supercritical RP-3 in tubes, in-house experimental data from an electrically heated circular tube at 3 MPa were used for comparison. The experimental system operates as an open-loop circuit, as shown in Fig. 5. The test section uses a 1Cr18Ni9Ti stainless steel circular tube with an outer diameter of 3 mm, a wall thickness of 0.5 mm (yielding an inner diameter of 2 mm), and an effective heating length of 1000 mm. During the experiments, RP-3 aviation kerosene was pressurized to 3 MPa using a plunger metering pump, and the tube wall was directly heated via a low-voltage, high-current DC power supply to simulate the wall heat flux. Considering the temperature measurement accuracy (fuel temperature uncertainty of ±5 K and wall temperature uncertainty of ±2 K) and the uncertainty in the heat-flux boundary condition introduced by electrical measurements and heat-loss calibration, these deviations are within an acceptable range. The net heating power was corrected using an ambient-temperature calibration.

A set of representative operating conditions is summarized in Table 3: the inlet mass flow rate was maintained at 1.25 g/s with an inlet temperature of approximately 304 K. By adjusting the heating power, the net heat-flux density ranged from 316.985–691.487 kW/m2, corresponding to an increase in the outlet temperature from 780.795 K to 1022.021 K. This range spans typical conditions from the near-pseudocritical temperature region to the regime where pyrolysis becomes significant. The numerical simulations employed inlet and thermal boundary conditions consistent with those used in the experiments. The mixed-mean outlet temperature was extracted and compared with the experimental data in Fig. 6a. The results show that the predicted trend of the increase in the outlet temperature with increasing heat-flux density is in agreement with the experimental results, with a maximum relative error of 5.69% and errors below 5.5% for all the other cases, which are within an acceptable range.

In addition to the outlet-temperature validation, the pressure-drop model was further validated using new measurements obtained with the redesigned test section at min = 1.0 g/s (Tin ≈ 307 K, p = 3 MPa, Tout = 323–1073 K). As shown in Fig. 6b, the simulated total pressure drop agrees well with the measurements over the tested temperature range, with relative deviations within ±15%. This additional validation supports the reliability of the pressure-drop prediction used in the subsequent maldistribution analysis.

Further validation was performed for the axial wall-temperature distribution. For brevity, two representative cases with outlet temperatures of approximately 700°C and 750°C are presented in Fig. 6c,d, which are consistent with the outlet-temperature range discussed in the subsequent sections. The pointwise relative error of the wall temperature at 16 axial measurement locations ranges from 0.521% to 11.318% (mean 4.07%) for the Tout ≈ 700°C case and from 0.261% to 7.663% (mean 2.74%) for the Tout ≈ 750°C case. Notably, this study primarily validates the predictive accuracy of the model for supercritical heat transfer, pyrolysis reactions, and property variations. Owing to the generality of the governing equations for mass, momentum, and energy conservation, the model is therefore considered suitable for qualitatively revealing the flow-maldistribution mechanism in parallel channels and the associated transverse crossflow characteristics.

Figure 5: Schematic of the experimental system.

Figure 6: Experimental validation of the numerical model: (a) outlet fuel temperature; (b) pressure drop; (c) axial wall temperature for T_out ≈ 700°C; (d) axial wall temperature for T_out ≈ 750°C.

3 Results and Discussion

3.1 Transient Flow Distribution under Nonuniform Heat Flux

In this section, the transient evolution of the flow distribution in a parallel-channel system under nonuniform heat flux is examined. To distinguish the two response types, namely, property-jump-dominated behavior and pyrolysis-coupled behavior, two groups of cases with average heat fluxes of 1.0 MW/m2 and 1.5 MW/m2 were compared. At 1.0 MW/m2, the maximum bulk temperature along the channel was 778 K, and the RP-3 conversion remained below 1.5% throughout; thus, maldistribution was induced mainly by pseudocritical property jumps. At 1.5 MW/m2, the fuel temperature entered a significant pyrolysis regime, and the conversion increased rapidly with time; the accumulation of pyrolysis products altered the mixture properties and flow resistance, leading to a stronger transient amplification of maldistribution. To evaluate the suppression effect of the interconnection structure, the results without interconnection were used as the baseline, and cases with interconnection under the same inlet and thermal boundary conditions were used for comparison, with a focus on how the interconnection location affected different stages of flow redistribution.

3.1.1 Monitoring Strategy and Evaluation Metrics

To quantitatively characterize the flow distribution, resistance evolution, and heat-transfer/pyrolysis response under nonuniform heating, monitoring quantities, including the mass flow rate, pressure drop, outlet temperature, and species mass fractions, were recorded for each channel. For the no-interconnection cases, channel mass flow rates were extracted at representative cross-sections of the heated section for comparison. Four representative times, t = 6 s, 10 s, 20 s, and 30 s, were selected. The physical quantities are defined as follows:

To describe the flow distribution among parallel channels, the flow deviation coefficient βi is introduced: βi=mi−m¯m¯(9) where mi is the mass flow rate in channel i and m ¯ is the ideal average mass flow rate of all channels.

The total pressure drop of channel i is as follows: ΔPi(t)=ΔPf,i(t)+ΔPa,i(t)(10) where ΔPf is the frictional pressure drop and ΔPa is the acceleration pressure drop.

The acceleration pressure drop ΔPa is calculated as follows:

The frictional pressure drop ΔPf is calculated as follows: ΔPf=∫0Lfi2Dhmi2ρiA2dx(12) where A is the flow area, Dh is the hydraulic diameter, and L is the local heated length.

Accordingly, the equivalent resistance coefficient can be obtained as follows:

Defining the term in braces as the equivalent resistance coefficient Cr,i gives:

The outlet RP-3 mass fraction YRP3,out serves as an indicator of pyrolysis intensity; RP-3 conversion is further defined as follows:

To compare the progression of pyrolysis among channels and its influence on flow resistance and heat transfer.

3.1.2 Formation and Amplification of the Transient Maldistribution without Interconnection

From a transient perspective, as shown in Fig. 7, the flow redistribution in the Q1.0 cases exhibited a consistent three-stage behavior. Stage I (early heating): The channel mass flow rates were nearly uniform. Stage II (9–11 s): A sudden divergence occurred; the high-heat-flux channel mass flow rate decreased rapidly, whereas that of the low-heat-flux channel increased correspondingly. Stage III (t ≥ 15 s): The system reached a new dynamic balance, and the maldistribution stabilized. The onset of maldistribution is strongly dependent on θ: a larger θ triggers earlier onset (from 9.22 s at θ = 0.05 to 8.19 s at θ = 0.20) and a greater magnitude of flow deviation.

The fundamental driving force of this flow redistribution originated from the pronounced divergence of the pressure-drop components. As shown in Fig. 8, within the first 10 s, the frictional pressure drops ΔPf in all channels decreased. This occurred because the initial temperature increase primarily caused a pronounced decrease in the dynamic viscosity μ, which dominated the resistance behavior and reduced the friction factor, as supported by the local Re and μ evolution in Fig. 9.

Figure 7: Transient evolution of mass flow rates in Q1.0 cases for different nonuniformity factors θ: (a) 0.05; (b) 0.10; (c) 0.15; (d) 0.20.

Figure 8: Transient evolution of ΔP_f and ΔP_a in Q1.0 cases for different nonuniformity factors θ: (a) 0.05; (b) 0.10; (c) 0.15; (d) 0.20.

Figure 9: Temporal evolution of the Reynolds number and dynamic viscosity at three axial locations: (a) x = 175 mm; (b) x = 350 mm; (c) x = 525 mm.

However, the actual trigger for maldistribution occurred when the high-heat-flux channel approached the pseudocritical temperature Tpc. The thermodynamic state underwent an abrupt change: the sharp density drop associated with the pseudocritical transition enhanced the velocity gradient and momentum change, causing a nonlinear surge in the acceleration pressure drop ΔPa. Because the pressure drop scales with u2, this inertial effect rapidly overwhelms the viscous decay shown in Fig. 10. Under the parallel constraint of equal total pressure drop, the high-heat-flux channel was forced to reduce its mass flow rate. The field characteristics (Fig. 11) confirm this: at t = 10 s, the stronger density gradient near the outlet of the high-heat-flux channel coincides with the onset of ΔPa deviation. These results indicate that the pseudocritical property jump alone is sufficient to trigger transient maldistribution through divergence in acceleration pressure drop among channels.

Figure 10: Transient evolution of C_r,i in Q1.0 cases for different θ values: (a) 0.05; (b) 0.10; (c) 0.15; (d) 0.20.

Figure 11: Density contours at θ = 0.20 in the Q1.0 case: (a) t = 6 s; (b) t = 10 s; (c) t = 20 s; (d) t = 30 s.

3.1.3 Amplification of Transient Maldistribution under Endothermic Pyrolysis

Building upon the thermophysical trigger mechanism established in Section 3.1.2, the flow redistribution in the Q1.5 cases (Fig. 12) demonstrated the superimposed impact of endothermic pyrolysis. While the macroscopic flow evolution retained a three-stage characteristic, the involvement of chemical reactions fundamentally altered the timeline and the ultimate severity of the maldistribution.

Figure 12: Transient evolution of mass flow rates in Q1.5 cases for different nonuniformity factors θ: (a) 0.05; (b) 0.10; (c) 0.15; (d) 0.20.

The primary distinction is the accelerated onset and the continuous amplification of the deviation. Triggering occurs significantly earlier (e.g., at 4.64 s for θ = 0.20, compared to 8.19 s in the Q1.0 case). Furthermore, instead of reaching a moderate quasi-steady balance, the maldistribution at high nonuniformity (θ > 0.15) became deeply entrenched. In the final stage for θ = 0.20, the high-heat-flux channel mass flow rate decreased by 38.06%, whereas the low-heat-flux channel mass flow rate increased to 2.927 g/s.

This severe deviation originates from the reaction–property coupling, which prevented the equivalent resistance Cr,i from stabilizing. As shown in Fig. 13, after the brief initial friction-reduction phase, Cr,H entered a prolonged rising stage. Unlike the baseline cases where density variation is purely temperature dependent, intensive pyrolysis in the high-heat-flux channel generates extensive low-molecular-weight products. The RP-3 mass–fraction contours (Fig. 14) and density contours (Fig. 15) illustrate this coupling: by t = 10 s, the pronounced RP-3 consumption coincides with a further reduction in density. By t = 40 s, the high-heat-flux channel reaches a local conversion rate of 90.10%. This deep pyrolysis drastically reduces the fluid density beyond the pseudocritical effect, forcing a continuous surge in the flow velocity and accelerating the pressure drop. Consequently, the high-heat-flux channel is trapped in a sustained flow-reduction state, culminating in a locked-in maldistribution that does not recover spontaneously within the simulated time window.

Figure 13: Transient evolution of C_r,i in Q1.5 cases for different θ values: (a) 0.05; (b) 0.10; (c) 0.15; (d) 0.20.

Figure 14: RP-3 mass–fraction contours at θ = 0.20 in the Q1.5 case: (a) t = 6 s; (b) t = 10 s; (c) t = 20 s; (d) t = 40 s.

Figure 15: Density contours at θ = 0.20 in the Q1.5 case: (a) t = 6 s; (b) t = 10 s; (c) t = 20 s; (d) t = 40 s.

3.2 Suppression Mechanism and Performance Evaluation of the Interconnection Structure

3.2.1 Evaluation Approach and Quantitative Indicators

The interconnection structure was installed in the middle of the heated section. It not only altered the axial pressure-drop distribution of each channel but also induced lateral crossflow between channels at the interconnection so that the channel mass flow rate became discontinuous along the axial direction across the interconnection. For clarity, the mass flow rates “before” and “after” the interconnection are defined as the monitored channel mass flow rates at the cross sections upstream and downstream of the interconnection, respectively. To avoid potential errors caused by local recirculation/backflow and strong mixing in the immediate vicinity of the interconnection, the two monitoring sections were placed at a fixed offset from the interconnection centerline, i.e., 175 mm upstream (before) and 175 mm downstream (after). If the flow distribution is evaluated only at a single cross-section, it would be difficult to distinguish the formation of maldistribution from the redistribution introduced by the interconnection. Therefore, the flow deviation coefficient βi and the equivalent resistance coefficient Cr,i are evaluated at both monitoring sections, and cases without interconnection (Case 0) and with interconnection (Case 1) are compared together with internal flow-field details.

3.2.2 Effect of Interconnection Location on Flow Equalization

Fig. 16 shows the transient evolution of fuel temperature in the three channels for different interconnection locations. All cases underwent rapid heating at the beginning and gradually approached a quasi-steady state, but the temperature difference among the three channels was highly sensitive to the interconnection location. With a midstream interconnection (e.g., x = 350 mm), the temperature curves of the three channels converged more strongly, and the temperature difference diminished more substantially. When the interconnection was positioned too far upstream (x = 100 mm) or too far downstream (x = 600 mm), the convergence of the temperature difference was much slower, and the residual difference was larger. This finding indicated that the interconnection is not necessarily more effective when placed further upstream or downstream; its effectiveness depends on whether the intervention timing and redistribution strength match the critical interval of maldistribution development during the transient process.

Fig. 17 shows the direct effect of the interconnection on mass-flow distribution. The interconnection causes a sudden change and redistribution of mass flow rates before and after the interconnection. In the upstream interconnection case (x = 100 mm), mixing is introduced too early before the thermal imbalance is fully established; although it reduces the flow-rate disparity in the inlet section, it fails to suppress the stronger downstream divergence of properties caused by cumulative heating effects. In the midstream interconnection case (x = 350 mm), the redistribution is milder, and the mass flow rate variations in the downstream section are relatively controllable. In the downstream interconnection case (x = 600 mm), redistribution still occurs to some extent, but the overall equalization efficiency decreases because the upstream maldistribution pattern has already formed, resulting in more residual differences in temperature and mass flow rate. These observations indicate that the interconnection does more than simply exchange fluid between channels; more importantly, it modulates the timing and intensity of intervention in the maldistribution-development process during the transient evolution.

Figure 16: Temporal variations in the fuel temperature in three channels for different interconnection locations: x = 100 mm; x = 200 mm; x = 350 mm; x = 500 mm; and x = 600 mm.

Figure 17: Temporal variations in the mass flow rate in three channels for different interconnection locations: x = 100 mm; x = 200 mm; x = 350 mm; x = 500 mm; and x = 600 mm.

3.2.3 Suppression of Flow Maldistribution by the Interconnection Structure

The results indicated that under nonuniform thermal loading, the interconnection region formed a typical cross-channel jet with recirculation vortices. The crossflow was predominantly directed from the high-heat-flux channel toward the low-heat-flux channel. This behavior did not simply represent a transfer of useful coolant away from the most severely heated channel. Instead, it provided a lateral release path for overheated, low-density, and strongly cracked fluid accumulated in the high-heat-flux channel. As a result, the downstream accumulation of thermal nonuniformity, compositional nonuniformity, and hydraulic resistance in the high-heat-flux channel was weakened, and the outlet-state disparities among channels were reduced, thereby improving the overall flow uniformity of the parallel system.

The transient regulation process induced by the interconnection can be divided into three stages. As shown in Fig. 18, the lateral pressure gradient induced by thermal imbalance becomes the driving force for mass transport. Stage I (trigger stage, t ≈ 0–6 s): The three channels undergo rapid heating and property variation; the interchannel pressure difference begins to increase, and a jet starts to form at the interconnection. Stage II (enhancement stage, t ≈ 6–20 s): As the temperature in the high-heat-flux channel increases more rapidly, the density decreases and the resistance divergence intensifies, producing a stronger lateral pressure difference across the interconnection. Consequently, stronger crossflow from the high-heat-flux channel into the low-heat-flux channel is established. Although the high-heat-flux channel locally loses part of its mass flow through the interconnection, this lateral leakage extracts the most overheated and strongly cracked fluid from the deteriorated region, thereby weakening the downstream reinforcement of thermal-state disparity, resistance divergence, and further amplification of maldistribution. Stage III (stabilization stage, t ≈ 20–40 s): Crossflow regulation drives the thermal states and segmented pressure drops of channels toward a new balance; the jet strength and recirculation structure become stable, and the maldistribution indicators approach a quasi-steady plateau.

Notably, while improving flow uniformity, the interconnection also introduced complex local flow structures. These secondary flows enhance turbulent mixing and heat exchange but also generate additional local hydraulic losses through vortex dissipation. The vortex-dissipation-induced additional loss, defined as the excess local pressure drop across the interconnection segment (x = 335–365 mm) relative to that in the no-interconnection case under the same conditions (Q1.5–0.2), is 71.9, 22.7, 168.6, and 538.0 Pa at t = 6, 10, 20, and 40 s, respectively. At t = 40 s, this additional loss accounts for 2.2% of the total pressure drop and remains 1–2 orders of magnitude smaller than the distributed frictional and accelerational components. Therefore, the net effectiveness depends on the competition between the resistance-optimization benefit from lateral mass transfer and the pressure-drop penalty from local vortices. The midstream layout performs best because it provides a sufficient lateral driving pressure difference while keeping the post-redistribution flow disturbance within an acceptable range, avoiding new maldistribution instability induced by excessive local losses.

Figure 18: Secondary flow and streamlines at the interconnection (x = 350 mm) at different times: t = 6 s; t = 10 s; t = 20 s; t = 40 s.

3.2.4 Comprehensive Performance Evaluation of the Interconnection Structure

To evaluate the overall performance of the interconnection structure, four aspects are considered: flow uniformity, thermal-state uniformity, pyrolysis uniformity, and pressure-drop cost. The baseline case without interconnection was compared with the midstream interconnection case. The quantitative metrics are summarized in Table 4. Overall, the interconnection effectively suppresses flow maldistribution during the amplification and quasi-steady stages by weakening the outlet thermal, compositional, and hydraulic disparities generated by the deteriorated high-heat-flux channel, leading to improved thermal-state and pyrolysis uniformity among parallel channels. This improvement is achieved with only a minor increase in the total system pressure drop, indicating a favorable performance trade-off under the Q1.5–0.2 condition.

4 Conclusions

In this study, the dynamic flow distribution characteristics of supercritical RP-3 in parallel regenerative cooling channels under nonuniform heat flux were investigated. By establishing a transient model coupled with endothermic pyrolysis kinetics, the present work revealed that the deterioration of flow distribution is fundamentally driven by a positive-feedback cycle involving pseudocritical property fluctuations and reaction-induced resistance amplification. A midstream interconnection structure can passively suppress maldistribution by inducing pressure-driven lateral leakage from the high-heat-flux channel to the low-heat-flux channel. Rather than directly increasing the mass flow rate of the high-heat-flux channel, this crossflow releases overheated, low-density, and strongly cracked fluid from the most deteriorated region, thereby weakening the downstream accumulation of thermal and compositional nonuniformities and the associated resistance amplification.

• 1.The average heat flux directly determines the dominant mechanism of maldistribution evolution. At lower heat fluxes, maldistribution is driven mainly by differences in acceleration pressure drop induced by variations in pseudocritical properties; at higher heat fluxes, the endothermicity of pyrolysis and the low-density characteristics of products dominate, triggering a strong positive feedback loop of mass flow rate reduction, temperature increase, intensified pyrolysis, and sharp resistance increase. Increasing heat flux significantly accelerated the intensification of this feedback loop, leading to a nonlinear increase in maldistribution magnitude during the quasi-steady stage.
• 2.The maldistribution does not increase linearly with the heat-flux nonuniformity θ but is governed by the exponential nature of chemical kinetics. As θ increases, differences in reaction rate and property gradients among channels are progressively amplified, causing the high-heat-flux channel to enter the deep-pyrolysis regime earlier. High nonuniformity not only accelerated the establishment of the initial maldistribution but also markedly lowered the thermal stability boundary, making it difficult to maintain flow stability when inlet-flow-rate adjustment alone was applied.
• 3.The suppression effectiveness of the interconnection strongly depends on its axial location. The midstream layout best matches the development timing of transient maldistribution and uses the lateral pressure difference between channels to induce adaptive crossflow. This crossflow is directed predominantly from the high-heat-flux channel to the low-heat-flux channel and acts by releasing overheated and strongly cracked fluid from the deteriorated high-heat-flux region, rather than by directly supplementing its mass flow rate. In this way, the downstream accumulation of resistance and outlet-state disparities is weakened, and the positive-feedback chain of maldistribution amplification is effectively interrupted. Compared with the poor equalization of upstream layouts and the delayed response of downstream layouts, the midstream interconnection reduced the flow deviation coefficient by approximately 48.3% in the deep-pyrolysis stage, achieving an optimal balance between the equalization benefit and pressure-drop cost.

The present work is limited to a single rectangular interconnection of fixed geometry; future studies should conduct systematic parametric investigations covering (i) cross-sectional area and area ratio, (ii) the number of interconnections and their spacing, and (iii) shape effects, to further optimize the design for practical regenerative cooling applications.