Conjugate Heat Transfer Model Based on Simple and Coupled Energy and Heat Equations

In this study, a numerical weak coupling strategy for the modeling of a conjugate heat transfer phenomenon is considered. Where the incompressible Navier Stokes equations are solved using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) as a first step, and then the heat conduction equation for solid is solved in a second step considering the convective velocity field resulting from the first step. A finite-difference approach is used for both discretized time and spatial operators. In this paper, a two-dimensional simulation case study of a steady uniform stream wise flow around heated rectangular and triangle solids is presented. The simulation is forward in time until the steady-state regime is reached as the residuals converge and tend to zero. The spatial analysis of the temperature is obtained through the numerical resolution of the incompressible Navier Stokes, energy equation, and the heat diffusion equation for the fluid and solid media, respectively. The results show the temperature, velocity, and pressure fields in the space domain. The coding is conducted in MATLAB ® , and the flow chart of the method is provided. It was noted that the convection was more dominant than the diffusion.


INTRODUCTION
The fluid mechanics and heat transfer fields of science are bonded closely in many engineering fields and scientific work [1,2]. Fluid mechanics is one of the oldest disciplines of physics, which is focused on the study of gases and liquids' motion and their interaction with the atmosphere. Heat transfer is a branch of engineering which studies the thermal energy transport from one point or object to another due to temperature gradient. Generally, heat transfers through conduction, convection, radiation, or a combination of these three mechanisms [3]. The combination of both physics and their coupling lead to Conjugate Heat Transfer (CHT), enabling the description of complex processes which involve variations of temperature within solids and fluids due to thermal interaction between the solids and fluids present in most of the engineering applications where heat transfer is involved [4][5][6][7][8].

Incompressible Navier Stokes equations
Computational fluid mechanics (CFD) is the study of the motion of a fluid, or its effects, by numerically solving the equations governing the fluid [9]. Depending on the approximations chosen, which are generally the result of a compromise in terms of physical representation needs in relation to the computational or modeling resources available, the equations solved may be Euler equations, Navier Stokes equations, etc. The solutions in CFD methods are obtained by providing discrete approximations of temporal and spatial operators; thus, they always contain a certain amount of error that is issued from the numerical approximation, the physical models, or a combination of both [10][11][12]. With the recent advanced developments in CFD, complex engineering problems can be simulated as a substitution of a complicated or impossible analytical solution or costly and hazardous experiments [13].
where x and y (m) denote the Cartesian coordinates, t (s) denotes the time, u and v (m/s) denote velocity components in spatial directions, (kg/m3), (Pa), = / (m2/s), (Pa.s), , (N), (K), (J/(K kg)), (W/(m⋅K)), = (m2/s), and ̇ ̇ denote the density, the pressure, the kinematic viscosity, the dynamic viscosity, the body force per unit mass acting on the fluid element on Cartesian coordinates, the temperature, the specific heat at constant pressure capacity, the thermal conductivity, the thermal diffusivity, and the heat generation term, in turn. These equations take several different forms based on the application and case study. Equations (1-4) are partial differential equations (PDE) that have complex analytical solutions; thus, they should be solved numerically. To solve the equations numerically by a finite difference method, they should be discretized in time and space domains on a finite discrete set of grid points in the two-dimensional directions. Physical variables such as pressure, velocity, and temperature are numerically evaluated at mesh grid points [27][28][29][30]. The size and number of mesh grid points define the extent of calculations. Considering the symmetry of the cases and specific physical assumptions, realistic and complex 3-D problems can be degenerated to less complex 2-D problems, therefore, providing an offering of a considerable CPU memory consumption gain [31].
The two-dimensional incompressible Navier Stokes equations (Eqs.1-4) are numerically solved using the Semi-Implicit method for Pressure-Linked Equations (SIMPLE [32,33]) iterative numerical technique for continuity and momentum equations (Eqs.1-3) and a finite difference discretization for energy equation (Eq.4). The SIMPLE algorithm and its derivative methods SIMPLER, SIMPLEC, SIMPLEST, and PISO can be implemented to solve turbulence, heat exchanger analysis, chemical reactions, thermal diffusion, multiphase flows, radiation, and several other applications. These methods are usually complex; hence sometimes they need supercomputers to reduce extensive run times. The main iterative steps in SIMPLE algorithms and developing temperature field are the following [33][34][35][36][37][38][39][40]. 1. Solving the discretized u-momentum equation (Eq.2) to produce initial u-velocity 2. Solving the discretized v-momentum equation (Eq.3) to produce initial v-velocity 3. Solving the pressure correction equation to produce the differential pressure dp [32,33] 4. Implementing dp to obtain corrected values ( and ) of initial u and v [32,33] 5. Correcting the initial pressure p to obtain the corrected pressure = p + dp 6. Solving the discretized temperature equation (Eq.3) for the temperature variable ( ) considering the previously computed corrected velocity fields ( and ) at step 4 in order to obtain the temperature field

