A Numerical Method for Conjugate Heat Transfer Problems in Multicomponent Flows

The 3D numerical model of thermal interaction of a multicomponent flow and a solid body is presented. The approach is based on computational fluid dynamic simulation in which the heat transfer in a fluid domain and in a solid are fully coupled. This is known as the conjugate heat transfer problem. According to the our numerical scheme, the computation of a single time step is split into the sequence of hyperbolic and parabolic stages. The energy equation for solid and fluid domains is solved as a unified equation with the explicit-iteration scheme. At the fluid-solid interfaces the matching conditions are the continuity of temperature and the normal components of the heat flux. The continuity conditions are natural for the unified heat conduction equation and there is no need separate consideration of these conditions. The proposed method is generalized for multiblock frame: we implement the conjugate heat transfer algorithm in the case of interaction of fluid domains with various solid domains using multiblock unstructured conformal grids. Our technique can be used to develop a numerical conjugate heat transfer algorithm for calculating magnetohydrodynamic flows in channels.


Introduction
We develop the numerical model of thermal interaction of compressible viscous heat conducting multicomponent flow with a solid body (aircraft, walls of a combustion chamber of a propulsion system, etc). Heat transfer problem in computational fluid dynamic (CFD) simulation of supersonic vehicles involving considerable heating of an aerodynamic construction have become of particular significance. The present paper describes our approach to the conjugate heat transfer (CHT) problem. The algorithm and computer code are implemented for multicomponent flow, but multicomponent nature does not introduce principal specificity to the CHT problem, only it radically increases the computational complexity of the entire problem especially in 3D.
The solver MCFL (MultiComponent FLows) [1,2] may run as a standalone CFD solver for thermal multicomponent fluid problems. Here we propose an approach to extend functionality of MCFL by coupling it with a solid thermal solver in order to solve the complete conjugate heat transfer problem. We develop in-house solid thermal solver in the same framework as the solver MCFL (unstructured 3D grid, finite-volume discretization, parallel implementation). The heat conduction code is strongly coupled with the MCFL code. This coupling is based on the usage of a unified scheme for heat conduction in the fluid and in the solid. The basis of such coupling is provided by the scheme MCFL. According to the approach [1], computation of a single time step is split into the sequence of hyperbolic and parabolic stages. In the first (hyperbolic) substep, we consider the convective fluxes; in the second (parabolic) substep the all dissipative fluxes are considered. In our splitting approach the hyperbolic subproblem is solved with the modified Godunov type scheme [3,4].
On parabolic stage we solve a few diffusion subproblems which take into account processes of viscosity, multicomponent diffusion and heat conduction. In the CHT problem the heat conduction equations in the solid and in the fluid are solved with the unified algorithm based on the explicit iterative time integration scheme called LINS (Local Iterations for Navier-Stokes equations) [1,2]. Due to this scheme the coupling between the solid and the fluid is very tight at the solid-fluid interface. On each time step of the explicit scheme, or on each iteration of the LINS scheme,the fluid and solid temperature distributions are computed separately then they couple on interface with the simple interface treatment of the energy grid function. The proposed method is implemented within the NOISEtte code framework [5,6].
Hence the energy equation for solid and fluid domains is solved with the LINS scheme, which is a generalization of the time integration scheme [7] for parabolic equations. At the fluid-solid interface the matching conditions are the continuity of temperature and normal components of the heat flux. These conditions are natural for the heat conduction equation and there is no need their special treatment [8]. The developed approach demonstrates here in application to solution of a model problem of compressible heat conducting flow, namely the problem of supersonic flows (fluid domain) over a finite thickness flat plate (solid domain).
The proposed technique is recommended for problems in which convective and dissipative processes are in close interaction, therefore this technique can be used to develop a simulation tool for conjugate heat transfer for magnetohydrodynamic flows in channels. The numerical implementation in this work can be considered as a preamble to introduce the new CHT methodology for CFD simulation on multiblock 3D unstructured grids.
The objective of this paper is to present the conjugate heat transfer algorithm in framework of the MCFL technique which essentially uses the explicit iterative LINS scheme for the time integration.
The paper is organized as follows. In Section 2 we present a brief description of the LINS scheme. The structure of this scheme is of key importance for constructing the conjugate heat transfer algorithm. In Section 3 the presentation of this algorithm is given. In Section 4 we discuss the proposed algorithm. In Section 5 a test problem is described and preliminary results are presented. Some conclusions are given in the last section.

