A Parametric Numerical Study of the Airﬂow and Thermal Performance in a Real Data Center for Improving Sustainability

: In recent years, reducing energy consumption has been relentlessly pursued by researchers and policy makers with the purpose of achieving a more sustainable future. The demand for data storage in data centers has been steadily increasing, leading to an increase in size and therefore to consume more energy. Consequently, the reduction of the energy consumption of data center rooms is required and it is with this perspective that this paper is proposed. By using Computational Fluid Dynamics (CFD), it is possible to model a three-dimensional model of the heat transfer and air ﬂow in data centers, which allows forecasting the air speed and temperature range under diverse conditions of operation. In this paper, a CFD study of the thermal performance and airﬂow in a real data center processing room with 208 racks under di ﬀ erent thermal loads and airﬂow velocities is proposed. The physical-mathematical model relies on the equations of mass, momentum and energy conservation. The ﬂuid in this study is air and it is modeled as an ideal gas with constant properties. The model of the e ﬀ ect of turbulence is made by employing a k– ε standard model. The results indicate that it is possible to reduce the thermal load of the server racks by improving the thermal performance and airﬂow of the data center room, without a ﬀ ecting the correct operation of the server racks located in the sensible regions of the room. results accuracy with the computational cost and simulation time. Three key parameters were analyzed to measure the mesh quality: aspect ratio, skewness and orthogonal quality. The aspect ratio is a measure of the stretching of a cell. It is computed as the ratio of the largest and smallest dimensions of the edges of the faces of the control volume. A good aspect ratio has a value close to one. Skewness is deﬁned as the di ﬀ erence between the shape of the cell and the shape of an equilateral cell of equivalent volume.


Introduction
In the last few decades, the increase in energy consumption due to a rising utilization of information and communication technologies (ICT) is a difficult challenge, from the energetic, ecological and economical points of view, which compromises sustainability [1]. The creation and rapid dissemination of grid and cloud computing [2], in recent years, has led to an increase in the number of data centers as well as in their size [3].
A data center is a physical site that comprises ICT installations with the main function being the storage and distribution of data through Internet access or an internal network. Data centers throughout

Methodology
In this section, the mathematical formulations describing the air flow, heat transfer and the model of turbulence used for this study are presented. When analyzing any air flow problem, it is necessary to establish the conservation laws that are represented mathematically and physically by the governing equations.
The governing equations are composed of the state equation, continuity equation, momentum equation, and energy equation. In this paper, a three-dimensional case is considered, where the working fluid is assumed to be an ideal gas. The air flow regime is assumed to be turbulent and the heat transfer process is assumed to be steady.
Appl. Sci. 2019, 9, 3850 4 of 30 The state equation is a constitutive relation of a fluid which relates several properties. In engineering, the substances of interest are mostly gases, with moderate pressures and temperatures. In this paper, as mentioned above, an ideal gas is used, and the equation of state is given as follows [37]: in which ρ is the specific mass and is given in kg·m −3 , the pressure is given by p and it is measured in Pa, M m represents the molar mass of the gas and it is measured in kg·kmol −1 , the ideal gas constant is represented by R u and it is measured in J·mol −1 ·K −1 , and the temperature T is given in K. Modeling the air as an ideal gas allows taking into account the buoyancy effect.
The continuity equation, also referred to as the mass conservation equation, expresses that the mass entering a control volume is equal to the mass exiting that control volume. Consequently, no losses occur in this process. The continuity equation is expressed as follows [3]: where t is the time in seconds, x, y and z are the Cartesian coordinates and u, v and w are the corresponding velocity components in m·s −1 .
The equation of momentum can be obtained through Newton's second law, which asserts that the outcomes of the forces applied to a particle is equal to the rate of change of its linear momentum. These equations are presented as follows [37]: where µ is the dynamic viscosity in kg·m −1 ·s −1 and g is the gravitational acceleration in m·s −2 . The gravitational force term is included in the z momentum equation. The energy equation is the mathematical expression based on the First Law of Thermodynamics which asserts that the degree of change of the internal energy of a system is equal to the difference between the heat exchange with the external medium and the work done by it, whose expression is given by the following equation [38]: where C a represents the specific heat in W·kg −1 ·K −1 , λ represents the thermal conductivity in W·m −1 ·K −1 , and S T represents the source term of heat generation, measured in W·m −3 .

Turbulence Model
According to Almoli et. al [28], the thermal air flows in data centers are generally complex, with a Reynolds number, Re, based on the air velocity at the inlet air vent openings at around 1 m·s −1 with server racks with a height of around 2.4 m, a Re = 10 5 is obtained, which indicates that the air flow regime is turbulent. When the studies performed in CFD in data centers are analyzed, it can be observed that the model that provides more accurate forecasts makes use of the Reynolds Averaged Navier-Stokes (RANS) turbulence model, as was the case in [39].
In this study, the turbulence k-ε standard model is based on the conservation equations of continuity, momentum and energy in conjunction with the turbulence k-ε standard model, where k and ε are the variables kinetic turbulent energy and turbulent dissipation, respectively. According to the authors of [39], this model is more suited to turbulent flow due to the way it calculates turbulent viscosity and conductivity. It is also the most widely validated and frequently used in commercial CFD codes. In addition, according to Alkharabsheh et al. [40], previous studies have shown that the k-ε turbulent model displays a better performance and results when compared to the SST, k-ω, RSM and RNG k-ε models.
The mathematical expressions that characterize the model proposed and used in this study are based on the ones presented in [41] and are given as follows: where k is kinetic turbulent energy m 2 ·s −2 , U represents the velocity vector in m·s −1 , µ T represents the turbulent viscosity in kg·m −1 ·s −1 , σ k gives the turbulent Prandtl number for turbulent kinetic energy, P t characterizes the turbulent energy production in kg·m −1 ·s −3 , ε gives the turbulent kinetic energy dissipation in m −2 ·s −3 , the Prandtl number for turbulent energy dissipation is given by σ ε , and finally C 1ε and C 2ε are the constants of the production term and dissipation term, respectively. The abovementioned constants σ k , σ ε , C 1ε , and C 2ε assume the values of 1, 1.3, 1.44, and 1.92, respectively, as in the work of Hoang et al. [41]. Finally, the turbulent viscosity is defined by Equation (9) [41]: where C µ is the constant of the turbulent viscosity term. In this study, C µ = 0.09.