Convection Heat Transfer
The thermal coupling between the fluid and the solid media occurs at the boundary through the definition of appropriate boundary conditions that define a well-posed problem and guaranties the continuity properties across the boundary. One may remark, that the fluid's energy equation (Eq.4) and the solid's thermal conduction equation (Eq.5) are of the same form without the convection term + for the latest solid equation [47].
The weak coupling from fluid to solid is considered through the use of heat flux boundary condition at the wall (fluid/solid interface) without radiative heat flux, and it is defined by Eq. 6 [43].
Where ℎ , and denote the fluid-side local heat transfer coefficient, local fluid temperature and solid wall surface temperature, respectively.
In this work, the velocity and pressure field has developed for a 2-D space domain using SIMPLE algorithm. The thermal effects are analyzed by developing the transient temperature field in the space domain through the consideration of that the energy equation (Eq.4) and thermal conduction equation (Eq.5) for the fluid and solid domains, respectively. In the following sections, the SIMPLE algorithm, the discretization approach, the numerical test cases, and the obtained results are discussed. In the "Methodology" section, the numerical test cases and the flowchart of the inhouse MATLAB® code are presented. The "Results and discussion" sections provide the pressure, velocity, and temperature fields results and their analysis, in particular their interaction with each other.

METHODOLOGY
In this study, we consider as a two-dimensional simulations of a steady uniform streamwise flow with a constant velocity around heated rectangular and triangle solids where a constant temperature and constant heat flux generation within the solid are considered for rectangular and triangle solids, respectively. The two configurations are shown in figures 1 and 2. The whole computational domain is split into two domains: the objects (rectangular and triangle solids) are considered as solid domain and the rest of the space domain as fluid. In the case of the rectangular-shaped solid, the computational domain was extended in the x-direction in order to avoid numerical interactions between the constant temperature imposed Dirichlet boundary condition and the solid.
In this section, we consider the resolution of the incompressible Navier Stokes continuity and momentum equations (Eqs. 2 and 3) with a SIMPLE approach and energy equation (Eq. 5) with a finite-difference discretization to the velocity and temperature fields, respectively. The velocity field is calculated by implementing the SIMPLE algorithm. The velocity values (u & v) from the SIMPLE method are used as input values in the energy equation (Eq. 5) to obtain the temperature field in the fluid domain, and in particular, around the objects for the final determination of the heat conjugate coupling boundary condition heat flux at the wall (Eq. 6) to be used as a heat flux boundary condition for solving the final heat conduction equation for solid (Eq. 5) using a finite discretization approach in time and space.
The methodology has three parts, first obtaining velocity field by SIMPLE algorithms, second solving temperature field for space domain with a rectangular object, and third solving temperature field for the space domain with a triangular object.  Incompressible flow is assumed for both cases, which means the density does not change with time throughout the space domain; hence will have a constant value. The code is written with arbitrary values for initial and boundary conditions in both cases. The objects are considered as solid walls with no-slip boundary condition, which means the velocities are zero on the objects. The space domain with a rectangular object has no heat generation inside it, while the triangular object is solved with having a constant heat generation within it.
Since all the equations should be solved numerically, they need to be discretized in time and space domain. In this work, the space domain is discretized in a staggered grid, as shown in Figure 3. The discrete form of the incompressible Navier Stokes equations, which are implemented in the SIMPLE algorithm, is given by Eqs. 7 and 8.
where ∆ and and ∆ are the constant grid size in x and y directions, respectively; and ∆ is the timestep value used in the SIMPLE algorithm with respect to the Courant Friedrichs Lewy criteria to guarantee the overall numerical stability. The CFL condition is given by Eq. 9 [40,48].
The SIMPLE method is explained in details in [40], and the MATLAB code is available from [49].
The convective terms of the energy equation are discretized using the second-order central difference method, and it is shown in Eqs. 13 and 14.
Finally, the complete discretized fluid's energy equation is given by Eq. 15.
( , ) For the rectangular solid shape scenario, the temperature is assumed to be constant, and thus the heat conduction equation (Eq.5) discrete resolution is obsolete and not required. However, for the triangular solid shape scenario, a constant heat source of Q =0.1 is considered and the heat conduction equation (Eq.5) needs to be solved numerically with the additional coupling boundary conditions (Eq.6) to incorporate the surrounding fluid's coupling effects. The transient term of the heat conduction equation (Eq.5) is discretized by the forward difference method, and the diffusion terms are discretized by the central difference method, as follows by Eq. 16.