Brief description of the LINS scheme
We assume that a gas mixture consists of a set of chemical components with mass fractions Y m , m = 1, ..., N sp . In this paper we consider fluid dynamic flows with multicomponent diffusion and heat transfer, described by the system of Reynolds averaged Navier-Stokes equations which is supplied by convection-diffusion equations, see, e.g., [2]. This system is discretized in the integral form with conservative variables. The finite-volume node-based algorithm [5] is employed for spatial discretization based on construction of dual grid cells. Integrating over each dual cell the semi-discrete scheme is written. The convective, diffusive and heat fluxes are determined at the dual cell interfaces. All the discrete constructions are based on the templates of the package NOISEtte [6]. As a result, the our new computer code MCFL, combined with solid thermal code, inherits the high parallel efficiency of the NOISEtte code.
We introduce the state vector U ≡ ρ (1, u 1 , u 2 , u 3 , E, {Y m , m = 1, ..., N sp }) of conservative variables and two discrete operators on a spatial grid with a characteristic cell diameter h, namely, the convective operator C h and the diffusive operator D h , associated with dissipative processes of viscosity, multicomponent diffusion and heat conduction. In general these operators are nonlinear but for ease of explanation we consider these operators to be linear and write down We omit some unimportant details (source of various origins, boundary conditions and others). The linear operator D h arises from approximation of an differential elliptic operator, this operator is a non-negative definite self adjoint operator in the grid function space. For the explicit scheme (1) the time step τ is restricted by the combined convective-diffusion stability condition. To overcame this restriction we propose the technique that preserves the time step restriction that correspond only to the hyperbolic Courant-Friedrichs-Lewy (CFL) condition.
To do this we use the method of splitting the system into convective and diffusion parts. At each time step we implement two stages. At the first stage we set D h = 0 and solve the intermediate discrete hyperbolic subproblem (corresponding to convective flows) On this stage the modified Godunov scheme with the exact solution of the Riemann problem is used [3]. The time step τ = τ conv is restricted by the CFL condition, aτ /h ≤ 1 , where a is the velocity of propagation of perturbations in a cell of diameter h. The state vectorŪ j is obtained at the hyperbolic stage as intermediate ones. At the second stage we set C h = 0 and solve the discrete parabolic subproblem (corresponding to diffusion flows) with initial dataŪ j The parabolic subtask combines three parabolic subsystems which are taking into account viscosity, multicomponent diffusion and heat conduction. General description of MCFL computational scheme is done for the whole system, but when we consider CHT problem we deal with the heat conduction equation only.
In the scheme (3) one can take U j instead ofÛ j , then the summation of (3) and (2) leads to the explicit scheme (1). This explicit scheme requires the diffusion restriction on the time step size, τ ∼ h 2 . For diffusion-dominated process the usage of the explicit scheme (1) with a such time step restriction is practically impossible. In this case, one can write an implicit scheme, which can be solved with an iterative algorithm. Usually such algorithms require empirical settings (accuracy tolerance, preconditioner parameters, etc). Instead of an implicit scheme we propose the explicit iterative LINS scheme. Interpretation of this scheme: first, in iterations, U j is prepared, and then the scheme (3) is implemented. This whole procedure is based on the construction of the Chebyshev polynomials. The number q of LINS iterations is determined by Here the function s is the smallest integer greater than or equal to s, λ max is the upper bound of the spectrum of the operator D h . This estimate of λ max is obtained by the known Gershgorin theorem. The resulting scheme provides stability and approximation under the hyperbolic time step restriction τ ∼ const · h. The locality and simple algorithmic structure of the LINS scheme makes it an ideal candidate for parallelization.
The maximum step τ dif for explicit diffusion computations is estimated by an inequality τ dif · λ max ≤ 1, which is stronger than the known stability condition | 1 − τ dif · λ max |≤ 1. To ensure the correct transition from the LINS to the explicit scheme, we use the following sufficient condition for the stability of the explicit scheme (1), see [8]: This condition hold if the step size for the explicit time integration of the entire system is determined by the condition τ ≤ 0.5 min(τ conv , τ dif ).

