Numerical Analysis of Supercritical Fuel Cracking in Trapezoidal Rib Channels

1 Introduction

With the development of space transport systems towards high specific impulse and reusability, the thrust chambers of liquid rocket engines face severe heat flux challenges under extreme operating conditions. Utilizing the physical sensible heat and chemical cracking latent heat of on-board hydrocarbon fuels for regenerative cooling has become the mainstream solution for ensuring the structural integrity of propulsion systems [1,2,3]. Although related studies have confirmed the significant contribution of the “chemical heat sink” to enhancing heat transfer capabilities [4,5,6], this process involves extremely complex mechanisms of strong thermal-hydraulic-chemical coupling, posing severe challenges to thermal protection design. On the one hand, the drastic nonlinear changes in physical properties within the pseudo-critical region can easily induce buoyancy effects and heat transfer deterioration [7,8]; on the other hand, the generation of cracking gas products alters flow resistance characteristics, significantly increasing flow complexity [4,5]. To overcome the heat-transfer limitation, structural enhancement strategies such as profiled cooling channels, local geometric disturbance, and other passive channel-design modifications are often introduced to improve near-wall heat transfer. However, these techniques are often accompanied by severe pressure loss or local hotspot accumulation, making them difficult to adapt to non-uniform high heat flux environments [9,10,11,12]. Recent numerical studies have further explored the impact of various rib configurations (e.g., rectangular and triangular ribs) on the hydrodynamic forces and thermal performance in heated channels [13,14]. Research on conjugate heat transfer accompanied by strong property variations and intense chemical reactions remains relatively scarce, and there is an urgent need to explore novel structures that balance efficient heat transfer with low flow resistance penalties [15]. However, how to find the optimal balance between increased flow resistance, improved cracking conversion efficiency and enhanced heat transfer via turbulence structure introduction remains to be further explored. Against this background, this paper aims to investigate the comprehensive effects of an improved microchannel structure on the cracking and conjugate heat transfer characteristics of supercritical n-decane. Based on the OpenFOAM open-source calculation framework [16], this study develops a high-fidelity numerical solver suitable for the multi-field coupling of supercritical fluids. While traditional CFD methods face inherent challenges in the pseudo-critical region, emerging techniques such as deep learning-based flow field reconstruction and dual-attention FlowNet models have shown great potential in accurately reconstructing complex flow fields of supercritical fluids [17,18], and such data-driven approaches provide an efficient new path for the analysis of supercritical fluid flow fields. Meanwhile, considering that the focus of this paper is to reveal the intrinsic coupling mechanisms between fuel cracking and heat transfer, it is necessary to adopt research methods with explicit physical interpretability. Therefore, to solve the problems of computational oscillation and low prediction accuracy caused by sharp property changes in the pseudo-critical region, this paper independently develops a high-precision piecewise variable property transport model (for dynamic viscosity and thermal conductivity) based on an 8-segment 7th-order polynomial. This model is strongly coupled with the Peng-Robinson (PR) real gas equation of state [19] and chemical reaction kinetic mechanisms. Meanwhile, the k-ω SST model [20,21] is introduced to close the turbulence terms, thereby establishing a complete numerical calculation method. With n-decane as the research object, this study uses the developed solver to compare and analyze the flow characteristics, cracking conversion rates, and wall temperature distributions of a baseline straight channel and an improved trapezoidal ribbed channel under different mass flow rates. Finally, the engineering application potential of the improved structure is quantitatively evaluated using the thermal-hydraulic Performance Evaluation Criterion (PEC) [22]. Compared with previous studies, the present work provides three main contributions. On the one hand, a high-fidelity numerical solver is developed to simulate supercritical hydrocarbon fuel cracking with strong thermal-hydraulic-chemical coupling and improved numerical stability under steady-state conditions. In addition, an 8-segment 7th-order piecewise transport-property model for dynamic viscosity and thermal conductivity is established to improve the prediction accuracy in the pseudo-critical region. On the other hand, the trapezoidal ribbed channel is designed to enhance heat transfer using a limited number of rib units, which differs from conventional full-bottom distributed disturbance layouts. Therefore, this study provides not only a reliable numerical tool, but also a more compact enhancement strategy for regenerative cooling channels.

2 Development of Numerical Solver for Supercritical Reacting Flows

To address the conjugate heat transfer challenges induced by drastic property variations coupled with chemical reactions during the flow of supercritical hydrocarbon fuels in microchannels, a high-fidelity numerical solver was independently developed based on the OpenFOAM open-source framework. The numerical methods constructed within the solver including the reconstruction of governing equations for multi-component reacting flows, the development and integration of a high-precision piecewise variable property transport model, the numerical implementation strategy for simplified cracking reaction kinetics, and the construction of a flow-thermal-chemical multi-physics strong coupling algorithm were proposed. Finally, the reliability of the solver is comprehensively verified through multi-level test cases ranging from thermophysical property validation to reacting flow comparisons.

2.1 Reconstruction and Closure of Governing Equations for Multi-Component Reacting Flows

To accurately characterize the strong gradient variations of supercritical fluids in the pseudo-critical region, the conservation equations for mass, momentum, energy, and species transport were reconstructed within the fluid domain solver module. On this basis, source terms for cracking reactions were explicitly introduced to close the chemical reaction process.

