A Non-Intrusive DSMC-FEM Coupling Method for Two-Dimensional Conjugate Heat Transfer in Rarefied Hypersonic Conditions

Abstract

Accurate conjugate heat transfer (CHT) analysis is critical to the thermal management of hypersonic vehicles operating in rarefied environments, where non-equilibrium gas dynamics dominate. While numerous sophisticated CHT solvers exist for continuum flows, they are physically invalidated by rarefaction effects. This paper presents a novel partitioned coupling framework that bridges this methodological gap by utilizing the preCICE library to non-intrusively integrate the Direct Simulation Monte Carlo (DSMC) method, implemented in SPARTA, with the finite element method (FEM) via FEniCS for high-fidelity simulations of rarefied hypersonic CHT. The robustness and accuracy of this approach are validated through three test cases: a quasi-1D flat plate benchmark confirms the fundamental coupling mechanism against a reference finite difference solution; a 2D flat-nosed cylinder demonstrates the capability of the framework to handle highly non-uniform heat flux distributions and resolve the ensuing transient thermal response within the solid; finally, a standard cylinder case confirms the compatibility with curved geometries and its stability and accuracy in long-duration simulations. This work establishes a validated and accessible pathway for high-fidelity aerothermal analysis in rarefied gas dynamics, effectively decoupling the complexities of multi-physics implementation from the focus on fundamental physics.

1. Introduction

Hypersonic vehicles encounter significant thermal challenges during atmospheric reentry [1], continuous operation in ultra-low Earth orbit or in the vicinity of the Karman line [2], and orbital maneuvers [3], all of which involve rarefied flow environments. Accurately predicting the coupled temperature fields between the fluid and solid domains through conjugate heat transfer analysis is therefore critical for thermal protection system design and mission safety. The criticality of CHT is profoundly amplified under rarefied hypersonic conditions due to the complex, bidirectional interplay between the gas dynamics and the thermal response of the vehicle structure.

The transition to rarefied conditions introduces distinct physical phenomena not captured by continuum models. Near the vehicle surface, the Knudsen layer thickens as intermolecular collisions become less frequent than gas–wall collisions. This leads to significant velocity slip and temperature jump at the wall boundary, fundamentally altering momentum and energy exchange. Furthermore, rarefaction effects contribute to the thickening of shock waves and modifications in flow structures [4], which directly impact the heat flux distribution on the structure. These non-equilibrium effects demand a departure from traditional fluid dynamics assumptions.

For continuum hypersonic flows, significant efforts have been made to solve CHT problems by coupling computational fluid dynamics (CFD) solvers with solid heat conduction solvers. Ferrero and D’Ambrosio [5] made an attempt to tightly couple the Navier–Stokes fluid solver and heat conduction solver in the finite-volume framework for hypersonic CHT analysis. The study introduced a tight coupling technique suitable for unsteady time-accurate calculations and an alternative “quasi-stationary” approach for steady problems. The method’s reliability and accuracy were demonstrated by comparing numerical results with experimental data, and the impact of threshold settings on solutions was investigated. Again, in the finite volume framework, Peetala and Kulkarni et al. [6] used a loose coupling method to bridge the fluid and solid solvers to investigate the CHT of the shock wave–boundary layer interactions on a double wedge. It shows that CHT accurately predicts separation bubble length, heat flux, skin friction coefficient, and pressure distributions, which are critical for designing hypersonic vehicles. The use of different wall materials significantly affects the simulation results, with CHT providing more accurate predictions than fluid flow solvers alone. Buck [7] implemented tight coupling of finite-volume fluid and solid solvers in the Navier–Stokes Multi-Block framework to conduct CHT analysis in both two-dimensional (2D) and three-dimensional (3D) hypersonic flows. This method was validated through multiple test cases and applied to analyze the aerodynamic heating of vehicle components, demonstrating significant improvements in predicting wall temperatures and heat flux distributions compared to fluid-only simulations. The study also highlighted the importance of including conductive heat flux from species diffusion and convective boundary layer in CHT simulations. To enhance coupling efficiency and flexibility, more advanced strategies have been explored. Zhang and Chen et al. [8] presented a hybrid numerical framework that strategically combines the finite volume method for fluid and the finite element method for solid with adaptive time-stepping. Zhou and Lu et al. [9] presented a loosely coupled gas-kinetic BGK scheme for simulating CHT problems in hypersonic flows. They employed a modified Bhatnagar–Gross–Krook (BGK) scheme to solve the Navier–Stokes equations. The BGK scheme is extended to model heat conduction in solids by constructing an appropriate equilibrium distribution function, allowing both fluid and solid domains to be simulated using the same kinetic framework. Validation results of both methods show significant improvements in predicting long-duration CHT problems in terms of both efficiency and accuracy. Despite these methodological advancements, all such CFD-based approaches are predicated on the Navier–Stokes equations and are therefore fundamentally limited by the continuum assumption, rendering them incapable of accurately modeling the rarefaction physics described above.

