Determination of Turbulent Prandtl Number for Thermal Fluid Dynamics Simulation of HVAC Unit by Data Assimilation

A heating ventilation and air conditioning (HVAC) unit is an essential unit to adjust temperature for passenger's comfortability in an automotive cabin as shown in Fig. 1, dehumidifying/cooling and heating the air with the evaporator and the heater core inside the HVAC unit, and then adjusting the temperature with the airflow control door and switching the air outlet. In recent years, due to high performance requirements of automotive air conditioning systems, the complexity of HVAC internal flow path is progressing, i.e., left/right independent process control type for automobile, rear air conditioning for the rear seat passengers, etc., understanding the heat transfer phenomenon inside the HVAC unit becomes more and more important in HVAC unit design [1–3]. To achieve the demand target in performance improvement at an early stage of development, it is indispensable to grasp the information on flow and temperature field using thermal fluid dynamics (CFD) simulation. However, if a default setting of the commercial CFD simulation software is utilized to predict temperature at duct outlet of HVAC unit, in which Prandtl number is usually set to about 0.7–0.9 as laminar flow, there is a big difference between simulation results and the measurement data, because the airflow inside HVAC unit has a very complex turbulence vortex structure and the structure of the turbulence airflow changes with the variation of the gap between the airflow control door and the associated shut rib, due to the various opening degrees of the control door. Therefore, as for high-accuracy prediction of HVAC performance at an initial stage of the development using CFD simulation, the determination of the turbulent Prandtl number $(Prt)$ in the CFD model is a key issue [4,5].

The $Prt$⁠, which represents the dissimilarity between turbulent momentum and heat transfer, is useful for solving the heat transfer problem of turbulent flow. So far, based on Reynolds' analogy, many experimental and theoretical investigations have been carried out for determination of the $Prt$⁠; Blom [6] experimentally determined the $Prt$ in a development temperature boundary layer, Reynolds [7] examined more than 30 ways of predicting the $Prt$⁠, Browne and Antonia [8] measured the $Prt$ in the self-preserving region of a slightly heated plane jet, Kays [9] examined the experimental data on the $Prt$ for the two-dimensional turbulent boundary layer and for fully developed flow in a circular tube or a flat duct, Morris et al. [10] investigated the role of four alternative turbulent Prandtl number functions on the predicted heat transfer coefficients obtained by numerical analysis, Churchill [11] presented a detailed interpretation of the $Prt$⁠, Redjem-Saad et al. [12] investigated the $Prt$ by a direct numerical simulation (DNS) of heat transfer in a fully developed turbulent pipe flow, Venayagamoorthy and Stretch [13] derived a general relationship for the $Prt$ for homogeneous stably stratified turbulence from the turbulent kinetic energy and scalar variance equations, Li [14] theoretically investigated the $Prt$ in the atmospheric boundary layer, Basu and Holtslag [15] proposed a simple analytical approach based on the variance and flux budget equations to quantitatively evaluate the stability dependence of the $Prt$⁠, Xu et al. [16] investigated the behavior of the $Prt$ for buoyancy-affected flows near a vertical surface using machine learning techniques based on DNS data base, etc. However, due to theoretical deficiencies and lack of measurement precision, the results showed that qualitative and quantitative differences were too large to draw conclusive conclusions about the general situation. Fortunately, several pieces of evidence have shown that the $Prt$ which should be a function of local turbulence properties is often assumed in many engineering calculations to be a function of velocity and temperature, although a clear variation characteristic of the $Prt$ in the practical range has not yet been clarified.

In this study, we endeavor to improve the accuracy of the temperature field prediction in turbulent mixing for the HVAC unit, tuning up model parameter (the $Prt$⁠) in the CFD simulation by functionalization and data assimilation with actual measurement data, and the effects are compared and discussed. First, to understand the mixing process of hot and cold air in HVAC unit, we carry out an experimental study on mixing field of a simple HVAC model, and the velocity and temperature fields of mixing zone are simultaneously measured using particle image velocimetry (PIV) and thermocouples. Second, to compensate for the deficiency that the spatial and temporal measurement resolution of the temperature field is lacking compared to velocity field measurements, data assimilation is performed based on these experiments that simulate the mixing field inside the simple HVAC model and the $Prt$ distribution that fits the experimental data which are predicted. Finally, to employ the predicted $Prt$ in the thermal fluid simulation of the actual HVAC unit, we use multidimensional analysis to derive two functions of the $Prt$ for high velocity and low velocity regions.

