JP2016091320A - Thermal fluid analysis method and apparatus for server decentralization/centralization - Google Patents

Abstract

PROBLEM TO BE SOLVED: To dynamically distribute and centralize a server by introducing an objective function related to server distribution and centralization while dynamically handling a thermo-fluid equation to minimize the power consumption of an air conditioner in a certain period. Provided are a thermo-fluid analysis method and an apparatus. The present invention includes a display unit having a UI for setting conditions, a user-defined unit for setting model parameters and conditions input through the UI, an arithmetic processing control unit, and a data storage unit for arithmetic processing control. The department calculates the thermo-fluid equation based on the set model parameters and conditions until a predetermined period of time elapses, and sums the weighted average of the distance between the server and the air conditioner and the exhaust temperature between adjacent servers in the rack. It is characterized in that the load amount in each server is visualized on the display unit by performing a calculation that minimizes the included objective function, selecting a combination of model parameters and conditions that minimizes the cost after a lapse of a predetermined period. [Selection diagram] Fig. 2

Description

  The present invention relates to a thermal fluid analysis method and apparatus for server distribution / concentration, and more specifically, introduces an objective function related to server distribution / concentration while dynamically handling a thermofluid equation to perform air conditioning over a certain period of time. The present invention relates to a thermal fluid analysis method and apparatus for dynamic server distribution / centralization that minimizes power consumption of a machine.

  In data centers that are showing rapid expansion year by year, efforts to save power have become an urgent issue. In data centers where servers and various electronic devices are installed, reducing power consumption in air conditioners that cool the heat generated by these devices is an important issue, but because a thermal fluid phenomenon is involved It is necessary to consider the efficiency of heat circulation.

  In order to increase the efficiency of heat circulation, the cold air blowing method from the floor or ceiling, the arrangement and spacing of racks, the introduction of natural cold air, and server load distribution called virtualization technology have been used. However, since the load on the server changes with time and a local flow is formed near the structural features in the data center and near the rack, further improvements in heat circulation efficiency are required. It was done.

