Computational Modeling of Vehicle Radiators Using Porous Medium Approach

A common tool for the determination of thermal characteristics of vehicle radiators is the experimental testing. However, experimental testing may not be feasible considering the cost and labor-time. Basic understanding of the past experimental data and analytical/ computational modeling can significantly enhance the effectiveness of the design and development phase. One such computational modeling technique is the utilization of computational fluid dynamics (CFD) analysis to predict the thermal characteristics of a vehicle radiator. However, CFD models are also not suitable to be used as a design tool since considerable amount of computational power and time is required due to the multiple length scales involved in the problem, especially the small-scale geometric details associated with the fins. Although fins introduce a significant complexity for the problem, the repetitive and/or regular structure of the fins enables the porous medium based modeling. By porous modeling, a memory and time efficient computational model can be developed and implemented as an efficient design tool for radiators. In this work, a computational methodology is described to obtain the hydrodynamic and thermal characteristics of a vehicle radiator. Although the proposed methodology is discussed in the context of a vehicle radiator, the proposed methodology can be implemented to any compact heat exchanger with repetitive fin structures which is an important problem for many industrial applications.


Introduction
One of the main components of a cooling system of an engine is a radiator.Vehicle radiators are typically fin-and-tube-type compact heat exchangers (HXs) and composed of inlet manifolds, outlet manifolds, tubes and fins as shown in Figure 1.Simply, a radiator works with two fluids which are air and anti-freeze water mixture.Hot antifreeze water mixture flows through the tubes, whereas cooling air flows through the fins resulting in heat exchange between both streams.Due to the strong competition in the automotive industry, radiators with better performance (higher cooling capacity, less hydrodynamic loss, less weight, etc.) have been desired.A common tool for the determination of thermal characteristics of vehicle radiators is the experimental testing.However, experimental testing may not be feasible considering the cost and labor-time.Basic understanding of the past experimental data and analytical/computational modeling can significantly enhance the effectiveness of the design and development phase.There are techniques available to analyze HXs such as log mean temperature difference (LMTD) and effectiveness-NTU (ε-NTU).However, these techniques require some parameters known a priori such as overall heat transfer coefficients and/or NTU relations for a given HX.There are no general expressions for overall heat transfer coefficients and/or ε-NTU relations valid for any HX.Therefore, these parameters need to be predicted either from analytical expressions [1], experimental data [2,3] and/or computational models [3][4][5][6].A priori knowledge of these parameters is required for the designer.Therefore, implementation of LMTD and/or ε-NTU is not feasible especially for vehicle radiators which may include customdesigned fin configurations.Alternatively, computational fluid dynamics (CFD) analysis can be applied to predict the thermal characteristics of a radiator.However, CFD analysis of a fullsize HX is not feasible due to extremely high number of cells required to resolve the complex nature of the HXs; especially the fin structures.This point is more problematic when the number of fins is high in the case of heavy-duty vehicle radiators.Although fins introduce a significant complexity for the problem, the repetitive and/or regular structure of the fins enables the porous medium based modeling.From computational point of view, this approach offers some unique advantages.The complex fluid flow occurring through fins can be introduced into the model through porous parameters.Although the determination of these porous parameters requires a rigorous, detailed computational model with very fine mesh structure especially within the regions mainly responsible for the fluid friction and heat transfer, this modeling can be performed on a representative unit cell due to the repetitive nature of the fins.Once these effects are included through the porous parameters, the mesh structure simplifies dramatically and considering the whole geometry, the number of degree of freedom of the system drops down to a feasible number (in the order of 10 millions).Besides, the porous modeling does not require any boundary layer meshing since the friction and heat transfer parameters are already included through the porous parameters.