Numerical Model of the Data Center
The geometric model was developed taking into account the characteristics of the computational processing in order to avoid complications during the developed simulations. Thus, the threedimensional geometry was simplified to avoid that the mesh generated by the software being very refined. To achieve this, the server racks are represented as parallelepipeds, and the cables, rails, access doors and fixtures have been concealed to facilitate the construction of the mesh. In the case of lamps, these were concealed not only by the complexity they would generate in the developed mesh, but also because they were disconnected during the operation of the ICT room, generating no heat production. Figure 1 shows the three-dimensional geometry developed for the present study.
The computational mesh is composed of control volumes so that the equations of continuity, momentum and energy, presented in Section 2, could be solved. A mesh that has high quality is very important for the accomplishment of the simulations proposed in this paper, because this is fundamental for the numerical solution to converge and so that the achieved results are consistent and trustworthy. The resulting mesh consists of tetrahedral (cold corridors) and hexahedral (hot corridors and server racks) elements. This solution was adopted because the mesh would have a high number of control volumes if a mesh composed of only one type of elements were used. Additionally, as the air flow changes from horizontal direction when crossing the racks to vertical direction when going through the air outlet vent areas, the use of tetrahedral elements is required to reduce false mathematical diffusion and thus improve the accuracy of the numerical predictions. Figure 2 shows the mesh resulting from the study in an isometric perspective. Table 1 presents the characteristics of the computational mesh. dimensional geometry was simplified to avoid that the mesh generated by the software being very refined. To achieve this, the server racks are represented as parallelepipeds, and the cables, rails, access doors and fixtures have been concealed to facilitate the construction of the mesh. In the case of lamps, these were concealed not only by the complexity they would generate in the developed mesh, but also because they were disconnected during the operation of the ICT room, generating no heat production. Figure 1 shows the three-dimensional geometry developed for the present study.   The computational mesh is composed of control volumes so that the equations of continuity, momentum and energy, presented in Section 2, could be solved. A mesh that has high quality is very important for the accomplishment of the simulations proposed in this paper, because this is fundamental for the numerical solution to converge and so that the achieved results are consistent and trustworthy. The resulting mesh consists of tetrahedral (cold corridors) and hexahedral (hot corridors and server racks) elements. This solution was adopted because the mesh would have a high number of control volumes if a mesh composed of only one type of elements were used. Additionally, as the air flow changes from horizontal direction when crossing the racks to vertical direction when going through the air outlet vent areas, the use of tetrahedral elements is required to reduce false mathematical diffusion and thus improve the accuracy of the numerical predictions. Figure 2 shows the mesh resulting from the study in an isometric perspective. Table 1 presents the characteristics of the computational mesh.  The mesh density is an important metric and it is utilized for controlling accuracy and a highdensity mesh delivers results with high accuracy. The meshes creation was an iterative process as its quality plays a significant role in the accuracy and stability of the numerical computation. The attributes associated with mesh quality are node point distribution, smoothness, and skewness. The quality of the mesh elements was assessed during the meshes creation in order to balance the results accuracy with the computational cost and simulation time. Three key parameters were analyzed to measure the mesh quality: aspect ratio, skewness and orthogonal quality. The aspect ratio is a measure of the stretching of a cell. It is computed as the ratio of the largest and smallest dimensions of the edges of the faces of the control volume. A good aspect ratio has a value close to one. Skewness is defined as the difference between the shape of the cell and the shape of an equilateral cell of equivalent volume. The mesh quality improves as skewness value approaches the null value. The orthogonal quality for cells is computed using the face normal vector, the vector from the cell centroid to the centroid of each of the adjacent cells, and the vector from the cell centroid to each of the faces. The orthogonal quality for a cell is computed as the minimum of these quantities. The mesh quality increases as the value of orthogonal quality is towards unity. Several computational meshes were  The mesh density is an important metric and it is utilized for controlling accuracy and a high-density mesh delivers results with high accuracy. The meshes creation was an iterative process as its quality plays a significant role in the accuracy and stability of the numerical computation. The attributes associated with mesh quality are node point distribution, smoothness, and skewness. The quality of the mesh elements was assessed during the meshes creation in order to balance the results accuracy with the computational cost and simulation time. Three key parameters were analyzed to measure the mesh quality: aspect ratio, skewness and orthogonal quality. The aspect ratio is a measure of the stretching of a cell. It is computed as the ratio of the largest and smallest dimensions of the edges of the faces of the control volume. A good aspect ratio has a value close to one. Skewness is defined as the difference between the shape of the cell and the shape of an equilateral cell of equivalent volume.
The mesh quality improves as skewness value approaches the null value. The orthogonal quality for cells is computed using the face normal vector, the vector from the cell centroid to the centroid of each of the adjacent cells, and the vector from the cell centroid to each of the faces. The orthogonal quality for a cell is computed as the minimum of these quantities. The mesh quality increases as the value of orthogonal quality is towards unity. Several computational meshes were created until the size of mesh elements provided a good aspect ratio, a small skewness and a good orthogonal quality, ensuring a good relation between the computational cost and the accuracy of the results. Although there are perceptible steps in the predictions of temperature and velocity fields, it is necessary to highlight that the selected meshes present a reasonable solution considering the results accuracy and the computational time required. Additionally, it must be noted that the numerical results simulate a real physical situation and their analysis and conclusions are in accordance with the computational resources required and computational time suitable for an engineering application.
The equations described in Section 2 were solved by a numerical process where the control volume technique was used. This technique is used when the differential equations have no analytical resolution and the numerical model must allow the discretization of the equations. To achieve a possible resolution of the numerical method, the computational domain was first divided into 2,341,172 control volumes. As stated above, a high-density mesh delivers results with high accuracy.
The second-order UpWind Scheme was applied to momentum, turbulent kinetic energy, turbulent dissipation rate, and energy to achieve more precise results. The accuracy of this method is reached on the faces of the control volumes by the expansion of the Taylor Series of the solution centered on the control volume compared to its centroid [42].
The PREssure STaggering Option (PRESTO!) interpolation method was used for pressure discretization. Its use is justified for controlling the sudden pressure variations between the centroids of the control volumes, considerable body forces, the presence of flows involving porous medium and other situations. It is based in the adoption of a staggered grid to calculate the pressure via a discrete continuity balance. In the staggered grid, the values of pressure in the center are known and these are the values at the interfaces in the normal grid. The PRESTO! scheme uses this technique, which resembles in characteristics the staggered-grid schemes used with structured meshes proposed by Patankar [42]. The Semi-Implicit Method for Pressure Linked Equations-Consistent (SIMPLEC) algorithm, published by Van Doormal and Raithby [43], was chosen in detriment to the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm for the pressure-velocity coupling [42], by allowing the convergence of the solution more quickly and ensuring the stability of the method.