R.I. Bourisli, “Numerical Investigation Of Optimal Spacing Between Heat-Generating Blocks In Laminar Forced Convection,” In Proc. IEEE ITHERM, p.1-8, February 2010 M. Iyengar, R. Schmidt, and J. Caricari, “Reducing Energy Usage In Data Centers Through Control Of Room Air Conditioning Units,” In Proc. IEEE ITHERM, p.1-11, 2010

  In order to improve the cooling efficiency of the server, it is also considered to actually experiment by changing the server rack layout and cooling method, but the actual experiment in the data center is limited and the actual operation Experiments in the state can be difficult. Therefore, as shown in Non-Patent Documents 1 and 2, numerical simulation is widely applied as an important evaluation means.

  However, the conventional numerical simulation uses a model that assumes a constant server load due to the complexity of mathematical problems, that is, a static state, and uses a thermal fluid equation without fully considering thermal fluid phenomena. There were many models that considered only the temperature change near the server. For this reason, it has not been possible to realize an optimized expression that solves an energy optimization problem to which a numerical simulation model that simulates dynamic load fluctuations of an actual server and a dynamic thermofluid equation is applied.

  Therefore, the present invention considers the dynamic heat generation change of the server while dynamically handling the thermofluid equation, and introduces an objective function related to server distribution / concentration to optimize the air conditioner for a certain period as an optimization problem. It is an object of the present invention to provide a thermal fluid analysis method and apparatus for dynamic server distribution / centralization that minimizes power consumption.

  In order to solve the above-described problem, the thermal fluid analysis apparatus for server distribution / concentration according to claim 1 is a server in a server room in which a plurality of racks in which servers are stacked are installed and the server is cooled by indoor air conditioning. A thermal fluid analysis apparatus for server distribution and centralization for operation, comprising a display unit having a user interface for setting conditions relating to mathematical equations, and model parameters and conditions input through the user interface A user definition unit to be set, an arithmetic processing control unit, and a data storage unit are provided, and the arithmetic processing control unit is based on model parameters and conditions set by the user definition unit until a predetermined period elapses. Calculate the thermal fluid equation and use the calculation result of the thermal fluid equation to calculate the weighted average and the rack Calculation for minimizing the objective function so that the objective function including the sum of exhaust temperatures between adjacent servers is minimized, storing the minimized objective function in the data storage unit, and the arithmetic processing control unit Selects a combination of model parameters and conditions that minimizes the cost based on a calculation result stored in the data storage unit after the predetermined period has elapsed, and selects the selected combination of model parameters and conditions as the data The load on each server calculated based on the combination of the selected model parameters and conditions during the predetermined time is stored on the storage unit and visualized on the display unit.

  The thermal fluid analysis apparatus for server distribution / centralization according to claim 2 is the thermal fluid analysis apparatus for server distribution / centralization according to claim 1, wherein the calculation of the thermofluid equation is as follows: Navier-Stokes equations and energy equations are used.

  The thermal fluid analysis method for server distribution / concentration according to claim 3 is a server distribution / concentration method for server operation in a server room in which a plurality of racks in which servers are stacked are installed and the server is cooled by indoor air conditioning. A thermal fluid analysis method for centralization, a display unit having a user interface for setting conditions relating to mathematical equations, a user definition unit for setting model parameters and conditions input through the user interface, and an operation Using the thermal fluid analysis device for server distribution / concentration provided with a processing control unit and a data storage unit, set in the arithmetic processing control unit by the user definition unit until a predetermined period elapses. The thermal fluid equation is calculated based on the model parameters and conditions, and the calculation result of the thermal fluid equation is used to calculate the distance between the server and the air conditioner. Calculation to minimize the objective function so that the objective function including the sum of the weighted average and the exhaust temperature between adjacent servers in the rack is minimized, and the minimized objective function is stored in the data storage unit Selecting a model parameter and condition combination that minimizes the cost based on a calculation result stored in the data storage unit after the predetermined period has elapsed, and selecting the selected model parameter and condition combination. Storing in the data storage unit, and visualizing the load amount in each server calculated based on the combination of the selected model parameter and condition during the predetermined time on the display unit, Features.

  The thermal fluid analysis method for server distribution / centralization according to claim 4 is the thermal fluid analysis method for server distribution / centralization according to claim 3, wherein the calculation of the thermofluid equation is as follows: Navier-Stokes equations and energy equations are used.

  According to the present invention, while dynamically handling the thermofluid equation, considering the dynamic heat generation change of the server, the objective function related to server distribution / concentration is introduced to optimize the air conditioner over a certain period as an optimization problem. A thermal fluid analysis method and apparatus for dynamic server distribution / concentration that minimizes power consumption can be provided, and its effect has been shown for numerical experiments. As a result, it has become possible to further reduce power consumption by taking into account the thermal fluid phenomenon in the conventional virtualization technology.