CHT under rarefied conditions has long been studied using the Navier–Stokes–Fourier equations with Maxwell velocity slip [10] and Smoluchowski temperature jump boundary [11] conditions to model slip-flow regimes. Alkhalidi et al. [12] conducted a comprehensive finite-volume-based numerical investigation of buoyancy-driven flow in an enclosed cavity, revealing the Knudsen number as the dominant parameter causing a substantial decrease in Nusselt number with increased rarefaction. Similarly, Kabar et al. [13] applied a finite-volume approach to parallel-plate microchannels and highlighted the diminishing effect of axial conduction at high Knudsen numbers. Samian et al. [14] employed a thermal lattice Boltzmann method for micro/nanochannels and analyzed interfacial resistance variations with thermal conductivity ratios. Van Heumen [15] adopted a moment method with finite elements, comparing monatomic and polyatomic gases and underscoring the significance of collision models and internal degrees of freedom in conjugate thermal transport. These methods restrict the study of CHT under rarefied conditions in the near-continuum and slip-flow regimes. For CHT problems with moderate to high levels of rarefaction, numerical methods based on the Navier–Stokes–Fourier equations become inaccurate.

The Direct Simulation Monte Carlo method is the established high-fidelity tool for such rarefied flows and it has been extensively applied in studies of rarefied flow physics [4], hypersonic vehicles [16], and even astrophysics [17]. Despite its suitability, research on DSMC-based CHT remains scarce. A representative work can be found in the work of Moghadam et al. [18], where the DSMC method is used to analyze heat transfer characteristics in thermal cavities under rarefied flow conditions, but with a specific focus on constant heat flux boundary conditions. A notable effort by Apper and Kumar [19] developed a coupled method using the Stochastic Parallel Rarefied-Gas Time-Accurate Analyzer (SPARTA) codes and an in-house thermal code to investigate the thermal environment effects of hypersonic reentry, including the influence of rarefied flow, pyrolysis ablation effect, thermal radiation effect, and flow control. It is found that the coupled method can enhance the understanding of the thermal environment of hypersonic reentry vehicles.

However, the aforementioned coupled frameworks are typically monolithic or tightly integrated, necessitating either the development of bespoke code from the ground up or the creation of custom adapters deeply embedded within specific solvers. This approach introduces significant implementation hurdles, including complexities in managing disparate computational grids, reconciling different programming languages, and ensuring consistent data exchange. Consequently, researchers are often encumbered by these numerical challenges, diverting focus from the underlying physics of rarefied CHT and hindering broader scientific progress [20,21,22,23].

This inherent limitation underscores a critical research gap. There is a pressing engineering demand for a robust, flexible, and non-intrusive partitioned coupling scheme. Such a framework is crucial not only for seamlessly integrating the DSMC solver with advanced solid-conduction solvers but also for dramatically accelerating the design cycle to address critical thermal challenges in scenarios such as atmospheric reentry, sustained hypersonic flight, and orbital maneuvers. By abstracting away implementation complexities, it empowers researchers to swiftly construct and iterate upon simulation tools, thereby enabling a dedicated focus on solving these foundational engineering problems.

In recent years, the preCICE library [24] provides a framework, enabling non-intrusive coupling of existing solvers through a standardized API. The coupling approach using preCICE conceptualizes specialized solvers as modular components that interact via a standardized interface, thereby promoting flexibility, solver agnosticism, and the reuse of existing codes. A thriving ecosystem has emerged around preCICE, supported by official adapters for a wide range of solvers in continuum mechanics and structures (e.g., OpenFOAM [25], FEniCS [26], SU2 [27], and CalculiX [28]).

Therefore, to bridge this gap, a coupling method based on preCICE v3.3.0 [29] is used to connect SPARTA and FEM method to solve CHT problems in two-dimensional rarefied hypersionic flows. This partitioned strategy marks a departure from conventional monolithic approaches by delivering greater flexibility and solver independence, as it avoids any modification to the source codes of the individual solvers. Built upon the standardized preCICE infrastructure, the framework promotes reusability and minimizes development effort, allowing researchers to concentrate on the analysis of rarefied flow physics and thermal response rather than implementation challenges. It is worth noting that while FEniCS-0.9.0 is used here for demonstration, the methodology is solver-agnostic and compatible with any preferred FEM software. We demonstrate the coupling method by simulating CHT phenomena for a flat plate, a flat-nosed cylinder, and a cylinder, coupling SPARTA with FEM-based heat conduction solvers.

2. Numerical Model Description

2.1. The DSMC Method for Gas Flows

The DSMC [30] method is the predominant numerical approach for solving rarefied flow problems in engineering and scientific applications. Originally pioneered by Bird [31] in the 1970s, this technique avoids the direct computation of the Boltzmann equation by adopting a stochastic simulation strategy. DSMC typically involves discretizing the computational domain into grid cells whose sizes are smaller than the local mean free path. As a particle-based methodology, DSMC uses computational particles that each model a large ensemble of actual gas molecules or atoms. The gas dynamics are captured through a decoupled treatment of particle movement and collisions. Collisions are statistically sampled within each cell, enabling the reconstruction of flow characteristics through probabilistic averaging.

The DSMC method is employed in this work to model rarefied gas flows, implemented through the open-source SPARTA solver [32,33]. The vibrational energy mode of the gas is not included in the test cases, as it is not the primary focus of this work. For comprehensive discussions of both the fundamental DSMC theory and specific implementation details of the SPARTA solver, readers are referred to the foundational works by Bird [30] and the SPARTA technical manual and literature on the validation and verification tests [32,34].

2.2. The Finite Element Method for Heat Conduction

Internal heat conduction within a structure, caused by aerothermal loads, is simulated using the FEM. The finite element method is a numerical technique for solving problems governed by partial differential equations or formulated as functional minimization [35]. It discretizes a computational domain into smaller, simpler subdomains called finite elements (e.g., triangles, quadrilaterals in 2D; tetrahedra, hexahedra in 3D). Within each element, the solution is approximated using shape functions. By enforcing continuity and boundary conditions, a global system of algebraic equations is assembled and solved numerically [36].