We experimentally create several hot and cold air mixing fields with different velocity and temperature conditions in a simple HVAC model, and measure the air velocity and temperature of these mixed fields using particle image velocimetry and thermocouples, respectively.

Figure 2 shows a schematic diagram of experimental apparatus: two dryers with multi-stage speed and temperature adjustment functions are installed between the smoke chamber and the simple HVAC model. The outside air drawn into the smoke chamber is split into two air streams, which are heated and cooled by each dryer. After passing through each rectifier, the hot and cold air streams are mixed inside the simple HVAC model. By adjusting the velocity and temperature of the blowing air with the two dryers, it is possible to reproduce the various mixing field that are expected in an actual HVAC unit.

The simple HVAC model is made of transparent acrylic plate (thickness 5 mm) as shown in Fig. 3, and two guide plates are installed instead of the control door in the actual HVAC unit. In this experiment, the smoke of incense sticks is used to visualize the velocity field by PIV, and the temperature data of 60 points in the air mixing zone of the simple HVAC are measured by thermocouples, where the diameter of thermocouple is 0.05 mm, and be fixed on the pins with a height of 25 mm. However, if the temperature data of 60 observation points are measured at the same time, it is concerned that the pin in the upstream disturbs the downstream and the physical quantity in the mixing field cannot be observed correctly. Therefore, under a predetermined experimental condition, only one row of 10 pins is installed in the x direction, and multipoint temperature measurement is performed six times with moving the position of pins from the front to the rear in the y direction.

Figure 4 displays an overview of the experimental apparatus: the PIV system consists of a high speed camera (Photron, Fastcam SA-X2) equipped with Nikkor 50 mm f/1.2 that has a narrow band bandpass filter to avoid outside light, a double cavity Nd:YLF laser (Photonics Industries, DM30-527DH), and an analysis software (Dantec Dynamics, Dynamic Studio 6.1). The pair of particle images is captured at 2 kHz by using the frame-straddling method and instantaneous velocity field is calculated with adaptive PIV algorithm. Mean velocity fields and statics are calculated from 3638 instantaneous velocity files.

In this experiment, to create as many variations in the velocity and temperature as possible, considering the range of velocity and temperature inside the actual HVAC unit (0–15 m/s and 5–85 °C), the experimental conditions are set using multi-stage adjustment function of the dryer, and these specific values are obtained by anemometer and thermocouple measurements at inlets of the simple HVAC model as shown in Table 1. The surroundings temperature and relative humidity were maintained at 20 °C and 40%, respectively.

As a simultaneous measurement of the velocity by PIV and the temperature by thermocouple, it is necessary to discuss three influences on the flow field, such as the particle size, the particle concentration, and the temperature detection part of thermocouple. First, since the particle size in incense smoke is small enough (approximately less 1 μm), the influence of the particle size can be omitted [17]. Next, the fans in the two dryers diffuse the incense smoke to produce uniform smoke in this experiment, and we also carry out a preliminary experiment to determine the number of incense sticks in the chamber so that the variation in smoke concentration have no effects on measurement, so the influence of the particle concentration can also be ignored. Last, in order to avoid the influence from the temperature detection part of thermocouple, with adjusting the installation position of laser light sheet, the temperature detection part of thermocouple is fixed as close as possible to the measurement surface of PIV without disturbing the flow field; Fig. 5 compares the measured data of the velocity component without thermocouples and that with thermocouples under experimental conditions of test 04; it shows that thermocouples have almost no effect on the flow field of PIV measurement.

Figure 6 shows the error of the measured value and the corrected value under the experimental conditions of test 02; the maximum relative error and the absolute error are less than 2.0% and below 0.8 °C, respectively. However, the corrected value of temperature is used in the research.

Hereafter, only the results in which remarkable differences are observed under each of the experimental conditions are described. We select the case of high velocity (test 03) and the case of high velocity difference (test 04) with large temperature difference between the hot air and the cold air.

Figure 7 shows the mean velocity field and the mean temperature field in case of test 03 and test 04. Although the temperature difference between the hot air and the cold air is almost the same, the temperature fields are mixed well in the case of test 03, the reason is that a stronger vortex distribution is generated in the high velocity field of test 03 and the streamwise vortex promotes mixing as shown in Fig. 8, where the dimensionless vorticity is defined as $ωzd/Ul$ using the distance between thermocouples d and the velocity of cold air $Ul$⁠. As a result, by damming the flow with guide plates, a strong exfoliation shear layer is formed, and it also appears that the mixing is facilitated by mechanisms such as Kelvin–Helmholtz instability.