Boundary Conditions
The boundary conditions defined in the numerical models were based on the data collected in the studied real data center. To the values collected was also added the value of pressure loss in the server racks, which are provided in [44].
It is important to point out that, from now on all, coordinate values and parameters are presented with dimensionless quantities to simplify the comparison between case studies. Dimensionless values are identified with a subscript "adm" in all equations of this study. The dimensionless quantities range from 0 to 1. It should be noted that, in the case of the temperature, the maximum value corresponds to a temperature value exceeding the operating specifications of the server racks. It is an air temperature value with the capacity to cause malfunctions or to damage to the equipment present in the ICT rooms, while this equipment is intended to operate uninterruptedly.
In this study, as previously mentioned, the PRESTO! interpolation method was used mainly due to the presence of porous medium [45]. The porous medium is an agglomerate of particles that are close, leaving small voids between them. Due to the existence of these voids, it is allowed to circulate liquid or gaseous elements that can fill them. Thus, the larger are these voids, the greater is the porosity of this material, that is, a greater amount of another element is required to fill them. In this study, only porosity was considered for the server racks, due to micro perforated doors that allow air to be introduced into the server rack. Besides the porosity present in the doors of the server rack, there are still the components inside that occupy the empty space inside the server rack, thus not allowing free air circulation. Therefore, for this study, the porosity of the server rack was considered as a whole (carcass and elements) and with the value of 80%, provided by the data center.
Due to the existence of porous media, it is necessary to calculate the coefficient of inertial loss "C 2 " applied to the server racks [45], which is given by the following expression: where K L' is the loss factor, ∆p is the loss variation in Pa, and v is the velocity, in m·s −1 . The coefficient of inertial loss, C 2 is expressed as follows: The coefficient of inertial loss (C 2 ) is measured in m −1 . For the resolution of Equation (10), the pressure variation value of 7.5 Pa was adjusted to the velocity of m·s −1 at the velocities present in [44].
The equations presented in Section 2 are solved in the fluid zone. The fluid zone was considered having an initial dimensionless temperature of 7.09 × 10 −4 , that is, the average dimensionless temperature measured within the ICT room.
The turbulence intensity and the hydraulic diameter are the parameters that allow the iterative calculation of the turbulent kinetic energy and the coefficient of dissipation of the turbulent kinetic energy. Regarding the intensity of turbulence, it indicates the direction of the flow. Generally it is considered 1% for low turbulence intensity and for high turbulence intensity it is considered 10% [45]. In this study, a turbulence intensity of 5% was considered in the extraction zone, as it is located in the containment zone of the hot aisles, which provides a well-defined direction of air flow. Concerning the intensity of turbulence in the air vent area, 5% was also considered since the air vent grids allow directing the air. These values were based on the experimental and numerical results provided in [46][47][48][49].
The boundary conditions are presented in Table 2 and are characterized by dimensionless values. The area of air vents presents values that are always equal, while in the extraction area different values occur, since one of the elements that compose this area is smaller. Table 2. Boundary conditions imposed in the air vent and extraction areas.

All Air Vents Extraction Smaller Extraction
Turbulence intensity

Solution Convergence
The convergence of the solution is obtained through the relaxation of the variables, decreasing the abrupt variations during the iteration. It is considered that a solution converges when all the residuals reach the stopping criterion, λ = 1 × 10 −4 for all variables, except for the energy variable that assumes a stop criterion of λ = 1 × 10 −6 . In order for a solution to be considered convergent. there are several aspects to consider, such as refinement of the computational mesh, established numerical methods, relaxation factors and the number of imposed iterations. There are cases where it is not possible to achieve the convergence established by the stopping criterion, because the residuals stabilize before this criterion, thus assuming that the solution obtained is the possible convergence for the numerical method under analysis. The relaxation factors have no influence on the numerical solution' they only change the time of the iterative process. Table 3 presents the various relaxation factors applied to the numerical methods, following the data provided in [46][47][48][49]. The iterative process of the various numerical models was defined for 1500 iterations, since the solution would not reach the stopping criteria, only if it could reach the possible convergence. The authors ran few scenarios with 7500 iterations, and no divergence was verified. However, due to the reason that these simulations took a high amount of time to compute, the authors decided to run the simulation up to 1500 iterations, since the variable residuals accomplished the possible convergence. Thus, this number of iterations was established to reduce the processing time of the simulation. Each case took about 14 h to complete all the iterations on a work station with two 3.47 GHz processors with 16 cores and with a memory of 48 GB of Random-access memory (RAM).