The governing equations for the fluid domain are established as follows, including the mass conservation equation, momentum conservation equation, energy conservation equation and species transport equation: ∇(ρv)=0(1) ∇·{ρvv}=−∇p+∇·[μ(∇v+∇vT)]−∇·(ρv′v′¯)+ρg(2) ∇(ρvh)+∇(ρv12|v|2)=∇(αh,eff ∇h)+Qdot+ρ(v·g)(3) ∇·(ρvYi)=∇·(αi,eff∇Yi)+ω˙i(4) where ρ refers to fluid density, kg/m3, v is fluid velocity vector, m/s, μ dynamic viscosity, Pa·s, g is gravitational acceleration, m/s2, p is fluid pressure, Pa, h is sensible enthalpy, J/kg, Qdot is heat source term due to chemical reaction per unit volume, W/m3, Yi is mass fraction of species, and ω ˙ i is chemical reaction source term of species i, kg/(m3·s). ρ v ′ v ′ ¯ denotes the turbulent Reynolds stress tensor introduced by Reynolds averaging.

For convenience in describing heat diffusion in the fluid region and maintaining a unified formulation for conjugate heat transfer between the fluid and solid domains, α h , eff is defined as the effective enthalpy diffusion coefficient, kg/(m·s). It includes both molecular and turbulent contributions: αh,eff=λfcp,f+μtPrt(5) where λ f denotes the fluid thermal conductivity, W/(m·K). c p , f denotes the specific heat capacity of the fluid at constant pressure, J/(kg·K). μ t the turbulent dynamic viscosity, kg/(m·s). P r t the turbulent Prandtl number.

Similarly, α i , eff denotes the effective diffusion coefficient of species i, kg/(m·s), including both molecular and turbulent diffusion contributions: αi,eff=ρDi,m+μtSct(6) where D i , m is the molecular diffusion coefficient of species i, m2/s. S c t is the turbulent Schmidt number.

For the solid computational domain, the heat conduction equation is expressed in terms of specific enthalpy, so that energy exchange at the fluid–solid interface can be treated in a unified manner within the enthalpy-based conjugate heat transfer formulation: ∇·(αs∇hs)=0(7) where h s is the solid enthalpy, J/kg, and α s is the solid enthalpy diffusion coefficient, kg/(m·s). It is defined as: αs=λscp,s(8) where λ s is the thermal conductivity of the solid, W/(m·K), and c p , s is the specific heat capacity of the solid at constant pressure, J/(kg·K).

Turbulence closure was achieved using the k-ω SST model [20,21]. The hydraulic-diameter-based Reynolds number (Re) for both inner and outer channels is greater than 4900 across the entire flow field (including the inlet liquid-like region), indicating that the flow is fully turbulent throughout the computational domain. This model combines the solution advantages of both the near-wall and far-field regions and has been verified to accurately capture the near-wall flow details and complex heat transfer phenomena (such as buoyancy-induced effects) [23,24] of supercritical fluids under conditions of drastic property variations.

2.2 Independent Development of High-Precision Piecewise Variable Property Transport Model

2.2.1 Equation of State and Thermodynamic Properties

Based on the constant pressure assumption, the Peng-Robinson (PR) equation of state [19] was adopted in this study to calculate thermophysical properties. Thermodynamic properties were obtained by superimposing departure functions onto ideal gas properties (calculated using the built-in JANAF polynomials in OpenFOAM) [7].

2.2.2 Variable Property Transport Model

In the vicinity of the pseudo-critical region, the transport properties (dynamic viscosity and thermal conductivity) of supercritical hydrocarbon fuels exhibit strong nonlinear variations, which significantly affect the stability and accuracy of numerical simulations for conjugate heat transfer and fluid flow. To overcome the limitations of conventional models in describing such drastic property variations, a high-fidelity variable-property transport model suitable for constant-pressure processes was independently developed. This model employs a refined piecewise fitting strategy, in which the transport properties over the entire temperature range are reconstructed using 8-segment high-order polynomials (up to the 7th order). The reason for adopting the 8-segment 7th-order polynomial formulation is twofold. First, in simulations of supercritical hydrocarbon fuel cracking, the thermophysical properties of the multi-component fluid must be evaluated repeatedly at the control-volume level during the iterative solution process. Compared with conventional table-lookup approaches, the polynomial formulation provides higher computational efficiency and is therefore more suitable for strongly coupled large-scale simulations. Second, the 8-segment piecewise formulation can more accurately capture the complex thermophysical property variations of supercritical hydrocarbons, especially in the pseudo-critical region, thereby avoiding the error accumulation associated with single-function fitting and improving the solver’s capability to resolve steep near-wall property gradients. Specific model parameters and coefficients are detailed in Table 1.

By implanting it into the OpenFOAM thermophysical property library, the aforementioned high-order segmented polynomial model was successfully applied. Thereby, the solver can smoothly capture the drastic variations in physical properties in the pseudo-critical region, significantly improving its numerical stability and prediction accuracy in handling extreme operating conditions.

The maximum relative error of the model’s predicted dynamic viscosity is within 10% across the full temperature range, and the relative error of thermal conductivity is less than 1%, which fully meets the engineering simulation accuracy requirements. Meanwhile, the maximum relative pressure drop Δp/p of the system is only ~0.068% under all working conditions. Such a tiny pressure fluctuation is sufficiently small to justify the constant pressure assumption in the present simulation, as its impact on thermophysical properties is negligible compared to the model’s inherent error margin.