The velocity field of turbulent flow can be split into a mean term and a fluctuating term using Reynolds decomposition; the Reynold stresses component $u′v′¯$ reflects the fluctuating velocity intensity distribution in the turbulent field as shown in Fig. 9; it shows that the fluctuation strength in the case of test 03 is larger than that in the case of test 04, and especially in the area between the two guide plates.

Since the heat transfer between the hot air and the cold air predominantly occurred in the direction perpendicular to the mainstream direction, a local coordinate $(xi*,yi*)$ is introduced as shown in Fig. 10 

Figures 11(a) and 11(b) plot the turbulent heat flux $v*′T′¯$ and the average temperature gradient in the direction perpendicular to the mainstream direction in case of test 03 and test 04, respectively. In comparison with two distributions, it is hard to find a clean relationship with turbulent heat transfer and temperature gradient, but both of these maximum intensity values appeared in the mixing zone between two guide plates.

To sum up, since the streamwise vortex and the fluctuating velocity intensity enhance the mixing process of the hot air and the cold air, and the maximum intensity value of those always appeared in the area behind the two guide plates, it meant that installing ribbed airflow control doors in the flow path of HVAC unit which perform the same function as the guide plates was an effective method to promote mixing.

The $Prt$ is the model parameter for data assimilation estimation. From Eq. (5), the heat conductivity of the air changes with respect to the variation of the $Prt$⁠, so the value of the $Prt$ is a parameter that controls the heat transfer within the mixing field. Normally, it is often assumed that the value of the $Prt$ is about 0.9 in the case of air, but the temperature field of simulation does not match the experiment. Then data assimilation modifies the value of $Prt$ at each location so that the simulation results of the mixing field match the observed data in the experiment.

In the study, CFD simulation is carried out by a commercial CFD software Siemens star-ccm+ 13.06.012. Regarding the simulation model, the traditional setting for previously validated CFD simulation is employed as follows, the semi-implicit method for pressure-linked equation is used to solve the governing equations discretized by a finite volume method, and the quadratic precision upwind method using the gradient by the cell-based minimum square method is applied to each term of the Navier–Stokes equation, the realizable k–ε turbulence model is chosen here. As the boundary conditions, air velocity and temperature are fixed at the inlet, which are taken from the experimental conditions in Table 1, and the atmospheric pressure and temperature are imposed at the outlet, the turbulence intensity and the turbulence length scale are set to 1% and 10 mm at the inlet and outlet, and the two layer all wall y+ model is adopted on the wall of the HVAC model with no-slip condition and a constant thermal resistance.

First in the flowchart, the method of generating the initial ensemble $x0|0(i)$ (i = 1, 2, …, N_ensemble = 927,594) is arbitrary, but here it was generated based on the $Prt$ function of Ito et al. [4]. First, a convergent solution of the flow field is obtained using the $Prt$ function under predetermined test conditions. The x-direction velocity u, y-direction velocity v, temperature T, and the $Prt$ obtained are stored in the initial ensemble $x0|0(i)$⁠. Finally, an initial ensemble with a different $Prt$ distribution is generated by adding random numbers to the $Prt$ part of each ensemble member. Since each ensemble has smooth $Prt$ distribution and a trend appears in the model parameter distribution, this method has an advantage that it is easier to construct the $Prt$ function from the estimated $Prt$ distribution, compared to an initial ensemble that was generated by simply throwing random numbers and adding them to a reference value.

In addition, w is the observation error vector, which is currently set as w ∼ Norm (0, R) with R = 0.1I (I is the identity matrix). The current value of the observation error is tentatively employed for successful data assimilation.

Data assimilation is carried out under each experimental condition, and the $Prt$ function [4] is used as a reference value for generating the initial ensemble; first obtain the convergence solution of the flow field by the reference value under the predetermined experimental conditions, since the $Prt$ is a function of the flow field, and it is a smooth distribution at the end of the fluid analysis. Here, we create one ensemble member by adding one random number from the uniform distribution to the entire distribution.