The general computational CHT scheme
To explain CHT algorithm we remark that parabolic stage contains the procedure of solving the energy equation where ε(T ) is the specific internal energy per unit volume of a mixture, ε = ρ · e and e is the specific internal energy per unit mass, ρ is density, κ = κ(x, y, z, T ) is the thermal conductivity coefficient of a gas mixture. For an ideal gas mixture we have ε(T ) = ρC V · T , where C V is the mixture heat capacity at constant volume. The right hand side f takes into account the dissipation of energy due to friction, etc. In a solid, in cartesian coordinates (x, y, z) the quasilinear heat conduction equation is written for temperature T (x, y, z, t) where C(T ) and κ are the coefficients of heat capacity per unit volume and thermal conductivity, respectively, f is the right-hand side, taking into account internal heat sources (by default in solid f ≡ 0). The inputs are known functions of spatial coordinates, temperature, and time.
In isotropic materials, the thermal conductivity k is a scalar; in anisotropic materials, it is a second-rank tensor. The function C(T ) is determined by the density and specific heat of a unit of mass: Instead of (5), we can write the general equation (4) for the internal energy per unit volume In the fluid domain, the equation of the total energy conservation law is written. By the adopted splitting scheme, this equation can be reduced to equation (4), with its input data and the righthand side, which takes into account energy dissipation due to friction, etc.
As an example we consider supersonic heat conducting flow (fluid domain) over a finite thickness flat plate (solid domain). The solid is bounded by a streamlined surface Γ 0 and, possibly, a surface Γ 1 that has no contact with the fluid. On Fig. 1   The surface Γ 0 is the fluid-solid ideal contact boundary and the heat exchange conditions are set on it: continuity of temperature and heat flux is on this interface. The flow domain R F ⊂ R 3 is bounded by a compound boundary Γ = Γ 0 ∪ Γ 2 , where Γ 2 is the boundary of the fluid domain determined by the formulation of the problem. Here on Fig.1 the left and right faces of the fluid is inlet and outlet boundaries respectively, the upper face is far field boundary. The domain R S of a solid, generally speaking, can be inhomogeneous, its boundary may consist of non-intersecting parts: ∂R S = Γ 0 ∪ Γ 1 . The boundary Γ 1 is not of interest, any external boundary conditions for the heat equation must set on it. For example, we set temperature T b on the lower solid shell (isothermal wall), and zero heat fluxes on the left and right faces (see Fig. 1). In the third direction which is perpendicular to the vertical section, i.e. drawing plane, we set the symmetry boundary conditions. In according to adopted scheme, we assume that in the fluid domain, at one time step, we implement all stages of the calculation, i.e. convective, viscous and multicomponent diffusion stages, except for the heat conduction one. At the final heat conduction stage, the equation of the form (4) The equation (4) can be considered as a unified one for fluid and solid media. With conservative approximation of this equation, the heat transfer conditions (6) are satisfied automatically [8].
The presence of boundary conditions for heat transfer on the fluid-solid contact, combines the fluid dynamic and thermal problems in the CHT problem. It is possible to build an unified grid for fluid and solid with grid nodes on the interface Γ 0 , but the equations for these blocks are different and the solvers might be different. Therefore, we define as separate objects two calculation blocks and the fluid-solid interface Γ 0 . At each time step for the explicit scheme, or at each iteration of the LINS scheme, each block is processed separately, then the interface values (specific interior energy and temperature) are corrected. The result is equivalent to the solution of the heat equation on the combined grid using homogeneous discretization. This requires a few explanations. We write down a homogeneous explicit time integration scheme for the equation (4) considering for simplicity the equation of state ε = ρC V ·T for the mixture of ideal gases. Integrating (4) over the dual volume associated with an interface node marked by 0, we obtain at this node 0 the following time integration scheme where Ω 0 is the volume of the dual cell, D(T n ) is the heat flux through the dual cell boundary. The volume of the dual cell can be represented as the sum of partial volumes: Ω 0 = Ω 1 0 + Ω 2 0 . The sum can contain more than two members if there is more different domains share this node.
We assume that the heat capacity and the right-hand side as grid functions are defined at the grid nodes. On the interface at a node 0, they are two-valued functions C 1 0 , C 2 0 and f 1 0 , f 2 0 , and their unique approximation at a node 0 is based on an equality of the form Therefore The expression for the heat capacity C 0 is written in the same way. With such definition of the grid functions, the homogeneous scheme (7) nevertheless masks an inhomogeneity namely the  (7), we compute the temperature distribution at the interface: T n+1 0 = ε n+1 0 C 0 , thereby ensuring the continuity of temperature and internal energy per unit volume. In the case of a general equation of state ε 0 = ε 0 (T ), the temperature is calculated by inverting this equation.
Instead of the explicit scheme (7) we propose to use the explicit iterative LINS scheme as an unified algorithm for solving the heat equation (4) in the fluid and in the solid. The only input parameter of this scheme is the estimate λ max of the upper bound for the discrete operator L h which arises from approximation of the differential operator: where C min is a minimal value of the grid function ∂ε/∂T that is calculated at the lower time layer.
Above we explained the discretization for two heat conducting domains with the fluid-solid interfaces. It is easy to generalize this approach for a multiblock case. The present multiblock algorithm uses point-to-point matching across the interface of neighboring blocks, but a block face might adjoin several neighbors. We explain the multiblock algorithm using the example of the explicit scheme for equation (4). At the interface node of multiplicity l ≥ 2 for the energy equation at the thermal stage it is necessary to obtain the scheme (7). On the right-hand side of formula (7) the member D(T n ) is the heat flux through the boundary of a dual volume. This heat flux can be represented in the additive form here D i (T n ) is the heat flux through the boundary ∂ Ω i of the partial volume Ω i . The total volume is Ω 0 = Ω 1 + ... + Ω l . Let us write the result of partial discretization in each block i for the interface node of multiplicity l ≥ 2: where ε n+1 i , ε n i are old and new values of the grid function ε at the interface node, T n is the grid function obtained at the previous step, or previous LINS iteration, by inverting the equation of state ε = ε(T ). Here n is the time layer number, i is the number of a block, adjacent to the considered interface node. This function T n is already unique value defined in the interface nodes at the lower time layer n, as well as the internal energy: ε n 1 = ... = ε n l = ε n , these are the averaged values after previous interface treatment. Summing the equations (10) over the index i = 1, ..., l, after a simple transformation we obtain Comparing formulas (7) and (11), we obtain an expression for a unique value ε n+1 in an interface node: Above we obtained the same formula for the calculation of the right hand side f in the interface nodes, see (8).
The discrete scheme is constructed in each block including the block boundaries. In the fluid domain, the time step τ conv is determined by the CFL condition. Let the stages of convective, viscous and multicomponent diffusion have already been implemented, then a transition is made to the calculation of the thermal stage. The time step algorithm for calculating the conjugate heat conduction problem is as follows.
1. Setup. Calculate a global estimate λ max of the upper bound of the operator L h (9). Calculate the number p of iterations of the LINS scheme according to λ max and τ conv ; calculate the Chebyshev parameters.
2. The cycle over k = 1, ..., 2p − 1 and the cycle over the domains. Calculate one iteration of the LINS scheme. After completing the cycle over the domains, make correction of the energy at the interfaces.Correct the interface temperature with inverting the linearized equation of state.
3. Calculate the internal energy ε n+1 at the interface according to formulas (12) and recalculate the temperature T n+1 according to the equation of state ε n+1 = ε(T n+1 ).
For multiblock structure, the calculation of the estimate λ max is determined in each heat conducting block by the Gershgorin theorem on the circles of the spectrum. The resulting estimate λ max is found as the maximum of such regional estimates for related heat conducting domains. This estimate may need to be adjusted using the template coefficients on the interface.
The explicit structure of the iterations provides automatic fulfillment of heat transfer conditions and an elegant additive form of the method in the multiblock case, when the discretization becomes non-standard since an interface discretization pattern might lay in several adjacent blocks.

