Numerical Simulation of the Melting of Solid Particles in Thermal Convection with a Modified Immersed Boundary Method

Abstract

A new immersed boundary method is proposed for the numerical simulation of the melting of solid particles in its own liquid at a high temperature. The main feature of the new method is the use of the modified direct-forcing immersed boundary method for the solution of the flow field and the sharp-interface immersed boundary method for the temperature field. The accuracy of the proposed method is validated via three problems: the sedimentation of a non-melting particle, the melting of a fixed particle under mixed thermal convection, and the sedimentation of a melting particle. The method is then applied to the investigation of the effects of various parameters, the particle interactions and the particle shape on the particle melting time. A correlation for the melting time of a circular particle in forced thermal convection is established as a function of the Reynolds, Prandtl, and Stefan numbers. The melting time of a particle in mixed thermal convection first increases and then decreases, as the Grashof number increases. The effects of the particle interactions on the melting time are complicated due to the natural convection between two particles. The sufficiently strong natural convection can even render the downstream particle melt faster than the single particle. For the same particle area, the elliptic particle with the aspect ratio being around 1.4 melts most slowly.

1. Introduction

Melting phenomenon is commonly encountered in nature and engineering applications such as in the melting of ice in the ocean, the particle melting in the plasma spraying, and the unwanted melting of catalyst particles in a fluidized bed. When a solid particle is immersed in a fluid medium at a temperature higher than the melting temperature of the solid material Tm, melting takes place. There have been some experimental and numerical studies on the melting of particles. Kranse and Schenk investigated experimentally the melting of solid benzene spheres in liquid benzene [1], and Schenk and Schenkels studied the melting of ice spheres in water [2]. Argyropoulos et al. [3,4] examined the melting time and the heat flux on the surface of melting and rotating spheres of steel and aluminum immersed in the liquid of the same material, respectively. Hao and Tao [5,6] investigated the melting and heat transfer characteristics of ice particles in flowing water. Gan et al. [7] studied numerically the sedimentation of melting solid particles with an arbitrary Lagrangian–Eulerian (ALE) finite-element method. In their simulations, the surface of a melting particle was tracked by Lagrangian points, and the melting rate was determined from the local heat flux. Melissari and Argyropoulos [8,9] studied both numerically and experimentally the influences of the various parameters on the melting of an immersed sphere in a liquid metal and argued that the sphere melting technique could be used to measure the velocity of liquid metal. Kumar and Roy [10,11] simulated numerically the melting of a steel spherical particle in its own liquid with a finite-volume method, and analyzed the melting rate, the shape of the particle, and the heat transfer characteristic. The total melting time of particles in thermal convection is of practical interest, but its parametrical dependence is not well understood. For example, to our knowledge, no works on the effects of the particle interaction and the particle initial shape on the melting time have been reported.

The purposes of the present study are twofold: the first is to present a new numerical method for modeling the melting of solid particles in liquids of the same material in the presence of thermal convection, and the second is to gain some new insights in the melting time of particles via the applications of the proposed method. In the previous numerical simulations of Argyropoulos et al. [3,4] and Kumar and Roy [10,11], the surfaces of melting particles were determined from the liquid volume fraction, which was modeled as the function of the mixture temperature or enthalpy. In the present study, we adopted the same melting model as Gan et al. [7], who tracked the melting surface with Lagrangian points so that the interface between the solid and liquid phases was kept sharp, and the melting rate was determined from the local heat flux. Gan et al. [7] used an arbitrary Lagrangian–Eulerian (ALE) finite-element method [12] to solve the flow and temperature fields, whereas here we developed a combined immersed-boundary (IB) method to solve the problems. The IB method was originally proposed by Peskin to simulate the motion of flexible membranes in the human heart [13]. It is highly efficient for problems with moving interfaces due to the fact that a non-body-fitted Cartesian grid is used to solve the flow field and has been extensively studied and applied to various problems [14,15,16,17,18,19,20,21,22,23]. In the present study, we employed a direct-forcing-type IB method for the fluid/particle motion, and a sharp-interface-type IB method for the fluid/particle heat transfer, in accordance with the melting model we adopted. The new method can be regarded as an extension of our previous combined direct-forcing fictitious domain (FD) [24] and sharp-interface IB method for the conjugate heat transfer between fluids and non-melting particles [25]. The main modification is that the FD method is replaced with the IB method for the solution of the flow fields. Although the FD method can be regarded as an IB method in a wide sense, here, we distinguish between them to emphasize the modification: in our original FD method, the Lagrangian points are distributed inside the particle boundary to enforce the rigid-body motion of the fictitious fluids over the particle domain, whereas in the present direct-forcing IB method, the Lagrangian points are distributed only on the solid/fluid interface. The reason for the modification is that for the melting problem, the shape of the particle is changing, and it is a difficult task to distribute the Lagrangian points evenly over the particle domain. The sharp-interface (SI) IB method [26] was developed to handle the Poisson equation with jump coefficients across the interface. Its main feature is that the jump condition on the interface is used to modify the discretization of the differential operators on the Cartesian grids in the immediate vicinity of the interface, and its advantages is that the interface is kept sharp and the jump condition on the interface is accurately captured [27]. The sharp-interface IB method has an advantage over the direct-forcing IB (or FD) method in that the strong discontinuity across the interface could be better captured, and it has a disadvantage in that the scheme is more complex and less efficient. The SI method has been applied to the droplet [27,28], dielectrophoresis problems [29], and the heat transfer between fluids and non-melting particles [25]. We apply it here to the melting problem.