The target of this data assimilation is to determine only the $Prt$⁠, since the $Prt$ only affects the temperature field and does not change the velocity field, and we will focus on the temperature field from now on. Figure 14 shows the calculation results of temperature using the $Prt$ of the data assimilation estimation at each observation point under the conditions of test 03 and test 04. The blue diamond points are the measured values, the short vertical blue lines represent 95% confidence interval error bars, the green square points are the calculation results using the $Prt$ estimated for data assimilation, and the red triangle points are the calculation results using the conventional function of the $Prt$ by Ito et al [4]. The average temperature absolute errors for each measured temperature are 4.59 °C and 2.55 °C for data assimilation and 7.34 °C and 6.77 °C for the conventional $Prt$ function, and this reduction in errors is confirmed under other experimental conditions as well, −1.82 °C for test 01, −2.88 °C for test 02, −0.47 °C for test 05, and −1.13 °C for test 06; it indicates that the prediction accuracy of the temperature field is improved by correcting the $Prt$ for data assimilation. Here note that the deviation between the experimental data and the calculation result obtained through data assimilation increases with the numbered measurement points in the case of test 03 but such scenario is not there in the case of test 04. This is because that the high-numbered measurement points are installed close to the exit of the simple HVAC model, where the air flow is more unstable and turbulent in the case of test 03 due to high air velocity. Since the velocity field in the turbulent region is hard to resolve without considering the more precise turbulence model, the deviation of velocity similarly becomes larger with the numbered measurement points as shown in Fig. 15.

This is an issue to be addressed in future work. It indicates that the data assimilation leads to the converged CFD solution, which is consistent with the actual measurements of temperature as shown in Fig. 14. Regarding velocity shown in Fig. 15, on the other hand, there is still a large deviation even by the data assimilation in the case of high velocity air flow, and the same trend is observed under other test conditions. This may be because the discrepancy mainly comes from the turbulent model itself, and the current used turbulent model is simple and not suitable for real turbulent flow. For further improvement in CFD simulation accuracy assisted by the data assimilation, therefore, we are thinking of a possibility to optimize parameters of the turbulent model with data assimilation.

Next, Fig. 16 plots the variation of the $Prt$ with air velocity under all test conditions. The triangles indicate the calculation results of Eq. (4), and several unrealistic data are excluded due to the lack of the measured data; the circles represent the results from data assimilation. It shows that the $Prt$ is reducing with decreasing mainstream velocity, and when the mainstream velocity exceeds a certain value (u_c = 5.0 m/s), the $Prt$ becomes larger when velocity increases. In comparison with data assimilation results and calculation results of Eq. (4), it seems that the data assimilation results have a clear relationship with the flow velocity, and the $Prt$ can be easily functionalized using the data assimilation result.

Functionalization is a convenient way to apply the $Prt$ data assimilation results to thermal fluid simulations of actual HVAC units, and it is generally thought of as an exponential function of Reynolds number $Re$ and molecular Prandtl number $Pr$ such as $Prt≡a×Reb×Prc$⁠. However, in this study, according to the variation characteristics of the $Prt$ with respect to velocity (see Fig. 16), two functions for low velocity $(U<uc)$ and high velocity $(U≥uc)$ regions are proposed, and then the key variables for the function of $Prt$ are determined from several physical quantities by sensitivity analysis on a commercial optimization analysis software (modefrontier); here these physical quantities include magnitude of velocity U, temperature T, velocity gradient $dU≡(∂U/∂x)2+(∂U/∂y)2$⁠, and temperature gradient $dT=(∂T/∂x)2+(∂T/∂y)2$⁠. Figure 17 shows the results of sensitivity analysis, and it shows that two important variables are found to be U and T dominated in the low velocity region and only one key variable U dominates in the high velocity region. However, the influence of second sensitivity variable $dT$ in the high velocity region is also considered in the study.

Hereafter, for the purpose of distinguishing, the $Prt$ function value and the $Prt$ data assimilation results are denoted as $Prt_Eq$ and $Prt_DA$⁠, respectively.

Figure 18 plots the variations of $Prt_Eq$ and $Prt_DA$ with respect to U, and the coefficient of determination $R2$ has a high value of 0.93 as shown in Fig. 19; it seems that the $Prt$ function could roughly reproduce the data assimilation results. However, an issue is that the relative error is large in the small value region, it is conceivable that the amount of experimental data was insufficient, or that the current proposed $Prt$ function is still not the optimal forms. Regarding improvement of the $Prt$ function, there are two ways which are increasing experimental data under small velocity conditions and establishing more appropriate functional forms using optimization analysis. Anyway, if enough experimental data can be collected in the future, the machine learning model might be a final way instead of the $Prt$ function.

