Simulation of mixed convection over horizontal plate

Mixed convection over horizontal heated plate was simulated with four numerical models based on Reynolds stress, large eddy simulation (LES) and eddy viscosity approximations. Temperature distributions over plate and in adjacent volume of fluid were the main criteria of results assessment. Three-dimensional computational domain was considered with symmetry boundary conditions. Simulation was performed with Code_Saturne software package in unsteady formulation. Three orthogonal meshes were evaluated to validate initial guess about optimal cell size. u2-f and Smagorinsky LES models appeared to yield the most adequate results. However, temperature distribution in high-buoyancy region, located in the middle of heated plate, was reproduced much more accurate with classical LES and elliptic blending Reynolds stress models. Obtained results are suitable for industrial applications (e.g. cooling jackets) and might be a base ground for further research.


Introduction
Efficient heat exchange is usually associated with highly forced convection, while natural convection is negligible. However, maintenance of forced convection might be too expensive and complex, especially with liquid heat carriers. Thus many industrial heat exchangers (e.g. cooling jackets) tend to operate in mixed convection mode, when pressure and buoyant forces contribute to convection at approximately the same rate. Prediction of flow patterns and temperature distribution is crucial for effective implementation of mixed convection in heat exchange equipment. Understanding of flow affected by mixed convection can minimize thermal stresses and eliminate stagnation zones.
Mixed convection can occur when a downward flow impinges on horizontal heated plate and stagnation point flow (Hiemenz flow) above the plate provides action of buoyant motion perpendicular to forced motion. This common case has been studied extensively, but with a little insight into flow instability and actual temperature distribution. Most of recent studies of mixed convection over heated horizontal plate [1][2][3] reveal Hiemenz flow structure in two-dimensional formulation. At the same time some researchers [4][5][6] show asymmetric nature of Hiemenz flow. The most recent experimental study of Hiemenz flow over horizontal plate [7] provides an important insight into flow patterns and temperature distribution. However, flow instability is described in steady state formulation only. Numerical simulation was conducted to supplement basic dependencies which are based on experimental data [7]. Three-dimensional computational domain was considered since preliminary simulations revealed dramatic effect of dimension reduction on results. Preliminary simulations also gave erroneous results when problem was formulated in steady state. Considering preliminary results and suggestions of previous researchers, numerical simulation of downward water flow impinging on horizontal heated plate was conducted for three-dimensional domain in unsteady formulation. Obtained numerical model can be used for engineering of equipment that utilizes mixed convection over horizontal plate.

Computational domain and mesh
Main part of experimental rig [7] is square-profiled duct with heated plate which is perpendicular to downward water flow. Smooth confuser along with calming meshes provide a uniform downward flow. Upstream part of tested duct was considered as computational domain ( fig. 1). Symmetry conditions were applied to both sides due to prolate form of heated plate. Symmetry allowed reduction of computational domain down to 250x125x150 mm box. Computational mesh was generated with SALOMEan open source software platform designated for pre-and post-processing in numerical simulations. Three hexahedral meshes with cell sizes of 1.25, 1.5 and 2 mm were considered. Meshes consisted of 2.9, 1.8 and 0.83 million cells respectively. Each mesh was inflated near heat transfer surface ensuring y + <1. Considering flow velocity, first row of cells was limited to be 0.12 mm thick. Total thickness of inflated rows (delta) was in a range of 10-15 mm (see fig. 1).