2.3 Numerical Implementation of Simplified Cracking Reaction Kinetics

The n-decane cracking kinetic model proposed by Zhang et al. [25] was adopted as the benchmark mechanism in this study. To reduce the computational cost associated with three-dimensional conjugate heat transfer simulations, the original 23-step mechanism was simplified based on the molar fractions of the cracking products. The rationale for this simplification is that, in simulations of supercritical hydrocarbon fuel cracking, each additional reaction step increases the number of source-term evaluations and species-coupling calculations, thereby significantly increasing the overall computational cost. Therefore, according to the molar fractions of the major cracking products, the key species and dominant reaction pathways were retained, and a skeletal mechanism containing 9 reaction steps was constructed (the specific kinetic parameters are listed in Table 2). All reaction rates were calculated according to the Arrhenius law [20,23]. Although this simplified mechanism was adapted from the existing literature, it exhibits good compatibility with the present numerical framework. It provides reliable predictions of the overall cracking trend and species evolution with substantially improved computational efficiency, making it suitable for investigating the coupled transport phenomena in the present three-dimensional ribbed channel. The feasibility of the simplified mechanism is further validated in Section 3.3.3.

To account for the effect of fuel pyrolysis on both species evolution and thermal behavior, the chemical kinetics is incorporated into the present solver through a finite-rate Arrhenius formulation. The local reaction rates are determined from the local temperature, species concentrations, and stoichiometric relations, and are then used to construct both the species source terms and the chemical heat source term.

For the j-th reaction, the reaction rate constant is evaluated as: kj(T)=Ajexp(−EajRT)(9) where A j is the pre-exponential factor. Eaj is the activation energy, J/mol. R is the universal gas constant, J/(mol·K). T and T is the local fluid temperature, K. The parameter denoted by k in Table 2 corresponds to A j , while the actual reaction rate constant k j ( T ) is updated during the iterative solution. The activation energies listed in Table 2 were converted from kcal/mol to J/mol before being used in Eq. (9).

The molar concentration of species i is calculated as: Ci=ρYiWi(10) where C i is the molar concentration of species i, mol/m3. ρ is the fluid density kg/m3, Y i is the mass fraction of species. W i is the molecular weight of species, kg/mol. The local progress rate of reaction j is written as: rj=kj(T)∏m=1NsCmνm,j′(11) where r j is the reaction progress rate, mol / ( m 3 · s ) . Ns is the total number of species, and ν m , j ′ is the stoichiometric coefficient of reactant species m in reaction j. For unimolecular cracking reactions, Eq. (11) reduces to a first-order form.

The net production rate of species i is then given by: ω˙i=Wi∑i=1Nrνijrj(12) where ω ˙ i is the net production rate of species i, kg / ( m 3 · s ) . Nr is the total number of reactions, and ν i j is the net stoichiometric coefficient of species i in reaction j.

The volumetric chemical heat source term coupled into the energy equation is evaluated as: Qdot=−∑i=1Nshiω˙i(13) where Q d o t is the volumetric chemical heat source term, W/m3. h i is the specific enthalpy of species I, J/kg. Since n-decane pyrolysis is endothermic, Q d o t < 0 in the high-temperature reacting region.

2.4 Program Implementation of Multi-Physics Strong Coupling Algorithm

Based on the OpenFOAM framework, this study developed an iterative solution algorithm for fluid-thermal-chemical multi-field coupling. The solution process is shown in Fig. 1.

The core logic of the algorithm is as follows:

• (1)Fluid domain solution: The SIMPLE algorithm [26,27,28] is adopted to handle pressure-velocity coupling, solving the momentum equation and pressure correction equation.
• (2)Chemical reaction update: The Arrhenius rate constants, reaction progress rates, species source terms, and the chemical heat source term are updated and coupled into the species and energy equations.
• (3)Conjugate heat transfer coupling: This is the key to the algorithm development. When solving the energy equations for the fluid and solid domains, coupled boundary conditions are applied at the fluid-solid interface to strictly ensure temperature continuity (Tf = Ts) and heat flux continuity (qf = qs) at the interface. A second-order linear conservative interpolation scheme is employed for physical quantity transfer at the fluid-solid interface, which is consistent with the discretization order of the governing equations. The maximum relative temperature difference between the fluid and solid sides at the conjugate interface is merely 0.027% under all operating conditions, verifying that the interpolation error is fully negligible and imposes no effect on the solution accuracy.

Figure 1: Flowchart of the numerical solver algorithm.

Through this alternating iterative mechanism, the solver developed in this paper can achieve a fully coupled solution of the flow field, temperature field, and chemical reaction field until the residuals of all physical quantities converge. For the spatial discretization of the governing equations, the gradient terms are discretized using the Gauss linear scheme, whereas the convective terms are discretized using a bounded second-order TVD-type scheme, which ensures boundedness while retaining second-order spatial accuracy. The diffusion terms are discretized using the Gauss linear orthogonal scheme.

All numerical simulations in this study employ a steady-state solution strategy. A systematic sensitivity analysis of the iterative step size for the implicit solver was performed, with no obvious deviations in core physical quantities (wall temperature, fuel cracking conversion rate) under different step sizes, verifying the rationality of the selected step size. Steady-state convergence is defined by a dual coupling criterion, requiring residuals of all variables in the mass, momentum, energy and species transport equations to drop below 1 × 10−6 and key monitored parameters (outlet fluid temperature, n-decane mass fraction) to remain stable without obvious fluctuations. Strict adherence to this criterion avoids the limitations of single residual judgment and ensures the reliability and accuracy of coupled solution results for the flow, temperature and chemical reaction fields.