The rest of this paper is organized as follows. In the next section, the numerical method is described. In Section 3, we verify our code via three problems: the sedimentation of a non-melting particle, the melting of a fixed particle under mixed thermal convection, and the sedimentation of a melting particle. As discussed in Section 4, the new method was employed to investigate the melting time of a single circular particle, two circular particles, and an elliptic particle, respectively. Concluding remarks are given in the final section.

2. Numerical Model

2.1. Direct-Forcing IB Method for the Flow Fields

As mentioned earlier, a direct-forcing immersed boundary (IB) method was used to solve the flow field. Figure 1 is a schematic of the IB method, in which P represents the particle domain; $\partial P$, its boundary; $\Omega$, the entire computational domain; and $\Gamma$, its boundary. The spirit of the IB method is that the solid domain is filled with the fictitious fluid and a pseudo-body force (i.e., Lagragian multiplier) is introduced over the particle boundary to enforce the no-slip condition there.

To reasonably simplify the computation, the Boussinesq approximation is used to deal with the effect of temperature variation in the flow field. The change in temperature is assumed not to influence the properties of fluid medium, except for the density in the gravitational term, which has the following form:

$${\rho }_f={\rho }_{f_0}\left[1-{\beta }_f\left(T_f-T_0\right)\right],$$

(1)

where ${\rho }_{f_0}$ represents the reference density of the fluid at the reference temperature $T_0$, and ${\beta }_f$ is the fluid heat expansion coefficient.

The characteristic scales used for the non-dimensionlization scheme are $L_c$ for length, $U_c$ for velocity, $L_c/U_c$ for time, ${\rho }_{f0}{U}_c^2$ for pressure p, and ${\rho }_{f0}{U}_c^2/L_c$ for the Lagrange multiplier $\lambda$. Then, the dimensionless governing equations for the flow field and the particle motion can be written as follows:

$${\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}\cdot \nabla \mathbf{u}={\displaystyle \frac{{\nabla }^2\mathbf{u}}{\mathit{Re}}}-\nabla p-{\displaystyle \frac{Gr}{{\mathit{Re}}^2}}{\overline{T}}_f{\displaystyle \frac{\mathbf{g}}{g}}+\mathit{\lambda }\begin{array}{cccc}& & in& \Omega \end{array}$$

(2)

$$\nabla \cdot \mathbf{u}=0\begin{array}{cccc}& & in& \Omega \end{array}$$

(3)

$$\mathbf{u}=\mathbf{U}+{\mathit{\omega }}_s\times \mathbf{r}\begin{array}{ccc}& & \mathit{on}\begin{array}{cc}& \partial P\end{array}\end{array}$$

(4)

$${\rho }_r{V}_p^{*}{\displaystyle \frac{d\mathbf{U}}{dt}}=-{\displaystyle {\int }_{\partial p}\mathit{\lambda }}dx+{\displaystyle \frac{d}{dt}}{\displaystyle {\int }_p\mathbf{u}}d\mathbf{x}+\left({\rho }_r-1\right)V_p^{*}Fr{\displaystyle \frac{\mathbf{g}}{g}}$$

(5)

$${\rho }_r{J}^{*}{\displaystyle \frac{d{\mathit{\omega }}_s}{dt}}=-{\displaystyle {\int }_{\partial p}\left(\mathbf{r}\times \mathit{\lambda }\right)dx}+{\displaystyle \frac{d}{dt}}{\displaystyle {\int }_p\left(\mathbf{r}\times \mathbf{u}\right)d\mathbf{x}}$$

(6)