Numerical model
Numerical model was developed with Code_Saturnefree proprietary software designated for computational fluid dynamics and based on finite volume approach.
All simulations were carried out in unsteady formulation. Water was used as working medium with temperature-dependent physical properties. Viscosity, specific heat and thermal conductivity were defined by interpolation of tabular data while density was defined as: where T is local temperature in degrees Celsius. Boundary conditions were defined as follows: top face of computational domain was set as inlet with uniform velocity U=0.02 m/s; bottom face was split into two equal outlets (87.5x250 mm each) and heated wall (75x250 mm, q=7 kW/m 2 ) between them; side faces that are normal to x-axis were set as adiabatic walls with no-slip condition; side faces that are normal to y-axis were set as symmetry planes. Inlet temperature was set to $19 o C. Initial temperature in computational domain was equal to inlet temperature. Iterative method of gradient calculation was used. Unlimited iterative reconstruction of non-orthogonalities was used to receive distinctive feedback on mesh quality. Continuity and momentum equations were not combined since the continuity equation does not have temporal progress term. Transposed gradients and divergence source terms were handled in momentum equation to ensure convergence on relatively coarse meshes. Pressure relaxation (R=0.9) was used to stabilize calculation. Since high pressure gradients were not expected a standard pressure interpolation over computational domain was used. Very robust geometric-algebraic multigrid solver (GAMG) was utilized to solve linearized pressure equation. Momentum and energy equations were solved by diagonal solver (Jacobi). All linear solvers were limited by maximum number of iterations (N max =10 5 ) and maximum precision (epsilon = 10 -5 ). To ensure reasonable accuracy at such strict limitations three test runs were carried out with epsilon 5, 6 and 7. These test runs showed reasonable deterioration of accuracy with an average rate of 3.8 % per an order of epsilon magnitude. Least stable, but most accurate second-order (centered) discretization scheme was used.
Corrected semi-implicit algorithm (SIMPLEC) was utilized to solve Navier-Stokes equations. This algorithm allowed to set relatively large time step (Delta t = 0.01 s), ensuring recommended value of Courant number (Cr max = 5). Providing constant time step, 50 s of experiment were simulated at each run. Residuals analysis along with graphical evaluation of stability of velocity and temperature profiles showed that solution was usually converged in 40 s. Navier-Stokes and conservation equations were closed with Reynolds stress models (RSM) and scale-resolving simulation (SRS) models. Eddy viscosity models showed inadequate results except u2-f model, which was considered for mesh resolution analysis. Elliptic blending RSM (EBRSM) was preferred to Speziale-Sarkar-Gatski (SSG) model since its ability to handle inflated meshes (y + < 1). Classical dynamic large eddy simulation (LES) model and Smagorinsky LES model were considered as part of SRS. Wall-adapting local eddy (WALE) model was unable to reproduce required temperature distribution, whereas other SRS models are not implemented in current version of Code_Saturne (v 4.0). All considered models included convection and diffusion terms in every unknown variable. Reconstruction of convective and diffusive fluxes at the faces was enabled since mesh was orthogonal. Despite second-order terms of convection, diffusion and source terms were expressed in first-order by setting theta = 1:

Data reduction
Models adequacy was assessed by flow patterns and temperature distribution over heated plate. Adopted experimental study [7] revealed distinctive crests of heated water being pushed by buoyancy forces against downward flow. Consistency and height of these crests and corresponding hot stripes on heated plate were major criteria of reproduction assessment of flow patterns. Crests height was assessed by filtering out volume regions with temperature more than T in +1 o C. Distribution of plate temperature was derived from a series of measurements along y-oriented lines, which were discretized by 100 data points. Arithmetic mean of each line data over five time steps (46-50 s) was considered. Spatial-temporal standard deviations were obtained for each dataset.

Results
In general all considered models yielded adequate results ( fig. 2). At the same time one can mention significant spatial and temporal deviations of calculated temperature of heated plate represented on fig. 2 as vertical bars. These deviations were a result of crests movement across the plate.

Fig. 2. Models comparison (error bars show temporal standard deviation).
U2-f and LES models reproduced slight oscillations of hot crests, which were distributed over heated plate in mostly stationary manner. Hot crests modeled with EBRSM (see fig. 3) were moving with approximately constant speed of 0.01 m/s from the center of the plate towards side walls of the duct. While Smagorinsky model provided the lowest crests, other models reproduced well proportioned crests, which were in a good agreement with experimental vitalizations [7]. These crests were rising and falling, forming wave-like structures.

Discussion
Four numerical models based on LES, RSM and Eddy viscosity turbulence approximations were developed to simulate mixed convection over horizontal plate. EBRSM and Smagorinsky models yielded the most accurate results, providing adequate temperature distribution on heat transfer surface and corresponding pattern of buoyancy driven flow over it. Heat transfer coefficient overestimation in the middle of the plate might be the result of large mesh size. As one can see on fig. 3 (b) higher mesh resolution provides more even crests distribution, though these crests are less prominent. Lower profile of crests means higher surface temperatures and consequently lower heat transfer intensity. Nevertheless, mesh size had relatively low impact on obtained results, namely a negligible variation of controlled parameters. It is a known fact that mesh inflation affects computation results to a certain extent, providing dimensionless distance y + <1 at heat transfer surface. Therefore, obtained results supported initial guess about optimal cell size (2 mm). None of tested models was able to predict distribution of local heat transfer coefficient in exact manner. However, all models gave mean error below 2 %. One should mention asymmetry of experimental data [7] at points with |x|<10 mm, which were adopted «as is». Temporal deviations of obtained results are represented on fig. 2 by vertical bars, which illustrate uncertainty rate. Clearly, u2-f and classical LES models adequately reproduced a drop of heat transfer efficiency that occurs in the middle of the plate, whereas EBRSM and Smagorinsky models has failed to make consistent prediction in that region. Closer to the plate edges the situation was an opposite. Considering overall error values EBRSM and Smagorinsky models have been selected as more adequate. Numerical modeling of mixed convection over heated horizontal plate revealed some caveats in terms of handling high-buoyancy driven flow. Nevertheless, EBRSM and Smagorinsky models considered fairly adequate for engineering purposes. Clearly, further research is required to bridge the gap between modeling natural-biased and forced-biased regimes of mixed convection.