2.5 Construction of Computational Domain Physical Model

2.5.1 Description of Topological Structure and Definition of Computational Domain

In this study, a single dual-channel symmetric unit from a four-channel structure was selected as the computational domain to comparatively analyze Plate A and Plate B. As shown in Fig. 2, the computational domain consists sequentially of a 100 mm inlet adiabatic section, a 415 mm central heating section, and a 100 mm outlet adiabatic section along the flow direction. The boundary conditions were set as follows: the mid-plane is a symmetry boundary, a constant heat flux is applied to the bottom heating surface, and the remaining walls are adiabatic and no-slip.

• (1)Baseline Plate A: The flow passage consists of two 2 mm × 2 mm rectangular channels, with both the inter-channel rib and the outer wall having a thickness of 1 mm.
• (2)Improved Plate B: Periodic trapezoidal disturbance ribs were introduced based on Structure A, as shown in Fig. 2B. The ribs span across the channel (width 2 mm) and feature a longitudinally symmetric trapezoidal structure with a height of 0.5 mm. A total of 3 ribs are arranged within the heating section with a spacing of 100 mm, aiming to enhance heat transfer by inducing flow separation.

Figure 2: Schematic diagram of physical model and computational domain. (A) Baseline channel and (B) trapezoidal ribbed channel.

2.5.2 Boundary Condition Application and Numerical Case Settings

The inlet was configured with mass flow rate and fixed temperature boundary conditions. The outlet back pressure was 3.5 MPa, and the walls were no-slip. For the thermal boundary, constant heat flux heating was applied to the bottom single side, and the fluid-solid interface was set as a coupled boundary. To investigate the influence of flow rate, two sets of comparative cases were established (see Table 3).

3 Results and Discussion

3.1 Grid Independence Verification

To ensure sufficient resolution of near-wall flow and the grid independence of the computational results, this study conducted a systematic verification for both channel structures under the 2.0 g/s operating condition using structured hexahedral meshes (Fig. 3). The mesh generation was performed with the prerequisite of ensuring that the dimensionless wall distance of the first layer grid, y+, was less than 1. As shown in Fig. 4, a comparison of the calculation results from four sets of grids with different densities revealed that key physical quantities tend to stabilize as the grid was refined. Based on this convergence analysis and comprehensively considering the computational accuracy and efficiency, the final grid sizes were determined to be 2,003,645 cells for structure A and 2,703,345 cells for structure B. This grid resolution effectively balances the resolution of flow field details with overall computational efficiency, while guaranteeing that y+ ≤ 1.

Figure 3: Schematic diagram of the cooling channel computational mesh.

Figure 4: Results of grid independence verification.

3.2 Calculation Methods for Characteristic Parameters

To quantitatively evaluate the coupled thermal-hydraulic-chemical performance of the improved structure under supercritical pressure, a series of key characteristic parameters are defined in this study.

Cross-sectional average parameters. Given the non-uniform distribution of the physical properties of supercritical fluids across the cross-section, the bulk fluid temperature (Tf) and velocity (v) are both calculated using the mass-weighted average method: Tf=∫AcpρvTdA∫AcpρvdA(14) where A represents channel area, m2.

Heat transfer and flow characteristics. The local convective heat transfer coefficient (h) and local Nusselt number (Nu) are defined as follows: h=qwTw−Tf(15) Nu=hdhλ(16) where qw is heat flux, W/m2, Tw is wall temperature, K, and Tf is bulk fluid temperature, K, dh is hydraulic diameter of the channel, m and λ is thermal conductivity, W/(m2·K).

The flow resistance characteristics are characterized by the friction factor (f): f=2·Δp·dhρ·v¯2·L(17) where Δp is the streamwise pressure drop of the computational section (Pa), dh is the hydraulic diameter of the channel (m), ρ is the density of the fluid (kg/m3), v ¯ is the mass-weighted average velocity magnitude of the fluid (m/s), and L is the length of the computational channel (m).

Chemical reaction characteristics. To quantify the cracking reaction progress of n-decane along the flow direction, the local cracking conversion rate (ηx) is defined. First, the bulk mass fraction of n-decane ( Y C 10 H 22 ) on any flow channel cross-section is calculated using the mass-weighted average method as follows:

In fact, Y C 10 H 22 is the local mass fraction of n-decane. Based on this, the local cracking conversion rate along the flow path is defined as the proportion of the component consumption at the current cross-section relative to the inlet cross-section. Given that the inlet consists of pure n-decane working fluid (Yin = 1), ηx can be calculated using the formula:

This parameter reflects the cumulative degree of cracking of the fuel as it flows through the cross-section.

Comprehensive performance evaluation [22]. To comprehensively weigh the benefits brought by enhanced heat transfer against the cost of increased flow resistance, the performance evaluation criteria (PEC) factor is introduced: PEC=Nu/Nu0(f/f0)1/3(20) where Nu and f are the parameters of the improved structure (Plate B), and Nu0 and f0 are the corresponding parameters of the baseline straight structure (Plate A) under the same operating conditions. 1 PEC value more than 1 indicates that the comprehensive thermal-hydraulic performance of the improved structure is superior to that of the baseline structure.

