- Email: [email protected]
Numerical study on a heat transfer model in a Lagrangian fluid dynamics simulation
International Journal of Heat and Mass Transfer 103 (2016) 635–645
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Numerical study on a heat transfer model in a Lagrangian fluid dynamics simulation Kazuya Takabatake a, Xiaosong Sun a, Mikio Sakai b,⇑, Dimitrios Pavlidis c, Jiansheng Xiang c, Christopher C. Pain c a b c
Department of Nuclear Engineering and Management, School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Resilience Engineering Research Center, School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Applied Modeling and Computation Group, Department of Earth Science and Engineering, Imperial College London, Prince Consort Road, London SW7 2BP, UK
a r t i c l e
i n f o
Article history: Received 14 April 2016 Received in revised form 5 July 2016 Accepted 20 July 2016
Keywords: Moving particle semi-implicit method Lagrangian particle method Heat transfer Phase transition
a b s t r a c t Phenomena related to phase change heat transfer are often encountered in engineering. These phenomena are regarded to be complex, since not only phase transition from solid to liquid occurs but also movement of fluid interface has to be taken into consideration. Detailed numerical modeling of these complex systems is essential to better understand them and optimize industrial designs. Lagrangian methods are promising for simulating such complex systems. The Moving Particle Semi-implicit (MPS) method, which is one of the Lagrangian methods, is employed here to simulate the free surface fluid flows involving heat transfer and phase change. On the other hand, the existing MPS method could not apply Neumann boundary condition such as heat flux in the heat transfer simulations. This is because the surface direction could not be readily defined on the surface of the spherical fluid particles in the MPS method. Hence, prescribing the heat fluxes becomes problematic in the existing MPS method. To solve this problem, a new heat flux model is developed, where the divergence operator is applied in the heat transfer simulation. Simple verification tests are performed to demonstrate the heat flux model, where the calculation results are compared against analytically derived solutions. In addition, application of the signed distance function is also investigated in the heat transfer simulation for arbitrary shaped boundary. In simple verification tests, the computation results are shown to agree well with the analytical solutions. Consequently, adequacy of the novel heat transfer model developed here is shown in the Lagrangian fluid dynamics simulation. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction Phase change heat transfer is often encountered in various engineering [1–3] such as chemical engineering, food engineering, mechanical engineering, nuclear engineering, resilience engineering, etc. At present, introduction of numerical technologies are desired for accurate product design and/or better understanding of the phenomena. When phase change heat transfer phenomena are numerically simulated, not only heat transfer model but also interfacial fluid flow modeling is required. For modeling the interfacial fluid flow, the Volume-of-Fluid (VOF) [4–8], Smoothed Particle Hydrodynamics (SPH) [9–12] and Moving Particle Semiimplicit (MPS) method [13–15] have been employed in the computations so far. Lagrangian methods such as the SPH and MPS method are regarded to be accurate in simulating free surface fluid ⇑ Corresponding author. Fax: +81 3 5841 6977. E-mail address: [email protected] (M. Sakai). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.07.073 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
flows, since convection term is not required in the modeling. Accordingly, the Lagrangian approach becomes important in various engineering. Recently, weakly compressible MPS method [16,17] was developed for efficient calculation of the interfacial flows, where the pressure computation was performed explicitly. Hence, the calculation efficiency of the MPS method has been improved substantially and is now capable of industrial application. The MPS method was applied to complex thermal hydraulic systems such as a simulation of a severe accident in nuclear engineering [18,19]. Although the MPS method seems to be established, there are some issues in the heat transfer simulations. Specifically, in the existing MPS method, Neumann boundary condition such as heat flux boundary could not be set. This is because the surface direction could not be evaluated easily at the surface of MPS particles, and hence heat flux could not be set at the boundary. Only Dirichlet boundary could be used in the existing MPS method [18–21], and thus only the temperature could be set at the boundary. In addition, in the existing MPS method, accuracy
636
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Nomenclature C c cp D g h k L l N n n0 P Q q r r re T Tm t u w wg
Greek symbols variable [–] c phase fraction [–] C constant for surface detection [–] k Laplacian constant in the MPS method [–] l viscosity [Pa s] m fluid kinematic viscosity [m2/s] s viscous tensor q fluid density [kg/m3]
empirical constant for viscosity [–] sound speed [m/s] specific heat [J/(kg K)] dimension number [–] gravitational acceleration [m/s2] enthalpy [J/m3] thermal conductivity [W/(m K)] latent heat [J/kg] length [m] unit void vector [–] particle number density [–] initial particle number density [–] pressure [N/m2] heat source [J/m3] heat flux vector [W/m2] coordinate vector [m] distance [m] cutoff radius [m] temperature [K] melting point [K] time [s] fluid velocity [m/s] basic kernel function [–] special kernel function for pressure gradient [–]
U
Subscript i, j L S Wall Max 1,2
particle number liquid phase solid phase property of wall maximum property type number
Superscript ⁄ auxiliary value
of the heat transfer calculation is known to be worse in arbitrary shape wall systems. This is because the wall boundary was created by stepwise location of the MPS particles. Consequently, the existing MPS method has problems in heat transfer simulation. To solve the above issues, novel heat transfer model is developed in the MPS method. Specifically, heat flux model is newly developed in the current study. The weakly compressive MPS method is used in this study. To apply the MPS simulation of heat transfer in an arbitrary shaped domain, the signed distance function (SDF) is introduced for the first time. In addition, the implicit calculation algorithm [22] is newly introduced into the MPS method to simulate highly viscous fluid flow involving heat transfer efficiently. Verification tests are performed in a simple system to prove adequacy of the new models. Calculation results agree well with theoretical ones. Moreover, the numerical model is shown to simulate phase change system. Thus, the new heat transfer model is verified in Lagrangian fluid dynamics. 2. Methodology
Dh ¼ r ðkrTÞ þ Q ; Dt
ð3Þ
where h, k, T and Q are the enthalpy, thermal conductivity, temperature and energy inflow, respectively. 2.2. The MPS method In the MPS method [13], the differential operator is modeled by the weighted average between neighbor particles. A basic kernel function is given [16,17] as
( wðrÞ ¼
r re
þ rre 2 r < r e
0
r P re
;
ð4Þ
where r and re are distance and cutoff radius of interaction, respectively. Assuming that virtual particles are set in a structurally aligned particle gird, the initial particle number density and Laplacian constant are given as
n0 ¼ 2.1. The governing equations
X
wðjrj ri jÞ;
ð5Þ
j–i
For the weakly compressible fluid motion, the mass and momentum conservation equations can be written as
Dq þ qr u ¼ 0; Dt
ð1Þ
Du 1 ¼ rP þ r s þ g; Dt q
ð2Þ
where q, t, u, P, s and g are the fluid density, time, fluid velocity, pressure, viscous tensor and gravity, respectively. The s is given by
s ¼ r u þ r uT : Considering heat transfer, the energy conservation equation is given as follows:
k¼
1X jr j r i j2 wðjrj ri jÞ: n0 j–i
ð6Þ
These values are constant during the calculation. The particle density can be approximated to be proportional to the initial particle number density as follows:
ni ¼
X wðjrj ri jÞ;
ð7Þ
j–i
qi ¼
ni q: n0 0
ð8Þ
In the weakly compressible approach [16,17], pressure can be derived as a function of density and is given by
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
8 < q 0 c 2 n i 7 1 ni P n0 0 7 n : P¼ :0 ni < n0
2.3. Heat transfer model
ð9Þ
And by using the MPS gradient model, the pressure gradient term in Eq. (2) is derived as
rP i ¼
D X P i þ Pj rj ri wg ðr ij Þ : ng j–i jr j ri j jrj ri j
ð10Þ
Here, wg is a special kernel function for pressure gradient and is given by
wg ðrÞ ¼
r re
0
þ rre
r < re r P re
:
637
ð11Þ
Compared with conventional mesh methods, the heat transfer model is very limited in the MPS method. It is modeled by the same manner in the SPH method [24]. The energy conservation equation can be discretized in the MPS method as follows:
Dh 2D X 2ki kj ¼ ðT j T i Þwðjr j ri jÞ þ Q : Dt kn0 j–i ki þ kj
ð15Þ
In this study, two different properties are used to simulate the phase transition based on the enthalpy according to the past studies [18]. The density variation between solid and liquid phase is assumed to be small enough to be negligible compared with the variation between the liquid and gas phases. The temperature is derived as a function of enthalpy and is given by:
Besides, ng is calculated by using Eq. (11) in the same manner as n0 and makes it possible to discretize the pressure gradient term while the momentum remains constant. Eq. (2) can be decomposed into the following equations:
8 > T þ hhs > < m qcp T ¼ Tm > > : T m þ hhl
u u k 1 ¼ g rP; Dt q
ð12Þ
ukþ1 u ¼ r skþ1 : Dt
ð13Þ
where Tm, hs, hl and cp are the melting point, enthalpy of solid, that of liquid and specific heat, respectively. As shown in Fig. 1, temperature during phase transition is assumed to be constant. Enthalpy of the solid phase hs and that of the liquid phase hl are expressed as follows:
This is the similar to the fractional time step method [22] except for the pressure term because the value of pressure gradient is already known in the weakly compressible MPS method. Here, Eq. (13) is modeled by the same manner in previous studies
rs¼
2d X 2mi mj ðu u Þwðr r Þ : j i j i kn0 j–i mi þ mj
ð14Þ
qcp
h < hs hs 6 h 6 hl ;
ð16Þ
h > hl
hs ¼ qcp T m ;
ð17Þ
hl ¼ hs þ qL;
ð18Þ
where L is latent heat. In order to simulate phase transition, the phase fraction coefficient g is introduced and is calculated as:
c¼
8 > <1 > :
hl h hl hs
0
h < hs hs 6 h 6 hl :
ð19Þ
h > hl
This modeling has been successfully applied in various SPH and MPS simulations [18,19,23,24]. The auxiliary velocity u⁄ in Eq. (12) is estimated explicitly from Eqs. (9) and (10). The updated velocity uk+t in Eq. (13) is then calculated implicitly by using conjugate gradient method.
As shown in Fig. 2, the phase is determined from enthalpy and is interpolated by the linear function during phase transition.
Fig. 1. Temperature model.
Fig. 2. Phase change model.
638
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
In the current study, the material was modeled by a fluid. When the material is in solid state, the fluid viscosity was set to be sufficiently high. The viscosity is a function of enthalpy. Hence, the modeling of viscosity is indispensable for the heat transfer simulation. In this study, the viscosity model is expressed by the Arrhenius viscosity formula as:
lr ¼ ll expð2:5C cÞ;
ð20Þ
where lr, ll and C are the relative viscosity, viscosity of liquid phase and empirical constant, respectively. This empirical constant varies between 4 and 8 [25]. Eq. (20) makes it possible to simulate the behavior of a rigid body as an extremely high viscous liquid. 2.4. Boundary conditions 2.4.1. Dirichlet boundary conditions Heat transfer is calculated by the Laplacian of the temperature between neighboring particles as shown in Eq. (3). When a Dirichlet boundary condition is applied, the temperature at the boundary is fixed. Thus, the source term (the second term in the right hand side of Eq. (3)) is given by:
Q ¼ 0:
ð21Þ
For Dirichlet boundary conditions, no special treatment of surface detection is required. 2.4.2. Neumann boundary conditions In the case of Neumann boundary conditions, the heat flux at the boundary is fixed. In the meshless particle method, it is difficult to evaluate the surface area and its normal direction. The heat flux can be set as a volumetric heat source by introducing the divergence theorem. As such the heat fluxes at boundaries are modeled by:
Q ¼ r q;
ð22Þ
where q is heat flux vector. In the Neumann boundary condition, surface must be detected depending on the direction of heat flux vector. First, an ordinary surface detection is applied for each particle as shown in Eq. (23):
ni < Cn0 ;
ð23Þ
where C is a constant for surface detection and is set to 0.97 in this study as in the previous study [13,14]. Then, the void vector N, which is perpendicular to the surface, is calculated for each surface particle detected by Eq. (24) as follows:
P j–i
ðr j ri Þwðr j ri Þ jr r j
j i : N ¼ P ðr r Þwðr r Þ j–i j jri j ri jj i
ð24Þ
Finally, assuming that a conceptual particle is set at the head of the void vector, the surface where heat flux enters can be detected by:
N q < 0:
ð25Þ
tance from walls). One can see the details of the SDF in our previous studies. The SDF could be applied to a simulation of a granular flow in a complex shape boundary system such as a twin-screw kneader [29], a die-filling [30] and a ribbon mixer [31]. Adequacy of the SDF based wall boundary was demonstrated through validation tests. In the MPS method, the distance from walls is calculated at all the sampling points in the calculation domain at the beginning of each step. Assuming that the virtual particles are aligned as the structural grid in the wall, the value of the kernel function from the wall at each sampling point can be calculated from the distance. The value of the kernel function from wall at each particle is then calculated by interpolation using information of neighboring sampling points. Application of the SDF to simulations of solid– liquid flows has already been presented in our past studies [16,17], where the MPS method was employed to simulate the liquid phase. In this study, the SDF was introduced into the heat transfer simulation of the MPS method for the first time. 3. Numerical simulations 3.1. Conditions Three kinds of verification tests were performed to show the adequacy of the heat flux model. Besides, the applicability of the SDF to MPS method was proved in the heat transfer simulation. In the verification tests, a simple shape system such as a rectangular bar was chosen because of comparison with the analytical solution. Usage of the simple shape system is of great importance because the new technology is shown to be applied to a great deal of systems in engineering. Our objective is to show the agreement between the calculation and analytical results. Accordingly, applicability of the SDF was investigated in a simple system. In these verification tests, the temperature distribution was compared with the theoretical solution. In addition, the phase transition simulation was also performed due to heat flux by using a Neumann boundary condition. In the verification tests, time step, speed of sound and the empirical constant related to viscosity were, respectively, 1.0 105 s, 100 m/s and 8.0. The fluid particle size, which is so-called initial particle distance, was 5.0 mm. Value of the parameters used in the verification tests are summarized in Table 1. 3.1.1. Case 1: Dirichlet boundary condition Case 1 was performed to show application of the SDF in the heat transfer simulation. Heat transfer in the rectangular bar was evaluated in the verification test. Temperature distributions were compared with the calculation results between a particle based wall (Case 1-1) and the SDF based wall (Case 1-2). Fig. 3 shows the calculated particle arrangement. Dirichlet boundary condition was set at both the top and the bottom wall and adiabatic boundary at the side wall. In Case 1-1, the wall boundary was given by wall particles shown in Fig. 3(a). In Case 1-2, the wall was created by the SDF
By using the algorithm outlined from Eqs. (23)–(25), entrance boundaries can be detected. 2.4.3. Signed distance function Originally the SDF [26,27] was developed for granular flow simulations using the discrete element method [28]. The SDF is defined by sign and minimal distance from the calculation domain. The SDF value becomes positive inside the calculation domain, and vice versa. In the DEM, SDF makes it possible to relate displacement for the contact force and the solid particle location (or dis-
Table 1 Simulation conditions. Item
Unit
Value
Time step Gravitational acceleration Sound speed Particle diameter Empirical constant for viscosity
s m/s2 m/s mm –
1.0 105 9.8 100 5.0 8.0
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
(a) Case 1-1: particle wall
639
(b) Case 1-2: SDF based wall
Fig. 3. Initial arrangement of Case 1.
Table 2 Physical properties (Case 1). Item
Unit
Material 1
Density Kinematic viscosity Specific heat Thermal conductivity
kg/m3 m2/s J/(kg K) W/(m K)
1000 1.0 103 10 50
as designated in Fig. 3(b). The temperature was given at 300 K and 400 K at the top and bottom walls. The initial temperature was set to be 300 K in the bar. The bar was composed of one material, where density, kinematic viscosity, specific heat and thermal conductivity are prescribed as 1000 kg/m3, 1.0 103 m2/s, 10 J/(kg K) and 50 W/(m K), respectively. Note that the value of kinematic viscosity of the solid state became 5.0 108 times larger than that of the liquid state by using Eq. (20). Table 2 shows the physical properties. Objective of this verification is to show the adequacy of the heat transfer modeling, and hence these properties were arbitrarily chosen. Especially, heat capacity was extremely low in order to reduce the computational time. 3.1.2. Case 2: Neumann boundary condition without phase transition Case 2 was performed to show adequacy of the heat flux model in the MPS method. Fig. 4 shows the calculation systems. The geometrical configuration of the bar was the same as that in Case 1. To investigate the effects of material properties, a single bar composed of material-1 and a combined bar composed of material-1 and material-2 were used in Cases 2-1 and 2-2. In both cases, the bottom wall was defined by the SDF. Neumann boundary condition, namely, heat flux, was
applied at the top surface. Dirichlet boundary condition, namely, 300 K of constant temperature, was given at the bottom SDF wall. Adiabatic boundary condition was given at the side SDF wall. The intensity of heat flux was 25,000 W/m2. The initial temperature was 300 K in the bar. In the materials, the density, kinematic viscosity and specific heat were 1000 kg/m3, 1.0 103 m2/s and 10 J/(kg K). On the other hand, the thermal conductivities were different in Cases 2-1 and 2-2. In Case 2-2, the different thermal conductivity was 50 W/(m K) and 500 W/ (m K) bottom half and top half, respectively. Table 3 summarizes the physical properties of Case 2.
3.1.3. Case 3: Neumann boundary condition with phase transition Case 3 was performed to simulate phase change material under the heat flux model by the MPS method. Fig. 5 illustrates the calculation systems. The bar was composed of two different materials as in Case 2. In both cases, bottom wall was expressed by the SDF. Neumann boundary condition, namely, heat flux, was applied at the top surface. The intensity of heat flux was 100,000 W/m2 in order to finish the calculation in short time. Adiabatic condition and constant temperature were given at both the side and the bottom in Cases 3-1 and 3-2. The initial temperature was 300 K in the bar. As far as the density, kinematic viscosity and latent heat were concerned, common values were used in Cases 3-1 and 3-2, namely, they were 1000 kg/m3, 1.0 103 m2/s and 1.0 104 J/ kg, respectively. On the other hand, specific heat, thermal conductivity and melting point were given by different values. The specific heat, thermal conductivity and melting point were 100 J/(kg K), 50 W/(m K) and 10,000 K in material 1. In material 2, the specific heat, different thermal conductivity and melting point were 10 J/ (kg K), 500 W/(m K) and 320 K, respectively. Table 4 shows the physical properties in Case 3.
640
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
(a) Case 2-1: one material
(b) Case 2-2: two materials
Fig. 4. Initial arrangement of Case 2. Table 3 Physical properties (Case 2).
Table 4 Physical properties (Case 3).
Item
Unit
Material 1
Material 2
Density Kinematic viscosity Specific heat Thermal conductivity
kg/m3 m2/s J/(kg K) W/(m K)
1000 1.0 103 10 50
1000 1.0 103 10 500
(a) Case 3-1: adiabatic
Density Kinematic viscosity Specific heat Thermal conductivity Melting point Latent heat
kg/m3 m2/s J/(kg K) W/(m K) K J/kg
(b) Case 3-2: Dirichlet boundary
Fig. 5. Initial arrangement of Case 3.
Material 1
Material 2
1000 1.0 103 100 50 10,000 1.0 104
1000 1.0 103 10 500 320 1.0 104
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
0s
1s
2s
641
3s
Fig. 6. Snapshot of heat transfer simulation in Case 1-1.
ature was gradually increased upside due to Dirichlet condition. The temperature distribution reached steady state after around 3 s. Fig. 9 shows the transient change of temperature distribution in Case 1-2. The maximum temperature of the bar became 399.9 K at the particle on the top boundary. The simulation results were compared with the theoretical solution. The theoretical solution can be easily obtained in this problem. This is because the temperature can be linearly interpolated between the walls. The temperature can be conducted theoretically to be 398.8 K in the calculation system, when the sampling point was same as that of the simulations. In spite of such a simple case, temperature distribution obtained in Case 1-1 was slightly different from that of Case 1-2. This is due to the deformation by the gravitational force when the particle wall boundary was applied. The deformation may increase or decrease the number of wall particles inside cutoff radius, because the coordinates of wall particles are fixed during the calculation. On the other hand, the value of kernel function was determined by the distance from wall when the SDF was used. Therefore, the deformation was not intrinsically an important problem. Nevertheless, the difference of calculation results was less than 1% from theoretical one. The error might be occurred due to the spatial resolution. Hence, the temperature between the simulations and theoretical solution was in good agreement. Besides, application of the SDF was also proved through the verification tests.
4. Results and discussion 4.1. Case 1: Dirichlet boundary condition In order to show application of the SDF in the heat transfer simulation, Case 1-1 and Case 1-2 were performed under Dirichlet condition. As mentioned above, in Case 1-1, the wall boundary was created by the MPS particles. Fig. 6 illustrates the snapshot of the heat transfer simulation in Case 1-1. The temperature was gradually increased upside due to Dirichlet condition. The temperature distribution became steady state in 3 s. Fig. 7 shows the transient change of temperature distribution in Case 1-1. The maximum temperature was 396.3 K at the particle on the top boundary of the bar. In Case 1-2, the wall was generated by the SDF. Fig. 8 designates the typical snapshot of the heat transfer simulation in Case 1-2. By the same manner as Case 1-1, the temper-
4.1.1. Case 2: Neumann boundary condition without phase transition In order to demonstrate the heat flux model by the MPS method, Cases 2-1 and 2-2 were performed under the Neumann boundary condition. In Case 2-1, the bar composed of a single material. Fig. 10 illustrates a snapshot of the heat transfer simulation in Case 2-1. The temperature was gradually increased from the top due to the heat flux. The temperature distribution became
Fig. 7. Transient change of temperature distribution in Case 1-1.
0s
1s
2s
3s
Fig. 8. Snapshot of heat transfer simulation in Case 1-2.
642
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
ture distribution in Case 2-2. The maximum temperature of the bar became 350.1 K. The simulation results were compared with the theoretical solution. The maximum temperature of each case can be calculated theoretically by using the following equations:
T wall T max ; l
ð26Þ
T wall T max ; l1 =k1 þ l2 =k2
ð27Þ
q ¼ k
q¼
Fig. 9. Transient change of temperature distribution in Case 1-2.
steady state in 20 s. Fig. 11 shows the transient change of temperature distribution in Case 2-1. The maximum temperature of the bar was 392.6 K at the particle on the top boundary of the bar. The linear temperature gradient was obtained. In Case 2-2, the bar composed of material-1 and material-2. Fig. 12 illustrates the snapshot of the heat transfer simulation in Case 2-2. As in Case 2-1, the temperature was gradually increased from the top due to the heat flux. The temperature distribution became steady state after around 20 s. Fig. 13 shows the transient change of tempera-
where l is length. Substituting the physical properties in Eqs. (26) and (27), the theoretical solution of the maximum temperature of the bar in Case 2-1 and Case 2-2 were 398.8 K and 354.9 K, respectively. The error became around 1.5 % in Cases 2-1 and 2-2. The error might be occurred due to the spatial resolution. In addition, the error of Case 2-2 was slightly larger than that of Case 2-1. This is because the heat input varied from place to place such as at the surface, at the side and the vertex. Through the verification tests, the application of the heat flux model was proved. 4.2. Case 3: Neumann boundary condition with phase transition Case 3 was performed to simulate phase change under the heat flux model by the MPS method. Case 3-1 was performed to simulate melting in an adiabatic environment. Case 3-2 was performed
0 sec
1 sec
2 sec
5 sec
10 sec
20 sec
Fig. 10. Snapshots of Case 2-1.
3 sec
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Fig. 11. Transient change of temperature distribution in Case 2-1.
643
Fig. 13. Transient change of temperature distribution in Case 2-2.
0 sec
1 sec
2 sec
5 sec
10 sec
20 sec
3 sec
Fig. 12. Snapshots of Case 2-2.
to simulate the resolidification of the molten material under Dirichlet condition. Fig. 14 illustrates a snapshot of the heat transfer simulation in Case 3-1. The temperature was gradually increased from the top due to heat flux. The bar began to melt in 2 s because the local temperature exceeded the melting temperature. Upper part was almost molten in 10 s. All the material on the top half was molten completely in 20 s. The molten fluid was distributed on the floor. Hence, melting phenomena due to heat flux could be simulated qualitatively by the proposed method. Fig. 15 illustrates the snapshot of the heat transfer simulation in Case 3-2. The temperature was grad-
ually increased from the top due to the heat flux. As with Case 3-1, the bar began to melt in 2 s because the local temperature exceeded the melting temperature. The molten fluid re-solidified at the floor due to rapid cooling. All the material on the top half of the bar melted completely in 20 s. The molten fluid was distributed near the bar unlike Case 3-1. Hence, resolidification could be simulated qualitatively by the proposed method. These results are valuable because the melting and solidification under heat flux could be simulated by the MPS method for the first time. At the early stage of the melting phenomenon in Figs. 14 and 15, the represented processes seem to be symmetric. On the other
644
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Fig. 14. Snapshots of Case 3-1.
0 sec
1 sec
2 sec
4 sec
5 sec
10 sec
3 sec
Fig. 15. Snapshots of Case 3-2.
hand, after 3 s, they seem to be asymmetric in both cases. It is well known that the pressure oscillation occurred in the MPS method [32]. This pressure oscillation might make the fluid flow asymmetric.
the finite difference method not only at the equilibrium state but also at the non-equilibrium state. The adequacy of our divergence model was shown through the investigation. 5. Conclusions
4.3. Convergence study regarding spatial resolution As already mentioned, the MPS method was in good agreement with the theoretical solution by using our divergence model. In this section, resolution effect on the divergence model was examined by comparing with the finite difference method. Fig. 16 shows the verification results where spatial resolution was investigated in the MPS method. As shown in this figure, the calculation error was improved due to the spatial resolution in the MPS method. Besides, it was shown that this model was in good agreement with
Lagrangian approaches such as the SPH and MPS methods are often employed in simulating free surface fluid flows. In recent years, the MPS method has been used for complex fluid flows such as free surface fluid flows involving heat transfer. In the existing Lagrangian approach, heat flux could not be set because of difficulty of surface direction definition. In addition, creation of arbitrary shape wall boundary was also problematic. This is because the MPS particles should be located in heat transfer simulation in the existing MPS method. In order to solve these problems, a
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Fig. 16. Comparison between the difference method and the MPS method in Case 2.
new heat transfer model was developed in the MPS method. Moreover, in this study, application of the SDF to simulate heat transfer was investigated for the first time. In order to introduce heat fluxes into the MPS method Neumann boundary condition was modeled by using a divergence model. Adequacy of this modeling was demonstrated using a simple system by comparing the calculation results against the analytical solutions. Employing the simple shape system is important because the new model is shown to be used to a great deal of systems in engineering. Application of the SDF to the MPS simulation of heat transfer was verified for the first time in this study. The verification test was performed in a simple system. The temperature obtained from the SDF was in good agreement with that of the particle based wall boundary. Finally, it was shown that melting material due to heat fluxes and resolidification could be simulated successfully by this new technologies. Consequently, innovative heat transfer model was developed in a Lagrangian approach. Acknowledgement This study was financially supported by the Initiatives for Atomic Energy Basic and Generic Strategic Research by the Ministry of Education, Culture, Sports, Science, and Technology of Japan. References [1] G. Quarini, Y. Chang, Heat transfer characteristics of ice melting in water and salt solutions, Chem. Eng. Res. Des. 80 (2002) 320–324. [2] K.W. Chua, Y.T. Makkawi, M.J. Hounslow, Time scale analysis for fluidized bed melt granulation III: binder solidification rate, Chem. Eng. Sci. 66 (3) (2011) 336–341. [3] Y. Allouche, S. Varga, C. Bouden, A.C. Oliveira, Validation of a CFD model for the simulation of heat transfer in a tubes-in-tank PCM storage unit, Renew. Energy 89 (2016) 371–379.
645
[4] K. Yokoi, Efficient implementation of THINC scheme: a simple and practical smoothed VOF algorithm, J. Comput. Phys. 226 (2) (2007) 1985–2002. [5] X. Sun, M. Sakai, Numerical simulation of two-phase fl ows in complex geometries by using the volume-of-fluid/immersed-boundary method, Chem. Eng. Sci. 139 (2016) 221–240. [6] X. Sun, M. Sakai, Three-dimensional simulation of gas–solid–liquid flows using the DEM–VOF method, Chem. Eng. Sci. 134 (2015) 531–548. [7] D. Pavlidis, Z. Xie, J.R. Percival, J.L.M.A. Gomes, C.C. Pain, O.K. Matar, Two- and three-phase horizontal slug flow simulations using an interface-capturing compositional approach, Int. J. Multiphase Flow 67 (S) (2014) 85–91. [8] Z. Xie, D. Pavlidis, J.R. Percival, J.L.M.A. Gomes, C.C. Pain, O.K. Matar, Adaptive unstructured mesh modelling of multiphase flows, Int. J. Multiphase Flow 67 (S) (2014) 104–110. [9] J.J. Monaghan, An introduction to SPH, Comput. Phys. Commun. 48 (1988) 89– 96. [10] X. Sun, M. Sakai, Y. Yamada, Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method, J. Comput. Phys. 248 (2013) 147–176. [11] A. Eitzlmayr, J. Khinast, Co-rotating twin-screw extruders: detailed analysis of conveying elements based on smoothed particle hydrodynamic. Part 1: Hydrodynamics, Chem. Eng. Sci. 134 (2015) 861–879. [12] A. Eitzlmayr, J. Khinast, Co-rotating twin-screw extruders: detailed analysis of conveying elements based on smoothed particle hydrodynamics. Part 2: Mixing, Chem. Eng. Sci. 134 (2015) 880–886. [13] S. Koshizuka, A. Nobe, Y. Oka, Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nucl. Sci. Eng. 123 (1996) 421–434. [14] X. Sun, M. Sakai, K. Shibata, Y. Tochigi, H. Fujiwara, Numerical modeling on the discharged fluid flow from a glass melter by a Lagrangian approach, Nucl. Eng. Des. 248 (2012) 14–21. [15] M. Sakai, Y. Shigeto, X. Sun, T. Aoki, T. Saito, J. Xiong, S. Koshizuka, Lagrangian– Lagrangian modeling for a solid–liquid flow in a cylindrical tank, Chem. Eng. J. 200–202 (2012) 663–672. [16] Y. Yamada, M. Sakai, Lagrangian–Lagrangian simulations of solid–liquid flows in a bead mill, Powder Technol. 239 (2013) 105–114. [17] X. Sun, M. Sakai, M.-T. Sakai, Y. Yamada, A Lagrangian–Lagrangian coupled method for three-dimensional solid–liquid flows involving free surfaces in a rotating cylindrical tank, Chem. Eng. J. 246 (2014) 122–141. [18] G. Li, Y. Oka, M. Furuya, Experimental and numerical study of stratification and solidification/melting behaviors, Nucl. Eng. Des. 272 (2014) 109–117. [19] R. Chen, Y. Oka, G. Li, T. Matsuura, Numerical investigation on melt freezing behavior in a tube by MPS method, Nucl. Eng. Des. 273 (2014) 440–448. [20] G. Li, M. Liu, G. Duan, D. Chong, J. Yan, Numerical investigation of erosion and heat transfer characteristics of molten jet impinging onto solid plate with MPS–LES method, Int. J. Heat Mass Transf. 99 (2016) 44–52. [21] Y. Liang, Z. Sun, G. Xi, L. Liu, Numerical models for heat conduction and natural convection with symmetry boundary condition based on particle method, Int. J. Heat Mass Transf. 88 (2015) 433–444. [22] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (2) (1985) 308–323. [23] P.W. Cleary, Extension of SPH to predict feeding, freezing and defect creation in low pressure die casting, Appl. Math. Model. 34 (11) (2010) 3189–3201. [24] M. Tong, D.J. Browne, An incompressible multi-phase smoothed particle hydrodynamics (SPH) method for modelling thermocapillary flow, Int. J. Heat Mass Transf. 73 (2014) 284–292. [25] M. Ramacciotti, C. Journeau, F. Sudreau, G. Cognet, Viscosity models for corium melts, Nucl. Eng. Des. 204 (1–3) (2001) 377–389. [26] Y. Shigeto, M. Sakai, Arbitrary-shaped wall boundary modeling based on signed distance functions for granular flow simulations, Chem. Eng. J. 231 (2013) 464–476. [27] M. Sakai, How should the discrete element method be applied in industrial systems?: a review, KONA Powder Part. J. 33 (33) (2016) 169–178. [28] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65. [29] M. Sakai, Y. Shigeto, G. Basinskas, A. Hosokawa, M. Fuji, Discrete element simulation for the evaluation of solid mixing in an industrial blender, Chem. Eng. J. 279 (2015) 821–839. [30] Y. Tsunazawa, Y. Shigeto, C. Tokoro, M. Sakai, Numerical simulation of industrial die filling using the discrete element method, Chem. Eng. Sci. 138 (2015) 791–809. [31] G. Basinskas, M. Sakai, Numerical study of the mixing efficiency of a ribbon mixer using the discrete element method, Powder Technol. 287 (2016) 380– 394. [32] M. Kondo, S. Koshizuka, Improvement of stability in moving particle semiimplicit method, Int. J. Numer. Methods Fluids 65 (6) (2011) 638–654.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Numerical study on a heat transfer model in a Lagrangian fluid dynamics simulation Kazuya Takabatake a, Xiaosong Sun a, Mikio Sakai b,⇑, Dimitrios Pavlidis c, Jiansheng Xiang c, Christopher C. Pain c a b c
Department of Nuclear Engineering and Management, School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Resilience Engineering Research Center, School of Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Applied Modeling and Computation Group, Department of Earth Science and Engineering, Imperial College London, Prince Consort Road, London SW7 2BP, UK
a r t i c l e
i n f o
Article history: Received 14 April 2016 Received in revised form 5 July 2016 Accepted 20 July 2016
Keywords: Moving particle semi-implicit method Lagrangian particle method Heat transfer Phase transition
a b s t r a c t Phenomena related to phase change heat transfer are often encountered in engineering. These phenomena are regarded to be complex, since not only phase transition from solid to liquid occurs but also movement of fluid interface has to be taken into consideration. Detailed numerical modeling of these complex systems is essential to better understand them and optimize industrial designs. Lagrangian methods are promising for simulating such complex systems. The Moving Particle Semi-implicit (MPS) method, which is one of the Lagrangian methods, is employed here to simulate the free surface fluid flows involving heat transfer and phase change. On the other hand, the existing MPS method could not apply Neumann boundary condition such as heat flux in the heat transfer simulations. This is because the surface direction could not be readily defined on the surface of the spherical fluid particles in the MPS method. Hence, prescribing the heat fluxes becomes problematic in the existing MPS method. To solve this problem, a new heat flux model is developed, where the divergence operator is applied in the heat transfer simulation. Simple verification tests are performed to demonstrate the heat flux model, where the calculation results are compared against analytically derived solutions. In addition, application of the signed distance function is also investigated in the heat transfer simulation for arbitrary shaped boundary. In simple verification tests, the computation results are shown to agree well with the analytical solutions. Consequently, adequacy of the novel heat transfer model developed here is shown in the Lagrangian fluid dynamics simulation. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction Phase change heat transfer is often encountered in various engineering [1–3] such as chemical engineering, food engineering, mechanical engineering, nuclear engineering, resilience engineering, etc. At present, introduction of numerical technologies are desired for accurate product design and/or better understanding of the phenomena. When phase change heat transfer phenomena are numerically simulated, not only heat transfer model but also interfacial fluid flow modeling is required. For modeling the interfacial fluid flow, the Volume-of-Fluid (VOF) [4–8], Smoothed Particle Hydrodynamics (SPH) [9–12] and Moving Particle Semiimplicit (MPS) method [13–15] have been employed in the computations so far. Lagrangian methods such as the SPH and MPS method are regarded to be accurate in simulating free surface fluid ⇑ Corresponding author. Fax: +81 3 5841 6977. E-mail address: [email protected] (M. Sakai). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.07.073 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
flows, since convection term is not required in the modeling. Accordingly, the Lagrangian approach becomes important in various engineering. Recently, weakly compressible MPS method [16,17] was developed for efficient calculation of the interfacial flows, where the pressure computation was performed explicitly. Hence, the calculation efficiency of the MPS method has been improved substantially and is now capable of industrial application. The MPS method was applied to complex thermal hydraulic systems such as a simulation of a severe accident in nuclear engineering [18,19]. Although the MPS method seems to be established, there are some issues in the heat transfer simulations. Specifically, in the existing MPS method, Neumann boundary condition such as heat flux boundary could not be set. This is because the surface direction could not be evaluated easily at the surface of MPS particles, and hence heat flux could not be set at the boundary. Only Dirichlet boundary could be used in the existing MPS method [18–21], and thus only the temperature could be set at the boundary. In addition, in the existing MPS method, accuracy
636
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Nomenclature C c cp D g h k L l N n n0 P Q q r r re T Tm t u w wg
Greek symbols variable [–] c phase fraction [–] C constant for surface detection [–] k Laplacian constant in the MPS method [–] l viscosity [Pa s] m fluid kinematic viscosity [m2/s] s viscous tensor q fluid density [kg/m3]
empirical constant for viscosity [–] sound speed [m/s] specific heat [J/(kg K)] dimension number [–] gravitational acceleration [m/s2] enthalpy [J/m3] thermal conductivity [W/(m K)] latent heat [J/kg] length [m] unit void vector [–] particle number density [–] initial particle number density [–] pressure [N/m2] heat source [J/m3] heat flux vector [W/m2] coordinate vector [m] distance [m] cutoff radius [m] temperature [K] melting point [K] time [s] fluid velocity [m/s] basic kernel function [–] special kernel function for pressure gradient [–]
U
Subscript i, j L S Wall Max 1,2
particle number liquid phase solid phase property of wall maximum property type number
Superscript ⁄ auxiliary value
of the heat transfer calculation is known to be worse in arbitrary shape wall systems. This is because the wall boundary was created by stepwise location of the MPS particles. Consequently, the existing MPS method has problems in heat transfer simulation. To solve the above issues, novel heat transfer model is developed in the MPS method. Specifically, heat flux model is newly developed in the current study. The weakly compressive MPS method is used in this study. To apply the MPS simulation of heat transfer in an arbitrary shaped domain, the signed distance function (SDF) is introduced for the first time. In addition, the implicit calculation algorithm [22] is newly introduced into the MPS method to simulate highly viscous fluid flow involving heat transfer efficiently. Verification tests are performed in a simple system to prove adequacy of the new models. Calculation results agree well with theoretical ones. Moreover, the numerical model is shown to simulate phase change system. Thus, the new heat transfer model is verified in Lagrangian fluid dynamics. 2. Methodology
Dh ¼ r ðkrTÞ þ Q ; Dt
ð3Þ
where h, k, T and Q are the enthalpy, thermal conductivity, temperature and energy inflow, respectively. 2.2. The MPS method In the MPS method [13], the differential operator is modeled by the weighted average between neighbor particles. A basic kernel function is given [16,17] as
( wðrÞ ¼
r re
þ rre 2 r < r e
0
r P re
;
ð4Þ
where r and re are distance and cutoff radius of interaction, respectively. Assuming that virtual particles are set in a structurally aligned particle gird, the initial particle number density and Laplacian constant are given as
n0 ¼ 2.1. The governing equations
X
wðjrj ri jÞ;
ð5Þ
j–i
For the weakly compressible fluid motion, the mass and momentum conservation equations can be written as
Dq þ qr u ¼ 0; Dt
ð1Þ
Du 1 ¼ rP þ r s þ g; Dt q
ð2Þ
where q, t, u, P, s and g are the fluid density, time, fluid velocity, pressure, viscous tensor and gravity, respectively. The s is given by
s ¼ r u þ r uT : Considering heat transfer, the energy conservation equation is given as follows:
k¼
1X jr j r i j2 wðjrj ri jÞ: n0 j–i
ð6Þ
These values are constant during the calculation. The particle density can be approximated to be proportional to the initial particle number density as follows:
ni ¼
X wðjrj ri jÞ;
ð7Þ
j–i
qi ¼
ni q: n0 0
ð8Þ
In the weakly compressible approach [16,17], pressure can be derived as a function of density and is given by
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
8 < q 0 c 2 n i 7 1 ni P n0 0 7 n : P¼ :0 ni < n0
2.3. Heat transfer model
ð9Þ
And by using the MPS gradient model, the pressure gradient term in Eq. (2) is derived as
rP i ¼
D X P i þ Pj rj ri wg ðr ij Þ : ng j–i jr j ri j jrj ri j
ð10Þ
Here, wg is a special kernel function for pressure gradient and is given by
wg ðrÞ ¼
r re
0
þ rre
r < re r P re
:
637
ð11Þ
Compared with conventional mesh methods, the heat transfer model is very limited in the MPS method. It is modeled by the same manner in the SPH method [24]. The energy conservation equation can be discretized in the MPS method as follows:
Dh 2D X 2ki kj ¼ ðT j T i Þwðjr j ri jÞ þ Q : Dt kn0 j–i ki þ kj
ð15Þ
In this study, two different properties are used to simulate the phase transition based on the enthalpy according to the past studies [18]. The density variation between solid and liquid phase is assumed to be small enough to be negligible compared with the variation between the liquid and gas phases. The temperature is derived as a function of enthalpy and is given by:
Besides, ng is calculated by using Eq. (11) in the same manner as n0 and makes it possible to discretize the pressure gradient term while the momentum remains constant. Eq. (2) can be decomposed into the following equations:
8 > T þ hhs > < m qcp T ¼ Tm > > : T m þ hhl
u u k 1 ¼ g rP; Dt q
ð12Þ
ukþ1 u ¼ r skþ1 : Dt
ð13Þ
where Tm, hs, hl and cp are the melting point, enthalpy of solid, that of liquid and specific heat, respectively. As shown in Fig. 1, temperature during phase transition is assumed to be constant. Enthalpy of the solid phase hs and that of the liquid phase hl are expressed as follows:
This is the similar to the fractional time step method [22] except for the pressure term because the value of pressure gradient is already known in the weakly compressible MPS method. Here, Eq. (13) is modeled by the same manner in previous studies
rs¼
2d X 2mi mj ðu u Þwðr r Þ : j i j i kn0 j–i mi þ mj
ð14Þ
qcp
h < hs hs 6 h 6 hl ;
ð16Þ
h > hl
hs ¼ qcp T m ;
ð17Þ
hl ¼ hs þ qL;
ð18Þ
where L is latent heat. In order to simulate phase transition, the phase fraction coefficient g is introduced and is calculated as:
c¼
8 > <1 > :
hl h hl hs
0
h < hs hs 6 h 6 hl :
ð19Þ
h > hl
This modeling has been successfully applied in various SPH and MPS simulations [18,19,23,24]. The auxiliary velocity u⁄ in Eq. (12) is estimated explicitly from Eqs. (9) and (10). The updated velocity uk+t in Eq. (13) is then calculated implicitly by using conjugate gradient method.
As shown in Fig. 2, the phase is determined from enthalpy and is interpolated by the linear function during phase transition.
Fig. 1. Temperature model.
Fig. 2. Phase change model.
638
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
In the current study, the material was modeled by a fluid. When the material is in solid state, the fluid viscosity was set to be sufficiently high. The viscosity is a function of enthalpy. Hence, the modeling of viscosity is indispensable for the heat transfer simulation. In this study, the viscosity model is expressed by the Arrhenius viscosity formula as:
lr ¼ ll expð2:5C cÞ;
ð20Þ
where lr, ll and C are the relative viscosity, viscosity of liquid phase and empirical constant, respectively. This empirical constant varies between 4 and 8 [25]. Eq. (20) makes it possible to simulate the behavior of a rigid body as an extremely high viscous liquid. 2.4. Boundary conditions 2.4.1. Dirichlet boundary conditions Heat transfer is calculated by the Laplacian of the temperature between neighboring particles as shown in Eq. (3). When a Dirichlet boundary condition is applied, the temperature at the boundary is fixed. Thus, the source term (the second term in the right hand side of Eq. (3)) is given by:
Q ¼ 0:
ð21Þ
For Dirichlet boundary conditions, no special treatment of surface detection is required. 2.4.2. Neumann boundary conditions In the case of Neumann boundary conditions, the heat flux at the boundary is fixed. In the meshless particle method, it is difficult to evaluate the surface area and its normal direction. The heat flux can be set as a volumetric heat source by introducing the divergence theorem. As such the heat fluxes at boundaries are modeled by:
Q ¼ r q;
ð22Þ
where q is heat flux vector. In the Neumann boundary condition, surface must be detected depending on the direction of heat flux vector. First, an ordinary surface detection is applied for each particle as shown in Eq. (23):
ni < Cn0 ;
ð23Þ
where C is a constant for surface detection and is set to 0.97 in this study as in the previous study [13,14]. Then, the void vector N, which is perpendicular to the surface, is calculated for each surface particle detected by Eq. (24) as follows:
P j–i
ðr j ri Þwðr j ri Þ jr r j
j i : N ¼ P ðr r Þwðr r Þ j–i j jri j ri jj i
ð24Þ
Finally, assuming that a conceptual particle is set at the head of the void vector, the surface where heat flux enters can be detected by:
N q < 0:
ð25Þ
tance from walls). One can see the details of the SDF in our previous studies. The SDF could be applied to a simulation of a granular flow in a complex shape boundary system such as a twin-screw kneader [29], a die-filling [30] and a ribbon mixer [31]. Adequacy of the SDF based wall boundary was demonstrated through validation tests. In the MPS method, the distance from walls is calculated at all the sampling points in the calculation domain at the beginning of each step. Assuming that the virtual particles are aligned as the structural grid in the wall, the value of the kernel function from the wall at each sampling point can be calculated from the distance. The value of the kernel function from wall at each particle is then calculated by interpolation using information of neighboring sampling points. Application of the SDF to simulations of solid– liquid flows has already been presented in our past studies [16,17], where the MPS method was employed to simulate the liquid phase. In this study, the SDF was introduced into the heat transfer simulation of the MPS method for the first time. 3. Numerical simulations 3.1. Conditions Three kinds of verification tests were performed to show the adequacy of the heat flux model. Besides, the applicability of the SDF to MPS method was proved in the heat transfer simulation. In the verification tests, a simple shape system such as a rectangular bar was chosen because of comparison with the analytical solution. Usage of the simple shape system is of great importance because the new technology is shown to be applied to a great deal of systems in engineering. Our objective is to show the agreement between the calculation and analytical results. Accordingly, applicability of the SDF was investigated in a simple system. In these verification tests, the temperature distribution was compared with the theoretical solution. In addition, the phase transition simulation was also performed due to heat flux by using a Neumann boundary condition. In the verification tests, time step, speed of sound and the empirical constant related to viscosity were, respectively, 1.0 105 s, 100 m/s and 8.0. The fluid particle size, which is so-called initial particle distance, was 5.0 mm. Value of the parameters used in the verification tests are summarized in Table 1. 3.1.1. Case 1: Dirichlet boundary condition Case 1 was performed to show application of the SDF in the heat transfer simulation. Heat transfer in the rectangular bar was evaluated in the verification test. Temperature distributions were compared with the calculation results between a particle based wall (Case 1-1) and the SDF based wall (Case 1-2). Fig. 3 shows the calculated particle arrangement. Dirichlet boundary condition was set at both the top and the bottom wall and adiabatic boundary at the side wall. In Case 1-1, the wall boundary was given by wall particles shown in Fig. 3(a). In Case 1-2, the wall was created by the SDF
By using the algorithm outlined from Eqs. (23)–(25), entrance boundaries can be detected. 2.4.3. Signed distance function Originally the SDF [26,27] was developed for granular flow simulations using the discrete element method [28]. The SDF is defined by sign and minimal distance from the calculation domain. The SDF value becomes positive inside the calculation domain, and vice versa. In the DEM, SDF makes it possible to relate displacement for the contact force and the solid particle location (or dis-
Table 1 Simulation conditions. Item
Unit
Value
Time step Gravitational acceleration Sound speed Particle diameter Empirical constant for viscosity
s m/s2 m/s mm –
1.0 105 9.8 100 5.0 8.0
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
(a) Case 1-1: particle wall
639
(b) Case 1-2: SDF based wall
Fig. 3. Initial arrangement of Case 1.
Table 2 Physical properties (Case 1). Item
Unit
Material 1
Density Kinematic viscosity Specific heat Thermal conductivity
kg/m3 m2/s J/(kg K) W/(m K)
1000 1.0 103 10 50
as designated in Fig. 3(b). The temperature was given at 300 K and 400 K at the top and bottom walls. The initial temperature was set to be 300 K in the bar. The bar was composed of one material, where density, kinematic viscosity, specific heat and thermal conductivity are prescribed as 1000 kg/m3, 1.0 103 m2/s, 10 J/(kg K) and 50 W/(m K), respectively. Note that the value of kinematic viscosity of the solid state became 5.0 108 times larger than that of the liquid state by using Eq. (20). Table 2 shows the physical properties. Objective of this verification is to show the adequacy of the heat transfer modeling, and hence these properties were arbitrarily chosen. Especially, heat capacity was extremely low in order to reduce the computational time. 3.1.2. Case 2: Neumann boundary condition without phase transition Case 2 was performed to show adequacy of the heat flux model in the MPS method. Fig. 4 shows the calculation systems. The geometrical configuration of the bar was the same as that in Case 1. To investigate the effects of material properties, a single bar composed of material-1 and a combined bar composed of material-1 and material-2 were used in Cases 2-1 and 2-2. In both cases, the bottom wall was defined by the SDF. Neumann boundary condition, namely, heat flux, was
applied at the top surface. Dirichlet boundary condition, namely, 300 K of constant temperature, was given at the bottom SDF wall. Adiabatic boundary condition was given at the side SDF wall. The intensity of heat flux was 25,000 W/m2. The initial temperature was 300 K in the bar. In the materials, the density, kinematic viscosity and specific heat were 1000 kg/m3, 1.0 103 m2/s and 10 J/(kg K). On the other hand, the thermal conductivities were different in Cases 2-1 and 2-2. In Case 2-2, the different thermal conductivity was 50 W/(m K) and 500 W/ (m K) bottom half and top half, respectively. Table 3 summarizes the physical properties of Case 2.
3.1.3. Case 3: Neumann boundary condition with phase transition Case 3 was performed to simulate phase change material under the heat flux model by the MPS method. Fig. 5 illustrates the calculation systems. The bar was composed of two different materials as in Case 2. In both cases, bottom wall was expressed by the SDF. Neumann boundary condition, namely, heat flux, was applied at the top surface. The intensity of heat flux was 100,000 W/m2 in order to finish the calculation in short time. Adiabatic condition and constant temperature were given at both the side and the bottom in Cases 3-1 and 3-2. The initial temperature was 300 K in the bar. As far as the density, kinematic viscosity and latent heat were concerned, common values were used in Cases 3-1 and 3-2, namely, they were 1000 kg/m3, 1.0 103 m2/s and 1.0 104 J/ kg, respectively. On the other hand, specific heat, thermal conductivity and melting point were given by different values. The specific heat, thermal conductivity and melting point were 100 J/(kg K), 50 W/(m K) and 10,000 K in material 1. In material 2, the specific heat, different thermal conductivity and melting point were 10 J/ (kg K), 500 W/(m K) and 320 K, respectively. Table 4 shows the physical properties in Case 3.
640
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
(a) Case 2-1: one material
(b) Case 2-2: two materials
Fig. 4. Initial arrangement of Case 2. Table 3 Physical properties (Case 2).
Table 4 Physical properties (Case 3).
Item
Unit
Material 1
Material 2
Density Kinematic viscosity Specific heat Thermal conductivity
kg/m3 m2/s J/(kg K) W/(m K)
1000 1.0 103 10 50
1000 1.0 103 10 500
(a) Case 3-1: adiabatic
Density Kinematic viscosity Specific heat Thermal conductivity Melting point Latent heat
kg/m3 m2/s J/(kg K) W/(m K) K J/kg
(b) Case 3-2: Dirichlet boundary
Fig. 5. Initial arrangement of Case 3.
Material 1
Material 2
1000 1.0 103 100 50 10,000 1.0 104
1000 1.0 103 10 500 320 1.0 104
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
0s
1s
2s
641
3s
Fig. 6. Snapshot of heat transfer simulation in Case 1-1.
ature was gradually increased upside due to Dirichlet condition. The temperature distribution reached steady state after around 3 s. Fig. 9 shows the transient change of temperature distribution in Case 1-2. The maximum temperature of the bar became 399.9 K at the particle on the top boundary. The simulation results were compared with the theoretical solution. The theoretical solution can be easily obtained in this problem. This is because the temperature can be linearly interpolated between the walls. The temperature can be conducted theoretically to be 398.8 K in the calculation system, when the sampling point was same as that of the simulations. In spite of such a simple case, temperature distribution obtained in Case 1-1 was slightly different from that of Case 1-2. This is due to the deformation by the gravitational force when the particle wall boundary was applied. The deformation may increase or decrease the number of wall particles inside cutoff radius, because the coordinates of wall particles are fixed during the calculation. On the other hand, the value of kernel function was determined by the distance from wall when the SDF was used. Therefore, the deformation was not intrinsically an important problem. Nevertheless, the difference of calculation results was less than 1% from theoretical one. The error might be occurred due to the spatial resolution. Hence, the temperature between the simulations and theoretical solution was in good agreement. Besides, application of the SDF was also proved through the verification tests.
4. Results and discussion 4.1. Case 1: Dirichlet boundary condition In order to show application of the SDF in the heat transfer simulation, Case 1-1 and Case 1-2 were performed under Dirichlet condition. As mentioned above, in Case 1-1, the wall boundary was created by the MPS particles. Fig. 6 illustrates the snapshot of the heat transfer simulation in Case 1-1. The temperature was gradually increased upside due to Dirichlet condition. The temperature distribution became steady state in 3 s. Fig. 7 shows the transient change of temperature distribution in Case 1-1. The maximum temperature was 396.3 K at the particle on the top boundary of the bar. In Case 1-2, the wall was generated by the SDF. Fig. 8 designates the typical snapshot of the heat transfer simulation in Case 1-2. By the same manner as Case 1-1, the temper-
4.1.1. Case 2: Neumann boundary condition without phase transition In order to demonstrate the heat flux model by the MPS method, Cases 2-1 and 2-2 were performed under the Neumann boundary condition. In Case 2-1, the bar composed of a single material. Fig. 10 illustrates a snapshot of the heat transfer simulation in Case 2-1. The temperature was gradually increased from the top due to the heat flux. The temperature distribution became
Fig. 7. Transient change of temperature distribution in Case 1-1.
0s
1s
2s
3s
Fig. 8. Snapshot of heat transfer simulation in Case 1-2.
642
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
ture distribution in Case 2-2. The maximum temperature of the bar became 350.1 K. The simulation results were compared with the theoretical solution. The maximum temperature of each case can be calculated theoretically by using the following equations:
T wall T max ; l
ð26Þ
T wall T max ; l1 =k1 þ l2 =k2
ð27Þ
q ¼ k
q¼
Fig. 9. Transient change of temperature distribution in Case 1-2.
steady state in 20 s. Fig. 11 shows the transient change of temperature distribution in Case 2-1. The maximum temperature of the bar was 392.6 K at the particle on the top boundary of the bar. The linear temperature gradient was obtained. In Case 2-2, the bar composed of material-1 and material-2. Fig. 12 illustrates the snapshot of the heat transfer simulation in Case 2-2. As in Case 2-1, the temperature was gradually increased from the top due to the heat flux. The temperature distribution became steady state after around 20 s. Fig. 13 shows the transient change of tempera-
where l is length. Substituting the physical properties in Eqs. (26) and (27), the theoretical solution of the maximum temperature of the bar in Case 2-1 and Case 2-2 were 398.8 K and 354.9 K, respectively. The error became around 1.5 % in Cases 2-1 and 2-2. The error might be occurred due to the spatial resolution. In addition, the error of Case 2-2 was slightly larger than that of Case 2-1. This is because the heat input varied from place to place such as at the surface, at the side and the vertex. Through the verification tests, the application of the heat flux model was proved. 4.2. Case 3: Neumann boundary condition with phase transition Case 3 was performed to simulate phase change under the heat flux model by the MPS method. Case 3-1 was performed to simulate melting in an adiabatic environment. Case 3-2 was performed
0 sec
1 sec
2 sec
5 sec
10 sec
20 sec
Fig. 10. Snapshots of Case 2-1.
3 sec
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Fig. 11. Transient change of temperature distribution in Case 2-1.
643
Fig. 13. Transient change of temperature distribution in Case 2-2.
0 sec
1 sec
2 sec
5 sec
10 sec
20 sec
3 sec
Fig. 12. Snapshots of Case 2-2.
to simulate the resolidification of the molten material under Dirichlet condition. Fig. 14 illustrates a snapshot of the heat transfer simulation in Case 3-1. The temperature was gradually increased from the top due to heat flux. The bar began to melt in 2 s because the local temperature exceeded the melting temperature. Upper part was almost molten in 10 s. All the material on the top half was molten completely in 20 s. The molten fluid was distributed on the floor. Hence, melting phenomena due to heat flux could be simulated qualitatively by the proposed method. Fig. 15 illustrates the snapshot of the heat transfer simulation in Case 3-2. The temperature was grad-
ually increased from the top due to the heat flux. As with Case 3-1, the bar began to melt in 2 s because the local temperature exceeded the melting temperature. The molten fluid re-solidified at the floor due to rapid cooling. All the material on the top half of the bar melted completely in 20 s. The molten fluid was distributed near the bar unlike Case 3-1. Hence, resolidification could be simulated qualitatively by the proposed method. These results are valuable because the melting and solidification under heat flux could be simulated by the MPS method for the first time. At the early stage of the melting phenomenon in Figs. 14 and 15, the represented processes seem to be symmetric. On the other
644
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Fig. 14. Snapshots of Case 3-1.
0 sec
1 sec
2 sec
4 sec
5 sec
10 sec
3 sec
Fig. 15. Snapshots of Case 3-2.
hand, after 3 s, they seem to be asymmetric in both cases. It is well known that the pressure oscillation occurred in the MPS method [32]. This pressure oscillation might make the fluid flow asymmetric.
the finite difference method not only at the equilibrium state but also at the non-equilibrium state. The adequacy of our divergence model was shown through the investigation. 5. Conclusions
4.3. Convergence study regarding spatial resolution As already mentioned, the MPS method was in good agreement with the theoretical solution by using our divergence model. In this section, resolution effect on the divergence model was examined by comparing with the finite difference method. Fig. 16 shows the verification results where spatial resolution was investigated in the MPS method. As shown in this figure, the calculation error was improved due to the spatial resolution in the MPS method. Besides, it was shown that this model was in good agreement with
Lagrangian approaches such as the SPH and MPS methods are often employed in simulating free surface fluid flows. In recent years, the MPS method has been used for complex fluid flows such as free surface fluid flows involving heat transfer. In the existing Lagrangian approach, heat flux could not be set because of difficulty of surface direction definition. In addition, creation of arbitrary shape wall boundary was also problematic. This is because the MPS particles should be located in heat transfer simulation in the existing MPS method. In order to solve these problems, a
K. Takabatake et al. / International Journal of Heat and Mass Transfer 103 (2016) 635–645
Fig. 16. Comparison between the difference method and the MPS method in Case 2.
new heat transfer model was developed in the MPS method. Moreover, in this study, application of the SDF to simulate heat transfer was investigated for the first time. In order to introduce heat fluxes into the MPS method Neumann boundary condition was modeled by using a divergence model. Adequacy of this modeling was demonstrated using a simple system by comparing the calculation results against the analytical solutions. Employing the simple shape system is important because the new model is shown to be used to a great deal of systems in engineering. Application of the SDF to the MPS simulation of heat transfer was verified for the first time in this study. The verification test was performed in a simple system. The temperature obtained from the SDF was in good agreement with that of the particle based wall boundary. Finally, it was shown that melting material due to heat fluxes and resolidification could be simulated successfully by this new technologies. Consequently, innovative heat transfer model was developed in a Lagrangian approach. Acknowledgement This study was financially supported by the Initiatives for Atomic Energy Basic and Generic Strategic Research by the Ministry of Education, Culture, Sports, Science, and Technology of Japan. References [1] G. Quarini, Y. Chang, Heat transfer characteristics of ice melting in water and salt solutions, Chem. Eng. Res. Des. 80 (2002) 320–324. [2] K.W. Chua, Y.T. Makkawi, M.J. Hounslow, Time scale analysis for fluidized bed melt granulation III: binder solidification rate, Chem. Eng. Sci. 66 (3) (2011) 336–341. [3] Y. Allouche, S. Varga, C. Bouden, A.C. Oliveira, Validation of a CFD model for the simulation of heat transfer in a tubes-in-tank PCM storage unit, Renew. Energy 89 (2016) 371–379.
645
[4] K. Yokoi, Efficient implementation of THINC scheme: a simple and practical smoothed VOF algorithm, J. Comput. Phys. 226 (2) (2007) 1985–2002. [5] X. Sun, M. Sakai, Numerical simulation of two-phase fl ows in complex geometries by using the volume-of-fluid/immersed-boundary method, Chem. Eng. Sci. 139 (2016) 221–240. [6] X. Sun, M. Sakai, Three-dimensional simulation of gas–solid–liquid flows using the DEM–VOF method, Chem. Eng. Sci. 134 (2015) 531–548. [7] D. Pavlidis, Z. Xie, J.R. Percival, J.L.M.A. Gomes, C.C. Pain, O.K. Matar, Two- and three-phase horizontal slug flow simulations using an interface-capturing compositional approach, Int. J. Multiphase Flow 67 (S) (2014) 85–91. [8] Z. Xie, D. Pavlidis, J.R. Percival, J.L.M.A. Gomes, C.C. Pain, O.K. Matar, Adaptive unstructured mesh modelling of multiphase flows, Int. J. Multiphase Flow 67 (S) (2014) 104–110. [9] J.J. Monaghan, An introduction to SPH, Comput. Phys. Commun. 48 (1988) 89– 96. [10] X. Sun, M. Sakai, Y. Yamada, Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method, J. Comput. Phys. 248 (2013) 147–176. [11] A. Eitzlmayr, J. Khinast, Co-rotating twin-screw extruders: detailed analysis of conveying elements based on smoothed particle hydrodynamic. Part 1: Hydrodynamics, Chem. Eng. Sci. 134 (2015) 861–879. [12] A. Eitzlmayr, J. Khinast, Co-rotating twin-screw extruders: detailed analysis of conveying elements based on smoothed particle hydrodynamics. Part 2: Mixing, Chem. Eng. Sci. 134 (2015) 880–886. [13] S. Koshizuka, A. Nobe, Y. Oka, Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nucl. Sci. Eng. 123 (1996) 421–434. [14] X. Sun, M. Sakai, K. Shibata, Y. Tochigi, H. Fujiwara, Numerical modeling on the discharged fluid flow from a glass melter by a Lagrangian approach, Nucl. Eng. Des. 248 (2012) 14–21. [15] M. Sakai, Y. Shigeto, X. Sun, T. Aoki, T. Saito, J. Xiong, S. Koshizuka, Lagrangian– Lagrangian modeling for a solid–liquid flow in a cylindrical tank, Chem. Eng. J. 200–202 (2012) 663–672. [16] Y. Yamada, M. Sakai, Lagrangian–Lagrangian simulations of solid–liquid flows in a bead mill, Powder Technol. 239 (2013) 105–114. [17] X. Sun, M. Sakai, M.-T. Sakai, Y. Yamada, A Lagrangian–Lagrangian coupled method for three-dimensional solid–liquid flows involving free surfaces in a rotating cylindrical tank, Chem. Eng. J. 246 (2014) 122–141. [18] G. Li, Y. Oka, M. Furuya, Experimental and numerical study of stratification and solidification/melting behaviors, Nucl. Eng. Des. 272 (2014) 109–117. [19] R. Chen, Y. Oka, G. Li, T. Matsuura, Numerical investigation on melt freezing behavior in a tube by MPS method, Nucl. Eng. Des. 273 (2014) 440–448. [20] G. Li, M. Liu, G. Duan, D. Chong, J. Yan, Numerical investigation of erosion and heat transfer characteristics of molten jet impinging onto solid plate with MPS–LES method, Int. J. Heat Mass Transf. 99 (2016) 44–52. [21] Y. Liang, Z. Sun, G. Xi, L. Liu, Numerical models for heat conduction and natural convection with symmetry boundary condition based on particle method, Int. J. Heat Mass Transf. 88 (2015) 433–444. [22] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (2) (1985) 308–323. [23] P.W. Cleary, Extension of SPH to predict feeding, freezing and defect creation in low pressure die casting, Appl. Math. Model. 34 (11) (2010) 3189–3201. [24] M. Tong, D.J. Browne, An incompressible multi-phase smoothed particle hydrodynamics (SPH) method for modelling thermocapillary flow, Int. J. Heat Mass Transf. 73 (2014) 284–292. [25] M. Ramacciotti, C. Journeau, F. Sudreau, G. Cognet, Viscosity models for corium melts, Nucl. Eng. Des. 204 (1–3) (2001) 377–389. [26] Y. Shigeto, M. Sakai, Arbitrary-shaped wall boundary modeling based on signed distance functions for granular flow simulations, Chem. Eng. J. 231 (2013) 464–476. [27] M. Sakai, How should the discrete element method be applied in industrial systems?: a review, KONA Powder Part. J. 33 (33) (2016) 169–178. [28] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65. [29] M. Sakai, Y. Shigeto, G. Basinskas, A. Hosokawa, M. Fuji, Discrete element simulation for the evaluation of solid mixing in an industrial blender, Chem. Eng. J. 279 (2015) 821–839. [30] Y. Tsunazawa, Y. Shigeto, C. Tokoro, M. Sakai, Numerical simulation of industrial die filling using the discrete element method, Chem. Eng. Sci. 138 (2015) 791–809. [31] G. Basinskas, M. Sakai, Numerical study of the mixing efficiency of a ribbon mixer using the discrete element method, Powder Technol. 287 (2016) 380– 394. [32] M. Kondo, S. Koshizuka, Improvement of stability in moving particle semiimplicit method, Int. J. Numer. Methods Fluids 65 (6) (2011) 638–654.