The PDE of the transient heat conduction with constant thermal properties on an element based on Fourier’s law [37] is

(1)$\rho c_p{\displaystyle \frac{\partial T}{\partial t}}=k{\nabla }^{2}T+\dot{S}$

where T is the temperature, k is the thermal conductivity, $\rho$ is the density, $c_p$ is the specific heat capacity [38], and $\dot{S}$ is the term of the heat source and is zero in this work.

The DSMC solver will send the heat flux distribution at the interface to the FEM solver, which will use this heat flux distribution as a von Neumann boundary condition. The von Neumann boundary condition [35] of the interface is expressed as

(2)$-k{\displaystyle \frac{\partial T}{\partial n}}=q$

where q is the surface heat flux derived from the DSMC solver. The weak form approach [36] used in this study to solve the heat conduction equation offers an integral formulation of the governing differential equations. This formulation is developed by first multiplying the differential equations by an arbitrary function and then integrating the resulting product over the spatial domain of the problem. Integration by parts is applied to reduce the order of the derivatives, after which boundary conditions are included where applicable [36]. The weak form of Equation (1) at the interface can be acquired by integrating both sides of the equation on the domain $\Omega$

(3)${\int }_{\Omega }\rho c_p{\displaystyle \frac{\partial T}{\partial t}}vdx={\int }_{\Omega }k{\nabla }^{2}Tvdx$

where v is the test function. According to Green’s function, the right-hand side is

(4)${\int }_{\Omega }k{\nabla }^{2}T=-{\int }_{\Omega }k\nabla T·\nabla vdx+{\int }_{\partial \Omega }k{\displaystyle \frac{\partial T}{\partial n}}vds$

At the interface, substituting Equation (2) into Equation (4), the weak form becomes

(5)${\int }_{\Omega }{\displaystyle \frac{\partial T}{\partial t}}vdx+{\int }_{\Omega }\alpha \nabla T·\nabla vdx={\int }_{\partial \Omega }{\displaystyle \frac{q}{\rho c_p}}vds$

where $\alpha$ is $k/\rho c_p$, which is known as the thermal diffusivity, and $\partial \Omega$ is the boundary surface. The solver uses a backward Euler scheme to discretize time. Hence, the final weak form is

(6)${\int }_{\Omega }{\displaystyle \frac{T^{n+1}}{\Delta t}}vdx+{\int }_{\Omega }\alpha \nabla T^{n+1}·\nabla vdx={\int }_{\Omega }{\displaystyle \frac{T^n}{\Delta t}}vdx+{\int }_{\partial \Omega }{\displaystyle \frac{q}{\rho c_p}}vds$

where the superscript n is the time and $\Delta t$ is the time-step. It is important to note that the time-step for the FEM must satisfy the Fourier number condition, $Fo$, to achieve numerical stability, which is [38]

(7)$Fo={\displaystyle \frac{\alpha \Delta t}{{(\Delta x)}^2}}\le {\displaystyle \frac{1}{2}}$

All heat conduction parts in the test simulations of the following sections are solved using the FEniCS finite element solver [39].

2.3. Coupling Method

In traditional CFD-FEM solutions for CHT, CFD passes the surface temperature to FEM, while FEM passes the heat flux to CFD. In SPARTA, the gas temperature is modified through collisions between particles and walls. The wall statistically calculates the energy transferred from particles to the wall. Therefore, SPARTA passes the surface heat flux to FEM as a von Neumann boundary condition, and FEM solves the heat conduction PDE based on the new boundary condition.

The flow chart of the CHT simulation is presented in Figure 1. A simple SPARTA adapter is built according to the tutorial of the preCICE documentation [40] to realize the communication between SPARTA and preCICE. Before starting any simulations, the SPARTA adapter initializes the mesh, the internal interpolator, the SPARTA solver, and the data containers. The surface temperature obtained by the FEM solver and mapped by the preCICE module is read by the adapter at each time window. The SPARTA solver will use the modified surface temperature for particle–surface collisions. The surface total heat flow distribution will then be computed and extracted to the preCICE module for updating the von Neumann boundary condition in the FEM solver after the SPARTA solver runs the internal “run” command using the most recent surface temperature. The adapter will then move on to the following time window.

The serial explicit coupling scheme is used to integrate the SPARTA and FEniCS solvers. In this framework, the preCICE coupling time step is set to match the convergent time step of the FEM used to solve the solid heat conduction equation. This setup allows the DSMC solver to statistically average macroscopic quantities while treating the solid surface as a constant temperature boundary. Meanwhile, the FEM solver updates the solid temperature based on the averaged heat flux received from SPARTA.

Data transfer at the interface is managed using a nearest-neighbor mapping method with a “consistent” constraint, which ensures accurate projection of heat flux from the SPARTA surface grids to the solid mesh and vice versa for temperature.

The numerical stability of this explicit scheme is fundamentally based on the significant difference in physical time scales between the fluid and solid domains. The characteristic time for solid heat conduction is several orders of magnitude longer than the time scale of molecular motion resolved by DSMC. By aligning the coupling time window with the thermal response time of the solid, rapid statistical fluctuations inherent in the DSMC solution are naturally averaged over an appropriate time interval for the slower thermal diffusion in the solid. This physically consistent approach ensures that the coupled system evolves stably toward a steady state. The thermal inertia of the solid effectively dampens high-frequency numerical perturbations, guaranteeing robust and convergent simulations.