Description of the Modeled Case Studies
In this section, we present the parametric studies performed by changing different parameters that allowed the analysis of a good operation of an ICT room. For this analysis, the temperatures of the server racks and the hot aisles were studied, as well as the airflow in the ICT room. The investigated parametric studies differed among them, due to the variation of the thermal load of the server racks. Four case studies with different thermal loads were studied. Cases 1 and 2 represent a usual thermal load distribution in the data center, especially Case 2. This objective was to have a study that closely resembles the reality of the operation of the server racks located in the data center. Meanwhile, Cases 3 and 4 represent a more academic approach. The thermal is not usually constant throughout the row since the servers do not usually operate at the same level of thermal load. However, the values of the thermal load in Cases 3 and 4 are close to the reality since these servers do not usually operate at maximum performance.
Each thermal load was simulated with the maximum and minimum air flow velocity. The dimensionless values are identified with a subscript "adm". Figure 3 shows the upper plan view of the TI room, where the numbering corresponds to the server rack number, ranging from 1 to 208. In this figure, it is also possible to observe the hot aisle/cold aisle configuration of the data center. The smaller extraction aisle, mentioned in Table 3 and seen in Figure 3, is the one of the ranging racks from 181 until 208. In all four case studies addressed in this study, a dimensionless thermal load was used, qadm. In Case 1, different thermal loads imposed at the boundary conditions of the heat flow were applied symmetrically to every row of all server racks, as shown in Figure 4a. The thermal load is as follows: at the end of the rows of the server racks, 90% of the maximum thermal load was applied, while 70% of the maximum thermal load and 50% of the central racks were imposed on the intermediate server racks. For each case, the maximum and minimum air flow velocity was applied, simulated and studied. The purpose was to observe the behavior of the temperature and air flow in the room. To better understand how the loads were applied along the rows, a single row of server racks is shown in Figure 4 for each of the four case studies, where each color represents a value of the thermal load.  In all four case studies addressed in this study, a dimensionless thermal load was used, q adm . In Case 1, different thermal loads imposed at the boundary conditions of the heat flow were applied symmetrically to every row of all server racks, as shown in Figure 4a. The thermal load is as follows: at the end of the rows of the server racks, 90% of the maximum thermal load was applied, while 70% of the maximum thermal load and 50% of the central racks were imposed on the intermediate server racks. For each case, the maximum and minimum air flow velocity was applied, simulated and studied. The purpose was to observe the behavior of the temperature and air flow in the room. To better understand how the loads were applied along the rows, a single row of server racks is shown in Figure 4 for each of the four case studies, where each color represents a value of the thermal load. In all four case studies addressed in this study, a dimensionless thermal load was used, qadm. In Case 1, different thermal loads imposed at the boundary conditions of the heat flow were applied symmetrically to every row of all server racks, as shown in Figure 4a. The thermal load is as follows: at the end of the rows of the server racks, 90% of the maximum thermal load was applied, while 70% of the maximum thermal load and 50% of the central racks were imposed on the intermediate server racks. For each case, the maximum and minimum air flow velocity was applied, simulated and studied. The purpose was to observe the behavior of the temperature and air flow in the room. To better understand how the loads were applied along the rows, a single row of server racks is shown in Figure 4 for each of the four case studies, where each color represents a value of the thermal load.  Case 2 is similar to Case 1. However, in this case study, the thermal load applied on the end of the server racks is exchanged with the one of the central server racks, that is, the central server racks had a maximum heat load of 90% and the end server racks dissipated 50% of the maximum heat load. For Case 3, the thermal load applied was 75% of the maximum thermal load in the heat flow frontier condition, q. It was applied equally in every server rack, for the maximum and minimum velocity of the air flow. Finally, Case 4 is similar to Case 3, the difference being only the applied thermal load, which in this case was equal to the measured and nominal load value of the data center room.

Case 1
In Case 1, several simulations were made with different thermal loads that were applied in the boundary conditions of the heat flow. The loads imposed were applied symmetrically along the rows, as shown in Figure 4a. The boundary conditions taken in account for these simulations are as follows: at the end of the rows of the server racks, 90% of the maximum thermal load was applied, while 70% of the maximum thermal load and 50% of the central racks were applied: q adm = 0.892 (90%), 0.676 (70%) and 0.460 (50%), respectively. The minimum air flow dimensionless velocity was set to be as follows: v adm = 0. By observing Figures 5 and 6, the temperature field predictions of the y and z planes can be seen, respectively. These temperature fields predict the set of temperature values at all points in a given space at a given instant. For these specific cases, in steady-state conditions, the temperature field is independent of time. Case 2 is similar to Case 1. However, in this case study, the thermal load applied on the end of the server racks is exchanged with the one of the central server racks, that is, the central server racks had a maximum heat load of 90% and the end server racks dissipated 50% of the maximum heat load. For Case 3, the thermal load applied was 75% of the maximum thermal load in the heat flow frontier condition, q. It was applied equally in every server rack, for the maximum and minimum velocity of the air flow. Finally, Case 4 is similar to Case 3, the difference being only the applied thermal load, which in this case was equal to the measured and nominal load value of the data center room.

