Abstract: - In recent years, CAE (Computer-aided engineering) using FEM (Finite element method) simulations
were generally used in the design field. FEM simulations were classified as static implicit and explicit dynamics
methods. Particularly FEM simulations using the static implicit method were very generally and usefully used in
the industrial world. The static implicit method FEM consists of static, buckling, thermal, vibration analyses, and
so on. This FEM thermal analysis can't calculate the phenomena of heat transfer and internal forced cooling in
some enclosures of a machine tool. On the other hand, when the temperature distributions in a structure such as
that are calculated, FVM fluid simulation was used for that. However, this simulation can't exactly calculate a
complex structure in a machine tool, needs very long calculation times and the calculation accuracy is very poor.
Therefore, in this research, the calculation technique regarding FEM thermal simulation using 4 virtual fluid
elements was developed and evaluated for the phenomena of heat build-up and internal forced cooling in some
enclosures of a machine tool. The algorithms and the calculation models regarding 4 virtual fluid elements were
developed, then the proposed method was evaluated in the experiment. It is concluded from the results that; (1)
the 4 virtual fluid elements were developed for FEM thermal simulation instead of FVM fluid simulation, (2)
FEM thermal simulation with the developed 4 virtual fluid elements has high calculation accuracy and a short
calculation time, (3) the proposed method was very effective in the design.
Key-Words: - FEM thermal simulation, Virtual fluid elements, FVM fluid simulation, Machine tool design, CAE,
Virtual technique
Received: May 22, 2022. Revised: January 6, 2023. Accepted: February 17, 2022. Published: March 20, 2023.
1 Introduction
Recently, in the design phase, it is often necessary to
perform fluid analysis of industrial products using
FVM fluid simulation, [1]. Particularly, in the field
of thermal design, [2], [3], [4], [5], [6], it's a very
important technology. However, FEM thermal
simulation is more prevalent than FVM technology
in the industrial world, [7]. In addition, FEM thermal
simulation is easier to operate than FVM fluid
simulation and can be used more quickly. However,
FEM thermal analysis can't calculate the phenomena
of heat build-up and internal forced cooling in some
enclosures of a machine tool. On the other hand,
when the temperature distributions in a structure such
as that are calculated, FVM fluid simulation was used
for that. However, this simulation can't exactly
calculate a complex structure in a machine tool,
needs very long calculation times and the calculation
accuracy is very poor.
In this study, 4 virtual fluid elements (hereinafter
referred to as "virtual fluid elements") were
developed and evaluated in order to substitute FEM
thermal simulation for FVM fluid simulation. Then
the algorithm for these virtual fluid elements was
constructed and evaluated its overall industrial
effectiveness by experimentally evaluating its
computational time and error.
2 Explanation Regarding 4 Virtual
Elements for FEM Thermal Simulation
2.1 Outline for 4 Virtual Fluid Elements
Figure 1 shows the four virtual fluid elements of the
structure to be analyzed. This represents the behavior
of the flowing medium in the structure. It consists of
four elements; I: Simulated element for inflow, II:
Simulated element for outflow, III: Simulated
element for the boundary layer, and IV: Simulated
element for convection.
Development of Calculation Technique for FEM Thermal Simulation
Using Virtual Fluid Elements
IKUO TANABE
School of Engineering
Sanjo City University
5002-5 Kamisugoro, Sanjo, Niigata
JAPAN
HIROMI ISOBE
Department of Mechanical Engineering
Nagaoka University of Technology
1603-1 Kamitomioka, Nagaoka, Niigata
JAPAN
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
8
Volume 18, 2023
2.2 I: Simulated Element for Inflow
As shown in Figure 2(a), at time t = 0 [s], the flowing
medium in the micro-region near the structure inlet
has a heat quantity Qinlet. Then, during Δt [s], a
flowing medium with a heat quantity Qin smaller than
Qinlet flows into this micro-region, and the air in the
micro-region is cooled. As a result, the heat content of
the air in a small region is (Qin+Qinlet)/2.
Next, in the case of the I: Simulated element for
inflow, as shown in Figure 2(b), there is no flowing
medium from the outside during the period from time
t = 0 [s] to Δt [s] (no state shown in Figure 2(a)),
during which heat is dissipated on the inlet surface of
the structure by the heat transfer with the outside
(Qinlet －Qin)/2. The heat content of the flowing
medium in the micro-region near the inlet of the
structure during Δt s becomes (Qin + Qinlet)/2, which
is thermally equivalent to the state shown in Figure
2(a).
This equivalence is shown in Equation (1).
(Qinlet－Qin)/2
＝ (ρairAinVinΔt cairTinlet－ρairAinVinΔt cairT∞) /2
＝ αinletAin (Tinlet－T∞ )Δt (1)
αinlet ＝ ρair Vin cair / 2 (2)
where ρair is density [kg/m3], Ain is the cross-sectional
area of the inlet of the structure [m2], Vin is the velocity
of the inflowing medium [m/s], cair is the specific heat
of the flowing medium [J/kg-K], T∞ is outside flowing
medium temperature [K], Tinlet is structure inlet
temperature [K] and αinlet is heat transfer coefficient
with the outside on the αinlet structure inlet surface
[W/m2-K]. By hypothetically making the structure
inlet have a heat transfer coefficient αinlet is possible
to simulate the entry of a flowing medium with the
heat quantity Qin from the structure inlet.
2.3 II: Simulated Element for Outflow
As shown in Figure 3, at the outlet of the structure, a
fluidized medium with a heat quantity Qout [J] is
evacuated in Δt [s] by forced evacuation by a fan. In
the case of the II: simulated element for outflow, the
actual forced exhaust is simulated by dissipating the
heat quantity Qout [J] in Δt [s] by heat transfer with the
outside of the structure outlet. This equivalence is
Outlet
Fan
Heat source Qh
Inlet
Structure
(a) Schematic view of a structure with a heat source
min
mout
min：Entering mass flow , mout：Exiting mass flow，
αair：Heat transfer coefficient, T∞：Ambient
temperature
T∞
ℓ
αair
Outlet
Boundary
layer
Structure
Imaginary
air model
Fan
Inlet
I: Simulated element for inflow
II: Simulated element for outflow
III: Simulated element for boundary layer
IV: Simulated element for convection
(b) Schematic view of analysis object for FEM thermal simulation
instead of FVM fluid simulation
・
・
Heat source Qh
Figure 1: Analysis object for FEM thermal simulation
instead of FVM fluid simulation. Four virtual
elements I, II, III and IV were used for the FEM
thermal simulation.
(a) Phenomenon of inflow on inlet ( Qin inlet in the structure
during
Δ
t )
Qin + Qinlet
2
Qin
Boundary layer
Structure
Imaginary
air model
Air
Minute
region on
inlet
Qinlet
Inlet
(b) Heat transfer of I: Simulated element for inflow
((Qinlet－Qin ) / 2 outlet from the structure during
Δ
t
)
Qin + Qinlet
2\
Boundary layer
Structure
Imaginary
air model
Minute region
on inlet Qinlet
Inlet
Air
Qinlet − Qin
2
Tinlet
Figure 2: Explanation of
I: Simulated element for
inflow. Heat capacity of the inflow was
supposed by heat transfer.
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
9
Volume 18, 2023
shown in Equations (3) and (4).
Qout = ρair Aout Vout Δt cair Tout
= αoutlet (Tout － T∞) AoutΔt (3)
αoutlet = ρair Vout cair / (1 － T∞ / Toutlet) (4)
where Aout is the cross-sectional area of the exit
surface of the structure [m2], Vout is the velocity of
the outflowing medium [m/s], and Tout is the average
temperature of the outflowing medium [K]. αoutlet
is the heat transfer coefficient on the exit surface of
the structure[W/m2-K]. Where the medium volume
flowing out from the structure in Δt [s] (AoutVoutΔt
=AinVinΔt) is small, the average temperature Tout of the
outflowed medium is taken to be equivalent to the
temperature Toutlet on the exit surface of the structure
(Tout = Toutlet). By using the heat transfer coefficient
αoutlet as a boundary condition on the structure outlet in
FEM thermal simulation, it is possible to simulate the
forced exhaust from the structure.
2.4 III: Simulated Element for Boundary
Layer
As shown in Figure 4, the heat quantity Qb is
transferred by heat transfer between the inner wall of
the structure and the flowing medium in the structure
in Δt [s]. In this case, the amount of heat transfer at
the wall surface of the structure by the heat transfer is
referred to as Qb [J]. In the case of the III: simulated
element for the boundary layer, the thermal
conductivity λb [W/m-K] is set to simulate the heat
transfer. The equivalent states are shown in Equation
(5).
Qb = αw (Tb1－Tb2) AwΔt=(λb / lb) (Tb1－Tb2) AwΔt,
λb = αw lb (5)
Where, αw is the heat transfer coefficient of the inner
wall of the structure [W/m2-K], Tb1 is the temperature
of the boundary layer in contact with the III: simulated
element for boundary layer [K], Tb2 is the temperature
of the inner wall of the structure in contact with the
boundary layer [K], Aw is the surface area of the inner
wall of the structure [m2], λb is the thermal
conductivity of the boundary layer [W/m-K] and lb is
the thickness of the boundary layer [m]. Although this
lb can be any value, we henceforth fix lb = 0.003 m in
this research because there are enough virtual fluid
elements in the structure and the FEM mesh model is
not too fine.
2.5 IV: Simulated Element for Convection
The energy balance of the flow medium in the
structure is considered by replacing it with the
elements shown in Figure 5. Considering the energy
balance, [8], the following assumptions are made. a:
The velocity and temperature fields are steady, and
the airflow is laminar. b: The physical properties
(density, viscosity, thermal conductivity, and specific
heat at constant pressure) are constant. c: Ignoring
heat generation in viscous work. As the thermal
conductivity of the IV: simulated element for
convection is presented in Figure 5, the pseudo-
thermal conductivity λcx (heat transfer coefficient that
simulates the convection phenomenon) is defined.
The heat quantity Qcx [J] moving through the elements
in Δt seconds is given by Equation (6). The heat
transfer Qcx' [J/s] per unit time due to the heat transfer
caused by the temperature difference between the two
ends of the elements in Figure 5 is given by Equation
(7). In the x-axis direction, the difference between the
energy exerted by convection from the right side of
the microelement and the energy inflowed by
convection from the left side of the microelement is
given by Equation (8). Since the sum of the heat
transferred by heat conduction in Equation (7) and the
Qout
Air
Fan
Figure 3: Explanation of
II: Simulated element for
outflow (Qout outlet in the structure during
Δ
t ). Heat capacity of the outflow was
supposed by heat transfer.
Structure Boundary layer
Imaginary air model
Tout
Toutlet
・
λb
Tb1\
\\\\\\\
Qb
xb
Tb2
Figure 4: Explanation of III: Simulated element for
boundary layer (Qh transfer to the structure
during Δt ). Heat transfer of the boundary layer
was supposed by heat conduction.
Imaginary air model Boundary layer
Structure
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
10
Volume 18, 2023
energy transferred by convection in Equation (8) is
the heat of Equation (6), Equation (9) is obtained. To
organize this Equation (9), the pseudo-thermal
conductivity λcx is shown in the equation from
Equation (10), the pseudo-thermal conductivity of a
virtual fluid element to simulate convection can be
obtained; we have proposed the pseudo-thermal
conductivity λcx in the x-direction, but the pseudo-
thermal conductivities λcy and λcz in the y-direction
and z-direction can also be obtained by the same equation.
Qcx ＝( λcx / lcx) (Tc1－Tc2 ) Ac
Δ
t (6)
Qcx’＝( λair / lc) ( Tc1 － Tc2) Ac
Δ
t (7)
ρair cair Vcx Ac Tc1
Δ
t
－
ρair cair Vcx Ac Tc2
Δ
t
＝ ρair cair Vcx Ac ( Tc1 － Tc2 )
Δ
t (8)
ρair cair Vcx Ac (Tc1－Tc2)
Δ
t ＋( λair / lc) (Tc1－Tc2)
Ac
Δ
t＝( λcx / lcx) (Tc1－Tc2) Ac
Δ
t (9)
λcx ＝ λair ＋ ρair cair Vcx lcx (10)
Where λcx is the pseudo-thermal conductivity of the
element in the x-direction [W/m-K], lcx is the length
of the element in the x-direction [m], Tc1 is the
temperature of the left end of the element [K], Tc2 is
the temperature of the right end of the element [K], Ac
is the cross-sectional area of the yz plane of the
element [m2], λair is pseudo-thermal conductivity of
the flow medium in the x-direction [W/m-K], Vcx is
the average velocity of air in the x-direction [m/s].
3 Evaluation Experiment Regarding 4
Virtual Fluid Element Models
3.1 Experimental Models and Experimental
Conditions for Evaluation
Experiments were conducted on three models to
evaluate the effectiveness of the proposed 4 virtual
fluid elements. First, a complex model for the
evaluation experiment is shown in Figure 6, with one
fan, two shielding walls, and ceramic heaters in two
locations. The fluid medium is air. For the simple
model, two shielding walls and one heater were
removed from the complex model. As shown in
Figure 6(b), thermocouples were installed inside the
enclosure and around the periphery of the enclosure
in all experiments (16 points for simple models and
26 points for complex models) and measured until the
temperature of the air or cutting fluid reached steady
state, respectively. In the forced cooling with cooling
oil for the complex structural model shown in Figure
7 (for the third model), the previous complex
enclosure was appropriated and cooling oil was used
instead of air. The intake and exhaust ports of the
enclosure in the complex model were modified to a
nipple structure so that the coolant oil can be easily
supplied, and the enclosure was erected so that the
inlet is downward and the outlet is upward. As shown
in Table 1, two types of exhaust velocity, two types
of the heating value of the heater, and two types of
cooling oil inflow were used as the parameters of the
experiment. The room temperature was 20±1°C.
3.2 Calculation Models and Its Analysis
Condition for Evaluation
The three experiments in the previous section were
calculated using SolidWorks 2017. Figure 8 shows
the computational models of FEM thermal simulation
and FVM fluid simulation for the complex model,
respectively, as representatives of each
computational model. The FVM fluid simulation is
performed using Flow simulation, a fluid analysis
software on SolidWorks 2017, and the discretization
method is the finite volume method. For comparison,
the number of nodes in the FEM thermal simulation
model and the number of elements in the FVM fluid
simulation model are almost the same. Table 2 and
λcx，λair
Tc1
Tc2
lcx
ρair cair Vcx Ac Tc1Δt
ρair cair Vcx Ac Tc2Δt
lcy
Qcx＝ (Tc1－Tc2) AcΔt
λair
lc
Heat conduction
Convection
Figure 5: Energy balance of an element for x direction in structure. Convection in the structure was supposed by heat
conduction.
x
y
Imaginary flowing medium model
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
11
Volume 18, 2023
20
20
35
・
Case of material：SS400, t = 10 mm
x
y
z
Ceramics heater
Air
Fan
(0.5 m/s or 1.0m/s)
70
Partition plate1
・
50
50
25
45
50
20
50
50
25
25
Partition plate 2 ( SS400 , t =10mm )
140
70
The model is originally covered
(a) Schematic view of experimental set-up for the complex model
20
Point
Coordinate
Point
Coordinate
Point
Coordinate
①
(90,11.75,57.5)
⑪
(45,45,0)
e
(130,45,120)
②
(45,10,40)
⑫
(0,45,35)
f
(140,45,100)
③
(45,10,90)
⑬
(90,45,130)
g
(70,45,65)
④
(115,10,40)
⑭
(130,45,65)
h
(15,45,55)
⑤
(45,45,10)
⑮
(90,0,65)
i
(125,35,65)
⑥
(10,45,35)
⑯
(65,90,65)
j
(80,60,80)
⑦
(170,45,65)
a
(130,11.75,100)
Thermocouples installed
in the same place as only
①～⑯ in the case simple
model.
⑧
(155,15,20)
b
(90,45,50)
⑨
(90,45,57.5)
c
(75,45,85)
⑩
(20,50,95)
d
(50,45,120)
(b) Schematic view of experimental set-up with several thermocouples
Figure 6: Schematic view of experimental set-up for evaluation of the proposed virtual elements using the complex
model. In case of the simple model, the green two parts with several holes were removed. In case of the complex
model with forced cooling oil, the inlet and the outlet of the model were remade using two nipples.
Origin point (0,0,0)
①
②
③
④
⑤
⑥
⑦
⑧
⑨
⑩
⑪
⑫
⑬
⑭
⑮
e
b
d
a
c
f
h
g
i
●
x
y
z
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
⑯
●
●
●
●
●
●
The model is
originally
covered by a
plate
● : Thermocouples
j
Oil
x
y
z
The model is
originally
coveredby a
plate (180 ×
130 × 10)
Oil (3ℓ/min or 8 ℓ/min)
Pump
Figure 7: Schematic view of experimental set-up using
the complex model with forced cooling oil
(See Table 1 for the experimental
conditions).
Table 1. Experimental conditions and equipment for
evaluation of proposed
imaginary air element
Air outlet velocity m/s (Voltage
applied to the fan V)
0.5 , 1.0 (3.5 , 5)
Input power to ceramics heater W
5 , 12
Oil inflow rate ℓ/min (Oil inlet
velocity m/s)
3 , 8 (2.6 , 6.9)
Temperature of inflowing air and oil
℃(Oil temperature is controlled by a
thermostat)
20
Ambient temperature ℃
20
Thermocouple
Type T
Distance between fan and
anemometer mm
100
Anemometer (KANOMAX)
6151
Data logger (YOKOGAWA)
MV-230
Slidac (MATSUNAGA)
SD-1310
Power supply (TAKASAGO)
GPV 035-20
Pump (KYOWA)
KYC-300-1
Flowmeter (HORIBA)
LW5-TTN
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
12
Volume 18, 2023
Figure 8: Schematic views of calculation models of complex model
(a) For FEM thermal simulation
(b) For FVM fluid simulation
Imaginary air
Boundary layer
・
・
Nodes
The model is originally
covered by a plate.
αinlet
αoutlet
・
・
□：Area with high heat transfer
coefficient
The model is covered by
a plate unlike the figure
on the left.
Simple model ：9127
Complex model ：18442
Complex model with
cooling oil：30726
Elements
Simple model ：9060
Complex model ：18007
Complex model with cooling
oil：30609
Table 2: Calculation conditions of FVM fluid simulation.
Condition of fluid simulation (ρ：Density，c：Specific heat，λ：Thermal conductivity，μ：Viscosity)
Wall condition
αair = 5 [W/m2・K]
Initial solid temperature
20 ℃
Input power to ceramics heater
5 W , 12 W
Boundary condition
Simple model and complex model
Condition of inlet : Static pressure
Outlet velocity：0.5 m/s , 1.0 m/s
Complex model with cooling oil
Oil inflow rate
：
3 ℓ/min (2.6m/s) , 8 ℓ/min(6.9m/s)
Condition of outlet：Static pressure
Oil temperature entered
Convergence condition of steady
simulation
End of calculation of temperature and heat transfer coefficient or all possible values
Unsteady simulation
All physical time：3600 s , Time step：10 s
Table 3: Calculation conditions of FEM thermal simulation and boundary conditions for proposed
imaginary air element; four virtual elements in Figure 1 were supposed.
Condition of thermal simulation (ρ：Density，c：Specific heat，λ：Thermal conductivity，μ：Viscosity)
Simulation solver
Intel direct sparse
Unsteady simulation
All physical time：3600 s , Time step：10 s
Wall condition
αair = 5 [W/m2・K]
Initial solid temperature
20 ℃，Oil :Oil temperature entered
Temperature of inlet
Input power to ceramics heater P
5 W , 12 W
Boundary condition
Inlet
αinlet
Outlet
αoutlet
Physical
properties of
parts
Imaginary air element
Simple model and complex model
ρ：1.293 kg/m3 , c：1006 J/(Kg・K) , λ：λc
Complex model with cooling oil
ρ：893 kg/m3 , c：2878 J/(Kg・K) , λ：λc
Boundary layer
element
Simple model and complex model
ρ：1.293 kg/m3 , c：1006 J/(Kg・K) , λ：λb
Complex model with cooling oil
ρ：893 kg/m3 , c：2878 J/(Kg・K) , λ：λb
Boundary condition for proposed imaginary air element
FEM
model
Velocity of
liquid V
m/s
αinlet
(Equation 4) W/m2・K
λb
(Equation
10)
W/m
・
K
Vc (Equation 9)
(Vx , Vy , Vz) m/s
λc (Equation 9)
W/m・K
αoutlet (Equation 7)
W/m2・K
Simple
model
Outlet
velocity：
0.5
217
30 ℃:18400,40 ℃:9510,
50 ℃:6500,60 ℃:5050
0.00480
0.0129
(0.0698, 0.0305,
0.0480)
(0.521, 0.228,
0.358)
(13.6, 2.61, 6.41)
(101, 19.4, 47.8)
Outlet
velocity：
1.0
434
30 ℃:36900,40 ℃:19000,
50 ℃:13100,60 ℃:10100
0.00679
0.0199
(0.140, 0.0611, 0.0960)
(1.04, 0.456, 0.717)
(27.1, 5.20, 12.8)
(202, 38.7, 95.5)
Complex
model
Outlet
velocity：
0.5
217
30 ℃:18400, 40 ℃:9510,
50 ℃:6500,60 ℃:5050
0.00480
0.0129
(0.0698, 0.0305,
0.0480)
(0.521, 0.228,
0.358)
(13.6, 2.61, 6.41)
(101, 19.4, 47.8)
Outlet
velocity：
1.0
434
30 ℃:36900,40 ℃:19000,
50 ℃:13100
，
60 ℃:10100
0.00679
0.0199
(0.140 , 0.0611 ,
0.0960)
(1.04, 0.456,
0.717)
(27.1, 5.20, 12.8)
(202, 38.7, 95.5)
Complex
model
with
cooling
oil
Inlet
velocity：
2.6
3.1×
106
30 ℃:1.92×107,40 ℃:9.93×107,
50 ℃:6.83×107, 60 ℃:5.27×107
0.0532
4.06
(0.00284, 0.00649,
0.00446)
(2.55, 5.82, 6.48)
(496, 2590, 1230)
(4.45×105, 2.32×106,
1.78×106)
Inlet
velocity：
6.9
8.4×
106
30 ℃:5.11×108,40 ℃:1.90×108,
50 ℃:1.31×108, 60 ℃:1.01×108
0.0869
7.31
(0.00758, 0.0173,
0.0119)
(6.79, 15.5, 10.7)
(1320, 6910, 3270)
(1.19×106, 6.20×106,
2.93×106)
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
13
Volume 18, 2023
Table 3 shows the analytical conditions of the FVM
fluid simulation and FEM thermal simulation,
respectively.
3.3 Comparison of Experimental and
Calculated Results
Figure 9(a) and Figure 9(b) show examples of the
calculated results (temperature distribution) of FEM
thermal and FVM fluid simulations, respectively. A
simple model is used to show the distribution of
steady-state temperature rise values when 5 W of
power is applied to a single ceramic heater and the air
is expelled from the outlet at a rate of 0.5 m/s.
Figure 10 shows the temperature distribution in
Figure 9. The experimental values corresponding to
the thermocouple positions and the calculated results
of FEM thermal simulation, FVM fluid simulation
(convergence condition for temperature only), and
FVM fluid simulation (convergence condition for all
possible values) are shown, respectively. These
experiments and three types of analyses resulted in a
total of 24 pairs of results (= 3 × 4 × 2), with three
models, four combinations of conditions, and two
types of steady-state and non-steady-state analyses;
Figure 9 and Figure 10 are examples.
Figure 11 shows the results of FEM thermal
simulation and FVM fluid simulation using the
virtual fluid element proposed in this paper,
comparing the analysis time and the error based on
the experimental temperature distribution in the
enclosure, respectively. Here, we evaluate three types
of analysis results based on the results of the previous
Figure 9: Calculation results for temperature change using the simple model with 5 W power and 0.5 m/s
velocity. In case of the FVM fluid simulation, there are the disproportionate temperature distribution
because of the real flow influences.
(b)Temperature change by FVM fluid simulation
(a)Temperature change by FEM thermal simulation
55.0
52.5
50.0
47.5
45.0
42.5
40.0
55.0
52.5
50.0
47.5
45.0
42.5
40.0
Temperature ℃
Temperature ℃
A-A’ Section A-A’ Section
A A’ A A’
-0.14
0.08
0.13 0.14
0.04
0.05
-1
3
7
11
15
19
23
27
31
35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Temperature change ΔT ℃
Experiment
FEM thermal simulation
FVM fluid simulation(Convergence factors：only temperature)
FVM fluid simulation(Convergence factors：all factors )
Inner wall Air Outer wall
Figure 10: Results for temperature change distribution using the simple model with 5 W power and 0.5 m/s velocity.
Calculated results without the adjustment of the boundary conditions using FEM thermal simulation
unexpectedly were suit to the experimental results.
Heater
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
14
Volume 18, 2023
0
50
100
150
200
250
300
0.5
5
0.5
12
1.0
5
1.0
12
0.5
5
0.5
12
1.0
5
1.0
12
2.6
5
2.6
12
6.9
5
6.9
12
0.5
5
0.5
12
1.0
5
1.0
12
0.5
5
0.5
12
1.0
5
1.0
12
2.6
5
2.6
12
6.9
5
6.9
12
FEM thermal simulation with FVM fluid simulation
FVM fluid simulation
Figure 12: Accuracy of analysis at three models focusing on the conditions of radiation and flow (Please see Table
4 and Figure11 regarding “Low”, “High”, “Small” and “Large” for V and P respectively). When the air
velocity was changed from 0.5 m/s to 1.0 m/s, the oil velocity was changed from 2.6 m/s to 6.9 m/s and
the input power to the ceramics heaters were changed from 5 W to 12 W, the calculated results without
the adjustment of the boundary conditions using FEM thermal simulation unexpectedly were nearly to
the experimental results. The FEM thermal simulation using the proposed 4 virtual elements was used
for a design instead of FVM fluid simulation.
FEM thermal simulation with
proposed imaginary air element
FVM fluid simulation (Convergence
values: Only temperature)
FVM fluid simulation (Convergence
values: All values)
A
B
C
A
B
C
A
B
C
A
B
C
V = Low, P = Small
V = Low, P = Large
V = High, P = Small
V = High, P = Large
0.6
1.2
1.2
0
100
200
300
Calculation error Ce %
(A ：Simple model ，B：
Complex model,
C：Complex model with cooling oil)
Figure 11: Calculation time and accuracy of simulation at three models (V：Velocity, P：Input power to ceramics
heater). Calculation times for the FEM thermal unsteady simulation were very longer than that of the FVM
fluid simulation. And calculated results without the adjustment of the boundary conditions using FEM
thermal simulation unexpectedly were nearly to the experimental results. Therefore, the FEM thermal
simulation using the proposed 4 virtual elements was used for a design instead of FVM fluid simulation.
0
100
200
300
400
500
600
Calculation time s
Calculation error Ce %
V [m/s]
P [W]
FEM model
Simple model
Complex model
Complex model
with cooling oil
Simple model
Complex model
Complex model
with cooling oil
Type
Steady simulation
Unsteady simulation
FEM thermal simulation with
proposed imaginary air element
FVM fluid simulation (Convergence
values: Only temperature)
FVM fluid simulation (Convergence
values: All values)
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
15
Volume 18, 2023
24 pairs of experiments. Both steady-state and
unsteady-state conditions are summarized. In the
non-steady state, the results are calculated in 10-
second steps over the 1 hour from the start. The
calculation error was calculated by Equation (11),
using 16 points for the simple model and 26 points
for the complex model, both inside and outside the
enclosure.
Ce = [ Σ (Tc-n－Te-n)/Te-n ]÷(16 or 26)× 100 (11)
where Ce is the calculation error [%], Te-n is the
experimental temperature rise at n points in the
enclosure [°C], and Tc-n is the calculated temperature
rise at n points in the enclosure [°C].
For steady-state conditions, regarding calculation
time, the FEM thermal simulation using the proposed
virtual fluid element was shorter than the FVM fluid
simulation by about 4/9 to 1/285. The fluid analysis
of a complex model using a cutting oil agent did not
converge after more than 24 hours of calculation, so
the number of elements was reduced from 30726 to
7667, resulting in convergence. The calculation errors
for the FEM thermal simulation with virtual fluid
elements are 7/8 to 1/9 of one of the FVM fluid
simulations. The reason why the steady-state analysis
takes longer for the FVM fluid simulation than for the
FEM thermal simulation is thought to be due to the
difference in the discretization of the analysis. For the
unsteady state, regarding calculation time, the FEM
thermal simulation with the proposed virtual fluid
element takes about 2.1 to 3.1 times longer than the
FVM fluid simulation. This may be due to the time
required to calculate the convergence of the four
virtual fluid elements. The computational errors of
FEM thermal simulations with virtual fluid elements
are 6/7 to 1/9 of one of the FVM fluid simulations. It
can be concluded that the proposed FEM thermal
simulation using virtual fluid elements can be
calculated accurately for both steady-state and
unsteady-state analyses, and especially for steady-
state analyses, the analysis time is shorter than that of
FVM fluid simulation. In general, a general
evaluation is often made by steady-state thermal
analysis, as exemplified by the design of a heat sink
inside a personal computer. Therefore, the proposed
FEM thermal simulation with virtual fluid elements is
considered to be very effective in terms of both
analysis time and calculation accuracy.
Figure 12 shows a comparison of the calculated
steady-state temperature rise values in Figure 11,
adjusted for the high and low inflow and outflow
velocities of the cooling medium and the large and
small input power of the ceramic heater; the
calculated results showed similar trends for three
types of experiments (A; simple model + air, B;
complex model + air, and C; complex model +
cutting oil) without being affected by the high and
low inflow and outflow velocities of the cooling
medium or the large and small input power of the
ceramic heater. It can be seen that the FEM thermal
simulation using virtual fluid elements is more
accurate than the FVM fluid simulation, about 7/8 to
1/9. In addition, the FEM thermal simulation
proposed in this paper can accurately reflect the
effect of the temperature change of the enclosure due
to the difference in the characteristics of laminar-
turbulent flow due to the increase in the heating value
of the heat source and the change in the velocity of
the medium.
4 Considerations for Applying the
Proposed Method to Actual Machine
Tools
In this chapter, the four proposed virtual fluid
elements will be used in the future to study the
phenomenon of heat buildup in the actual machine
tool structure and the FEM thermal simulation of
forced air intake and exhaust into and out of the
machine structure.
As shown in Table 4, a machine tool consists of a
main structure such as a spindle head and bed, as well
as thermal-volumetric spaces (TVS) such as safety
enclosures and piping. Heat sources in machine tools
include spindle bearings, ball screws and ball nuts,
linear guides, and motors, as well as heat generated
by cutting tools and workpieces, accumulated chips,
and cutting fluid. These heat sources mainly conduct
heat within the machine structure, causing thermal
deformation and loss of machining accuracy. Apart
from that, these heat sources also transfer heat to the
fluid (mainly air) in the enclosure or thermal volume
space TVS, causing the fluid to rise in temperature,
which again transfers heat to the machine structure,
resulting in complex thermal deformation and
reduced machining accuracy. This also causes heat
buildup, which is counteracted by forced air intake
and exhaust into the machine structure.
In this research, the four proposed virtual fluid
elements were evaluated using simple models. When
applying the elements to actual machine tools with
complex specifications (shape, size, material, number
of machine elements, etc.), the FEM models can be
easily created using current CAD software, and
complex thermal conditions can be easily set.
Therefore, it is easy to set up I: Simulated element for
16 or 26
n=1
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
16
Volume 18, 2023
Table 4: Considerations for applying the proposed four virtual models to real machine tools.
Structures containing fluids in machine tools. Sources of heat generated during
machine tool operation
Thermal effects of the sources
・Main structures as a machine tool element (Head
stock, Bed and so on), and complex ribbed structures
in main structures
・Accessories (NC controller, oil tank and so on)
・Enclosures for safety
・Thermal-Volimetric Spaces (Other spaces such as
plumbing)
・Spindle bearings, ball screws and
ball nuts, linear guides and motors
for operation
・Workpieces, Tools and Chips for
heat generated during cutting
・
Cutting fluids for forced cooling (as
a heat source here)
・Heat conduction to machine structure →
Thermal deformation → Reduced
accuracy
・Heating the fluid in the enclosures and
the thermal-volimetric spaces → The
deformation for the fluid→ Reduced
accuracy
inflow and II: Simulated element for outflow out of
the four proposed virtual fluid elements. In addition,
as mentioned above, III: Simulated element for the
boundary layer can be easily created using both the
cavity and shell functions of CAD, and IV: Simulated
element for convection can be easily created using
the cavity function of CAD. To perform highly
accurate FEM thermal simulations, it is necessary to
know the characteristic values of the four virtual
models in advance.
5 Conclusion
(1) A method to calculate the temperature distribution due
to fluid behavior in a structure was proposed by FEM
thermal simulation using four virtual fluid elements; I:
inflow, II: outflow, III: boundary layer, and IV:
convection.
(2) In the steady-state analysis, the FEM thermal
simulation using the proposed virtual fluid element
was shorter than the FVM fluid simulation by about
4/9 to 1/285. The computational errors are 7/8 to 1/9
of those of the FVM fluid simulation, which is more
accurate. For the unsteady state, the FEM thermal
simulation with the proposed virtual fluid element
takes about 2.1 to 3.1 times longer than the FVM
fluid simulation. The computational errors are about
6/7 to 1/9 of those of the FVM fluid simulation.
References:
[1] A. Itou, T. Nakanishi, T. Saburi, S. Kubota, Y.
Ogata, High Performance Parallel Computing for
Computational Fluid Dynamics (CFD) -Second
Report, Komatsu technical report, Vol.51,
No.156, 2007, pp. 21-26 (in Japanese).
[2] J. Jedrzejewski, W. Modrzycki, Compensation
of Thermal Displacement of High-speed Precision
Machine Tools, Journal of Mechanical
Engineering, Vol. 7, No. 1, 2007, pp. 108-114
[3] Z. Winiarski, Z. Kowal, J. Jedrzejewski, Precise
Modelling of Machine Tool Drives with Ball
Screw Thermal Behaviour, Journal of
Mechanical Engineerig, Vol. 17, No. 1, 2017, pp.
31-45.
[4] RENESAS，“Mechanism of Heat radiation”
https://www.renesas.com/jajp/support/
technicalresources/package/characteristic/heat-
01.html [Accessed on 25 February, 2023].
[5] J. Jedrzejewski, W. Kwasny, Z. Kowal, Z.
Winiarski, Development of the Modelling and
Numerical Simulation of Thermal Properties of
Machine Tools”, Journal of Machine
Engineering, Vol.14, No. 3, 2014, pp.5-20.
[6] J. Glanzel, T. S. Kumar, C. Nauman, M. Puntz,
Parameterization of Environmental Influences by
Automated Characteristic Diagrams for the
Decoupled Fluid and Structural-machine
Simulations”, Journal of Machine Engineering,
Vol.19, No. 1, 2019, pp.98-113.
[7] Z. Lei, Y. Nagata, Development for a parallelized
CFD code-ADCS, JAXA Research and
Development Report, JAXA-RR-09-006, 2010,
pp. 1-5.
[8] H. Sugiyama, M. Sano, Y. Nagahasi, N. Kato,
Transfer Phenomenology to Learn for the First
time - to Understand Momentum, Heat, Mass
Transfer Integratedly – , Morikita Publication,
2014, pp. 110-155 (in Japanese).
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.2
Ikuo Tanabe, Hiromi Isobe
E-ISSN: 2224-3461
17
Volume 18, 2023
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The authors equally contributed in the present
research, at all stages from the formulation of the
problem to the final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
No funding was received for conducting this study.
Conflict of Interest
The authors have no conflicts of interest to declare
that are relevant to the content of this article.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US
