Modeling mesoscopic solidification using dissipative particle dynamics

International Journal of Thermal Sciences 101 (2016) 207e216

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Modeling mesoscopic solidification using dissipative particle dynamics
n a, Jinliang Yuan a Erik O. Johansson a, Toru Yamada a, b, *, Bengt Sunde a b

Department of Energy Sciences, Faculty of Engineering, Lund University, Box 118, SE-221 00 Lund, Sweden Graduate School of Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan

a r t i c l e i n f o

a b s t r a c t

Article history: Received 17 July 2015 Received in revised form 5 November 2015 Accepted 6 November 2015 Available online 9 December 2015

Dissipative particle dynamics with energy conversion (DPDe) is a simulation technique that has been used to model thermal transport characteristics and heat transfer at mesoscale. This study shows further development of the DPDe method capturing solid/liquid phase-change phenomena and its application to water freezing in a parallel-plate straight channel. In this work, the weighting functions of the random and dissipative forces are modeled as functions of temperature in order to correctly predict the temperature dependent properties of the fluid in a two-dimensional domain. An equation of state is incorporated in the model in order to model the solidification of water. Careful consideration is taken to couple the latent heat of the system to real world units, and the solidification predicted using this model is compared to a well known analytical solution. The developed model is employed to simulate the thermally developing flow in a parallel-plate channel with constant wall temperatures below the freezing point. A liquid pump is introduced along with a region initiating the liquid temperature in order to create the thermally developing flow. The investigations show that the fluid velocity has a small effect on the time it takes for the channel to freeze completely. The dominating factor will be the temperature of the solid walls in the domain. The simulations also show that when a higher wall temperature is applied, the solid/liquid interface will be rougher due to mesoscopic fluctuations of heat and momentum. © 2015 Elsevier Masson SAS. All rights reserved.

Keywords: Dissipative particle dynamics Phase-change Variable properties Liquid pump

1. Introduction Dissipative Particle Dynamics (DPD) is a particle-based, meshfree simulation method introduced by Hoogerbrugge and Koelman [1]. The method has been used extensively in simulations at the mesoscale, where the domains investigated are too large for molecular dynamics (MD) simulations to be efficient, but the materials are still not sufficiently homogeneous to justify a macro-scale continuum description. Polymers used in the membranes of fuel cells [2e4], red blood cells [5,6], and surfactants [7,8] are representative examples of DPD applications. In every time step of a DPD simulation, three types of interactions occur: (i) conservative repulsion caused by spatial arrangement and energetic interaction between the different particles, (ii) dissipative interaction caused by energy lost due to friction or viscosity within a particle, and (iii) random interactions stemming from the thermal motion of the

* Corresponding author. Graduate School of Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan. E-mail address: [email protected] (T. Yamada). http://dx.doi.org/10.1016/j.ijthermalsci.2015.11.002 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.

molecules within a DPD particle [9]. The dissipative and random interactions are linked together in order to satisfy the fluctuationedissipation theorem and act as a thermostat to maintain a constant simulation temperature. The result is that although momentum is conserved, energy is not [10]. To be able to model thermal energy transport in mesoscopic systems, the dissipative particle dynamics with energy conserva~ ol [10] and Avalos and Mackie tion (DPDe) was introduced by Espan [11]. The DPDe method has been used to study heat transfer by convection in several configurations [12e14], thermal conduction in nanoparticle suspensions [15e18], and ballistic-diffusive thermal conduction in thin films [19]. A comprehensive review of the current status of the DPDe method was recently published by Chaudhri and Lukes [20]. A challenge that should be taken into account for further enhancement of the DPDe method is to capture the phase-change phenomenon, a topic of significance at mesoscale. An example showing the importance of phase-change is the ice formation in the porous materials used in proton exchange membrane fuel cells (PEMFCs). Recent studies have shown crucial performance

208

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

Nomenclature A a Cv D DPD DPDe e f hsf K kb L m MD ndim P PEMFC Pr Q r rc R Sc T Tfr t tcl v

Strength of conductive heat transfer Strength of the conservative force Specific heat at constant volume, J kg1 K1 Diffusivity, m2 s1 Dissipative particle dynamics Dissipative particle dynamics with energy conservation Unit vector Force, N Latent heat, J kg1 Proportionality constant Boltzmann constant Spatial length, m Mass, kg Molecular dynamics Number of dimensions Pressure, Pa Proton exchange membrane fuel cell Prandtl number (na1) Heat flux, W m2 Inter-particle distance Cutoff-radius Radius, m Schmidt number (nD1) Temperature, K Freezing temperature, K Time, s Closing time, s Velocity vector