3. Validation Tests

Three validation tests are considered in this work, all of which are simulated under hypersonic conditions. The first case is a quasi-one-dimensional conjugate heat transfer problem. Due to the rarity of conjugate heat transfer simulations under rarefied hypersonic conditions, the decoupling verification strategy is adopted to progressively verify the entire part of conjugate heat transfer in the first case. The test case will then be expanded to include a two-dimensional axisymmetric problem. Finally, a verification will be conducted for a two-dimensional problem involving a curved surface geometry.

3.1. Rarefied Hypersonic Flow Around a Finite Rectangular Flat Plate

The first test case is a hypersonic rarefied $N_2$ gas flow around a 0.1 m rectangular flat plate, which derives from the benchmarking case of SPARTA [32]. The computational domain is shown in Figure 2. The thickness of the flat plate is 0.005 m. The flow velocity is (1504, 0, 0) m/s. The inlet gas number density and temperature are 3.716 × ${10}^{20}$ m⁻³ and 13.32 K, respectively, corresponding to a Knudsen number of 0.016. In SPARTA, the computational domain is discretized into 250 × 250 cells to guarantee the cell size is smaller than the local mean free path. The time-step is 3.102 $\times {10}^{-7}$ s.

For pure DSMC simulation, the surface temperature is 290 K with full accommodation. The simulation is performed for 100,000 steps after reaching the steady state to reduce the statistical errors, as shown in Figure 3a. When the number of DSMC steps reaches 50,000, the statistical fluctuation is significantly reduced compared with results using 10,000 and 20,000 samples, and a small discrepancy is observed between the results using 50,000 and 100,000 samples. The field of the translational temperature is presented in Figure 4. The comparison of the pressure distributions of the upper surface between the results in Refs. [41,42] and the current simulation is presented in Figure 3b. Excellent agreement can be seen between the two DSMC simulations, proving the accuracy and effectiveness of the DSMC solver.

For the CHT test, the structure is decomposed into 500 elements in the x direction with only 1 element in the y direction. The heat conduction of the flat plate becomes a one-dimensional problem. The left surface of the flat plate is the interface between SPARTA and FEniCS. The heat flux is sampled and written by the SPARTA–preCICE adapter. Then, it is read by the FEniCS–preCICE adapter and passed to FEniCS. To validate the accuracy of the data transmission, the heat flux on the interface written by the SPARTA–preCICE adapter is set to 100,000 W/m² rather than passing the calculated heat flux from SPARTA, because the heat flux from SPARTA is too small to generate a noticeable temperature increase. The preCICE time window is 1 × 10⁻⁴ s and the coupling simulation ends at 2 ms.

The material of the flat plate is copper whose thermal conductivity, density, and specific heat capacity are 400 W/(m · K), 8960 kg/m³, and 385 J/(kg · K) [43]. To demonstrate the precision of the CHT calculation, the result of the temperature distribution inside the copper flat plate is compared to that of a finite difference method (FDM) [44]. The Python code for the FDM used in this test can be found in Ref. [45]. The space and time are discretized using the centered differencing scheme and implicit backward Euler scheme, respectively. The equations are solved with Newton’s method.

The result from the CHT test is compared with that from the FDM, as shown in Figure 5. Because of insufficient time for the constant heat flux to penetrate into the flat plate, the temperature increase is limited to the region close to the stagnation point. However, excellent agreement between the results from the coupling method and the FDM can be seen and the difference between the two results is 0.0005%, proving the correctness of the mathematical model in FEniCS.

3.2. Hypersonic Flow over a Flat-Nosed Cylinder

The second test case is a 2D hypersonic flow over a flat-nosed cylinder, which derives from the validation test of the DSMC solver [30,46]. Bird [30] conducted a 2D axisymmetric problem while Scanlon and Roohi et al. simulated a 3D case. Due to limited computational resources, this paper considers a 2D axisymmetric condition.

According to our tests, it is found that the surface heat flux generated by using the velocity in the literature [30,46] is small, and it takes a very long time to observe obvious temperature changes inside the structure under the condition of using real materials. The purpose of the test is to verify the correct mapping of the non-uniform heat flux and the correct conduction process. Therefore, the flow Mach number at the inlet is increased to 10.74 to obtain a higher surface heat flux, and a virtual material with ultra-high thermal conductivity, low density, and low specific heat capacity is used.

The computational domain of the DSMC simulation is discretized into 200 × 150 cells. The free-stream gas flow is assumed to be argon flow to cancel the influence of chemical reactions. The number density, temperature, and velocity at the inlet boundary are 1 × ${10}^{21}$/m³, 100 K, and 2000 m/s, respectively. A diffusive cylinder wall with full accommodation is considered and the time-step is 2 × 10⁻⁸ s.

According to the mesh independence test for the flat-nosed cylinder with a constant heat flux of 100,000 W/m² at the left and upper boundaries, as seen in Figure 6, the negligible discrepancy can be found when the spatial resolution reaches 100 × 50, so this mesh level is used in the formal simulation. The material density, thermal conductivity, and specific heat capacity are 1000 kg/m³, 10,000 W/(m·K), and 200 J/(kg·K), respectively. The time-step in FEniCS is 2 × 10⁻⁴ s.