It is a figure which illustrates the composition of the thermal fluid analysis device for server distribution and concentration concerning the present invention. It is a figure which shows the flow of the main processes of the thermofluid analysis method for server distribution and concentration based on this invention. It is a figure for demonstrating the control method of the dynamic server heat generation amount in this invention. The figure which shows the server room which has arrange | positioned two racks displayed by UI in this invention and three servers per rack is shown. It is a figure which illustrates the number and definition in a some rack and server in a server room. Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). Is a diagram showing a state of temperature of the air conditioner in the server room in the case of variously changing the size of the weighting factor sigma ₁ and sigma ₂ in Equation (8). The time change of the power consumption for every six servers by the dynamic distribution and aggregation of the servers in the experiments shown in FIGS.

  The thermal fluid analysis method for dynamic server distribution / concentration in the present invention relates to server load distribution / concentration by dynamically using the thermofluid equation and considering the heat generation amount of the server that changes with time. The objective function is optimized so as to minimize a weighted average related to the distance between the server and the air conditioner and exhaust heat between adjacent servers in the rack over a certain period. This makes it possible to perform numerical simulations that simulate actual server load fluctuations and to take into account dynamic changes in the amount of heat generated by the server, making it possible to calculate server load distribution and concentration with higher accuracy than conventional simulations. It becomes.

  FIG. 1 illustrates the configuration of a thermal fluid analysis apparatus for dynamic server distribution / centralization according to the present invention. In FIG. 1, a display unit 100 having a user interface (UI), a user definition unit 101 for setting conditions regarding a mathematical equation input through the UI, an arithmetic processing control unit 102, and a data storage unit 103 are provided. Prepare.

  FIG. 2 shows a main processing flow of the thermal fluid analysis method for dynamic server distribution / centralization according to the present invention. As shown in FIG. 2, in the thermal fluid analysis method for dynamic server distribution / concentration in the present invention, in step 201, the user definition unit 101 includes model parameters related to mathematical equations input by the user through the UI, and Set conditions. In step 202, the arithmetic processing control unit 102 calculates a two-dimensional or three-dimensional thermal fluid equation within the range of the model parameters and conditions determined by the user definition unit 101. In step 203, the arithmetic processing control unit 102 uses the calculation result of the thermal fluid equation to obtain the distribution and aggregation of the server, and within the range of model parameters and conditions determined by the user definition unit 101, The objective function is minimized so that the objective function including the weighted average of the air conditioning distance and the exhaust temperature between adjacent servers in the rack is minimized. In step 204, the minimized objective function is stored in the data storage unit 103. In step 205, the arithmetic processing control unit 102 determines whether or not the predetermined time has elapsed, and determines whether or not the predetermined time within the range of the model parameters and conditions determined by the user definition unit 101. The thermal fluid equation calculation in step 202 and the objective function minimization calculation in step 203 are continued until the time elapses, the calculation result is stored in the data storage unit 103, and the calculation is terminated when a predetermined time elapses. In step 206, the arithmetic processing control unit 102 selects a combination of model parameters and conditions that minimizes the cost, based on the calculation result stored in the data storage unit 103. In step 207, the arithmetic processing control unit 102 stores the selected combination of model parameters and conditions as final data in the data storage unit 103, and sets the selected combination of model parameters and conditions during the predetermined time. The display unit 100 visualizes the load amount calculated on the basis of each server.

  The arithmetic processing control unit 102 uses a Navier-Stokes equation or an energy equation capable of time evolution calculation when calculating the thermal fluid equation in step 202. Further, when the arithmetic processing control unit 102 performs the objective function minimization calculation in step 203, the sum of the weighted average regarding the distance between the server and the air conditioner and the exhaust temperature between the adjacent servers in the rack becomes the minimum. Calculate as follows.

(Example)
The basic mathematical equation according to the present invention will be described. Since the object to be handled is an air flow phenomenon, the viscosity coefficient is set to zero in this embodiment. This eliminates the calculation of the diffusion term in Navier-Stokes (NS), which is the governing equation related to velocity and pressure, and enables a significant high-speed solution. A group of equations relating to unknowns such as the NS formula is discretized by a regular Cartesian grid, and variables relating to atmospheric pressure and velocity are MAC grids arranged on the center and grid points of the mesh, respectively. The NS expression is efficiently and accurately solved by the MacCorack algorithm and the MICCG (0) method under boundary conditions and initial conditions. The MICCG (0) method was applied to the projection term in order to be coupled with the NS equation under the constraint that the velocity divergence is zero (continuous), that is, ∇u = 0, assuming the incompressibility of the fluid. .

The total external force f in the NS type, f _Coldw wind from the air conditioner, the wind of the f _Hotw from the server, the buoyancy of the heat flow to f _tb, when the f _co and external force by the eddy particle method (Vortex Particle Method), f = f _coldw + f _hotw + f _tb + f _co Further, when ρ is the fluid density, T is the temperature of the fluid, T _amb is the environmental temperature, and α and β are adjustment coefficients, f _tb = −αρz + β (T−T _amb ) z. z = (0, 0, 1) is a vertical unit vector in three dimensions.

When u is a velocity, p is an atmospheric pressure, _ck is a thermal conductivity, and T _sour is a temperature of the indoor space, the NS equation u _t and the energy equation T _t regarding the time change of the temperature are _expressed by the following (Equation 1) and ( It is shown by Formula 2). ∇ and ∇ ² are operators related to the spatial first derivative and spatial second derivative, respectively.
u _t = − (u∇) u− (∇p / ρ) + u∇ ² u + f (Equation 1)
T _t = − (u∇) T + c _k ∇ ² T + T _sour (Equation 2)