degradation and mechanical damage of the materials caused by the ice formation during PEMFC operations [21]. The phase-change also has a significant impact for the manufacturing society. Threedimensional additive manufacturing devices have been rapidly developed. The resolution with good surface finish of these devices is highly related to the liquid transport phenomena involving liquidesolid phase transition [22]. Therefore, further development of the phase-change model in the DPDe method will enhance its applicability to address the current scientific and engineering challenges. Another challenge in using the DPDe method is to correctly predict the behavior of simple fluids such as water or glycerin. In the standard DPDe method, the viscosity of the simulated fluid is directly proportional to the temperature of the fluid. However, for simple fluids, the opposite is true. Recently, Li et al. [23] proposed an improvement to the DPDe-model where the temperature dependent properties of simple fluids can be correctly predicted by incorporating a temperature dependence in the weighting function of the dissipative and random forces. Their improvement is also applied in this paper, together with an equation of state proposed by Willemsen et al. [24] that allows the modeling of phase-change between solid and liquid in the simulation domain. The objective of this paper is to combine the proposed temperature dependent DPDe model for simple fluids [23] with the equation of state for solid/liquid interaction [24] into a model that can not only correctly predict the thermophysical properties of a simple fluid, but also the phase-change between liquid and solid state. Water is chosen as the working fluid, and careful consideration is taken to validate that the model behaves similar to water in

Y

a b g ε

z, ze k k0 l n r s u

Position of interface, m Thermal diffusivity, m2 s1 Ratio between fluid and solid density Dissipative strength Internal energy, J kg1 Random numbers with zero mean and unit variance Strength of conductive heat flux Heat friction coefficient Thermal conductivity, W m1 K1 Kinematic viscosity, m2 s1 DPDe Density Strength of random interactions Weighting function

Subscripts C Cold ext External i, j Indices s, f Solid, fluid x, y Direction Superscripts C Conservative D Dissipative HC Conductive heat HR Random heat HV Viscous heat ne Exponential of dissipative weighting function R Random < Real

the investigated range of temperatures with regard to diffusivity, viscosity, heat transfer, and rate of phase-change to the solid state. The length, heat capacity, and latent heat of the simulation are linked to the real world units via physically justified equations. The paper also covers an application of the proposed model: the simultaneous fluid flow and freezing of water in a parallel plate channel with constant sub-zero wall temperatures. This set of simulations can be useful in studies of, e.g., the behavior of microfluidic phase-change valves [25]. In order to study forced convection using the DPDe method, a liquid pump [26,27] is introduced into the simulation domain.

2. Methodology 2.1. Governing equations The DPD method is a mesh-free, particle-based method where the different particles in the simulation domain interact with one another through a set of distance and velocity dependent forces [9]. The time evolution of the simulated system is governed by Newton's second law:

mi

 dvi X C R f ij þ f D ¼ ij þ f ij þ f ext dt jsi

(1)

fC, fD, and fR are the forces acting between the particles in the simulation, namely conservative, dissipative, and random forces, respectively. fext is an external force that can be used to induce a

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

fluid flow in the simulation domain. The conservative, dissipative, and random forces are given by:

  f Cij ¼ aij uC rij eij    D fD ij ¼ gij u rij eij $vij eij   f Rij ¼ sij uR rij zij Dt 1=2 eij

(2)

(3)

 dεi X HV ¼ qij þ qHC þ qHR ij ij dt jsi

(4)



:

rij < rc



  rij  rc

:

(5)

The fluctuationedissipation theorem states that in order to achieve thermodynamic equilibrium from this method, the following relationships must hold [10]:

gij ¼

  s2ij Ti þ Tj 4kB Ti Tj

;

    2 uD rij ¼ uR rij

and change due to random heat flux caused by thermal fluctuation (qHR ij ). These three terms are expressed as

qHV ij ¼

  rij < rc

:

  rij  rc

(7)

ne is the exponent that controls the viscosity of the DPDe fluid. A challenge in using the DPDe method is to present a model that correctly describes the temperature dependent physical properties of a liquid. Many common liquids such as water and glycerin have viscosities decreasing, but diffusivities increasing with increasing temperature. In order to capture this phenomenon, Li et al. [23] introduced a model where the exponent ne is a function of the particle temperature. Their approach is incorporated in this work, and ne is given as

    ne Tij ¼ 0:41 þ 1:9 Tij  1

0 < ne < 1:5

(8)

where Tij is the average temperature between particles i and j. The temperature of the particles in the simulation is allowed to vary between 0.91 (corresponding to 273 K), and 1.2433 (corresponding to 373 K).

" ) # (    2 s2ij    1 uD rij gij eij $vij   sij uR rij eij $vij zij 2CV mi (10)

  1 1 HC qHC  rij ij ¼ kij u Ti Tj