Here, u is the fluid velocity. U and ω_s are the particle translational and angular velocities, respectively. r is the position vector with respect to the particle center. ${\rho }_r$ represents the solid–fluid density ratio defined by ${\rho }_r={\rho }_s/{\rho }_{f_0}$, with ${\rho }_s$ being the solid density. $V_p^{\ast }$ is the dimensionless particle volume (or area in two dimensions) defined by $V_p^{\ast }=M/\left({\rho }_s{L}_c^d\right)$, with M being the particle mass and d being the dimensionality of the problem. ${\mathbf{J}}^{\ast }$ is the dimensionless moment of inertia tensor defined by ${\mathbf{J}}^{\ast }=\mathbf{J}/\left({\rho }_s{L}_c^{d+2}\right)$. Re denotes the Reynolds number, defined by $Re={\rho }_{f0}{U}_c{L}_c/\mu$, with $\mu$ being the fluid viscosity, and Fr represents the Froude number, defined by $Fr=g{L}_c/U_c^2$. Gr is the Grashof number, defined by $Gr={\rho }_{f0}^2{\beta }_f{L}_c^{3}g\left(T_2-T_1\right)/{\mu }^2$, where T₁ and T₂ are two reference temperatures, and for the melting problem, we let T₂ be the initial temperature of the fluid medium T₀ and T₁ be the melting temperature of the particle T_m. The dimensionless temperature can then be defined by $\overline{T}=\left(T-T_m\right)/\left(T_0-T_m\right)$. Note that we do not assume the Boussinesq approximation to the fictitious fluid and the solid inside the particle boundary, so that the buoyant term arising from the temperature variation in Equation (2) should not applied to the particle domain and is explicitly not included in Equation (5), unlike our previous FD formulation [25,30]. Our numerical tests show that the two schemes yield roughly the same results.

A fractional-step time scheme is used to decouple the problem of Equations (2)–(6) into a fluid flow sub-problem and a particle motion sub-problem as follows:

Fluid-flow sub-problem for ${\mathbf{u}}^{\ast }$ and p:

$${\displaystyle \frac{{\mathbf{u}}^{*}-{\mathbf{u}}^n}{\Delta t}}-{\displaystyle \frac{{\nabla }^2{\mathbf{u}}^{*}}{2\mathit{Re}}}=-\nabla p-{\displaystyle \frac{1}{2}}[3{(\mathbf{u}\cdot \nabla \mathbf{u}+{\displaystyle \frac{Gr}{{\mathit{Re}}^2}}{\overline{T}}_f{\displaystyle \frac{\mathbf{g}}{g}})}^n-{(\mathbf{u}\cdot \nabla \mathbf{u}+{\displaystyle \frac{Gr}{{\mathit{Re}}^2}}{\overline{T}}_f{\displaystyle \frac{\mathbf{g}}{g}})}^{n-1}]+{\displaystyle \frac{{\nabla }^2{\mathbf{u}}^n}{2\mathit{Re}}}+{\mathit{\lambda }}^n,$$

(7)

$$\nabla \cdot {\mathbf{u}}^{*}=0.$$

(8)

Particle sub-problem for ${\mathbf{U}}^{n+1}$, ${\mathit{\omega }}_s^{n+1}$, ${\mathbf{u}}^{n+1}$, and ${\mathit{\lambda }}^{n+1}$:

$${\displaystyle \frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^{*}}{\Delta t}}={\mathit{\lambda }}^{n+1}-{\mathit{\lambda }}^n,$$

(9)

$${\mathbf{u}}^{n+1}={\mathbf{U}}^{n+1}+{\mathit{\omega }}_{\mathrm{s}}^{n+1}\times \mathbf{r},$$

(10)

$${\rho }_{\mathrm{r}}{V}_{\mathrm{p}}^{*}{\displaystyle \frac{{\mathbf{U}}^{n+1}-{\mathbf{U}}^n}{\Delta t}}=-{\displaystyle {\int }_{\partial p}{\mathit{\lambda }}^n}\mathrm{d}\mathbf{x}+{\displaystyle \frac{{\displaystyle {\int }_p{\mathbf{u}}^{*}}d\mathbf{x}-{{\displaystyle {\int }_p\mathbf{u}}}^{n}d\mathbf{x}}{\Delta t}}+({\rho }_{\mathrm{r}}-1)V_{\mathrm{p}}^{*}Fr{\displaystyle \frac{\mathbf{g}}{g}},$$

(11)

$${\rho }_{\mathrm{r}}{J}^{*}{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{s}}^{n+1}-{\mathit{\omega }}_{\mathrm{s}}^n}{\Delta t}}=-{\displaystyle {\int }_{\partial p}\mathbf{r}\times {\mathit{\lambda }}^n}\mathrm{d}\mathbf{x}+{\displaystyle \frac{{\displaystyle {\int }_p\mathbf{r}\times {\mathbf{u}}^{*}}d\mathbf{x}-{\displaystyle {\int }_p\mathbf{r}\times {\mathbf{u}}^n}d\mathbf{x}}{\Delta t}}.$$

(12)