Discussion
Still generally accepted that direct solution of three-dimensional CHT problem is computationally expensive, therefore algorithms with successive iterations over domains are used to achieve the continuity conditions (6), see, e.g., [9]. This means the balance equations (6) at the solid-fluid interface is solved iteratively. The question of the stability of such an iterative process remains open. For example, one can consider 1D two-domain boundary value problem for the heat equation u t = a u xx on a grid which is uniform in each domain. In [9] the conditions (6) are achieved iteratively: the heat flux transfers from the second domain to the first domain as a boundary condition, and temperature transfers from the first domain to the second domain and so on. The necessary stability condition for such algorithm has the form where a 1 , a 2 , h 1 , h 2 are the thermal coefficients and the grid steps in the first and second domains, respectively [10]. This result was obtained by V. Zaguskin and V. Kudryashov in the early 1960s in connection with calculations of two-dimensional gas dynamics with thermal conductivity. For a 3D conjugate heat transfer problem, it is practically impossible meet such stability conditions. Therefore we have proposed the direct time integration LINS scheme which is not an iterative solver for the balance interface equations. In general this approach is known. The novelty of our approach is threefold: (1) developing the LINS scheme whose properties are the locality and simple algorithmic structure, lack of tuning parameters makes it an ideal candidate for parallelization, (2) developing a very elegant additive technique for implementing computation on multiblock conformal unstructured grids, (3) developing multiblock technique for solving CHT problems with various fluid/solid domains.