! (11)

  e 1=2 HR rij zij Dt qHR ij ¼ Aij u

(12)

kij and Aij are the strengths of the conductive and the random heat fluxes, respectively. They are given by: kij ¼ Cv2 k Aij ¼

 2 Ti  Tj 4kb

(13)

qffiffiffiffiffiffiffiffiffiffiffiffiffi 2kB kij

(14)

where k in eq. (13) is a mesoscale heat friction coefficient. An expression for k as a function of Prandtl number and viscosity was derived by Li et al. [23]. Following their approach in a twodimensional domain one obtains:

k¼

:

(9)

The change in internal energy consists of three parts: change HC due to viscous heating (qHV ij ), change due to heat conduction (qij )

(6)

Ti is the temperature of DPDe particle i. The dissipative weight function is given by:

8  rij ne > < 1    rc uD rij ¼ > : 0

2.2. Conservation of energy The conservation of energy requires the tracking of an additional variable, the internal energy of the DPDe particle ε, and an additional governing equation [10]:

where rij is the distance between particles i and j, and eij is a unit vector between them. a in eq. (2) is the maximum repulsion between two particles, commonly chosen to be aij ¼ 75.0kBT/r where kB is the Boltzmann constant and the density used is r ¼ 4. In eq. (4), zij is a random number with zero mean and unit variance. Each interacting pair of particles have the property zij ¼ zji to ensure that the momentum of interacting particles is conserved. The factor Dt1/2 in eq. (4) arises from the integration of the random equations of motion [28]. uC, uD, and uR are weight functions for the different forces. These functions decide the strength of the inter-particle forces between the particles i and j as a function of the distance between them. The two particles i and j have no interaction beyond a certain cutoff radius, rc. The conservative weight function is given as:

 8 rij > < 1    rc uC rij ¼ > : 0

209

120 n prc4 Cv r Pr

(15)

The kinematic viscosity, n, and Prandtl number, Pr, are temperature dependent properties, and the other terms in eq. (15) are given as input in the simulation. The weight functions uHC and uHR in eqs. (11) and (12) are related as uHC ¼ (uHR)2 [10]. zeij is an antisymmetric random number with unit variance, zeij ¼ zeji . To be able to model phase change in the simulation domain, an equation of state presented by Willemsen et al. [24] is introduced. The temperature of each DPDe particle is determined by the internal energy, the heat capacity at constant volume, Cv, and the latent heat, hsf, as follows:

. 8 > Cv ; εi  hsf > > < Ti ¼ Tfr ; > > > : ε =C ; v i

 

εi > Cv Tfr þ hsf



Cv Tfr < εi < Cv Tfr þ hsf   εi < Cv Tfr



(16)

The equation of state is divided into three regions: a liquid region (εi > CvTfr þ hsf), a solid region (εi < CvTfr) and a transition region (CvTfr < εi < CvTfr þ hsf). Particles that have an internal energy corresponding to the transition region will not change their state, i.e., liquid particles remain liquid and solid particles remain solid. The velocity of a solid particle will be set to zero. The exponent of the weight function ne ¼ 1.0 is used for the solid liquid interactions.

210

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

2.3. Temperature dependence of the thermophysical parameters Three thermophysical properties of the working fluid are investigated: mass diffusivity (D), kinematic viscosity (n), and thermal diffusivity (a). The mass diffusivity is investigated by calculating the mean square displacement of the particles in the simulation, a method that has been extensively used in DPD, as well as other particle-based simulations [3,29]. The mass diffusivity can be calculated as:

D ¼ lim

t/∞

2 o 1 n rðtÞ  rð0Þ 2ndim

(17)

ndim is the number of dimensions in the simulation domain, i.e., two in this study. The kinematic viscosity and the thermal diffusivity are investigated in a fashion similar to one another, by using the periodic Poiseuille flow method proposed by Backer et al. [30] and the analogous periodic heat source method proposed by Li et al. [23]. The periodic Poiseuille flow method works as follows, a simulation domain with periodic boundary conditions is divided into two parts: a positive body force is added to all particles in one part of the domain and a negative body force of equal magnitude is added to all particles in the other part of the domain. In the periodic heat conduction method, the external body force is replaced by a periodic heat source (Qext), such as microwave irradiation. The kinematic viscosity and the thermal diffusivity can then be determined as [30,23]:

fext L2x n¼ 8vmax  2 Qext a¼

Lx 2

16



1 1  Tmax  T0 Tmin  T0

(22)

The latent heat of the simulation also needs a physical justification. To determine the magnitude of the latent heat in DPD units, the Clapeyron equation is used. The Clapeyron equation can be applied to any phase change process occurring at constant temperature and pressure [34]. It can be expressed as:



dP dT

 ¼ sat

hsf rf Tm 1  b

(23)

(19)

hsf ¼ h< sf

(20)

K is a constant given by solving the following equation [31].

!

Cv rk< B ¼ 11 nm Cv<

(18)

The DPDe phase-change was evaluated by conducting a simulation of a one-dimensional unsteady solidification problem in a semi-infinite plane. The fluid in the domain is initially set to be the same as the walls in the domain. At the beginning of the simulation, the temperature at the surface is instantaneously lowered to a temperature below the freezing point, TC < Tfr and the fluid particles near the surface begin to freeze. The time evolution of the ice thickness, Y(t), can be expressed as [31]:

  2 T  Tc ls exp K 2 b2 fr 4aS rl hsf b   ¼4 K 2 pffiffiffiffiffiffiffiffiffi pffiffiffiffi paS erf 2Kb aS

L< ¼

where hsf is the change in enthalpy between the solid and liquid phases. b is the ratio between density in the solid and liquid phases: b ¼ rl/rs. Using eq. (23), and the real world length unit from eq. (22), the latent heat in DPD units can be determined as



pffiffi YðtÞ ¼ K t

units long. A much longer length is used in the Y-direction in order to minimize the influence of the top boundary condition on the heat transfer. A time step dt ¼ 0.01 is used. Using a larger time step can cause a slight drift in the total energy of the system [23]. A total of 10,000 time steps is used. The initial temperature is set to T0 ¼ 0.95 and the lower wall temperature is set to TC ¼ (0.80, 0.85, 0.90), i.e., three values below the freezing point Tfr ¼ 0.91. The modeling parameters used to study mean square displacement, periodic Poiseuille flow, periodic heat conduction, and solidification are summarized in Table 1. The physical length scale of the simulation domain can be indirectly determined by the heat capacity of the simulation domain following the approach employed by Li et al. [23]. Setting the DPD heat capacity to Cv ¼ 50,000 one obtains:

   3 2 TH  Tfr lL exp K 4aL  5   pffiffiffiffiffiffiffiffiffi paL erfc 2pKffiffiffiffi a

 3 < r< f Tm L ¼ 2:68  104 < k< rf Tm B

(24)

In summary, the length scale of the simulation is determined by selecting the heat capacity, and the latent heat is determined based on the length scale. The validity of this assumption will be investigated in the results section. 2.4. Fluid flow and freezing in a parallel plate channel The freezing model is applied to analyse the combined fluid flow and the solidification in a parallel plate channel. This is a useful topic, for instance in the study of microfluidic phase change valves [25]. A fluid flow is induced in the parallel plate channel, and after the flow is equilibrated, the temperature of



L

(21) where l is the thermal conductivity, given as l ¼ arCv. The conductivity of the solid is approximated to be the same as that of the liquid. There are two unknowns in eq. (21), the proportionality constant K and the latent heat hsf. The latent heat in the simulation is chosen to correspond to the real world units of latent heat of water, as will be described below. K can be found by solving eq. (21) numerically. A freezing simulation is carried out in order to investigate how well the proposed method agrees with the analytical results obtained by combining eq. (20) with eq. (21) and applying physical parameters corresponding to liquid water [32,33]. The simulation domain is defined to be 12  100(Lx,Ly) DPD

Table 1 DPDe parameters used for mean square displacement (a), periodic Poiseuille flow (b) periodic heat conduction (c), and 1D solidification (d) simulations. The desired outputs are the mass diffusivity (D), the kinematic viscosity (n) the thermal diffusivity (a), and the solidification proportionality constant K. Parameters

a

b

c

d

(Lx, Ly) Cv rc ttot dt fext Qext/Cv T T < ðKÞ Desired output

(30, 30) 1.0  105 1.81 50,000 0.01 e e 0.91e1.2433 273e373 D(T)

(20, 30) 1.0  105 1.81 110,000 0.01 0.2 e 0.91e1.2433 273e373 n(T)

(20, 30) 1.0  105 1.81 110,000 0.01 e 0.001 0.91e1.2433 273e373 a(T)

(12, 100) 1.0  105 1.81 10,000 0.01 e e 0.80e0.90a 240e270a K

a Represents the cold wall temperature. The initial temperature isT0 ¼ 0.95, T0< ¼ 285 K for all cases.

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