Note that from Equations (5)–(11), we have used the approximation $-{\displaystyle {\int }_{\partial p}{\mathit{\lambda }}^n}\mathrm{d}\mathbf{x}+{\displaystyle \frac{1}{\Delta t}}{\displaystyle {\int }_p{\mathbf{u}}^{*}}d\mathbf{x}={\displaystyle {\int }_{\partial p}{\mathit{\lambda }}^{n+1}}\mathrm{d}\mathbf{x}+{\displaystyle \frac{1}{\Delta t}}{\displaystyle {\int }_p{\mathbf{u}}^{n+1}}d\mathbf{x}$, with is derived from Equation (9). A similar technique was used in the derivation of our direct-forcing FD formulation, and thus, the above IB and our previous FD formulas for the particle velocities are consistent. We are not aware of any other IB works in which this technique has been used.

For the particle sub-problem, the particle velocities ${\mathbf{U}}^{n+1}$ and ${\mathit{\omega }}_{\mathrm{s}}^{n+1}$ are first obtained from Equations (11) and (21), and then the Lagrange multiplier is calculated from Equations (9) and (10) as

$${\mathit{\lambda }}^{n+1}={\displaystyle \frac{{\mathbf{U}}^{n+1}+{\mathit{\omega }}_{\mathrm{s}}^{n+1}\times \mathbf{r}-{\mathbf{u}}^{*}}{\Delta t}}+{\mathit{\lambda }}^n.$$

(13)

Finally, the fluid velocity ${\mathbf{u}}^{n+1}$ is updated from Equation (9), which is a step that can also be skipped since ${\mathbf{u}}^{*}$ is a good approximation to ${\mathbf{u}}^{n+1}$, and in the previous direct-forcing IB method [20], ${\mathbf{u}}^{*}$ was directly taken as ${\mathbf{u}}^{n+1}$. The fluid flow sub-problem of Equations (7) and (8) is solved with a finite-difference projection method on a half-staggered grid, as in our direct-forcing FD method [24]. The present direct-forcing IB scheme differs from the previous direct-forcing FD method in that the Lagrange multiplier is defined on the particle boundary instead of over the particle domain, and the inertial force and torque of the fictitious fluid are obtained from the direct integration over the particle domain instead of using the rigid-body motion condition.

2.2. Sharp-Interface IB Method for the Temperature Field

The dimensionless governing equation for the temperature field is

$${\rho }^{\ast }{c_p}^{\ast }({\displaystyle \frac{\partial \overline{T}}{\partial t}}+u\cdot \nabla \overline{T})={\displaystyle \frac{1}{Pe}}\nabla \cdot \left(k^{\ast }\nabla \overline{T}\right)\begin{array}{cccc}& & in& \Omega \end{array}.$$

(14)

For the melting problem, one can assume the following boundary condition:

$$\begin{array}{ccc}\overline{T}={\overline{T}}_m& & on\begin{array}{cc}& \partial P\end{array}\end{array}.$$

(15)

Here, ${\rho }^{\ast }$, ${c_p}^{\ast }$, and $k^{\ast }$ are the dimensionless density, heat capacity, and heat conductivity in the individual medium, with the reference of ${\rho }_{f0}$, $c_{pf}$, and $k_f$, respectively. The subscript ‘f’ represents the fluid medium. From our definition, $k^{\ast }$ is unity in the fluid medium and equals the solid–fluid heat conductivity ratio $k_r$ (defined by $k_r=k_s/k_f$) in the particle domain. Pe denotes the Peclet number defined by $Pe={\rho }_{f0}{c}_{pf}{U}_c{L}_c/k_f$. Note that $Pe=RePr$, where Pr is the Prandtl number defined by $Pr=\mu c_{pf}/k_f$. ${\overline{T}}_m$, is the dimensionless melting temperature and equals zero from our definition of the dimensionless temperature.

In the present study, we assume that the thermal conductivity of the solid is large, and thus, the inner temperature of the particle is almost homogeneous and equals the melting temperature during the melting process. Since we are mainly concerned with the particle total melting time, we set the initial particle temperature to be the melting temperature for simplicity. The inner temperature then remains constant.

For the sharp-interface method in two-dimensions, the jump condition for the temperature and the temperature gradients in the coordinate directions on the solid–fluid interface need to be used instead.

$${\left[\overline{T}\right]}_I={\overline{T}}_I^{+}-{\overline{T}}_I^{-}=0\begin{array}{ccc}& & on\begin{array}{cc}& \partial P\end{array}\end{array},$$

(16)

$${\left[k^{\ast }{\overline{T}}_x\right]}_I={\left(k^{\ast }{\overline{T}}_x\right)}_I^{+}-{\left(k^{\ast }{\overline{T}}_x\right)}_I^{-}=a_I\left(t\right)\begin{array}{ccc}& & on\begin{array}{cc}& \partial P\end{array}\end{array},$$

(17)