Case 1
In Case 1, several simulations were made with different thermal loads that were applied in the boundary conditions of the heat flow. The loads imposed were applied symmetrically along the rows, as shown in Figure 4a. The boundary conditions taken in account for these simulations are as follows: at the end of the rows of the server racks, 90% of the maximum thermal load was applied, while 70% of the maximum thermal load and 50% of the central racks were applied: qadm = 0.892 (90%), 0.676 (70%) and 0.460 (50%), respectively. The minimum air flow dimensionless velocity was set to be as follows: vadm = 0. By observing Figures 5 and 6, the temperature field predictions of the y and z planes can be seen, respectively. These temperature fields predict the set of temperature values at all points in a given space at a given instant. For these specific cases, in steady-state conditions, the temperature field is independent of time.   In the case of the maximum air flow dimensionless velocity, vadm = 1, the temperature field predictions can be seen in Figures 7 and 8 for the y and z planes, respectively.  In the case of the maximum air flow dimensionless velocity, v adm = 1, the temperature field predictions can be seen in Figures 7 and 8 for the y and z planes, respectively. In the case of the maximum air flow dimensionless velocity, vadm = 1, the temperature field predictions can be seen in Figures 7 and 8 for the y and z planes, respectively.   More than a few hot spots are foreseen to arise at the termination of every server rack line, as illustrated in Figures 5 and 6. These spots are easily spotted if noticing the colors related to the values which are higher than the dimensionless temperature, Tadm = 0.350, visible on the left scale, except for the hot spots in the extraction zone, which is Tadm = 0.125. Through a careful analysis of Figure 6a, a greater heat intensity is foreseen to occur at the bottom the racks located at the end ow every row. It can be deduced that the airflow is inadequate in this area. The justification behind this inadequacy is that the airflow intensity is higher in the center of the air vents than at the bottom and the top. Figure  6a-c also shows that for this case study hot spots can be predicted in the extraction zone and at the termination of all the rows.
In Figure 7a,c, which have server racks at the end of the rows (yadm = 0.100 and yadm = 0.900), the temperatures are expected to be very close to the hot spot limit, with Rack 2 exceeding that limit at the base. As far as the plan representing the central server racks is concerned (Figure 7b), it is expected that the temperature does not vary much from the inlet air temperature.
By observing Figure 8, it can be seen that there is a possibility of hot spots in a few server racks, such as racks located at the extremity of each row. As far as the extraction zone is concerned, the probability of hot spots is very high, since the mean temperatures shown in Figure 8d are very close to the inlet air temperature.

Case 2
In Case 2, the applied loads were applied symmetrically along the rows, as shown in Figure 4b. The boundary conditions set for the simulations performed in this studies are as follows: at the end of the rows of the server racks, 50% of the maximum thermal load was applied, while 70% of the maximum thermal load and 90% of the central racks were applied: qadm = 0.460 (50%), 0.676 (70%) and More than a few hot spots are foreseen to arise at the termination of every server rack line, as illustrated in Figures 5 and 6. These spots are easily spotted if noticing the colors related to the values which are higher than the dimensionless temperature, T adm = 0.350, visible on the left scale, except for the hot spots in the extraction zone, which is T adm = 0.125. Through a careful analysis of Figure 6a, a greater heat intensity is foreseen to occur at the bottom the racks located at the end ow every row. It can be deduced that the airflow is inadequate in this area. The justification behind this inadequacy is that the airflow intensity is higher in the center of the air vents than at the bottom and the top. Figure 6a-c also shows that for this case study hot spots can be predicted in the extraction zone and at the termination of all the rows.
In Figure 7a,c, which have server racks at the end of the rows (y adm = 0.100 and y adm = 0.900), the temperatures are expected to be very close to the hot spot limit, with Rack 2 exceeding that limit at the base. As far as the plan representing the central server racks is concerned (Figure 7b), it is expected that the temperature does not vary much from the inlet air temperature.
By observing Figure 8, it can be seen that there is a possibility of hot spots in a few server racks, such as racks located at the extremity of each row. As far as the extraction zone is concerned, the probability of hot spots is very high, since the mean temperatures shown in Figure 8d are very close to the inlet air temperature.

Case 2
In Case 2, the applied loads were applied symmetrically along the rows, as shown in Figure 4b. The boundary conditions set for the simulations performed in this studies are as follows: at the end of the rows of the server racks, 50% of the maximum thermal load was applied, while 70% of the maximum thermal load and 90% of the central racks were applied: q adm = 0.460 (50%), 0.676 (70%) and 0.892 (90%), respectively. The minimum air flow dimensionless velocity was set as follows: v adm = 0. The results from the simulations of the temperature field can be seen, respectively, in Figures 9 and 10 for the y and z planes. Figures 9 and 10 show the temperature field predictions for the case study. The temperature fields are presented with the same planes as in Case 1.
By observing Figure 9, the probability of hot spots on the server racks can be accurately assessed. In the plans of Figure 9, it can also be easily noticed that the average of the non-dimensional temperature values in the extraction zone exceeds T adm = 0.125. Therefore, the minimum air flow dimensionless velocity is not sufficient to successfully cool all the racks and, thus, avoid the occurrence of any type of hotspots.
By analyzing Figure 10, it can be seen that a few server racks can reach temperatures that exceed the minimum limit considered for a hot spot. In the extraction zone, the mean value of the dimensionless temperature in this zone is expected to exceed T adm = 0.125 (the hot spot for this area), which means that hot spots in the extraction zone are very likely to occur.
In the case of the maximum level of air flow dimensionless velocity, v adm = 1, the results from the simulations of the temperature field can be seen, respectively, in Figures 11 and 12 for the y and z planes.    By observing Figure 9, the probability of hot spots on the server racks can be accurately assessed. In the plans of Figure 9, it can also be easily noticed that the average of the non-dimensional temperature values in the extraction zone exceeds Tadm = 0.125. Therefore, the minimum air flow dimensionless velocity is not sufficient to successfully cool all the racks and, thus, avoid the occurrence of any type of hotspots.
By analyzing Figure 10, it can be seen that a few server racks can reach temperatures that exceed the minimum limit considered for a hot spot. In the extraction zone, the mean value of the dimensionless temperature in this zone is expected to exceed Tadm = 0.125 (the hot spot for this area), which means that hot spots in the extraction zone are very likely to occur.
In the case of the maximum level of air flow dimensionless velocity, vadm = 1, the results from the simulations of the temperature field can be seen, respectively, in Figures 11 and 12 for the y and z planes.     As shown in Figure 11b, the predicted dimensionless temperature does not exceed T adm = 0.150. The zone that can reach this temperature is very small and very close to the base. In the other planes, the dimensionless temperature does not exceed T adm = 0.200, which is also only reached at or near the base of the server racks.
Due to the increase in inlet air flow velocity, the temperature fields shown in Figure 12 predict that the average temperature in the data center room is similar to the inlet air temperature and that the formation of hot spots is almost non-existent, since the predicted maximum dimensionless temperature is T adm = 0.250. As mentioned for Case 1, the hot spots can be easily spotted by noticing the colors related to the values, which are higher than the dimensionless temperature, T adm = 0.350. In addition, no hot spots are expected in the extraction zone.