the walls is lowered below the freezing temperature. A temperature gradient is formed in the channel, and the fluid starts to solidify. A challenge in using particle based simulations when a temperature gradient is applied is implementing boundary conditions in the flow direction. The most common boundary condition used for DPD simulations including a fluid flow is the periodic boundary condition. This boundary condition has been successfully used in isothermal simulations, but when a temperature gradient is present, the conservation of energy between the inlet and the outlet does not hold. In order to circumvent this problem, a liquid pump [26,27] is introduced in the fluid domain. The simulation domain is divided into three parts: a liquid pump region, a temperature reset region and the simulation domain, as shown schematically in Fig. 1. In the liquid pump region, an external force is added to all particles. The external force will induce a flow and due to viscous effects all particles in the simulation will be set to motion. The temperature reset region is placed directly after the pump region. Here, the temperature of all particles is set to the inlet temperature of the fluid domain. The walls of the solution domain are kept at a constant temperature below the freezing point. The fluid flow simulations are carried out in a 2-D simulation domain with (Lx,Ly) ¼ (150, 50), where the liquid pump and temperature reset regions are Lpump ¼ 37 and Ltr ¼ 3 DPD units long, respectively. A time step dt ¼ 0.01 is used. In all cases, two layers of equally spaced solid particles, paired with a bounce-back boundary condition [35], are used to represent the channel walls. How the freezing in the channel affects the fluid flow is investigated by applying different external forces to the liquid pump region, and different constant temperatures on the solid walls of the simulation domain. For all investigated cases, the simulation is allowed to equilibrate for 10,000 time steps. Then, the temperatures of the walls are instantaneously lowered while the temperature of the inlet is kept the same to be able to model freezing in a thermally developing region. The inlet temperature T0 ¼ 1.0 is used for all cases and the lower wall temperature is allowed to vary between 0.80 < Tc < 0.90. An additional 110,000 time steps are run, and the freezing behavior is analysed. All simulations are averaged over 1000 consecutive runs, employing different random seeds. All simulations are carried out by using an in-house FORTRAN code.

211

3. Results and discussion 3.1. Physical properties of the DPD fluid The results from the investigations of mass diffusivity, kinematic viscosity, and thermal diffusivity are summarized as two dimensionless numbers, the Schmidt number (Sc ¼ n/D), and the Prandtl number (Pr ¼ n/a) as displayed in Fig. 2. It can be seen that both the Schmidt, and Prandtl numbers show good agreement with the experimental data over the entire range of investigated temperatures, which gives a good indication that the investigated fluid behaves as liquid water in the simulation domain. 3.2. Modeling of phase change The phase change model is evaluated by comparing the predicted rate of solid growth in a domain having a sub-zero wall temperature to that of an analytical solution given by eqs. (20) and (21). Fig. 3 shows the time-evolution of the solid/liquid interface together with the corresponding temperature distribution in the domain with TC ¼ 0.80, at three different times, t ¼ 50, 75, and 100. The temperature distribution (bottom half of the figure) shows that even for the last time step, and the lowest wall temperature investigated, the simulation temperature is affected in less than half of the domain (only the bottom half of the domain is shown in the figure). This leads to the conclusion that the artificial upper boundary has little impact on the simulation, which can justify the assumption that the domain behaves like a semi-infinite plane. The upper portion of the figure depicts the positions of the particles in the simulation, solid particles being black and liquid particles being white. Note that only the bottom part of the domain is shown so that the interface can be seen more clearly. An analytical solution would show a perfectly even surface between solid and liquid phases. However, some surface roughness is seen in the simulations as a result of mesoscopic fluctuations of both heat and momentum. Fig. 4 shows a comparison between the analytical growth rate of the solid interface as represented by eqs. (20) and (21), and the interface obtained from the current simulations for three different temperatures of the cold surface: TC ¼ 0.80, 0.85, 0.90. The dashed lines shown in Fig. 4 present analytical solutions for the respective cold surface temperatures, with origin shifted to coincide with the numerical results. It can be seen that the predicted rate of the phase change shows a very good agreement with the analytical results, but while the analytical solution predicts the freezing to start instantaneously, the freezing takes a longer time to start in the simulations. This is due to the latent heat in the equation of state, eq. (16). An amount of internal energy corresponding to the latent heat must first be subtracted from the fluid particles before they can freeze, which does not happen instantaneously. 3.3. Fluid flow in a parallel plate channel

Fig. 1. Schematic image depicting the simulation setup including a liquid pump. The three different regions in the simulation are: A, external applied force (pump), B, temperature reset region, and C, solution domain (figure not to scale). The solution domain is cooled by applying a constant a temperature Tc to the walls in the solution domain.

The simultaneous fluid flow and heat transfer in a channel with sub-zero constant wall temperatures have been investigated. In this section, results regarding velocity profile, temperature distribution, and the mechanisms related to when the channel freezes shut, i.e., when the solid formed at the top and bottom walls come in contact with each other making it impossible for a fluid to flow through the channel, are discussed. When the temperature at the solid wall is set to a constant temperature lower than the freezing point, a layer of ice will start to form at the solid surface. Fig. 5 shows a representative example of the temperature distribution in the domain at three different