$${\left[k^{\ast }{\overline{T}}_y\right]}_I={\left(k^{\ast }{\overline{T}}_y\right)}_I^{+}-{\left(k^{\ast }{\overline{T}}_y\right)}_I^{-}=b_I\left(t\right)\begin{array}{ccc}& & on\begin{array}{cc}& \partial P\end{array}\end{array}.$$

(18)

In the above Equations, ${\left[\right]}_I$ represents a jump across the interface.

We assume that the jump values of the temperature gradients in the normal and tangential directions are c_I and d_I, respectively, i.e.,

$${\left[k^{\ast }{\overline{T}}_n\right]}_I=c_I\left(t\right);{\left[k^{\ast }{\overline{T}}_t\right]}_I=d_I\left(t\right)\begin{array}{ccc}& & on\begin{array}{cc}& \partial P\end{array}\end{array}.$$

(19)

Then, a_I and b_I can be calculated from c_I and d_I:

$$\begin{array}{l}{a}_I=c_I{n}_x+d_I{n}_y\\ b_I=c_I{n}_y-d_I{n}_x\end{array},$$

(20)

where n_x and n_y are two components of the outward unit normal on the particle boundary. Since the temperature on the interface is a constant, the temperature gradient along the tangential direction is zero on each side of the interface, and thus, d_I in Equation (20) is zero. Only $c_I$ is unknown a priori at any time step, and iteration is used to circumvent the problem. We set the temperature on the grids inside the particle boundary to be ${\overline{T}}_m$, and $k_r$ = 100, which ensures that the boundary conditions of Equations (16)–(18) are essentially equivalent to Equation (15).

A second-order semi-implicit Adam–Bashforth/Crank–Nicolson scheme is employed to discretize Equation (14) in time:

$${\rho }^{\ast }{c_p}^{\ast }({\displaystyle \frac{{\overline{T}}^{n+1}-{\overline{T}}^n}{\Delta t}})={\displaystyle \frac{1}{2}}\left[\nabla \cdot \left({\displaystyle \frac{k^{\ast }}{Pe}}\nabla {\overline{T}}^{n+1}\right)+\nabla \cdot \left({\displaystyle \frac{k^{\ast }}{Pe}}\nabla {\overline{T}}^n\right)\right]-{\displaystyle \frac{1}{2}}\left[3{(u\cdot \nabla \overline{T})}^n-{(u\cdot \nabla \overline{T})}^{n-1}\right]$$

(21)

The temperature Equation (21) with the boundary condition of Equations (16)–(18) can be solved with the sharp-interface method. The reader is referred to [25,29] for the principle and discretization schemes of the sharp-interface method.

2.3. Melting Model

The rate of melting is determined from the local heat flux on the particle surface [7]:

$${\rho }_s{h}_f{\displaystyle \frac{d{r}_n}{dt}}=-k{\displaystyle \frac{\partial T}{\partial n}},$$

(22)

where h_f is the latent heat of fusion, and r_n is the coordinate of the particle surface along the local normal n.

The dimensionless form of Equation (22) is

$${\displaystyle \frac{{\rho }_r}{St}}{\displaystyle \frac{d{r}_n}{dt}}=-{\displaystyle \frac{1}{RePr}}{\displaystyle \frac{\partial \overline{T}}{\partial n}},$$

(23)

where St represents the Stefan number, defined by $St=c_p\Delta T/h_f$, and $\partial \overline{T}/\partial n$ represents the dimensionless temperature gradient, which is actually the local Nusselt number (Nu).

In the melting process, the Lagrangian points on the surface of the particle retract along the local normal n and become closer to each other. The computation becomes unstable when two points become too close. The following mergence technique is used to overcome the problem. We define $d_{s,s+1}$ as the distance between the consecutive points s and s + 1. Once $d_{s,s+1}<d_{\min }$ at any time step, these two points are merged to a new point s + 1/2, which has the position and Lagrange multiplier as follows:

$$x_{s+1/2}=\left(x_s+x_{s+1}\right)/2,$$

(24)

$${\lambda }_{s+1/2}=\left({\lambda }_s+{\lambda }_{s+1}\right)/2.$$

(25)

For the freely moving particle in two dimensions, its area S_p, mass center (x_p, y_p) and moment of inertia J_p need to be computed. Since we have the Lagrangian points along the particle boundary, it is convenient to calculate these quantities using the boundary integration from the Gaussian law:

$$S_p={\displaystyle {\int }_{P}dA}={\displaystyle {\int }_{\partial P}\frac{1}{2}\left(x{n}_x+y{n}_y\right)}dl,$$

(26)