To avoid transient statistic noise of the surface heat flux distribution before the simulation reaches its steady state, the CHT will not start until the DSMC simulation reaches a steady state. The flow temperature field in steady state is presented in Figure 7. High-speed airflow stagnates on the plane, forming a boundary layer, and a strong shock wave is generated at the leading edge of the boundary layer. This will generate a massive heat flux on the solid surface, as shown in Figure 8. The main goal of this test is to validate that the structure heat conduction is physical under the non-uniform surface heat flux distribution. It can be found that the heat flux at x = 0 m is significantly higher than that of y = 0.01 m due to the strong main shock wave. In Figure 8, an increase in the heat flux can be found between 0 and 0.01 m (flat nose), and a slight decrease in the waist of the cylinder (y = 0.01 m) can be seen. This is caused by an oblique shock resulting from the further compression after the primary shock wave along the waist of the cylinder, which is demonstrated in the pseudo-schlieren image in Figure 9.

Therefore, during the CHT, the temperature in the upper left corner of the structure will be the highest because of the maximum heat flux. The results of the evolution of the temperature distribution of the structure are shown in Figure 10. Due to extremely high thermal conductivity, the heat passes rapidly in the cylinder. The internal temperature in the upper left corner raises 0.6 K in 0.01 s and the temperature of the elements close to the flat nose is higher than the other places because of higher surface heat flux, leading to an expected non-uniform temperature distribution in the cylinder.

3.3. 2D Hypersonic Flow Past a Cylinder

The last test case is a Mach 10 argon flow impacting a copper cylinder with a diameter of 0.3048 m. The thermal properties of the copper are the same as those in Section 3.1. This case has been widely used in the benchmarking of DSMC-related codes [46,47,48,49,50,51]. Again, the pressure $p_{\infty }$, temperature $T_{\infty }$ and velocity at the inlet are 1.173 Pa, 200 K, and (2634, 0, 0) m/s, corresponding to a gas number density of 4.237 × 10²⁰/m³ and a Knudsen number of 0.01 based on the cylinder diameter. The variable hard sphere collision model is used in the simulation. The temperature exponent, $\omega$, reference temperature, $T_{ref}$, and the atom diameter are 0.734, 1000 K, and 3.595 × 10⁻¹⁰ m according to Ref. [50]. The particle–surface interaction is assumed to be diffusive with full accommodation. The initial temperature of the structure is 500 K and no temperature gradient exists in the structure at the beginning. The computational domain is shown in Figure 11 and is discretized into 600 × 800 cells. The time-step used in the DSMC solver is 1 × 10⁻⁷ s, which is smaller than the mean collision time (9.6 × 10⁻⁶). The simulation using the coupled method is ended at 0.744 s here because of the expensive cost of the DSMC method.

A mesh independence test is conducted with a constant heat flux of 40,000 W/m² around the cylinder, see Figure 12. It is found that when the spatial resolution reaches 25,288 cells, the numerical difference is too small to be noticed, so this level of resolution is considered in the simulation. The material of the cylinder is copper again and the properties are the same as that mentioned in Section 3.1.

The fields of the velocity component in the x direction and the thermal temperature are demonstrated in Figure 13. As expected, the hypersonic gas flow stagnates at the point (−0.1524, 0) m and the temperature increases significantly, leading to the highest heat flux in the domain. As the gas flows over the cylinder, the heat flux gradually decreases. The distribution of heat flux on the upper surface is presented in Figure 14, which shows good agreement with the results in Refs. [46,48].

At the beginning, the high-temperature region is located around the stagnation point and the temperature at the back of the cylinder is relatively low, as shown in Figure 15. To prove the ability of the long-duration thermal analysis using the non-intrusive framework, a conventional ’quasi-steady’ method is used by substituting the surface heat flux distribution at the gas steady state into the FEniCS solver and a comparison is made between the ’quasi-steady’ result and that from the coupled method. An excellent agreement can be found in Figure 16a between the results from the two methods. The heat from the hypersonic flow gradually penetrates about 41.4 mm (0.1524 m minus 0.11 m) of the cylinder at 0.744 s and the maximum temperature increase is approximately 1.12 K. It can be found that the shape of the surface temperature distribution is consistent with that of the heat flux. As time passes, heat gradually penetrates the interior of the cylinder, causing a significant increase in the cylinder temperature, as seen in Figure 16b and Figure 17. The maximum temperature in the copper cylinder at 10 min is about 538.64 K. The comparison between the short-term temperature rise predicted by our coupled scheme and the initial phase of the long-term quasi-steady solution in Figure 16 demonstrates the accuracy of the convective heat transfer prediction. This strong agreement strongly suggests the ability of the framework to reliably project the thermal response of the solid over extended durations.

4. Conclusions and Future Works

This study has provided a flexible, partitioned numerical framework for conjugate heat transfer analysis in rarefied gas flows. The robustness of this framework was systematically demonstrated through a series of validation cases. The quasi-1D flat plate simulation verified the fundamental coupling mechanism, showing excellent agreement with a reference finite difference solution. The axisymmetric flat-nosed cylinder case proved the ability of the framework to accurately handle complex, non-uniform heating patterns and manage the transient thermal response. Finally, the standard cylinder case validated the compatibility of the adapter with curved geometries and its stability for long-duration coupling performance using realistic material properties.