212

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

Fig. 2. Schmidt number and Prandtl number: comparison between simulations and experimental results.

Fig. 3. Visualisation of the moving solid liquid interface, and the corresponding temperature profile for three different instants of time: t ¼ 50, 70, 100, and a solid wall temperature TC ¼ 0.80.

instants of time for a set external force (fext ¼ 0.2) and a set constant wall temperature (TC ¼ 0.85). It can be seen that the temperature distribution is very even in the entire domain: the heat transfer is not heavily influenced by the random mesoscopic interactions in

the simulations. The temperature only shows some fluctuations at the interface between solid and liquid (denoted by a thick black line), which will be discussed in more detail in the following sections.

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

Fig. 4. Comparison between analytical [31] and simulated time evolution of the solid/ liquid interface in a semi-infinite domain. The solid lines represent the analytical solution, and the dashed lines have the same slope, but are transposed so that the intersection coincides with the numerical results. Three temperatures of the lower wall are investigated: TC ¼ (0.80, 0.85, 0.90), as represented by blue, red, and black, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In order to evaluate the effects of the flow velocity on the water freezing, three different external forces were simulated with a wall temperature Tc ¼ 0.85. The velocity profile for each external force in the middle of the channel is shown in Fig. 6 after the simulation has run for 35,000 time steps. It is seen that for all three cases, the velocity is zero close to the walls where the liquid has frozen. The same portion of the channel is frozen for all three external forces that are investigated even though the maximum velocity in the center increases with more than 150% when the external force is increased from fext ¼ 0.1 to fext ¼ 0.4, leading to the conclusion that the freezing is primarily dependent on the wall temperature, and not on the velocity of the flow. In order to further clarify the impact of the fluid velocity on the ice growth, simulations were carried out using two different pumping forces and three different wall temperatures, all chosen below the freezing point. How the solid/liquid interface evolves

213

Fig. 6. Velocity profile in the middle of the parallel plate channel for three different added external forces. The constant wall temperature is TC ¼ 0.85 and the fluids have undergone freezing for 35,000 time steps.

close to the fluid inlet is shown in Fig. 7. Three different instants in time are shown: t ¼ 350, where the channel has not closed for any of the simulated wall temperatures, t ¼ 650, shortly after the channel closes for the lowest studied wall temperature (TC ¼ 0.80), and t ¼ 1100, shortly after an intermediate temperature (TC ¼ 0.85) starts to close. The highest wall temperature investigated, TC ¼ 0.90, does not close in the time period simulated. It can be seen in Fig. 7 that the applied external force has little to no effect on the time evolution of the solid interface for any instant in time, regardless of configuration. Instead the freezing is dictated by the temperature of the walls. This agrees with the results of Myers and Low [25], who conclude from a mathematical investigation of the solidification of fluid flow in a microchannel that the closing time does not depend on the flow rate, if the flow rate is sufficiently small. They propose a model for the closing time as a function of the fluid temperature, wall temperature, and channel height/radius only:

tcl f

R2 Tf  Tw

(25)

Fig. 5. Temperature distribution in a channel with external force fext ¼ 0.2, and constant wall temperature TC ¼ 0.85, shown at three different instants of time. The thick line represents the solid/liquid interface.

214

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

follow the prediction of Myers and Low [25] to some extent, but that the freezing takes a longer time than predicted when the wall temperature increases. The longer closing times can be attributed to the roughness of the interface between the solid and the liquid, caused by the mesoscopic fluctuations in the domain, whose effects are larger when the temperature difference between the solid walls and the fluid is smaller. Fig. 9 emphasizes how the temperature difference between the fluid inlet and the solid walls affects the behavior of the solid/liquid interface when the channel is about to close. Two different wall temperatures are investigated, and the solid/liquid interface is shown at three different instants in time for each: when the two solid fronts first touch each other, 20 DPD time units later, and 40 DPD time units later. Interactions between the warmer particles and the colder ones that are close to the solid interface can cause particles to freeze/unfreeze rapidly, which gives rise to the surface roughness seen in Fig. 9. For the lower wall temperature, shown on the left hand side of Fig. 9, the freezing is much faster and the effect from the warmer particles is smaller, causing a more even interface. On the right hand side of Fig. 9, where the higher wall temperature is applied, it can be clearly seen that parts of the channel that are frozen at t ¼ 990 melt and liquid appears again at t ¼ 1010, only to freeze again at t ¼ 1030. This freezing/unfreezing behavior prolongs the time it takes for the parallel plate channel to freeze shut as compared to what would be expected if mesoscopic transfer of heat and momentum was not accounted for.