$$\left.\begin{array}{l}{x}_p={\displaystyle \frac{{\displaystyle {\int }_{P}xdA}}{S_p}}={\displaystyle \frac{{\displaystyle {\int }_{\partial P}\frac{1}{2}{x}^2{n}_x}dl}{S_p}}\\ y_p={\displaystyle \frac{{\displaystyle {\int }_{P}ydA}}{S_p}}={\displaystyle \frac{{\displaystyle {\int }_{\partial P}\frac{1}{2}{y}^2{n}_y}dl}{S_p}}\end{array}\right\},$$

(27)

$$J_p={\displaystyle {\int }_P({x_r}^2+{y_r}^2)dA}={\displaystyle {\int }_{\partial P}\frac{1}{3}\left({x_r}^3{n}_x+{y_r}^3{n}_y\right)}dl.$$

(28)

Here, $\left(x_r,y_r\right)$ represents the position vector with respect to the particle mass center.

3. Validation of Our Code

The purpose of this section is to verify our numerical method. The following three test problems are considered: the sedimentation of a non-melting particle, the melting of a fixed particle under mixed thermal convection, and the sedimentation of a melting particle.

3.1. Sedimentation of a Non-Melting Particle

The accuracy of the previous direct-forcing FD method has been fully validated [24], and here, we intend to verify the direct-forcing IB code by comparing the results of the sedimentation of a non-melting particle in a box to those from the direct-forcing FD method. The characteristic length is taken as the particle diameter d (which means that d = 1 in the code) and the characteristic velocity as $U_c=\sqrt{{\displaystyle \frac{\pi d}{2}}\left({\rho }_r-1\right)g}$ (Yu and Shao [24]). The size of the enclosure is $4\times 8$. The velocities on the walls are zero. Constant temperatures are enforced on the bottom and top walls: $T_b=1$ and $T_t=0$. Side walls are adiabatic. The physical parameters are ${\rho }_r$ = 1.1, $k_r$ = 1, $c_{pr}$ = 0.09, Re = 40, Pr = 1, and $G_r=1\times {10}^5$. The mesh size is h = d/16, and the dimensionless time step is Δt = 0.005. The particle is released from the position of (2, 6).

Figure 2 shows the fluid and temperature fields obtained using the FD-SI method and the IB-SI method, respectively. The flow fields are solved with the FD and IB methods, respectively, and the temperature fields are solved with the SI method for the non-melting case [25]. We see that two results on the fluid fields outside the particles and the particle position agree well with each other.

3.2. Melting of a Fixed Particle Under Mixed Thermal Convections

The computational domain is a box of $4d\times 16d$, d being the particle initial diameter. The uniform flow condition is imposed on the bottom inlet and side walls, and the outflow condition is imposed on the top outlet. The particle initial diameter and the upward mainstream velocity are taken as the characteristic length and velocity, respectively. The particle is fixed, and its center is located at (2, 7.5). The dimensionless particle temperature is zero, and the dimensionless initial fluid temperature is unity. The temperatures on the bottom inlet and side walls are unity, and the vanishing heat-flux condition is imposed on the top outlet. The control parameters are the same as for Gan et al. [7]: Re = 20, ${\rho }_r=1$, Pr = 0.7, Ra = 800 (Ra = GrPr), and St = 0.025125.

Figure 3 shows the comparison of our shapes of the melting particle at t = 6.57 and t = 39.83 to those obtained by Gan et al. [7], and Figure 4 shows the comparison of the local Nusselt number on the particle surface at t = 6.57. In our simulations, two mesh sizes ($h=d/16$ and $h=d/32$) were used for the grid convergence check. The numbers of corresponding Lagrangian points on the particle surface were 52 and 100, respectively. From Figure 3 and Figure 4, our results obtained from two grids are in good agreement with each other and those of Gan et al. [7]. In the simulations below, the finer mesh was used.

As mentioned earlier, the consecutive Lagrangian points are merged when they approach very closely to void the computational instability. The effect of this mergence technique on the accuracy is examined. We set the critical distance as $d_{\min }=0.5h$. Figure 5 shows a comparison of the particle shapes obtained by merging and not merging the points, respectively. The results are in good agreement, even at t = 64 when the particle area decreases below 10 percent of its original value and the number of Lagrangian points reduces to 29. The computation becomes unstable after t = 64 without merging points.

3.3. Sedimentation of a Melting Particle

In the preceding problem, the melting particle is fixed, and now we consider the sedimentation of a melting particle. The computational domain is a box of $4d\times 16d$. To mimic an infinitely long channel, we shift the flow fields and the particle position upward one mesh distance once the particle falls below a vertical position that is 6d higher than the bottom inlet, so that the computational domain looks like moving with the particle. The particle initial diameter d and the velocity $U_c=\sqrt{{\displaystyle \frac{\pi d}{2}}\left({\rho }_r-1\right)g}$ are taken as the characteristic length and velocity, respectively. The bottom and side walls have constant temperatures. The particle is released from (2, 12). The control parameters are ${\rho }_r=1.00232$, Re = 38.7, Pr = 0.7, Gr = 100, and St = 0.025125. The particle Reynolds number is defined by ${\mathrm{Re}}_p=\mathrm{Re}\left|U_p\right|$. Without the consideration of melting, the parameters above lead to a steady Re_p of about 21, as with Gan et al. [7].