Porous modeling
Porous modeling is governed by three models.The simplest model is the Darcy's model which is suggested by Henry Darcy (1856) during his investigations on hydrology of the water supplies of Dijon [7].Darcy's equation is expressed as: where, Δp is the pressure drop, l is the pipe length, V is the average velocity, μ is the dynamic viscosity and α is permeability of porous domain.Permeability depends on the fluid properties and the geometrical properties of the medium.The dependence of the pressure drop on velocity in the Darcy's equation is linear; therefore, Darcy's equation is applicable when the flow is laminar.As the velocity increases, the dependence of the pressure drop on velocity becomes non-linear due to drag caused by solid obstacles.At this point, there are two extended models proposed in the literature namely Forchheimer and Forchheimer-Brinkman model.For moderate Reynolds numbers, including nonlinear effects, pressure drop is defined as Forchheimer's equation [7]: where C F is the dimensionless form-drag constant and ρ is the density of the fluid.The first term denotes the viscous characteristics of porous flow and the second term (also called Forchheimer term) denotes the inertial characteristics.Lastly, Forchheimer-Brinkman model includes additional Laplacian term in addition to Forchheimer's equation.Forchheimer-Brinkman model is expressed as [7]: ( where is the effective viscosity.In general, added Laplacian term (also known as Brinkman term) resolves effects of the flow characteristics in a thin boundary layer at the near wall regions.Strictly speaking, the last term becomes important for large porosity (ratio of the fluid volume to the solid volume in a porous medium) values which means the effect is negligible for many practical applications where typically porosity value is relatively small.Eq. ( 3) without the quadratic term is known as extended Darcy (or Brinkman) model.Therefore, Forchheimer-Brinkman model is the most general model, but the inclusion of the Brinkman and Forchheimer term on the left-hand side can be questionable since the Brinkman term is appropriate for large porosity values, yet there exists uncertainty about the validity of the Forchheimer term at larger porosity values [7].
Velocity definition in porous modeling is specified by using two different descriptions: superficial formulation and physical velocity formulation.Superficial velocity formulation does not take the porosity into account during the evaluation of the continuity, momentum and energy equations.On the other hand, physical velocity formulation includes porosity during the calculation of transport equations [8].The continuity and momentum transport equation for a porous domain using Forcheimer's model can be written as [2]: where γ is the porosity, C 2 is the inertial coefficient for porous domain and is the body force term.
Besides flow modeling, heat transfer modeling for porous flow is described by using two models which are (i) equilibrium model and (ii) nonequilibrium model.Equilibrium model (one-equation energy model) is used when the porous medium and fluid phase are in thermal equilibrium.However, in most cases, fluid phase and porous medium are not in thermal equilibrium.For such cases, nonequilibrium thermal model is more realistic.In the case of the radiator, this issue is important since the temperature difference between the solid (fins) and the fluid (air flowing through fins) is the driving mechanism for the heat transfer [4].Therefore, the nonequilibrium model includes two energy equations (known as also two-equation energy model): one is for the fluid domain and the other is for the solid domain.The coupling of these two models is via the term which represents the heat transfer between the fluid and the solid domains.The conservation equations for the two energy model can be written as [2]: where subscripts `s' and `f' stand for solid and fluid, respectively.E is the total energy, T is the temperature, k is the thermal conductivity, S is the energy source term and ( ) stands for the effect of enthalpy transport due to the diffusion of species.The last term in both of the equations is the coupling term which models heat transfer between the fluid and solid domains.In this coupling term,ℎ denotes heat transfer coefficient for the fluid/solid interface and denotes the interfacial area density that is the ratio of the area of the fluid/solid interface and the volume of the porous zone.Through Eqs. ( 3)-( 7), there are many parameters which are material's property and fixed once the materials for the fluid and the solid are selected.On the other hand, there are some parameters (i.e.porous parameters) which are functions of material, geometry and the flow condition.These parameters are γ, α, C 2 , ℎ and .Among these, γ and are purely geometric parameters and can be determined once the geometry of the porous structure is known.In the case of a radiator modeling, once the geometry of the fins is set, these two parameters can be determined beforehand.The other parameters are flow-dependent, meaning that they need to be determined for a specific flow condition.At this point, these parameters can be determined through some analytical expressions [9][10][11], experimental results (e.g.wind tunnel testing) [4,12], empirical correlations [13] and/or computational models [14][15][16][17][18][19][20][21] typically valid for a representative unit cell.All these approaches were implemented in the literature for different studies for the analysis of micro/macro heat sinks and HXs.

Computational modeling of heat exchangers
Porous modeling can be implemented to any geometry which would resemble a porous structure.Moreover, if the porous structure has repetitive nature, the porous coefficients can be obtained through a detailed modeling of representative unit cell through analytical, experimental or computational means.Heat sinks are very good examples of this case and porous modeling approach has been implemented for the analysis of micro/macro heat sinks [9][10][11][12].A two-equation energy model has been implemented to analyze a straight-finned heat sink [9] together with Darcy's model and implemented to perform an optimization for an internally finned tube [10] and to discuss the effect of aspect ratio and effective thermal conductivity on the thermal performance of a micro-heat sink [11] together with extended Darcy's model.The heat sinks contain a regular structure; therefore, there is a chance to derive analytical expressions to estimate the porous parameters [9][10][11].
Considering the HXs with complex fin structures, computational modeling is even more challenging; therefore, the computational models typically focus on specific subcomponents of HXs such as a representative unit cell for the fin structure [14][15][16][17][18][19][20][21], radiator fan [22] and inlet manifold [23,24].The thermal performance of a HX can be achieved by simply increasing the performance of the fin structure alone.A fin structure with higher heat transfer together with less pressure drop can significantly enhance the performance of the entire system.To investigate the thermal performance of a fin structure, experimental [14][15][16][17][18] and/or computational models [14,17,[19][20][21] can be realized for different fin geometries.Moreover, improving the flow maldistribution at the inlet manifold may also increase the thermal performance.Computational modeling of the flow maldistribution may lead to performance enhancement for HXs [23,24].
Analyzing subcomponents may lead to qualitative conclusion for the thermal performance of an HX, however to estimate the thermal performance quantitatively, a rigorous 3-D modeling of the entire HX is required.Since a rigorous modeling is not computationally feasible, a 2D model [4], hydraulic and thermal resistances-based models [12,25] and 3D mesoscale models (considering macro control volumes) have been introduced in the literature to predict the thermal performance quantitatively [26][27][28][29].A 2D model was developed to compare the equilibrium model (one-equation thermal model) and nonequilibrium thermal model (twoequation energy model) for a relatively small size matrix type HX [4].A resistance-based model was implemented to predict the hydrodynamic and thermal performance of a carbon-foamfinned HX which combined many different correlations from the literature to predict the hydrodynamic and thermal resistances [25].The success of the model strongly depends on the accuracy of the porous parameters.For this particular example, the model was proven to predict the hydrodynamic and the thermal performance within ±15% of the experimental data.A Compact Heat Exchanger Simulation Software (CHESS) has been developed [26][27][28] as a rating and design tool for industrial use based on the empirical correlation of the porous parameters to analyze the fin-and-tube part of a vehicle radiators (excluding inlet and outlet manifolds).It was demonstrated that by using CHESS, the thermal performance of different vehicle radiators was predicted within ±15% of the experimental values.Alternatively, a porous modeling-based CFD model for fluid flow and meso-scale ε-NTU-based modeling for thermal characteristic was utilized for an air-to-air cross-flow HX [29] to investigate the effect of the maldistribution on the thermal performance.A 3D CFD model coupled with porous medium approach has been developed to investigate the hydrodynamic performance of a plate-fin HX in which the porous parameters were also determined using a detailed CFD model on the unit cell [30].
A full-size 3D thermal modeling of a relatively small compact HX was conducted with different fin configurations and the heat transfer and friction factor parameters which can be used in conjunction with the LMTD or ε-NTU method [5].Since the size of the HX was small, the meshing was not a problem and the computational model was utilized for the design of an inlet manifold for a better performance.Considering the size of the vehicle radiator, this approach is not an option.Thermal and structural analysis of a heavy-duty truck radiator which had finned structure both on the liquid and air side has been performed using Commercial CFD software, FLUENT ® [31].Forchheimer's relation was used for the porous modeling together with the experimental data.One-equation energy model was used together with the averaged equivalent thermal conductivity.The local heat transfer coefficients and pressure distribution gathered from the thermal analysis were used as a boundary condition for finite element structural analysis through which the thermal stresses and strains were obtained.
One alternative to all these approaches can be the modeling of the vehicle radiator with porous medium approach where the porous parameters are also deduced from a rigorous CFD modeling on a unit cell.Moreover, this procedure may be performed with a commercial CFD software which would have very strong meshing, solving and post-processing capabilities.However, implementation of the two-temperature energy equation is crucial for an accurate prediction especially for the vehicle radiators.This may not be straightforward with a commercial software.At this point, FLUENT ® may be a viable solution since the two-temperature model capability has been included in version 14.5.More recently, a computational modeling of a fin-and-tube-type vehicle radiator has been conducted based on two-temperature model and the cooling capacity of a heavy-duty vehicle radiator has been estimated without any need for empirical and/or experimental data [32].In the upcoming section, the computational methodology of such a computational model is outlined.This approach may allow CFD modeling to be an efficient rating and design tool for vehicle radiators.Although the proposed computational methodology is discussed for a vehicle radiator, it may also be implemented to any compact HX with repetitive fin structures which is an important problem for many industrial applications.

Computational modeling
The proposed computational methodology is implemented for a 4-row 39-column commercial available heavy-duty vehicle (more specifically tractor) radiator as shown in Figure 1.Tractor that uses the manufactured radiator has a 64 HP Perkins engine which requires a minimum cooling capacity of 55 kW according to the catalog data.The cooling capacity of this radiator was reported as 55.8 kW by the tractor company as a result of in-house experiments following the SAE-J1393 protocol [33].Catalog data are tabulated in Table 1.The fin structure used on this radiator is a wavy fin (WF) structure which is a typical structure used in vehicle radiators due to its superior thermal performance.The selected wavy fin configuration is 84 mm in length.The computational procedure starts with the determination of the porous parameters for a given mesh configuration.The geometric parameters are determined using the CAD model.On the other hand, to determine the flow-based parameters, a parametric study needs to be performed on a unit cell with high resolution which consists of one repeating section of the fin structure.For the determination of the porous medium coefficients, the flow field should be analyzed only for the section with the finned structure (physical fin simulations).To verify the extracted porous medium coefficients, the flow field within the unit cell together with included upstream and downstream fluid domain needs to be modeled both using actual fin geometry and porous modeling.This analysis needs to be performed only once for each mesh configuration of interest.

Determination of the porous parameters
Fin analysis is progressed under three main steps: a. Simulating the unit cell straight fin model by using different air inlet velocities and obtaining the resultant pressure drop across the fin.
b. Fitting a second-order curve to the collected pressure versus velocity data gives the Darcy-Forchheimer's relation as: where a and b are the coefficients characterizing the flow.Obtaining the flow-based porous medium coefficients is followed by the determination of input parameters for heat transfer modeling.The necessary input parameters are the average heat transfer coefficient (HTC) and interfacial area density (IAD) for the two-equation energy model.Average heat transfer coefficient is obtained from FLUENT ® post-processing which can be calculated by using the following relation: The reference temperature in the above equation is the average temperature of the air between the inlet and outlet of the finned channel.

Physical fin simulations
Unit cell of a wavy fin model, Model-A shown in Figure 2(a), is analyzed in order to obtain the porous medium parameters.Flow parameters are obtained by using the Forchheimer's relation.Model-A is simulated using different Reynolds numbers to obtain Forchheimer's curve.Once the parameters are obtained, Model-B, which is a unit cell of the wavy fin with additional upstream and downstream domains as shown in Figure 2(b), is analyzed.Since the air domains (without fins) are attached at the inlet and the exit of the porous domain, the flow area contracts (at the inlet of the porous domain) and expands (at the exit of the porous domain) at the interfaces of these domains.To capture the physics, porous-jump boundary conditions are introduced to match the results of the two models [8].Boundary conditions for Model-A are set as follows: the velocity inlet and pressure outlet boundary conditions are assigned for the fin inlet and outlet, respectively.Wall boundary condition is applied for the upper and lower walls.Constant wall temperature boundary condition is assigned to the walls as the thermal boundary condition (which is close to the real situation, since the temperature variation in the z-direction is small).Periodic boundary condition is used for the right and left sides.For Model-B, additional upstream symmetry and downstream symmetry are assigned for upstream and downstream domains.For both simulations, SIMPLE method is used with a least square-based approach for gradient reconstruction.In addition, standard scheme for pressure and the second-order up-winding schemes for momentum, turbulent kinetic energy and turbulent dissipation rate are employed.Relaxation factors are set to their default values.
For both simulations, a minimum convergence of 1 × 10 −5 is obtained for all residuals.One important step is the determination of the appropriate turbulence model.At this point, some benchmark solutions, empirical/experimental results can be employed for the determination of the appropriate turbulence model.For fin structure under consideration, k-ε realizable turbulence model with standard wall function is used (the detailed discussion on the justification of the use of this turbulence model can be found elsewhere [32]).Afterward, mesh independence analysis needs to be performed to ensure the mesh-independent solutions.It is observed that approximately 4,900,000 number of cells with 30 layers of boundary mesh for the fins generates a mesh-independent result for this particular fin configuration.Table 2 contains the input parameters for the Model-A.In Figure 3, the pressure drop across the fin structure is plotted against velocity and a second-order curve is fitted to the simulation data.The corresponding inertial and viscous coefficients are determined as 17.3 and 4.01 × 10 6 , respectively.Heat transfer parameters are obtained from the simulations of Model-B.For the simulation of Model-B, the input parameters are defined as 7.0 m/s for the inlet velocity, 304.2K for the inlet temperature and 359.7 K for the temperature of fin walls in accordance the tabulated catalog data.Average surface heat transfer coefficient and tuned porous jump coefficients for the unit cell of a wavy fin are presented in Tables 3 and 4, respectively.

Fin simulations with porous modeling
Once the porous coefficients are obtained, the flow field of the air can be modeled using the porous modeling.Upstream and downstream domains are attached for this analysis as shown in Figure 4(a).To verify the porous modeling, the results are compared with the physical fin simulations.For the porous fin model, hexa-sweep meshing is used.The mesh of the porous model (Figure 4(b)) consists of 5320 cells.After completing the meshing process, boundary conditions are assigned.Besides the physical fin boundary condition configurations, additional porous-jump boundary conditions are introduced to match the physical fin simulations.All solver settings are taken to be the same as the physical fin simulations.After porous medium flow coefficients, porous-jump coefficients and heat transfer parameters are obtained from the simulation of a unit cell of a wavy physical fin and porous medium simulations are carried out with the same input parameters to verify the porous modeling.together with a reasonable accuracy.According to the presented results, pressure and temperature drop characteristics are coherent between the physical fin and porous medium.Contour representations for y+, velocity and temperature distribution across the fin are presented in Figure 6.It is seen from the Model-B results that y+ values are acceptable with respect to analysis results (for SST turbulence model maximum y+ value should be smaller than 1.0) [8] and velocity and temperature distributions have convenient characteristics.Cross-sectional velocity and temperature distributions for the air-side and the water-side streamlines are presented in Figures 7 and 8, respectively.Temperature gradients are successfully achieved in z-and y-directions as expected.Air-side temperature is increasing in the flow direction as a result of the heat transfer from the water-side, while the water-side temperature is decreasing in the flow direction.The flow is not distributed uniformly among the tubes as shown in Figure 8(a).However, to improve the performance of a radiator, the maldistribution of the flow at the header needs to be reduced [22][23][24].Therefore, one can clearly state that there is a room for improvement for the design at hand.This nonuniform distribution of the flow among the tubes contributes also to the temperature in the x-direction.According to the simulation, the average outlet water temperature is found to be 354.3K and total temperature drop of water through the radiator is calculated as 5.4 K which leads to a total heat capacity of: p Q mC T 2.41 x 4208 x 5.36 54.4 kW The pressure drop for water which is also an important performance parameter for radiators is found to be 6.5 kPa.According to the catalog data, the outlet water temperature, temperature drop across the radiator and the cooling capacity are 354.2K, 5.5 K and 55.8 kW.The same parameters are found to be 354.3K, 5.4 K and 54.4 kW with the proposed CFD analysis.The deviation of the CFD results with the catalog is within 2.5% which is quite acceptable for a thermal analysis.Moreover, the proposed model solves the problem within a reasonable computational time.Considering the accuracy of the result and computational cost, the proposed methodology can be used as a rating and design tool for vehicle radiators.

Concluding remarks
Although the repetitive fin structures introduce a challenge for the computational modeling of a radiator, the repetitive nature also enables an efficient porous medium modeling.Moreover, again due to the repetitive nature, the porous parameters can be obtained by CFD modeling of a representative unit-cell with high resolution.A successful implementation of porous modeling can lead a dramatic reduction in computational cost and time.The implementation of the computational methodology through a commercial software also benefits from the powerful meshing, solving and post-processing capabilities.As demonstrated, CFD analysis of a radiator by using porous medium approach gives reasonable and reliable results.By using CFD analysis, design cost may be decreased dramatically by easing the experimental testing process.The porous parameters of a given fin geometry can be obtained within a couple of hours which may enable the hydrodynamic and thermal optimization of a radiator.
Optimization of radiators in terms of size and weight is desired to keep up with the constraints within competitive automotive industry.An efficient computational model enables the optimization process to be performed computationally for a range of different design parameters.Furthermore, more realistic computational models may be developed such as the inclusion of the radiator fan into model or the inclusion of the under hood equipment together with the increasing computational power of the computers.On top of these, the coupling of the flow and temperature field with the structural analysis may lead to far more efficient and robust radiator designs.

c.
Obtaining the inertial coefficient and viscous coefficient using the extracted coefficients in step (b) as:

Figure 2 .Table 2 .
Figure 2. (a) Model-A: WF unit cell domain and (b) Model-B: WF unit cell with inlet and exit domains.

Figure 3 .
Figure 3. Unit cell physical WF simulation pressure versus velocity.

Figure 5
(a) compares the sectional-averaged pressure drop for the physical fin and porous medium simulations.Figure5(b)shows the same comparison for the sectional mass-flow-averaged temperature drop.As seen from Figure5, an acceptable consistency is achieved with the porous modeling.One should note that the porous medium requires only 5320 cells per unit cell mesh; on the other hand, physical fin requires 4,900,713 cells.If a full-sized radiator is modeled with physical fins, the required cell number will be approximately 20 billion which is not feasible to analyze even with today's computing technology; therefore, by using porous modeling approach, full-sized model can be analyzed within reasonable computing time

Figure 4 .
Figure 4. (a) Unit cell porous model with inlet and exit domains and (b) mesh configuration.

Figure 5 .
Figure 5.Comparison of physical WF and porous model: (a) sectional average pressure drop and (b) sectional massflow averaged temperature drop.

2. 4 .
Radiator modeling3-D CAD model of the 4-row 39-column radiator is prepared by using CAD software.After forming the 3-D model, the meshing process is progressed.Fin, upstream, downstream and tube domains are meshed with hexa-type elements, while the upper and lower tanks are meshed with tetra elements.Tubes are meshed with a boundary layer mesh having two layers with 0.1 mm first layer height.The generated mesh consists of 53,355,356 cells with an average skewness value of 0.178.Mass flow inlet, pressure outlet, velocity inlet, pressure outlet, upstream wall and downstream wall boundary conditions are assigned for water inlet, water outlet, air inlet, air outlet and the outer surface boundary of the upstream and the downstream domains, respectively.The air inlet velocity is taken as 7.0 m/s with an inlet temperature of 304.2 K temperature, while the mass flow rate of water is 2.41 kg/s with an inlet temperature of 359.7 K in accordance with the catalog data.Second-order upwind scheme is used for momentum, turbulent kinetic energy (TKE) and turbulent dissipation rate (TDR).Relaxation factors are selected as 0.05 for momentum, 0.3 for TKE and TDR and 0.4 for turbulent viscosity in order to obtain optimized convergence rate and solution time.The heat transfer coefficient between the fins and air is taken as 170 W/m 2 K referring to the previous unit cell simulations.A converged solution is obtained after 472 iterations when the minimum residual is smaller than 1 × 10 −4 .The simulations are performed on a DELL T5600 Workstation (Intel ® Xeon ® , 3.30 GHz, 2 processors, 16 cores, 128 GB RAM).The overall solution time is observed to be approximately 12 h and 40 min.

Figure 8 .
Figure 8. Water-side streamlines (a) colored according to the velocity and (b) colored according to the temperature.

Table 1 .
Catalog data for the four-row radiator.

Table 3 .
Porous parameters of a WF.

Table 4 .
Porous jump coefficients for a unit cell of a WF.