Case 3
For Case 3, the conditions of the boundary were set as follows: 75% of the maximum heat load in the heat flow boundary condition was imposed symmetrically, meaning that q adm = 0.729. As in previous cases, the minimum and maximum air flow velocity were applied. In the case of the minimum air flow velocity, the temperature field predictions can be seen, respectively, in Figures 13 and 14 for the y and z planes.

Case 3
For Case 3, the conditions of the boundary were set as follows: 75% of the maximum heat load in the heat flow boundary condition was imposed symmetrically, meaning that qadm = 0.729. As in previous cases, the minimum and maximum air flow velocity were applied. In the case of the minimum air flow velocity, the temperature field predictions can be seen, respectively, in Figures 13  and 14 for the y and z planes.  In Figure 13, it is possible to observe the temperature of the server racks above the established value for hot spots and average dimensionless temperatures in the hot aisles of Tadm = 0.350. Several hot spots occur with the minimum air flow velocity. It is possible to observe in Figure 14 that, for a minimum air flow velocity, when 75% is set for the maximum load, several hot spots are foreseen, as In Figure 13, it is possible to observe the temperature of the server racks above the established value for hot spots and average dimensionless temperatures in the hot aisles of T adm = 0.350. Several hot spots occur with the minimum air flow velocity. It is possible to observe in Figure 14 that, for a minimum air flow velocity, when 75% is set for the maximum load, several hot spots are foreseen, as confirmed by Figure 13, in the server racks at the end of all the rows. Regarding the extraction temperature, an average dimensionless temperature around T adm = 0.200 is expected, and, since it is higher than T adm = 0.125, it leads to the prediction of hot spots.
When the air flow velocity is simulated at the maximum level, the temperature field predictions show a different result, as can be seen, respectively, in Figures 15 and 16 for the y and z planes.
By observing Figure 15, it is possible to assess that some of the racks in the data center room exhibit temperature values similar to the inlet air temperature, especially in the intermediate server racks of each row.
As shown in the plane temperature fields of Figure 16, when applying the maximum level of air flow velocity, the thermal behavior of the data center room does not provide hot spots, neither for the server racks nor for the extraction zone, which has an average temperature similar to the inlet air temperature. Both figures confirm that no hot spots occur in this scenario. By observing Figure 15, it is possible to assess that some of the racks in the data center room exhibit temperature values similar to the inlet air temperature, especially in the intermediate server racks of each row.  By observing Figure 15, it is possible to assess that some of the racks in the data center room exhibit temperature values similar to the inlet air temperature, especially in the intermediate server racks of each row.

Case 4
In Case 4, the nominal thermal load in the heat flow frontier condition, as mentioned in Section 4, and provided by the data center room, is q adm = 0.340. As in previous cases, the maximum and minimum air flow velocity were applied. The scale was changed for Case 4 and it does not present the maximum value existing in the other three cases. The reason is that this case study has lower temperatures and small differences would not be perceptible. By focusing on a narrower scale, it is possible to differentiate between slight changes in temperature.
The result of simulations for the minimum level of air flow velocity and the obtained temperature field predictions can be seen, respectively, in Figures 17 and 18 for the y and z planes.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 20 of 30 As shown in the plane temperature fields of Figure 16, when applying the maximum level of air flow velocity, the thermal behavior of the data center room does not provide hot spots, neither for the server racks nor for the extraction zone, which has an average temperature similar to the inlet air temperature. Both figures confirm that no hot spots occur in this scenario.

Case 4
In Case 4, the nominal thermal load in the heat flow frontier condition, as mentioned in Section 4, and provided by the data center room, is qadm = 0.340. As in previous cases, the maximum and minimum air flow velocity were applied. The scale was changed for Case 4 and it does not present the maximum value existing in the other three cases. The reason is that this case study has lower temperatures and small differences would not be perceptible. By focusing on a narrower scale, it is possible to differentiate between slight changes in temperature.
The result of simulations for the minimum level of air flow velocity and the obtained temperature field predictions can be seen, respectively, in Figures 17 and 18 for the y and z planes.     In the planes shown in Figure 17, it can be observed that the greater development of hot spots occurs at the terminations of the rows of the server racks. The temperature fields of the x-y planes ( Figure 18) represent the temperature prediction for a usual thermal load with minimum air flow velocity. Due to the normal operation of the server racks, the air flow velocity is insufficient to avoid the formation of hot spots in several server racks. In the extraction zone, hot spots are not expected since the average dimensionless temperature does not exceed Tadm = 0.125.
By analyzing both Figures 17 and 18, it can be assessed that the server racks at the end of the rows reach higher temperatures, meeting the dimensionless temperature value of Tadm = 0.400. The hot aisles are expected to reach the maximum average mean temperatures, Tadm = 0.120.
The result of simulations for the maximum level of air flow velocity and the obtained temperature field predictions can be seen, respectively, in Figures 19 and 20 for the y and z planes.
By observing Figure 19, it is possible to assess that the values obtained in the dimensionless temperature range of the server racks do not exceed the value, Tadm = 0.160. In the extraction zone, the simulations show that hot spots are not foreseen to occur, since the value of the dimensionless temperature field does not exceed the value of Tadm = 0.06.  Figure 20 shows the absence of hot spots, since the inlet air flow velocity proves to be oversized for this case study, which naturally leads to the prediction of low values in the extraction zone. It can be seen that server racks at the end of the rows are very close to the values of the central server racks, since the scale does not present the maximum value existing in the other three cases. Thus, it can be seen that there are no hot spots in any server rack. The value of the temperature range in the extraction zone is almost identical to that of the inlet air temperature. In the planes shown in Figure 17, it can be observed that the greater development of hot spots occurs at the terminations of the rows of the server racks. The temperature fields of the x-y planes ( Figure 18) represent the temperature prediction for a usual thermal load with minimum air flow velocity. Due to the normal operation of the server racks, the air flow velocity is insufficient to avoid the formation of hot spots in several server racks. In the extraction zone, hot spots are not expected since the average dimensionless temperature does not exceed T adm = 0.125. Figures 17 and 18, it can be assessed that the server racks at the end of the rows reach higher temperatures, meeting the dimensionless temperature value of T adm = 0.400. The hot aisles are expected to reach the maximum average mean temperatures, T adm = 0.120.