In successfully demonstrating this capability, this work provides the rarefied flow community with a validated tool for aerothermal analysis. It significantly lowers the barrier to entry for performing complex CHT simulations by allowing researchers to leverage existing, powerful solvers without invasive code modifications. This method is directly applicable to the high-fidelity conjugate heat transfer analysis of ultra-low orbit vehicles and micro-electromechanical systems. Future work will focus on extending this framework to three-dimensional problems and can incorporate additional physics, such as surface ablation and radiation, to create a comprehensive tool for hypersonic vehicle design and analysis.

Author Contributions

Conceptualization, Z.C.; methodology, Z.C.; software, Z.C.; validation, Z.C.; formal analysis, Z.C.; investigation, Z.C.; resources, Z.C.; data curation, Z.C.; writing—original draft preparation, Z.C.; writing—review and editing, C.M.; visualization, Z.C. funding acquisition, Z.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data can be available if requred.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

$\rho$	density
u	velocity conponet in the x direction
T	temperature
q	heat flux
k	thermal conductivity
$c_p$	specific heat capacity
$\dot{S}$	term of heat source
x	spatial cartesian coordinate
t	time
$Fo$	Fourier number
Ma	Mach number
Kn	Knudsen number
∞	free stream properties
CHT	Conjugate Heat Transfer
FEM	Finite Element Method
FDM	Finite Difference Method
PDE	Partial Differential Equation
DSMC	Direct Simulation Monte Carlo

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 Flow chart for conjugate heat transfer using the partitioned coupling method. Y: yes; N: no.

Figure 2 The computational domain of a hypersonic flow around a rectangular flat plate.

Figure 3 (a) Comparison of the pressure distribution of the upper surface with different numbers of samples in DSMC. (b) Comparison of pressure distributions of the upper surface of the flat plate [41,42].

Figure 4 Field of temperature at the steady state around the flat plate using SPARTA.

Figure 5 Comparison of the internal temperature distribution of the flat plate at 2 ms.

Figure 6 Temperature distribution along the line between point (0, 0.01) and point (0.02, 0) with different spatial resolutions.

Figure 7 Gas temperature field at the steady state.

Figure 8 Surface heat flux distribution along the surface path s from (0, 0) to (0.02, 0.01) in clockwise direction. For 0 < s < 0.01 m, this part belongs to line x = 0 m. The part where s > 0.01 m belongs to line y = 0.01 m.

Figure 9 Pseudo-schlieren at the steady state.

Figure 10 Temperature variation of the flat-nosed cylinder with time.

Figure 11 Computational domain of a hypersonic flow past a cylinder. $p_{\infty }$ is the freestream pressure, $u_{\infty }$ is the freestream velocity in the x direction, and $T_{\infty }$ is the freestream temperature.

Figure 12 Mesh independence test of the cylinder under different spatial resolution.

Figure 13 Field of the velocity component in the x direction (a), field of the Mach number (b), and the field of translational temperature (c) at the steady state around the cylinder using SPARTA.

Figure 14 Heat flux distribution of the upper spherical surface at the steady state [46,48].

Figure 15 Field of the temperature of the cylinder at 0.744 s using SPARTA–FEniCS coupled method. The temperature range is modified for the purpose of visualization.

Figure 16 The temperature distribution along the x axis (a) in the short duration; (b) in the long duration using SPARTA-FEniCS and the quasi-steady method.

Figure 17 The variation of cylinder temperature in long duration. The maximum temperature is limited to 520 K for visualization.

References