In this study, data assimilation was performed based on the experiments that simulate the mixing field by hot and cold air inside a simple HVAC model, and the turbulent Prandtl number distribution that fits the experimental data were predicted. It was confirmed that the accuracy of the temperature field prediction is improved under all experimental conditions compared to the results using the conventional turbulent Prandtl number function by Ito et al. [4], as the average deviation of temperature between the experiment and simulation reduced by 2.2 °C and was less than 3.5 °C. It was demonstrated that the data assimilation is most effective for CFD simulation to simulate the temperature field consistent with the actual measurements. Consequently, the data assimilation is promising to enhance high fidelity as well as low cost in real\minus world product design and development. In addition, in order to apply the turbulent Prandtl number determined by data assimilation to CFD simulations of actual HVAC unit, two new turbulent Prandtl number functions are proposed, which are constructed from a power function of the dimensionless velocity $U*$ and an exponential function of the dimensionless temperature $T*$ or temperature gradient $dT*$⁠, such as $0.016×U*−0.187×e1.498×T*$ for the low velocity region and $0.040×U*0.274×e0.006×dT*$ for the high velocity region.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

d =

distance between thermocouples

h =

heat transfer coefficient in the Kramers equation

w =

observation error vector of data assimilation

$J$ =

cost function of data assimilation

$T¯$ =

Reynolds averaged temperature

U =

magnitude of velocity

H =

mapping matrix in data assimilation

I =

identity matrix of data assimilation

$cw$ =

heat capacity of wire in the Kramers equation

$kf$ =

thermal conductivity of fluid

$rw$ =

radius of wire in the Kramers equation

$uc$ =

certain value of velocity

$u¯j$ =

Reynolds averaged velocity component

$x^n+1|n$ =

ensemble state vector of data assimilation

$yactual$ =

vector of actual measurement data

N_flow =

mesh nodes in the CFD model

N_ensemble =

total ensemble member

$Nl$ =

number of measurement points in the low velocity region

$Nh$ =

number of measurement points in the high velocity region

$Tf$ =

temperature of fluid in the Kramers equation

$Th$ =

high temperature at inlet of a simple HVAC model

$Tl$ =

low temperature at inlet of a simple HVAC model

$Tw$ =

temperature of sensor in the Kramers equation

$Uh$ =

high velocity at inlet of a simple HVAC model

$Ul$ =

low velocity at inlet of a simple HVAC model

$Kn+1$ =

Kalman gain of data assimilation

$V^n+1|n$ =

ensemble covariance matrix of data assimilation

$T*$ =

dimensionless temperature of the $Prt$ functions

$U*$ =

dimensionless velocity of $Prt$ functions

$xn|n(i)$ =

state vector of data assimilation

$yn+1(i)$ =

observation vector of data assimilation

$alorh,alorh*,blorh,clorh$ =

coefficients of $Prt$ functions

$dT$ =

temperature gradient

$dT*$ =

dimensionless gradient temperature of $Prt$ functions

$dU$ =

velocity gradient

$u′v′¯$ =

Reynolds shear component (turbulent shear stress)

$v′T′¯$ =

turbulent heat flux

$v*′T′¯$ =

turbulent heat flux at a local coordinate system

$(xi,yi)$ =

Cartesian coordinate system

$(xi*,yi*)$ =

local coordinate system

$Pr$ =

molecular Prandtl number

$Prt$ =

turbulent Prandtl number

$Prt_DA$ =

$Prt$ data assimilation results

$Prt_Eq$ =

$Prt$ function value

$Prf$ =

Prandlt number in the Kramers equation

$Re$ =

Reynolds number

$Ref$ =

Reynolds number in the Kramers equation

ε_h =

heat transfer eddy diffusivity

ε_m =

momentum eddy diffusivity

λ =

thermal conductivity

$ρw$ =

linear density of wire in the Kramers equation

$τ$ =

thermal time constant in the Kramers equation

$υt$ =

eddy viscosity

$ϕ$ =

optimal square difference

$ωz$ =

vorticity

cal =

calculated value

exp =

measured value

n|n =

current step|next step of assimilation

(i) =

the i-th ensemble member