3.3 Multi-Level Verification of Solver Reliability

3.3.1 Verification of Static Thermophysical Property Calculation Accuracy

To verify the reliability of the physical property model, the static physical properties of n-decane at 3.5 MPa were first calculated and compared. As shown in Fig. 5, the density, specific heat capacity, dynamic viscosity, and thermal conductivity predicted by the present model were in high agreement with the CoolProp benchmark data [29,30], with all relative errors within 10%, confirming that the model possessed a good calculation accuracy and engineering applicability.

Figure 5: Comparison of thermophysical parameters of n-decane at 3.5 MPa with CoolProp data. (a) density; (b) specific heat capacity at constant pressure; (c) dynamic viscosity; (d) thermal conductivity.

To quantify the calculation accuracy of the proposed 8-segment 7th-order polynomial transport property model, we calculated the relative errors between model predictions and CoolProp benchmark data for supercritical n-decane at 3.5 MPa over 273–973 K, with results shown in Fig. 6. The relative error of thermal conductivity is consistently below 1% across the full temperature range. For dynamic viscosity, the error stays within 3% far from the pseudo-critical region, peaks strictly below 10% near the pseudo-critical point (~635 K), and rapidly drops back to ~1% beyond this region. This model fully meets the accuracy requirements of engineering numerical simulation, and provides a reliable physical property basis for subsequent high-fidelity simulation of supercritical fluid multi-physics coupling processes.

Figure 6: Relative error distribution of the proposed 8-segment 7th-order polynomial model for transport properties of supercritical n-decane at 3.5 MPa.

3.3.2 Verification of Variable Property Flow and Heat Transfer

A square duct with a cross-section of 2 mm × 2 mm and a wall thickness of 0.5 mm was selected as the verification model. Along the flow direction, it was configured with 100 mm inlet/outlet adiabatic sections and a 600 mm central heating section. The computational operating conditions were set as follows: Inlet mass flow rate of 1 g/s (298 K), outlet back pressure of 3.5 MPa, and bottom outer wall heat flux of 1 MW/m2. The k-ω SST model is adopted for the solution. The exploited geometric model and boundary conditions of the verification case were exhibited in Fig. 7.

Figure 7: Geometric model and boundary conditions of the verification case (square duct).

The computational domain employed a structured mesh with near-wall refinement. Grid independence verification indicated that as the number of grid cells increased from 1 million to 3.54 million, the relative error of outlet parameters was controlled within 2%. Considering both computational accuracy and efficiency, 1,971,132 grid cells were finally selected.

Subsequently, the calculation results based on the CoolProp database were compared with those from the PR model above. As shown in Fig. 8 and Fig. 9, the two streamwise centerline temperature curves highly coincided, and the distribution characteristics of temperature contours at different cross-sections were essentially consistent. This fully confirmed that the constructed thermophysical property calculation model in this paper possessed an extremely high reliability in coupled flow and heat transfer simulations.

Figure 8: Comparison of streamwise centerline temperatures of the cooling channel under different physical property models.

Figure 9: Comparison of streamwise temperature field distributions predicted by different physical property models.

3.3.3 Comparative Verification by Experimental Reacting Flow

The circular-tube cracking experimental data reported by Lei et al. [10] were selected for validation. The computational domain was a circular tube with an inner diameter of 2 mm and an outer diameter of 3 mm, consisting of a 250 mm heated section and 150 mm adiabatic sections at both the inlet and outlet. The boundary conditions were specified as follows: an inlet mass flow rate of 0.8 g/s at 617.5 K, an outlet back pressure of 3.5 MPa, and a wall heat flux of 0.48 MW/m2 applied over the heated section. A structured mesh with near-wall refinement was employed. According to the grid-independence study, with the total cell number ranging from 1.5 million to 4.5 million, the error of the medium-density grid was already below 1%. Therefore, approximately 2.7 million grid cells were finally adopted. Based on this mesh, the fluid temperature and the mass fraction of unreacted n-decane at the cross-section of x = 400 mm were extracted for comparison, as listed in Table 4.

The simulated results show that the relative errors of the unreacted n-decane mass fraction and fluid temperature at the selected cross-section are 1.05% and 2.04%, respectively. In addition, the streamwise fluid-temperature distribution shown in Fig. 10 agrees well with the experimental trend, and the maximum relative deviation does not exceed 2.4%. These results indicate that the present solver can reasonably reproduce the cracking behavior of supercritical hydrocarbon fuels under the present conditions.

Figure 10: Comparison of simulated and experimental streamwise fluid temperatures.

To further examine the applicability of the developed multi-physics coupled solver under different operating conditions, an additional validation was carried out using the published experimental data of Zhao et al. [11]. In this case, the experimental conditions were reproduced using a straight tube with an inner diameter of 2 mm and a length of 70 cm. The inlet temperature, inlet mass flow rate, outlet back pressure, and wall heat flux were 773 K, 1 g/s, 3.0 MPa, and 0.22 MW/m2, respectively. Key physical quantities, including outlet temperature, gas yield, and heat transfer coefficient, were compared with the corresponding experimental measurements, as summarized in Table 5.