Here, T _Hotsour the server exhaust heat temperature, when blown up cold air temperature of the T _Coldsour from the floor of the air conditioner (CRAC), a _{T sour = T hotsour + T coldsour} . T _coldsour and c _k are defined by the user definition unit 101 through the UI by the user. The MacCormack algorithm and the MICCG (0) method are also used when solving equations related to temperature.

T _in is the intake temperature of the server, P is the load power (load amount) of the server, C _air is the heat capacity of the air, Δt is the time interval, S is the cross-sectional area of the server outlet, and V is the airflow exhausted from the server If the speed is M and M is the total mass of the airflow, the exhaust temperature T _{hotsour of the} server is _expressed by the following (Equation 3). In the present specification, for simplification, description will be made assuming that the load amount = the heat generation amount.

Rate of air flow is exhausted to _ser V from the server, if the fluid density of the air current is exhausted to [rho _air from a server, the variables and V _ser defined above, the total mass M of the airflow about the server, the following (formula 4).
M = Δt · S · V _ser · ρ _air (Formula 4)

In the present embodiment, it is assumed that the exhaust temperature for each server is determined by the heat generation amount, heat capacity, and intake air temperature T _in of the server. _Assuming that P _ser is the load amount of the server, the exhaust temperature T _{ser_out of the} server is expressed by (Expression 5) below from (Expression 3). P _ser and V _ser are defined by the user definition unit 101 through the UI by the user.

If V _acw is the cool air blowing speed from CRAC, P _cr is the load of CRAC, T _{cr_in} is the intake temperature of CRAC, T _{cr_out} is the exhaust temperature of CRAC, and S _in is the cross-sectional area of the blowing part of CRAC, then P _cr is It is shown by the following (Formula 6). V _acw and T _{cr_out} are defined by the user definition unit 101 through the UI by the user.
P _cr = C _air M (T _{cr_in} −T _{cr_out} )
= C _air S _in V _ser ρ _air V _acw (T _{cr_in} −T _{cr_out} )
= 1.003 × 1.29 · S _in V _acw (T _{cr_in} −T _{cr_out} ) (Formula 6)

The CRAC set temperature T _cr was varied in three steps as follows. T _th1 and T _th2 are upper limit temperature and the lower limit temperature of the intake air temperature of the CRAC, T _a is the average temperature in the server room (ambient temperature: ambient temperature), T _s is the temperature in the vicinity of the suction port to the CRAC in the ceiling of the room It is.

In this example, T _th1 = 30 ° C., T _th2 = 20 ° C., and T _a = 17 ° C. were set empirically. T _cr is assumed to vary between 10 ° C. and 20 ° C. Where appropriate, the average temperature is used from the temperature on the calculation grid. While T _cr is higher than T _a , the average indoor temperature is suggested to be low and stable, and the CRAC is turned off. The longer the CRAC is off for a certain period, the more power is saved in the server room.

  Next, in order to minimize the load power of the CRAC in a certain time, the distribution and aggregation of the server load position and the load amount are considered under the condition that the total server load amount is constant. This was done within the framework of the energy minimization problem. In order to consider the optimization, minκ was calculated by defining the objective function κ shown in (Equation 12) and minimizing the objective function κ. Here, d indicates the distance between the server and the air conditioner.

Regarding (Equation 8), the first term on the right side represents a weighted average of the position of a certain server and the position of the air conditioner, and the second term on the right side represents that the exhaust temperature T _out of adjacent servers in the rack is continuous. Show. The first and second terms in (Equation 8) are given by weighting coefficients σ ₁ and σ ₂ (σ ₁ + σ ₂ = 1). Regarding the continuous temperature distribution, there is little heat circulation due to local vortices and the like, and a smooth air flow is expected. The standard for optimization is that the average temperature in the server room decreases.

The minimum distance is searched for the distance d between the server and the air conditioner according to the first term on the right side of (Expression 8). The minimum distance d _i between the server and the air conditioner is expressed by the following (Equation 9).

P _tot is the total load amount (corresponding to the heat generation amount) of the server, and is represented by the following (Equation 10).

  A dynamic server heat generation control method according to the present invention will be described with reference to FIG. In general, the server load arrangement is highly uncertain. Therefore, in the numerical simulation in the present embodiment, the temporal load variation characteristic is generated in a pseudo manner according to the sine wave function and the probability distribution. The user can freely select an average value and a standard deviation for the normal distribution and Poisson distribution via the UI, and can change the amplitude and period (frequency) for the sine function. A numerical simulation was performed by giving each server an amount corresponding to the size of the function given via the UI.

  3A shows a sine function in which the horizontal axis represents time and the vertical axis represents amplitude, and FIG. 3B shows the Gaussian function in which the horizontal axis represents server load and the vertical axis represents probability, and FIG. ) Indicates the average value of the server load (because the Poisson function uses the average value as a variable) on the horizontal axis, and the Poisson function indicating the probability on the vertical axis.

  FIG. 4 is a diagram illustrating an arrangement in the server room displayed by the UI according to the present invention. FIG. 4 shows an arrangement example of a server room in which two racks displayed on the UI in the present invention and three servers per rack are arranged. Here, an air conditioner intake air temperature measurement unit, a temperature measurement unit in a blowout unit of the air conditioner, six servers, and two racks are shown.

  In the server room shown in FIG. 4, the cold air from the air conditioner shows a method in which the airflow flows from the bottom floor to the upper right ceiling. Also, in the server room shown in FIG. 4, the server intake section is provided on the opposite side of the two racks, and the exhaust heat of the server is released from the outside of the two racks. The same applies to the server room shown below.