By analyzing both
The result of simulations for the maximum level of air flow velocity and the obtained temperature field predictions can be seen, respectively, in Figures 19 and 20 for the y and z planes.
By observing Figure 19, it is possible to assess that the values obtained in the dimensionless temperature range of the server racks do not exceed the value, T adm = 0.160. In the extraction zone, the simulations show that hot spots are not foreseen to occur, since the value of the dimensionless temperature field does not exceed the value of T adm = 0.06. Figure 20 shows the absence of hot spots, since the inlet air flow velocity proves to be oversized for this case study, which naturally leads to the prediction of low values in the extraction zone. It can be seen that server racks at the end of the rows are very close to the values of the central server racks, since the scale does not present the maximum value existing in the other three cases. Thus, it can be seen that there are no hot spots in any server rack. The value of the temperature range in the extraction zone is almost identical to that of the inlet air temperature.

The Velocity Vectors
A study of the velocity vectors is also presented in this subsection. In the four case studies addressed in this section, the minimum and maximum air flow velocity were applied, thus the behavior of the velocity vectors was equal in all of them. The only noticeable difference is in the minimum and maximum air flow velocity. Thus, only for these two cases are the velocity vectors illustrated. A representation of three-dimensional geometry of velocity vectors in the data center

The Velocity Vectors
A study of the velocity vectors is also presented in this subsection. In the four case studies addressed in this section, the minimum and maximum air flow velocity were applied, thus the behavior of the velocity vectors was equal in all of them. The only noticeable difference is in the minimum and maximum air flow velocity. Thus, only for these two cases are the velocity vectors illustrated. A representation of three-dimensional geometry of velocity vectors in the data center room can be observed in Figure 21.

Minimum Air Flow Velocity
In the case of the minimum air flow velocity, the velocity vectors predictions can be seen in Figure 22-24 for the x, y and z planes, respectively.

Minimum Air Flow Velocity
In the case of the minimum air flow velocity, the velocity vectors predictions can be seen in Figures 22-24 for the x, y and z planes, respectively.
It is possible to observe in Figure 22b,d for x adm = 0.463 and x adm = 0.716 that near the inlet air flow vent the air collides with the plates at the extremity of the server racks. Thus, it becomes noticeable that the air flow gains speed at the end of each row.
By observing Figure 23a-c, it is expected that air recirculation will occur in some racks. This is due to the difficulty that air has in entering them, since the injected inlet air collides with the plates at the extremity of the server racks as mentioned above. In Figure 23b, it is predicted that the air flow will follow the desirable path. In this plan, it is possible to understand the behavior of the air flow in most of the server racks present in the data center room. It can be seen in the hot aisles that the prediction of the flow direction of the air exiting the server racks follows the desirable path to the extraction zone, with the exception of the air adjacent to the ends of the aisles, which tends to enter the exit of the server racks.
By observing Figure 24a,b, it is predicted that, in the mid-height of the inlet air flow vent area, more pronounced values of air velocity can be noticed. In the inlet air flow vent area, vortices and small air recirculation are expected to occur due to the amount of movement blown in by the air vents. Air is still expected to be diverted from the desirable path towards the cold aisles as it collides with the plates at the extremity of the server racks and with the containment doors. Due to deviations in the trajectory, the air increases its velocity at the beginning of the cold aisles, gradually losing it as it moves towards the center of the cold aisle. It is also predicted that the floor and ceiling of the data center room will experience slower air flow velocities.

Minimum Air Flow Velocity
In the case of the minimum air flow velocity, the velocity vectors predictions can be seen in Figure 22-24 for the x, y and z planes, respectively.

Maximum Air Flow Velocity
In the case of the maximum air flow velocity, the velocity vectors predictions can be seen, respectively, in Figures 25-27 for the x, y and z planes.
In this case, the prediction is very similar to that of the minimum air flow velocity, since the air flow behaves in the same way. However, in this case, the velocity of the air flow is higher than the minimum air flow velocity scenario. Figures 25 and 27 indicate that the behavior of the air flow is similar to the minimum air flow velocity scenario. Figure 26 shows the prediction of airflow in the hot aisles. It can be seen that in these planes the air in the center of the hot aisles is correctly directed to the extraction zone, as it can be seen in Figure 26b. However, at the end of the hot aisles, as seen in Figure 26a,c, the air has a horizontal direction, since in the server racks at the end, the flow is made in reverse because the problem mentioned in Section 5.5.1 occurs in the cold aisles. Figure 25b,d and Figure 27b demonstrate the prediction of the air flow velocity increases since near the inlet air flow vent the air collides with the plates at the extremity of the server racks. Due to the increase of this velocity, the air near the server racks located at the end experiences circulation difficulties.
The results show that the airflow could be better close to the end of the rows of the server racks. Such a phenomenon occurs because the layout of the data center room is imperfect. This imperfect cooling can be noticed by observing the enhanced section of Figure 26a, where the air inadequately flows, from the hot aisle, in the direction of the cold aisle. Nevertheless, as reinforced by the information conveyed by the planes in Figures 25 and 27, the cooling is significantly more efficient in the middle of the room, as demonstrated by the air flow represented in the enhanced section of Figure 26b, in which the air adequately flows in the direction of the hot aisle from the cold aisle.
it moves towards the center of the cold aisle. It is also predicted that the floor and ceiling of the data center room will experience slower air flow velocities.