The simulated outlet temperature is 921.5 K, corresponding to a relative error of 1.2% compared with the experimental value of 933 K. The predicted gas yield is 14.2%, with a relative error of 9.2% relative to the experimental value of 13%. The simulated heat transfer coefficient is 3.73 kW/(m2·K), corresponding to a relative error of 5.09% compared with the experimental value of 3.93 kW/(m2·K). The maximum relative error among these key parameters remains within 10%, which is acceptable for the present engineering-oriented numerical simulations. This cross-validation under lower pressure, heat flux, and mass flow rate further supports the applicability of the present solver to the prediction of supercritical n-decane cracking and heat transfer under different operating conditions.

Overall, the validation results are not only qualitatively consistent in trend, but also quantitatively acceptable. For the first experimental case, the relative errors of the key variables are 1.05% and 2.04%, respectively, and the maximum relative deviation along the streamwise direction shown in Fig. 10 is 2.4%. For the second experimental case, the relative errors of outlet temperature, gas yield, and heat transfer coefficient are 1.2%, 9.2%, and 5.09%, respectively. These results indicate that the deviations between the numerical predictions and the experimental data remain within a reasonable range under the present conditions, thereby supporting the applicability of the present numerical model.

3.4 Analysis of Flow Resistance Characteristics and Flow Field Structure

Fig. 11 displays the streamwise pressure drop distributions of the two structures at mass flow rates of 2.5 g/s (Fig. 11a) and 2.0 g/s (Fig. 11b), respectively. From Fig. 11 it can be seen that under both 2.5 g/s and 2.0 g/s operating conditions, the pressure drop curves of the inner and outer channels highly overlapped, which indicated that no significant flow mal-distribution occurred in the system, demonstrating excellent flow distribution uniformity. And the streamwise pressure drop of structure B (trapezoidal fins) was obviously higher than that of structure A (baseline straight channel) under all operating conditions, with the pressure difference increase being most pronounced at the inlet section (x = 0.1 m), which indicating that while the trapezoidal fins enhanced heat transfer, they also introduced significant form drag, constituting a quantifiable penalty in flow resistance. In addition, comparing Fig. 11a and Fig. 11b, it was evident that when the mass flow rate decreased from 2.5 g/s to 2.0 g/s, the overall system pressure drops decreased substantially. This is mainly attributed to the reduction in flow velocity, which significantly weakens the wall friction resistance. From the perspective of hydraulic characteristics, this reveals the potential advantage of low flow conditions in reducing pump power consumption and improving system propulsion efficiency.

Figure 11: Comparison of streamwise pressure drop distributions in mass flow rate. (a) 2.5 g/s and (b) 2.0 g/s.

Fig. 12 illustrates the streamwise flow velocity distributions for the two structures at mass flow rates of 2.5 g/s (Fig. 12a) and 2.0 g/s (Fig. 12b). Overall, the working fluid exhibited significant non-linear acceleration characteristics in the heating section. This stemmed from the strong thermal-chemical coupling mechanism driven by the combined effects of heat absorption expansion and cracking reactions of the supercritical fluid. On the one hand, physically sensible heat absorption leads to a decrease in density and thermal expansion. On the other hand, the cracking of large molecules into smaller ones causes a significant increase in the number of moles (i.e., chemical expansion). The superposition of these two effects jointly drives a sharp reduction in density and a drastic volumetric expansion.

Through further analysis of the flow field non-uniformity and flow rate effects, it can be concluded that the velocity differences between channels originate from geometric topology and heat load distribution as shown in Fig. 12 and Table 6 because the flow velocities in the outer channels of both structures are significantly higher than those in the inner channels. The mechanism lies in the fact that the outer channel is subjected to single-sided heating (with the outer side being adiabatic) and must bear the entire thermal load from the outer wall. This results in its average fluid temperature (e.g., 734.7 K at 2.5 g/s) being significantly higher than that of the inner channel, which is cooled by fluid on both sides (716.9 K). The higher fluid temperature activates stronger Arrhenius cracking reactions, inducing more intense chemical expansion effects, thereby significantly increasing the flow velocity in the outer channel. At the same time, by comparing Fig. 12a,b it was observed that when the mass flow rate decreased from 2.5 g/s to 2.0 g/s, the outlet flow velocity did not decrease proportionally, suggesting a “residence time compensation” effect under low flow conditions. This is because the lower flow velocity extends the residence time of the fluid in the high-temperature zone, allowing the cracking reaction to proceed more sufficiently. The momentum gain generated by the reaction-induced chemical expansion partially compensates for the momentum loss caused by the reduction in mass flow rate, thereby maintaining a relatively high outlet flow velocity.

Figure 12: Streamwise fluid-velocity distributions at (a) 2.5 g/s and (b) 2.0 g/s.

3.5 Analysis of Convective Heat Transfer Enhancement Characteristics