1. Girija, A.P. Comparative Study of Planetary Atmospheres and Implications for Atmospheric Entry Missions. arXiv; 2023; [DOI: https://dx.doi.org/10.48550/arXiv.2307.16277] arXiv: 2307.16277

2. Lysenko, O.I.; Sparavalo, M.K.; Tachinina, O.M.; Yavisya, V.S.; Ponomarenko, S.O. Feasibility reasoning of creating ultra-low orbit communication systems based on small satellites and method of their orbits designing. Inf. Telecommun. Sci.; 2020; 11, pp. 59-70. [DOI: https://dx.doi.org/10.20535/2411-2976.12020.59-70]

3. Leomanni, M.; Bianchini, G.; Garulli, A.; Quartullo, R.; Scortecci, F. Optimal low-thrust orbit transfers made easy: A direct approach. J. Spacecr. Rockets; 2021; 58, pp. 1904-1914. [DOI: https://dx.doi.org/10.2514/1.A34949]

4. Cao, Z.; White, C.; Kontis, K. Numerical investigation of rarefied vortex loop formation due to shock wave diffraction with the use of rorticity. Phys. Fluids; 2021; 33, 067112. [DOI: https://dx.doi.org/10.1063/5.0054289]

5. Ferrero, P.; D’Ambrosio, D. A numerical method for conjugate heat transfer problems in hypersonic flows. Proceedings of the 40th Thermophysics Conference; Seattle, WA, USA, 23–26 June 2008; 4247.

6. Peetala, R.K.; Kulkarni, V.; Sahoo, N. Shock wave boundary layer interactions in hypersonic flows over a double wedge geometry by using conjugate heat transfer. Heat Transf.; 2021; 1, pp. 801-817. [DOI: https://dx.doi.org/10.1002/htj.21905]

7. Buck, A. Conjugate Heat Transfer Simulations for Hypersonic Flight. Ph.D. Thesis; Universität der Bundeswehr München: Neubiberg, Germany, 2024.

8. Zhang, S.; Chen, F.; Liu, H. Time-adaptive, loosely coupled strategy for conjugate heat transfer problems in hypersonic flows. J. Thermophys. Heat Transf.; 2014; 28, pp. 635-646. [DOI: https://dx.doi.org/10.2514/1.T4278]

9. Zhou, D.; Lu, Z.; Guo, T.; Chen, G. A loosely-coupled gas-kinetic BGK scheme for conjugate heat transfer in hypersonic flows. Int. J. Heat Mass Transf.; 2020; 147, 119016. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.119016]

10. Maxwell, J.C. On stresses in Rarefied Gases Arising from Inequalities of Temperature. Philos. Trans. R. Soc. Lond.; 1879; 170, pp. 231-256.

11. von Smoluchowski, M. Ueber wärmeleitung in verdünnten Gasen. Ann. Phys. Chem.; 1898; 64, pp. 101-130. [DOI: https://dx.doi.org/10.1002/andp.18983000110]

12. Alkhalidi, A.; Kiwan, S.; Al-Kouz, W.; Alshare, A. Conjugate heat transfer in rarefied gas in enclosed cavities. Vacuum; 2016; 130, pp. 137-145. [DOI: https://dx.doi.org/10.1016/j.vacuum.2016.05.013]

13. Kabar, Y.; Bessaïh, R.; Rebay, M. Conjugate heat transfer with rarefaction in parallel plates microchannel. Superlattices Microstruct.; 2013; 60, pp. 370-388. [DOI: https://dx.doi.org/10.1016/j.spmi.2013.05.016]

14. Samiana, R.; Abbassia, A.; Ghazanfarian, J. High-Knudsen Conjugate Heat Transfer in Micro/Nano-Channels for use in microelectronic devices: Thermal Lattice Boltzmann Study. Proceedings of the 6th International Conference on Nanostructures; Kish Island, Iran, 7–10 March 2016.

15. Van Heumen, J.M.W. Numerical Study on Conjugate Heat Transfer Between Rarefied Gases and a Solid Using the Method of Moments. Master’s Thesis; Eindhoven University of Technology: Eindhoven, The Netherlands, 2021.

16. Sohn, I.; Li, Z.; Levin, D.A.; Modest, M.F. Coupled DSMC-PMC radiation simulations of a hypersonic reentry. J. Thermophys. Heat Transf.; 2012; 25, pp. 22-35. [DOI: https://dx.doi.org/10.2514/1.T3633]

17. Weinberg, M.D. Direct simulation Monte Carlo for astrophysical flows–I. Motivation and methodology. Mon. Not. R. Astron. Soc.; 2014; 438, pp. 2995-3006. [DOI: https://dx.doi.org/10.1093/mnras/stt2406]

18. Moghadam, E.; Roohi, E.; Esfahani, J. Heat transfer and fluid characteristics of rarefied flow in thermal cavities. Vacuum; 2014; 109, pp. 333-340. [DOI: https://dx.doi.org/10.1016/j.vacuum.2014.06.009]

19. Appar, A.; Kumar, R.; Naspoori, S.K. Conjugate flow-thermal analysis of a hypersonic reentry vehicle in the rarefied flow regime. Phys. Fluids; 2022; 34, 026107. [DOI: https://dx.doi.org/10.1063/5.0082783]

20. Idelsohn, S.R. Numerical Simulations of Coupled Problems in Engineering; Springer: Berlin/Heidelberg, Germany, 2014.

21. Uekermann, B.W. Partitioned Fluid-Structure Interaction on Massively Parallel Systems. Ph.D. Thesis; Technische Universität München: München, Germany, 2016.

22. Besseron, X.; Adhav, P.; Peters, B. Parallel multi-physics coupled simulation of a midrex blast furnace. Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region Workshops; Nagoya, Japan, 25–27 January 2024; pp. 87-98.

23. Adhav, P.; Besseron, X.; Peters, B. Development of 6-way CFD-DEM-FEM momentum coupling interface using partitioned coupling approach. Results Eng.; 2024; 22, 102214. [DOI: https://dx.doi.org/10.1016/j.rineng.2024.102214]

24. Chourdakis, G.; Davis, K.; Rodenberg, B.; Schulte, M.; Simonis, F.; Uekermann, B.; Abrams, G.; Bungartz, H.; Yau, L.C.; Desai, I. . preCICE v2: A sustainable and user-friendly coupling library. Open Res. Eur.; 2022; 2, 51. [DOI: https://dx.doi.org/10.12688/openreseurope.14445.1]

25. Chourdakis, G.; Schneider, D.; Uekermann, B. OpenFOAM-preCICE: Coupling OpenFOAM with external solvers for multi-physics simulations. OpenFOAM J.; 2023; 3, pp. 1-25. [DOI: https://dx.doi.org/10.51560/ofj.v3.88]

26. Rodenberg, B.; Desai, I.; Hertrich, R.; Jaust, A.; Uekermann, B. FEniCS–preCICE: Coupling FEniCS to other simulation software. SoftwareX; 2021; 16, 100807. [DOI: https://dx.doi.org/10.1016/j.softx.2021.100807]

27. Economon, T.D.; Palacios, F.; Copeland, S.R.; Lukaczyk, T.W.; Alonso, J.J. SU2: An open-source suite for multiphysics simulation and design. AIAA J.; 2016; 54, pp. 828-846. [DOI: https://dx.doi.org/10.2514/1.J053813]

28. Uekermann, B.; Bungartz, H.; Yau, L.C.; Chourdakis, G.; Rusch, A. Official preCICE adapters for standard open-source solvers. Proceedings of the 7th GACM Colloquium on Computational Mechanics for Young Scientists from Academia; Stuttgart, Germany, 11–13 October 2017.

29. preCICE Release. 2025; Available online: https://github.com/precice/precice/releases/tag/v3.3.0 (accessed on 16 November 2025).

30. Bird, G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows; Oxford University Press: Oxford, UK, 1994.

31. Bird, G.A. Monte Carlo simulation of gas flows. Annu. Rev. Fluid Mech.; 1978; 10, pp. 11-31. [DOI: https://dx.doi.org/10.1146/annurev.fl.10.010178.000303]

32. Plimpton, S.J.; Moore, S.G.; Borner, A.; Stagg, A.K.; Koehler, T.P.; Torczynski, J.R.; Gallis, M.A. Direct simulation Monte Carlo on petaflop supercomputers and beyond. Phys. Fluids; 2019; 31, 086101. [DOI: https://dx.doi.org/10.1063/1.5108534]

33. SPARTA Official Website. 2025; Available online: https://sparta.github.io/ (accessed on 16 November 2025).

34. Klothakis, A.; Nikolos, I.K. Comprehensive Evaluation of the Massively Parallel Direct Simulation Monte Carlo Kernel “Stochastic Parallel Rarefied-Gas Time-Accurate Analyzer” in Rarefied Hypersonic Flows—Part A: Fundamentals. Computation; 2024; 12, 198. [DOI: https://dx.doi.org/10.3390/computation12100198]

35. Nikishkov, G.P. Introduction to the Finite Element Method; University of Aizu: Aizuwakamatsu, Japan, 2004.

36. Zienkiewicz, O.C.; Taylor, R.L.; Zhu, J.Z. The Finite Element Method: Its Basis and Fundamentals; Elsevier: Amsterdam, The Netherlands, 2005.

37. Cengel, Y.A.; Ghajar, A.J. Heat and Mass Transfer (in SI Units); Mcgraw-Hill Education-Europe: London, UK, 2014.

38. Whitaker, S. Fundamental Principles of Heat Transfer; Elsevier: Amsterdam, The Netherlands, 2013.

39. Alnæs, M.; Blechta, J.; Hake, J.; Johansson, A.; Kehlet, B.; Logg, A.; Richardson, C.; Ring, J.; Rognes, M.E.; Wells, G.N. The FEniCS project version 1.5. Arch. Numer. Softw.; 2015; 3, 100.

40. The preCICE Documentation. 2025; Available online: https://precice.org/docs.html (accessed on 16 November 2025).

41. Klothakis, A.; Nikolos, I.K. Modeling of rarefied hypersonic flows using the massively parallel DSMC kernel SPARTA. Proceedings of the 8th International Congress on Computational Mechanics; Athens, Greece, 12–15 July 2015.

42. Allegre, J.; Raffin, M.; Chpoun, A.; Gottesdiener, L. Rarefied hypersonic flow over a flat plate with truncated leading edge. Prog. Astronaut. Aeronaut.; 1994; 160, 285.

43. Technical Data for Copper. Available online: https://periodictable.com/Elements/029/data.html (accessed on 16 November 2025).

44. Fu, R.; Weng, H.; Wenk, J.F.; Martin, A. Thermomechanical coupling for charring ablators. J. Thermophys. Heat Transf.; 2018; 32, pp. 369-379. [DOI: https://dx.doi.org/10.2514/1.T5194]

45. 1D Heat Conduction Solver. 2019; Available online: https://github.com/rickfu415/heatConduction (accessed on 16 November 2025).

46. Scanlon, T.J.; Roohi, E.; White, C.; Darbandi, M.; Reese, J.M. An open source, parallel DSMC code for rarefied gas flows in arbitrary geometries. Comput. Fluids; 2010; 39, pp. 2078-2089. [DOI: https://dx.doi.org/10.1016/j.compfluid.2010.07.014]

47. Lofthouse, A.J.; Boyd, I.D.; Wright, M.J. Effects of continuum breakdown on hypersonic aerothermodynamics. Phys. Fluids; 2007; 19, 027105. [DOI: https://dx.doi.org/10.1063/1.2710289]

48. Goshayeshi, B.; Roohi, E.; Stefanov, S. A novel simplified Bernoulli trials collision scheme in the direct simulation Monte Carlo with intelligence over particle distances. Phys. Fluids; 2015; 27, 107104. [DOI: https://dx.doi.org/10.1063/1.4933251]

49. Darbandi, M.; Roohi, E. A hybrid DSMC/Navier–Stokes frame to solve mixed rarefied/nonrarefied hypersonic flows over nano-plate and micro-cylinder. Int. J. Numer. Methods Fluids; 2013; 72, pp. 937-966. [DOI: https://dx.doi.org/10.1002/fld.3769]

50. Lofthouse, A.J. Nonequilibrium Hypersonic Aerothermodynamics Using the Direct Simulation Monte Carlo and Navier-Stokes Models. Ph.D. Thesis; University of Michigan: Ann Arbor, MI, USA, 2008.

51. Vasileiadis, N.; White, C. hybridDCFoam: A coupled DSMC/Navier–Stokes–Fourier solver for steady-state multiscale rarefied gas flows. Adv. Eng. Softw.; 2024; 193, 103669. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2024.103669]