Maximum Air Flow Velocity
In the case of the maximum air flow velocity, the velocity vectors predictions can be seen, respectively, in Figure 25-27 for the x, y and z planes.     In this case, the prediction is very similar to that of the minimum air flow velocity, since the air flow behaves in the same way. However, in this case, the velocity of the air flow is higher than the minimum air flow velocity scenario. Figures 25 and 27 indicate that the behavior of the air flow is similar to the minimum air flow velocity scenario. Figure 26 shows the prediction of airflow in the hot aisles. It can be seen that in these planes the air in the center of the hot aisles is correctly directed

Result Discussion and Analysis
Due to changes in the boundary conditions imposed on the heat flow and the different applied air flow velocities, the case studies differ with respect to the temperature fields, which can cause a few hot spots in the server racks. Figure 28 represents the maximum dimensionless temperature predicted in each case study. In this figure, the temperature limit of hot spots in the server racks can be observed. When analyzing the case studies and Figure 28, it is predicted that not all generate hot spots in server racks, especially in the cases with the maximum air flow velocity. to the extraction zone, as it can be seen in Figure 26b. However, at the end of the hot aisles, as seen in Figure 26a,c, the air has a horizontal direction, since in the server racks at the end, the flow is made in reverse because the problem mentioned in Section 5.5.1 occurs in the cold aisles. Figures 25b,d and 27b demonstrate the prediction of the air flow velocity increases since near the inlet air flow vent the air collides with the plates at the extremity of the server racks. Due to the increase of this velocity, the air near the server racks located at the end experiences circulation difficulties.
The results show that the airflow could be better close to the end of the rows of the server racks. Such a phenomenon occurs because the layout of the data center room is imperfect. This imperfect cooling can be noticed by observing the enhanced section of Figure 26a, where the air inadequately flows, from the hot aisle, in the direction of the cold aisle. Nevertheless, as reinforced by the information conveyed by the planes in Figures 25 and 27, the cooling is significantly more efficient in the middle of the room, as demonstrated by the air flow represented in the enhanced section of Figure  26b, in which the air adequately flows in the direction of the hot aisle from the cold aisle.

Result Discussion and Analysis
Due to changes in the boundary conditions imposed on the heat flow and the different applied air flow velocities, the case studies differ with respect to the temperature fields, which can cause a few hot spots in the server racks. Figure 28 represents the maximum dimensionless temperature predicted in each case study. In this figure, the temperature limit of hot spots in the server racks can be observed. When analyzing the case studies and Figure 28, it is predicted that not all generate hot spots in server racks, especially in the cases with the maximum air flow velocity. In this sense, it is possible to reduce the thermal load of the server racks by improving the overall thermal performance of the processing room, without jeopardizing the correct operation of the server racks located in the sensible regions of the data center room and increasing the sustainability of the entire data center. It can also be observed that Case 4 is the one that displays the best behavior and in which there is a lower probability of hot spot occurrence.
As stated in Section 4, the maximum airflow dimensionless velocity, vadm, is 1, a velocity that is In this sense, it is possible to reduce the thermal load of the server racks by improving the overall thermal performance of the processing room, without jeopardizing the correct operation of the server racks located in the sensible regions of the data center room and increasing the sustainability of the entire data center. It can also be observed that Case 4 is the one that displays the best behavior and in which there is a lower probability of hot spot occurrence.
As stated in Section 4, the maximum airflow dimensionless velocity, v adm , is 1, a velocity that is not a common practice in the studied data center. The velocity is set only in cases when the demand for the rack services and operation is greater than usual and the larger part of the racks function at full capacity. By setting the CRAC units to a maximum airflow velocity, the data center room can be successfully cooled, yet this refrigeration is achieved at a high cost since it is accomplished through the consumption of a high quantity of energy. Thus, as an alternative to upgrading the existing CRAC units, which could signify high expenses, a possible and more sustainable resolution can be altering the way the racks are arranged inside the data center room. The purpose is optimizing the existing layout with the aim of efficiently refrigerating the server racks. One of these alterations could be the removal of a few server racks located at the extremity of the line and consequently avoiding hot spots by improving the air flow.

Conclusions
A parametric study was made by employing CFD modeling for a real data center room comprising more than 200 server racks with the purpose of studying the thermal performance and airflow. Four case studies with distinct thermal loads were addressed in this study. Several scenarios were simulated for the heat flow limit condition and the results for four case studies were obtained and then analyzed for the maximum and minimum air flow velocity. The main results are as follows: • For all four case studies, the results predict the occurrence of several hot spots in the case of the minimum air flow velocity. • By analyzing the results, no hot spots materialize in the case of the maximum air flow velocity except in Case 1.

•
Case 4 shows the best solution. Conversely, refrigerating through these means is unsustainable and consumes a great amount of energy.

•
From the results, it can also be deduced that every predicted hot spot is situated at the end of the aisles. The reason behind this occurrence is the ineffective refrigeration of the data center server racks as a result of an inappropriate position of the CRAC units and this happens in every case study.

•
Having to always use the maximum air flow velocity, as in Cases 1-3, is not the best solution. The airflow is suboptimal close to the end of the rows since the layout of the data center room is imperfect. This challenge can be overcome by withdrawing a few server racks located at the end of the aisles. Such a solution could signify a much better airflow and could even allow a reduction of the air flow velocity, thus diminishing the overall cost with entire cooling operation.