FIG. 5 illustrates numbers and definitions in multiple racks and servers. As shown in FIG. 5, the server numbers are defined so as to increase from the upper left to the lower and toward the right. S ₃ C ₀ and S ₃ C _{1 in} FIG. 5 are used, for example, when obtaining d _i in (Equation 9), and the positions of the temperature measuring units C ₀ and C ₁ in the server S ₃ and the blowing unit of the air conditioner. Represents the distance to each other.

6 to 12 illustrate the state of the temperature of the air conditioner in the server room when the magnitudes of the weighting coefficients σ ₁ and σ ₂ (σ ₁ + σ ₂ = 1) in (Equation 8) are variously changed. . In FIGS. 6 to 12, a numerical simulation experiment of changes in the blowout temperature of one CRAC is performed with the total load of the server set to 8000 W, and temporal changes of four typical thermal fluid calculations are compared with differences in temperature distribution. . The time on the horizontal axis of the graphs of FIGS. 6E to 12E is a dimensionless quantity and indicates the number of times of discretization on the computer.

Here, FIG. 6 shows a case where the CRAC blowing temperature is fixed at 10 ° C., and FIGS. 7 to 12 show a case where the CRAC is automatically controlled. 8 shows the case where σ ₁ = 1 and σ ₂ = 0, FIG. 9 shows the case where σ ₁ = 0 and σ ₂ = 1, and FIG. 10 shows the case where σ ₁ = 0.2 and σ ₂ = 0.8. 11 shows the case where σ ₁ = 0.5 and σ ₂ = 0.5, and FIG. 12 shows the case where σ ₁ = 0.8 and σ ₂ = 0.2, respectively. 8 and 12, since the time during which CRAC is OFF is long, it has been analyzed in (Equation 8) that setting the weighting coefficient as σ ₁ > σ ₂ leads to power saving. . This is because priority is given to σ _{1 so} that a high load is given to a server close to the air conditioner and a low load is given to a server far from the air conditioner, and the influence of the intake air temperature due to the positional relationship between the air conditioner and the server is considered. Due to Further, as shown in FIG. 8, it has been obtained through dynamic numerical simulation that the heat circulation efficiency and the power saving can be achieved by reducing the allocation of the server load amount toward the upper part of the rack. Such a result is an unprecedented analysis result, and it can be said that the present invention is an optimization method related to dynamic server distribution and aggregation.