Fig. 13 compares the outer wall temperature distributions under mass flow rates of 2.5 g/s (Fig. 13a,b) and 2.0 g/s (Fig. 13c,d). It was evident that the improved structure B consistently exhibited superior thermal protection performance. The wall temperature of the baseline structure A increased smoothly along the flow direction and reached a peak near the outlet (e.g., approximately 1226 K at 2.5 g/s). In contrast, due to the introduction of periodic trapezoidal ribs, the wall temperature distribution of structure B presented a unique “sawtooth” fluctuation characteristic. This fluctuation originates from the periodic flow separation and reattachment induced by the ribs, which effectively disrupts the continuous development of the near-wall thermal boundary layer and inhibits the formation and accumulation of local hot spots, thereby significantly reducing the peak wall temperature by 30 K to 40 K compared to structure A. When the mass flow rate was reduced to 2.0 g/s, although the prolonged residence time of the fluid within the channel led to an increase in the overall wall temperature level, the peak wall temperature of structure B remained consistently lower than that of structure A (Fig. 13c,d). This indicates that the flow disturbance and heat transfer enhancement mechanisms introduced by the trapezoidal fins function effectively under different flow rates, demonstrating that the structural design possesses good adaptability across a wide range of operating conditions.

Figure 13: Streamwise wall-temperature distributions of the upper and lower surfaces: (a) outer channel at 2.5 g/s, (b) inner channel at 2.5 g/s, (c) outer channel at 2.0 g/s, and (d) inner channel at 2.0 g/s.

Fig. 14 illustrates the distribution characteristics of the streamwise Nusselt number (Nu). Driven by the fluid expansion and acceleration effect induced by thermal cracking (Eq. (8)), the Nu numbers of both structures increase significantly along the flow direction, and the overall heat transfer intensity of the inner channel is higher than that of the outer channel. Unlike the smooth rising trend of structure A (with an inner channel Nu peak of approximately 73.9 at 2.5 g/s), structure B presents distinct periodic peak-valley oscillations, with its Nu peak increased to 80.4. Further combining this with topological structure analysis, it can be known that the heat transfer enhancement of structure B presented an oscillatory characteristic and followed a flow field mechanism.

Figure 14: Streamwise local Nusselt number distributions at (a) 2.5 g/s and (b) 2.0 g/s.

To visually reveal the flow evolution characteristics induced by the trapezoidal ribs and elucidate the underlying heat transfer enhancement mechanism, Fig. 15 presents the streamline distributions around the rib structure in the outer channel at a mass flow rate of 2.5 g/s. For the baseline straight channel (Structure A), the streamlines maintain a smooth and parallel distribution along the flow direction, with no obvious flow separation or secondary flow generated. Under this flow pattern, the thermal boundary layer develops continuously along the streamwise direction, leading to a gradual thickening of the boundary layer and a corresponding weakening of the local heat transfer capacity. In contrast, for the improved Structure B with trapezoidal ribs, the streamlines present a significant deflection and disturbance effect when flowing through the rib structure. The rib acts as a flow blockage, forcing the mainstream fluid to accelerate along the top surface of the rib, while a distinct flow separation region is formed in the wake area behind the rib. The recirculation vortex structure generated in the separation region continuously entrains the mainstream fluid, promoting the momentum and energy exchange between the near-wall low-temperature fluid and the core high-temperature fluid.

Figure 15: Comparison of streamline distributions around the rib structures in the outer channel at a mass flow rate of 2.5 g/s. ((Top): Improved Structure B with trapezoidal ribs; (Bottom): Baseline straight Structure A).

To further quantify the flow disturbance intensity and correlate it with the heat transfer enhancement effect, Fig. 16 presents the local turbulent kinetic energy (TKE) distribution of the outer channel at the end of the heating section (x = 0.409 m–0.416 m) at 2.5 g/s and 2.0 g/s. Here, the outer channel is selected as representative because it bears higher thermal load and flow velocity. The figure clearly shows that Structure B forms a significant high TKE region in the wake of the trapezoidal ribs, which is highly consistent with the flow separation and recirculation zone observed in the streamline distribution (Fig. 15). This confirms that the trapezoidal ribs successfully induce periodic flow separation and reattachment processes, and the resulting high-shear turbulence not only disrupts the continuous development of the near-wall thermal boundary layer but also significantly enhances radial fluid mixing. The periodic peaks of the Nu number in Fig. 14 are a direct manifestation of this intense turbulent mixing effect. In addition, although the rib structure introduces certain form drag, the heat transfer gain it brings is sufficient to offset any potential heat transfer attenuation caused by local flow structure adjustments, implying a balance between heat transfer and flow resistance. Even under the low flow condition of 2.0 g/s (Fig. 16b), Structure B still maintained significant Nu number enhancement and oscillation characteristics, indicating that this enhancement design possessed stable enhancement effects under different operating conditions.

Figure 16: Comparison of local turbulent kinetic energy (TKE) distributions near the trapezoidal fins at the end of the heating section (outer channel, x = 0.409 m~0.416 m). (a) 2.5 g/s; (b) 2.0 g/s.

3.6 Evaluation of Cracking Reaction Performance

Based on the streamwise average fluid temperatures (Fig. 17) and cracking conversion rate curves (Fig. 18), the reaction kinetic characteristics within the coupled thermal-fluid-chemical field can be further analyzed. Under both 2.5 g/s (Fig. 18a) and 2.0 g/s (Fig. 18b) operating conditions, the conversion rates exhibited strong non-linear dependence on temperature. When the fluid temperature Tf < 660 K, the reaction was in the induction period, and the conversion rate was close to zero. Once Tf exceeded ~ 700 K, the conversion rate increased steeply in a typical Arrhenius exponential manner, indicating that chemical heat absorption was mainly concentrated in the high-temperature region downstream of the heating section (Arrhenius-type temperature sensitivity). By comparing the two structures, it was found that the baseline structure A possessed a higher average fluid temperature due to weaker fluid mixing, which promoted more intense cracking reactions, thereby resulting in a slightly higher outlet conversion rate than structure B. However, although structure B suffered a minor loss in chemical heat sink, it gained significant benefits in thermal protection and coking inhibition. Therefore, structure B significantly reduced the peak wall temperature. This not only helps alleviate material thermal stress but, more critically, keeps the wall temperature away from the temperature zone where n-decane is prone to severe coking. Given that the coking rate increases exponentially with wall temperature, this design effectively suppresses the risk of carbon deposition within the channel by slightly reducing the conversion rate, thereby enhancing the long-term operational reliability of the cooling system. Obviously, the structural difference possessed a significant influence on reaction and coking characteristics.