RESULTS AND DISCUSSION
In this section, the results from the implementation of the SIMPLE algorithm with energy and heat equations are discussed. The temperature and velocity fields (u & v) are obtained for the 2-D space domain with rectangular and triangular objects in the center. The graphs are plotted with MATLAB®. It is to be noted that the solution is achieved for the steady-state by running the simulation until the temperatures stop deviation, which can be seen in the related graphs. It is to be noted that the temperature field is achieved for the steady-state by running the simulation for the timestep of 3000 in both cases. Figure 5 shows the temperature field for the 2-D space domain with a rectangular object in its center. It can be seen that how heat dissipates through the space domain. The temperature decreases by getting further from the object edges with 100°C boundary condition. The temperature field diffuses toward the right side of the space domain until it reaches the steadystate of the boundary with 25°C. The temperatures change faster near the object, and it decreases more gradually by getting further from the objects. The temperature contours show that the heat flow is more affected by convention rather than conduction. Figure 6 shows the pressure field of the space domain with a rectangular object developed by the SIMPLE algorithm. A symmetrical pattern can be seen in the pressure contour. It can be seen that pressure values increase towards the edges with different directions indicated by minus sign. Figures 7 and 8 show the velocity field in x and y directions in turn. By looking at figure 7, it can be seen that the velocity pattern is symmetrical. The x-velocity values are bigger above and below the object while it reduces towards the edges and corners. By looking at figure 8, it can be seen that the y-velocity values have bigger values around the corners of the object with opposite directions.      Figure 10 shows the temperature field for the 2-D space domain with a triangular object. The conduction thermal diffusion in the object and convection thermal diffusion out of the object can be seen in the space domain. The temperature changes from 25°C to around 90°C until it reaches the steady-state. Figure 11 shows the pressure field of the space domain. It can be seen that the higher values of pressure are closer to the edges of the object; however, the minus sign in the pressure contour indicates the change of direction. Figures 12 and 13 show the velocity field in x and y directions, respectively. In figure 12, the x-velocity field is symmetrical around the x-axis. The highest values of the velocity are close to the object's upper and lower corners. In figure 13, the y-velocity field can be seen. A symmetrical pattern exists in the field for the values of the velocity about the x-axis; however, the directions change. The highest values of the velocity can be seen near the object edges and corners. Figure 14 shows the steady-state convergence plot through the temperature-time-rate decrease over the number of time steps for the 2-D space domain with a triangular object. The temperature difference reduces with a higher number of simulation runs until the difference meets zero, and a steady-state solution is achieved. The temperature stops deviating around the time step of 1000 in this case.

. Limitations
The grid for this model is coarse and cannot capture the boundary layer. The fluid is modeled based on arbitrary values; thus, flow regime and different regions of the flow, such as separation, stagnation point, or wake, are not of focus in this study. This paper demonstrates a method for the coupling of heat and energy equation with fluid mechanics governing equations. The current model can be used to analyze different problems with various shapes and geometries in relation to conjugate heat transfer.

CONCLUSION
In this presented work, the SIMPLE CFD algorithm is implemented to obtain the velocity field and pressure field in a 2-D space domain with rectangular and triangular objects in its center. The pressure field showed that the pressure values are higher at the sides of the objects in different directions. The velocity fields showed that the maximum velocity exists in the regions near the edges of the objects while having opposite directions. The temperature field is developed for the space domain by discretization and coupling of heat equation and energy equation. The coupled energy and heat equations are integrated with SIMPLE method to achieve the solution. The solution is achieved for a steady-state by an explicit time advance approach until the converged steady state is reached. The obtained temperature field showed that heat diffusion is convection dominated throughout the space domain. The solution stopped deviating at high values of the timestep. The method of the solution and the written code in MATLAB® is provided in the form of a flowchart.