FIG. 13 shows the time change of the power consumption for every six servers due to the dynamic distribution and aggregation of the servers in the experiments shown in FIGS. # 1 to # 6 in FIG. 13 correspond to the server numbers shown in FIG. Here, FIG. 13B shows a case where σ ₁ = 1 and σ ₂ = 0, FIG. 13C shows a case where σ ₁ = 0 and σ ₂ = 1, and FIG. 13D shows a case where σ ₁ = 0. 2. When σ ₂ = 0.8, FIG. 13E shows σ ₁ = 0.5, and when σ ₂ = 0.5, FIG. 13F shows σ ₁ = 0.8 and σ ₂ = 0. .2 cases are shown respectively.

  In FIG. 13B and FIG. 13F, the load on the server # 1 is set low. As shown in FIG. 4 and FIG. 5, the return airflow port on the ceiling of the air conditioner is racked in the server room. This is because it is installed asymmetrically with respect to the arrangement.

  As shown in FIG. 13, it was confirmed that the load amount of the dynamically changing server was non-linearly assigned according to the arrangement in the rack by optimal calculation. In other words, since a stable temperature has been reached from the initial state, it is optimized by assigning servers with high load to the bottom of the rack and servers with low load to the top of the rack. Indicated. However, this is the case of the circulation method of the cold air blowing on the floor and the ceiling return, but the other methods can be analyzed in the same manner.

  The present invention can also be applied to industrial fields related to monitoring and image sensing in a real environment in the electric power field, energy field, communication field, and sensing field.

Display unit 100
User definition part 101
Arithmetic processing control unit 102
Data storage unit 103

Claims (4)

A thermal fluid analysis device for server distribution / concentration for server operation of a server room in which a plurality of racks in which servers are stacked are installed and the server is cooled by indoor air conditioning,
A display unit having a user interface for setting conditions relating to mathematical equations;
A user definition unit for setting model parameters and conditions input through the user interface;
An arithmetic processing control unit;
A data storage unit;
The arithmetic processing control unit, until a predetermined period of time,
Calculate the thermal fluid equation based on the model parameters and conditions set by the user-defined part,
Calculation that minimizes the objective function using the calculation result of the thermal fluid equation so that the objective function including the sum of the weighted average of the distance between the server and the air conditioning and the exhaust temperature between adjacent servers in the rack is minimized. To store the minimized objective function in the data storage unit,
The arithmetic processing control unit selects a combination of model parameters and conditions that minimizes the cost based on a calculation result stored in the data storage unit after the predetermined period has elapsed, and selects the selected model parameters and conditions. Is stored in the data storage unit, and the load amount in each server calculated based on the combination of the selected model parameters and conditions during the predetermined time is visualized on the display unit. Thermo-fluid analyzer for dynamic server distribution / centralization.   2. The thermal fluid analysis apparatus for dynamic server distribution / concentration according to claim 1, wherein the calculation of the thermal fluid equation uses a Navier-Stokes equation and an energy equation. A thermal fluid analysis method for server distribution / centralization for server operation of a server room in which a plurality of racks in which servers are stacked are installed and the server is cooled by indoor air conditioning,
A display unit having a user interface for setting conditions relating to mathematical equations;
A user definition unit for setting model parameters and conditions input through the user interface;
An arithmetic processing control unit;
Using a thermal fluid analysis device for server distribution / concentration provided with a data storage unit, the arithmetic processing control unit,
The thermal fluid equation is calculated based on the model parameters and conditions set by the user definition unit until a predetermined period elapses, and the weighted average and rack for the distance between the server and the air conditioning are calculated using the calculation result of the thermal fluid equation. Performing a calculation for minimizing the objective function so as to minimize the objective function including the sum of exhaust temperatures between adjacent servers, and storing the minimized objective function in the data storage unit;
After the predetermined period, based on a calculation result stored in the data storage unit, a combination of model parameters and conditions that minimizes the cost is selected, and the selected combination of model parameters and conditions is selected in the data storage unit. And storing the load amount in each server calculated based on the combination of the selected model parameters and conditions during the predetermined time on the display unit;
A thermal fluid analysis method for dynamic server distribution / concentration, characterized in that   4. The thermal fluid analysis method for dynamic server distribution / concentration according to claim 3, wherein the calculation of the thermal fluid equation uses a Navier-Stokes equation and an energy equation.