Flow rate also has an important influence on the reaction progress. Comparing Fig. 18a and Fig. 18b it was seen that at a lower flow rate (2.0 g/s), the conversion rate curve shifted upward overall. This is because the lower flow velocity extends the residence time of the fluid in the high-temperature zone, allowing the cracking reaction to proceed more sufficiently, thereby compensating for the potential impact on reaction progress caused by the relative weakening of convective heat transfer.

Figure 17: Streamwise mainstream fluid temperature distributions. (a) 2.5 g/s; (b) 2.0 g/s.

Figure 18: Variation of cracking conversion rate with fluid temperature. (a) 2.5 g/s; (b) 2.0 g/s.

3.7 Comprehensive Thermal-Hydraulic Performance Evaluation

Based on the calculation for the comprehensive performance evaluation criterion (PEC) defined in Section 3.2, the comprehensive thermal-hydraulic performance of the improved cooling channel structure (Structure B) relative to the baseline structure (Structure A) was quantitatively evaluated. The specific calculation results under different operating conditions were presented in Table 7.

The results of the comprehensive performance evaluation criterion (PEC) calculated based on Eq. (20) are listed in Table 7. The data show that the improved structure B yields PEC values close to or slightly above 1 under all investigated operating conditions, indicating a modest overall thermal-hydraulic performance improvement. At the higher flow rate of 2.5 g/s, structure B enhances turbulent mixing and increases the Nu number by 7.32% (to 64.58), accompanied by a 23.71% increase in the friction factor. However, the additional resistance offsets most of the heat transfer gain, resulting in a PEC of approximately 1.00. At the lower flow rate of 2.0 g/s, the Nu number increases by 6.16%, while the friction factor rises by only 1.3%, leading to a PEC value of 1.057. Overall, the trapezoidal turbulator structure provides a limited but positive improvement in comprehensive performance under the present conditions, suggesting that it may serve as a feasible enhancement strategy for balancing heat transfer enhancement and hydraulic resistance.

4 Conclusions

Addressing the conjugate heat transfer problem of supercritical hydrocarbon fuels, this paper’s core innovation is the independent development of a high-fidelity multi-physics coupled solver based on OpenFOAM, integrating a high-precision piecewise variable property transport model with 8 segments and 7th-order polynomials, the PR equation of state, and a simplified cracking mechanism. It overcomes computational oscillations and accuracy loss caused by drastic pseudo-critical property variations. After full model validation, the thermal-hydraulic-chemical coupling characteristics of trapezoidal ribbed structure (B) and baseline structure (A) are compared, with key conclusions as follows:

Trapezoidal fins induce periodic flow separation and high-TKE regions, disrupting the thermal boundary layer and enhancing radial mixing, reducing peak wall temperature by 30 K–40 K to alleviate hot spots. This accurate quantification benefits from the model’s high-precision prediction of thermal conductivity (error < 1%) and dynamic viscosity (max error < 10%).

The trade-off between cracking and wall thermal protection is revealed. Structure B slightly reduces cracking conversion but significantly lowers wall temperature to suppress coking. The solver’s strong coupling ensures accurate capture of this synergy, enhancing long-term system safety while guaranteeing chemical heat sink utilization.

At low flow rates, Structure B introduces only limited additional resistance and shows a relatively favorable thermal-hydraulic response. The calculated PEC values are close to or slightly above 1 under all investigated flow rates, indicating a modest but positive overall performance improvement and suggesting a certain engineering reference value.

Through the development of the 8-segment 7th-order variable-property model and the multi-field strongly coupled solver, this work provides a high-fidelity numerical framework for studying the conjugate heat transfer of supercritical hydrocarbon fuels with drastic property variations and intense reactions in microchannels.

It should be noted that this study has two main sources of uncertainty, which are also common challenges in the field of supercritical hydrocarbon fuel cracking research. First, the lack of DNS analytical solution data for supercritical hydrocarbon fuels restricts the more precise calibration of the developed multi-physics coupled solver, which is a major limitation in the current model validation process. Second, consistent with the mainstream research in this field, the finite-rate model is adopted to control the cracking reaction rate in this study, and the influence of turbulence on the chemical reaction kinetics is not considered, which is the key uncertainty factor affecting the prediction accuracy of the cracking reaction process. Future research will focus on the above uncertainty factors, carry out comprehensive global uncertainty quantification and multi-parameter sensitivity analysis, combine DNS analytical solution data to optimize the calibration of the multi-physics coupled solver, and establish a turbulence-chemistry interaction model suitable for supercritical hydrocarbon fuel cracking, so as to further improve the prediction accuracy and engineering application value of the numerical model.