Numerical test.
The hypersonic flow over the finite thickness flat plate of a high thermal conductivity material is considered as a verification test case for the presented solver. The test statement and thermophysical data are mainly taken from [11]. Schematic presentation of the problem is given in Fig. 1. The fluid parameters. The height of the fluid domain is taken to be same as of its length L = 1 m. Hypersonic flow conditions of free stream Mach number 3 are considered to simulate the steady state flow. The mixture is O 2 − N 2 , oxygen and nitrogen, the mass fractions are Y 1 = 0.2, Y 2 = 1−Y 1 . Freestream parameters of the mixture: temperature T ∞ = 223 K, velocity u ∞ = 900 m/s, pressure p ∞ = 10 5 P a, density ρ ∞ = 1.584 kg/m 3 , viscosity 9.54 · 10 −6 m 2 /s, thermal conductivity λ ∞ = 0.0204 W/(m · K), heat capacity c V = 723 J/(kg · K).
The solid plate parameters. The plate is b = 0.01 m thick. Density ρ = 7500 kg/m 3 , thermal conductivity λ = 2.52 W/(m · K), heat capacity c = 500 J/(kg · K). The boundary conditions: zero heat flux is applied at the left and right plate surfaces, the lower plate surface behaves isothermally T b = 280 K, (i.e. the Dirichlet boundary condition is considered). The initial temperature field at the beginning of the simulation inside the solid is T = T b = 280 K, in the fluid domain T = T ∞ = 223 K. The upper plate surface is the solid-fluid interface Γ 0 with the continuity conditions (6).
A sequence of structured grids are used to discretize the fluid and solid domains for present preliminary investigations. To demonstrate the result, we take a grid (without grid independence research) with 210 × 232 grid nodes for the liquid domain and 210 × 11 grid nodes for the solid domain. The grid consists of the very fine cells near the flat plate (∆x min = 10 −3 , ∆y min = 10 −4 ). We do not use any efforts to refine grid near the leading edge to capture the transition and to resolve singularities at point x = 0. Only we construct an additional block [−0.1; 0]×[0; 1] with 11 × 232 grid nodes in free stream domain.
On Fig. 2 we present preliminary results of computations, wall temperature in Kelvin degrees is shown.  The solid line corresponds to calculation of the flow without consideration of the interface heat transfer (boundary condition is adiabatic wall). The dashed line is the result of solving the conjugate heat transfer problem. The obtained interface temperature is in qualitative agreement with the approximate analytical solution [11]. The temperature peak in the vicinity of the point x = 0 is the result of interaction of the supersonic flow and the leading edge of the plate with the formation of a shock wave. This result is obtained by application of the explicit schemes on both hyperbolic and parabolic stages. Application of the LINS scheme is under development.
The three-dimensional code MCFL uses a parallel environment via message passing interfaces (MPI). The proposed method is generalized for multiblock frame: we implement the conjugate heat transfer algorithm in the case of interaction of the fluid domains with various solid domains using multiblock unstructured conformal grids. Our technique can be used to develop numerical conjugate heat transfer algorithm for calculating magnetohydrodynamic flows in channels.

Conclusion
The present paper shows a new feature of the multicomponent flow solver MCFL. This solver is developing to provide a tight coupling technique for solving the conjugate heat transfer problems in fluid-solid domain. Inside the solid, the heat equation is solved with the time integration finite-volume LINS scheme which is also used for the energy equation in the fluid. The locality and simple algorithmic structure of this scheme makes it an ideal candidate for parallelization. The proposed approach is subject to the testing. The calculation results presented in Section 5, as well as the results of solving simple model problems that are not included in this paper, show the prospects of the proposed approach. Further developments will be done: this technique is required to carry out comprehensive verification on known complicated problems. We plan to quantify the accuracy and computational cost of the LINS scheme on a sequence of refined grids.