4. Conclusions

Fig. 7. The position of the solid/liquid interface displayed for fext ¼ 0.1 (solid lines) and fext ¼ 0.2 (dashed lines). Three different wall temperatures are shown, TC ¼ (0.80, 0.85, 0.90) represented in blue, red and black, respectively. Three different instants in time are shown, t ¼ (350, 650, 1100). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In order to investigate how well this relationship applies to the model presented in this paper, a series of simulations were carried out using different wall temperatures ranging from TC ¼ 0.80 to TC ¼ 0.85. The resulting closing time as a function of wall temperature is presented in Fig. 8. It can be seen that the simulated results

Fig. 8. The closing time for the channel for six given constant wall temperatures.

The dissipative particle dynamics with energy conservation is used in order to successfully model the temperature dependent thermophysical properties of simple fluids such as water in a twodimensional domain. An equation of state consisting of solid, liquid, and a transition region between the two is incorporated into the model to enable modeling of freezing. The physical parameters, in the form of length scale, heat capacity and latent heat are linked to the real world units by using scaling analysis, and the model is validated with regards to dimensionless numbers, in the terms of Schmidt number and Prandtl number. The applied equation of state is validated by comparing the simulated rate of solidification to that of an analytical solution. It is found that the rate of phase-change follows the analytical solution for a wide range of lower wall temperatures, leading to the conclusion that the proposed model is viable. However, while the analytical solution predicts the liquid to freeze instantaneously, the simulations take some time to freeze. We attribute this behavior to the latent heat of the model: before the freezing can commence, an internal energy corresponding to the latent heat must be subtracted from the fluid particles, which does not happen instantaneously. The proposed model is then applied to study the simultaneous fluid flow and solidification of liquid water in a straight parallel plate channel with constant sub-zero wall temperatures. The fluid flow is driven by an external liquid pump. It is found that the fluid flow velocity has little to no impact on the rate of the solidification, instead the temperature of the walls dictate the rate of the solidification, which is in line with previous studies in the field. The time it takes for the liquid in a channel to freeze shut is analysed and compared to a correlation previously presented in the literature. It is found that when the constant wall temperature is closer to the freezing temperature, the closing of the channel becomes irregular: certain areas in the channel may freeze and unfreeze rapidly due to the mesoscopic random fluctuations in the domain, and the closing of the channel takes a longer time than what has been previously estimated.

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

215

Fig. 9. The interface between solid and liquid (red lines), for two different wall temperatures is shown for various times as the channel is about to close due to freezing. One temperature contour in the liquid (black) and one in the solid (blue) are shown for comparison. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Acknowledgment Financial supports from the Swedish Institute Scholarship, European Research Council (ERC) (MMFCS-226238) and the National Swedish Research Council are kindly acknowledged. References [1] P. Hoogerbrugge, J. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, EPL Europhys. Lett. 19 (1992) 155e160. [2] K. Malek, M. Eikerling, Q. Wang, Z. Liu, S. Otsuka, K. Akizuki, M. Abe, Nanophase segregation and water dynamics in hydrated nafion: molecular modeling and experimental validation, J. Chem. Phys. 129 (2008) 204702. [3] G. Dorenbos, K. Morohoshi, Chain architecture dependence of pore morphologies and water diffusion in grafted and block polymer electrolyte fuel cell membranes, Energy Environ. Sci. 3 (2010) 1326e1338. [4] S. Yamamoto, S.-A. Hyodo, A computer simulation study of the mesoscopic structure of the polyelectrolyte membrane nafion, Polym. J. 35 (2003) 519e527. [5] I. Pivkin, G. Karniadakis, Accurate coarse-grained modeling of red blood cells, Phys. Rev. Lett. 101 (2008) 118105. [6] X. Li, A.S. Popel, G.E. Karniadakis, Bloodeplasma separation in y-shaped bifurcating microfluidic channels: a dissipative particle dynamics simulation study, Phys. Biol. 9 (2012) 026010.