Figure 6 shows the evolutions of the particle mass and Reynolds number. The results are in excellent agreement with those of Gan et al. [7]. The particle accelerates at initial stage of the sedimentation and then decelerates as it melts and its size is reduced.

4. Applications of the Method

We are interested in the effects of various parameters on the particle melting time. The effects of the Reynolds, Prandtl, and Stefan numbers on the melting time of a single circular particle in forced thermal convection, and the effects of the Grashof number in case of mixed thermal convection, the effects of the particle interactions, and the effects of initial particle shape on the melting time were investigated.

4.1. Melting of a Single Circular Particle in Forced Thermal Convection

A particle with d = 1, ${\rho }_r=1$, and $T_p=0$ is fixed at the center of a box of $\left(0,4\right)\times \left(0,16\right)$. The boundary conditions are the same as those in Section 3.2. Without the effect of the natural convection, the dimensionless melting time of the particle $t_m$ is dependent on Re, Pr, and St. The parameters studied in the range of $40\le Re\le 200$, $0.1\le Pr\le 1$, and $0.01\le St\le 1$, and the calculated melting times are listed in

Table 1. Particle melting times at the specified Re, Pr, and St in the case of forced thermal convection.

St	Re = 40	Re = 80	Re = 120	Re = 160	Re = 200	Pr
0.01	69.9	108.5	139.2	164.5	186	0.1
0.01	145.4	222.9	282.1	329.6	369.9	0.3
0.01	257.1	389.5	486.3	558.7	620.6	0.7
0.01	327.4	493.6	611	694.9	766.7	1
0.02	35.1	54.5	70.2	83.2	94.4	0.1
0.02	73.1	111.9	142.1	166.3	186.7	0.3
0.02	129.1	195.7	244.2	282	311.7	0.7
0.02	164.3	247.9	307.2	349	385.3	1
0.1	7.2	11.4	14.7	17.5	19.9	0.1
0.1	15.2	23.3	29.6	35	39.7	0.3
0.1	26.8	40.6	51	59.7	66.4	0.7
0.1	34.1	51.2	63.7	73.6	81.3	1
0.5	1.41	2.5	3.4	4.1	4.7	0.1
0.5	3.5	5.5	7	8.2	9.3	0.3
0.5	6.3	9.6	12	14	15.7	0.7
0.5	8.1	11.9	14.8	17.3	19.3	1
1	0.72	1.34	1.86	2.27	2.64	0.1
1	1.89	3.1	4.03	4.77	5.41	0.3
1	3.6	5.53	6.91	8.05	9.1	0.7
1	4.63	6.99	8.56	9.96	11.24	1

. Using the least squares method, we obtain the optimal fitting expression for the melting time: $t_m=0.4589R{e}^{0.5976}P{r}^{0.6561}S{t}^{-0.9332}$. The computational results and the fitting line are plotted in Figure 7, and one can see that the fitting is remarkably good.

We now study the effect of the natural convection on the melting of the particle. The value of the Grashof number changed from 0 to 6000, while the other parameters are fixed: Re = 20, Pr = 0.7, and St = 0.025125. The particle melting times at different Grashof numbers are plotted in Figure 8. The melting time first increases, as Gr increases from 0 to 3500, and then decreases with further increasing Gr. The reason for the delay of the melting at Gr ≤ 3500 is that the cold fluid around the particle tends to induce downward natural convection, which weakens the upward forced flow and thereby the heat transfer between the particle and the fluid. As shown in Figure 9, the natural and forced convections have comparable intensity and almost cancel each other at the region around the particle at Gr = 3500, resulting in the longest melting process. Figure 8 shows that the melting time decreases, as Gr increases beyond 3500, which is obviously caused by the fact that the downward natural convection predominates over the upward forced convection for Gr > 3500 (Figure 9).

4.2. Melting of Particle Pairs

The computational domain is enlarged to be $\left(0,4\right)\times \left(0,32\right)$. One particle is fixed at (2, 10), and the position of the other particle is varied to examine the effects of the particle inter-distance, defined by $s=\sqrt{{\left(x_{p1}-x_{p2}\right)}^2+{\left(y_{p1}-y_{p2}\right)}^2}$. We set Re = 20, Pr = 0.7, and St = 0.025125. Gr is set to be 1000 and 4000, corresponding to the forced and natural convection domination regimes, respectively.

Figure 10 shows the particle melting time as a function of the distance s between the particle pair at Gr = 1000. Figure 11 shows the temperature and flow fields at t = 40 for s = 2, 6, 9, and 12. From Figure 10, the particle interaction plays a negative role in the melting of both particles, which is not surprising given both particles’ behavior as heat sinks and colder fluid fills the space between two particles. As the inter-distance increases beyond 7, the effect of the downstream (upper) particle on the upstream (lower) particle becomes negligible. By contrast, the effect of the upstream particle on the downstream particle is more significant, since the downstream particle is in the long wake of the colder fluid. An interesting observation is that the melting time of the downstream particle does not decrease monotonically with the inter-distance; the melting time at s = 9 is larger than those at s = 6 and s = 12, as shown in Figure 10. The reason for this phenomenon is that at s = 6, the colder fluid between the particle pair induces a relatively strong downward natural convection, and at s = 12, the downstream particle faces an upward forced convection, whereas at s = 9, the forced and natural convections cancel each other around the downstream particle (Figure 11).

The results on the melting time, and the flow and temperature fields at Gr = 4000 are shown in Figure 12 and Figure 13, respectively. Because the natural convection dominates the flow at Gr = 4000, the ‘upstream’ (lower) particle is immersed in the downward moving colder fluid (Figure 13); thus, its melting time becomes larger compared to the case of a single particle and decreases with increasing distance between the particles. By contrast, the melting of the upper particle can even be accelerated when two particles are spaced neither too close nor too far away, say, 2 ≤ s ≤ 4, as shown in Figure 12. For the moderate distance between two particles, the downward natural convection in this space becomes stronger compared to the single-particle case (as indicated by the comparison between Figure 13b,c), resulting in the accelerated melting of the upper particle.

4.3. Melting of an Elliptic Particle

The dimension of the computational domain is $\left(0,16\right)\times \left(0,16\right)$. An elliptic particle is fixed at the center. The boundary conditions are the same as the above problem. The characteristic length is taken as the effective diameter of the particle with the initial area of S_p, namely, $d_{\mathrm{e}}=2\sqrt{{\mathrm{S}}_{\mathrm{P}}/\pi }$. The parameters are Re = 10, Pr = 0.7, St = 0.025125, ${\rho }_{\mathrm{r}}=1$, and Gr = 0. Here, we intend to investigate the effects of the particle shape and orientation on the melting rate. For simplicity, we assume that the particle axis is aligned horizontally or vertically and define the particle aspect ratio as A = a_x/a_y, where a_x and a_y being the length of the axis in the x and y directions, respectively.

The results on the melting time, and the flow and temperature fields for different aspect ratios are shown in Figure 14 and Figure 15, respectively. Figure 14 shows that the melting time is longest for the particle with A ≈ 1.4 and decreases as the particle becomes slenderer. The observations can be explained with two facts. First, the elliptic particle has a larger circumference for the heat transfer with the fluid, and thus generally melts faster, compared to the circular particle with the same area. The heat transfer efficiency is much higher at the leading and side parts of the particle boundary than at the rear part, as indicated in the flow and temperature fields in Figure 15. Therefore, the elliptic particle with the same shape melts faster when its longer axis is aligned with the streamwise direction, and the particle that melts most slowly is not the circular particle but the elliptic particle of A ≈ 1.4, whose longer axis is aligned with the cross-stream direction.

5. Conclusions

We presented a new immersed boundary method for the numerical simulation of the melting of solid particles in its own liquid. The new method is the combination of the modified direct-forcing immersed boundary method for the flow field and the sharp-interface immersed boundary method for the temperature field. The accuracy of the proposed method was verified via three problems: the sedimentation of a non-melting particle, the melting of a fixed particle under mixed thermal convection, and the sedimentation of a melting particle. The method was employed to investigate the effects of various parameters, the particle interactions, and the particle shape on the particle melting time. The following main conclusions were drawn:

(1)

The melting time of a circular particle in forced thermal convection can be approximated as $t_m=0.4589R{e}^{0.5976}P{r}^{0.6561}S{t}^{-0.9332}$.

(2)

The melting time of a circular particle in mixed thermal convection first increases and then decreases as Gr increases, as a result of the competition between the forced and natural convections.

(3)

The effects of the particle interactions on the melting time are complicated due to the natural convection between two particles. The sufficiently strong natural convection can even render the downstream particle melt faster than the single particle.

(4)

For the same particle area, the elliptic particle with the aspect ratio being around 1.4 melts most slowly, as its longer axis is aligned with the cross-stream direction.

The contributions of the present study are that we not only developed a new numerical method for simulations of the melting of solid particles in the presence of thermal convection but also provide new results that significantly enhance the understanding of the melting rate of particles in thermal convection. Our method may be applied to the problems on the melting of ice in the ocean, particle melting in plasma spraying, and the unwanted melting of catalyst particles in a fluidized bed.