[7] R.D. Groot, Electrostatic interactions in dissipative particle dynamicsdsimulation of polyelectrolytes and anionic surfactants, J. Chem. Phys. 118 (2003) 11265e11277. [8] L. Rekvig, M. Kranenburg, J. Vreede, B. Hafskjold, B. Smit, Investigation of surfactant efficiency using dissipative particle dynamics, Langmuir 19 (2003) 8195e8205. [9] R. Groot, P. Warren, Dissipative particle dynamics: bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (1997) 4423e4435. ~ ol, Dissipative particle dynamics with energy conservation, EPL [10] P. Espan Europhys. Lett. 40 (1997) 631e636. [11] J. Avalos, A. Mackie, Dissipative particle dynamics with energy conservation, EPL Europhys. Lett. 40 (1997) 141e146. [12] T. Yamada, A. Kumar, Y. Asako, O.J. Gregory, M. Faghri, Forced convection heat transfer simulation using dissipative particle dynamics, Numer. Heat Transf. Part A Appl. 60 (2011) 651e665. [13] E. Abu-Nada, Dissipative particle dynamics simulation of combined convection in a vertical lid driven cavity with a corner heater, Int. J. Therm. Sci. 92 (2015) 72e84. [14] Z.-H. Cao, K. Luo, H.-L. Yi, H.-P. Tan, Energy conservative dissipative particle dynamics simulation of mixed convection in eccentric annulus, Int. J. Heat Mass Transf. 74 (2014) 60e76. [15] T. Yamada, Y. Asako, O.J. Gregory, M. Faghri, Simulation of thermal conductivity of nanofluids using dissipative particle dynamics, Numer. Heat Transf. Part A Appl. 61 (2012) 323e337. [16] R. Qiao, P. He, Simulation of heat conduction in nanocomposite using energy-conserving dissipative particle dynamics, Mol. Simul. 33 (2007) 677e683. [17] R. Qiao, P. He, Mapping of dissipative particle dynamics in fluctuating hydrodynamics simulations, J. Chem. Phys. 128 (2008) 126101.

216

E.O. Johansson et al. / International Journal of Thermal Sciences 101 (2016) 207e216

[18] P. He, R. Qiao, Self-consistent fluctuating hydrodynamics simulations of thermal transport in nanoparticle suspensions, J. Appl. Phys. 103 (2008) 094305.
n, K. Park, M. Faghri, Diffusive-ballistic heat [19] T. Yamada, S. Hamian, B. Sunde transport in thin films using energy conserving dissipative particle dynamics, Int. J. Heat Mass Transf. 61 (2013) 287e292. [20] A. Chaudhri, J.R. Lukes, Energy-conserving dissipative particle dynamics for mesoscopic heat transfer simulations, Ann. Rev. Heat Transf. 17 (2014) 267e302. [21] Z. Wan, H. Chang, S. Shu, Y. Wang, H. Tang, A review on cold start of proton exchange membrane fuel cells, Energies 7 (2014) 3179e3203. [22] M. Vaezi, H. Seitz, S. Yang, A review on 3d micro-additive manufacturing technologies, Int. J. Adv. Manuf. Technol. 67 (2013) 1721e1754. [23] Z. Li, Y.-H. Tang, H. Lei, B. Caswell, G.E. Karniadakis, Energy-conserving dissipative particle dynamics with temperature-dependent properties, J. Comput. Phys. 265 (2014) 113e127. [24] S. Willemsen, H. Hoefsloot, D. Visser, P. Hamersma, P. Iedema, Modelling phase change with dissipative particle dynamics using a consistent boundary condition, J. Comput. Phys. 162 (2000) 385e394. [25] T. Myers, J. Low, An approximate mathematical model for solidification of a flowing liquid in a microchannel, Microfluid. Nanofluid. 11 (2011) 417e428.

[26] S. Ge, Y. Gu, M. Chen, A molecular dynamics simulation on the convective heat transfer in nanochannels, Mol. Phys. 113 (2015) 703e710. [27] A.J. Markvoort, P.A.J. Hilbers, S.V. Nedea, Molecular dynamics study of the influence of wall-gas interactions on heat flow in nanochannels, Phys. Rev. E 71 (2005) 066702. ~ ol, P. Warren, Statistical mechanics of dissipative particle dynamics, [28] P. Espan EPL Europhys. Lett. 30 (1995) 191e196. [29] A.A. Milischuk, B.M. Ladanyi, Structure and dynamics of water confined in silica nanopores, J. Chem. Phys. 135 (2011) 174709. [30] J. Backer, C. Lowe, H. Hoefsloot, P. Iedema, Poiseuille flow to measure the viscosity of particle model fluids, J. Chem. Phys. 122 (2005) 154503.
n, Introduction to Heat Transfer, WIT Press, 2012. [31] B. Sunde € ckh, T. Wetzel, Heat Transfer: Basics and Practice, Springer, 2012. [32] P. Bo [33] M. Holz, S.R. Heil, A. Sacco, Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1H NMR PFG measurements, Phys. Chem. Chem. Phys. (Incorporating Faraday Transactions) 2 (2000) 4740e4742. [34] Y.A. Çengel, M.A. Boles, Thermodynamics e an Engineering Approach, McGraw-Hill, 2011. ~ ol, Boundary conditions in dissipative particle [35] M. Revenga, I. Zuniga, P. Espan dynamics, Comput. Phys. Commun. 121 (1999) 309e311.