Finite volume method for multiphase flows with radiation and phase change

International Journal of Thermal Sciences 149 (2020) 106201

Contents lists available at ScienceDirect

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

Finite volume method for multiphase flows with radiation and phase change A. Hasečić a , S. Muzaferija b ,∗, I. Demirdžić a a b

Mašinski Fakultet, Vilsonovo šetalište 9, 71000 Sarajevo, Bosnia and Herzegovina Siemens Digital Industries Software, Nordostpark 3-5, 90411 Nürnberg, Germany

ARTICLE Keywords: Finite volume method Multiphase flow VOF Radiative heat transfer Phase change

INFO

ABSTRACT A mathematical model which can describe flow and phase change of a number of phases at high temperatures is presented. It combines an interface capturing multiphase model, the P–1 radiation model, and a melting/solidification model. The resulting equations are solved by employing the finite volume discretization, a segregated solution procedure and the SIMPLE algorithm. The melting/solidification model is a finite rate model which in the limiting case behaves like a thermodynamic equilibrium model and it can also be used in situations where the phase change occurs within a range of temperatures as well as for problems where the phase change occurs at a constant temperature. The method is verified on a number of problems. The results obtained show a good agreement with exact solutions or results which can be found in literature.

1. Introduction Multiphase flows at high temperatures, where the radiative heat transfer and phase change play an important role and cannot be neglected, are often found in nature and engineering practice. Typical examples can be found in different devices and manufacturing processes such as casting, laser material processing, thermal energy storage, additive manufacturing or nuclear engineering. Numerical simulations of such flows require a mathematical model which can describe the motion of multiphase systems (typically three-phase systems made of liquid, gas, and solid phases), energy transfer which includes all modes of heat transfer including radiation, and the phase-change processes (typically melting and solidification). In the current research we focus on the problems characterized by large scale interfaces between phases which can efficiently be resolved with numerical meshes and where the interactions between phases can be directly modeled. A number of techniques which belong to the family of interfacecapturing methods are available for this type of simulations [1–6]. A discretization which is based on the finite volume method can easily be made conservative if the transport of phases also represents the mass conservation of phases. The methods described in [3,4] meet this requirement and in this respect have an advantage over the level set method [5,6]. The method described in [4] is successfully used in [7,8] to describe multiphase flows with phase change. It is computationally very efficient, can deal with arbitrarily shaped interfaces, and can be extended to systems with an arbitrary number of phases and phase changes [8].

The work found in the literature, which studies the radiative heat transfer in multiphase systems, mainly investigates melting and solidification processes and radiative heat transfer in two-phase solid–liquid semitransparent systems where the interfaces between liquid and gas or solid and gas phases are considered as fixed boundaries and the modeling focuses on the development of the interface between the solid and the liquid phase as a consequence of melting and solidification. The radiative heat transfer through multiple layers in contact was first analyzed by Gardon [9] and Tarshis et al. [10]. They analyzed simple one-dimensional, steady state diffusion between two partially transparent media in contact with fixed and known in advance interface. They showed that optical properties of phases influence the temperature profiles, and that radiative heat transfer in general cannot be neglected. Spherical harmonics methods for calculation of radiative heat transfer was used by Matshushima and Viskanta [11] who used P–1 radiation model in their analysis of the effect of internal radiative transfer on natural convection and heat transfer in the crystal growth, and by Larsen et al. [12] who used simplified P–N radiation model in their investigation of radiative heat transfer in glass where they considered a discontinuity in the opacity, corresponding to two different materials in the left and right halves of the slab in the case of radiative cooling. In all the above researches, the interface between phases was known in advance. Habib [13,14] analytically solved the problem where the interface between phases was not known in advance. He solved a phase-change problem of a semi-infinite semitransparent material by conduction

∗ Corresponding author. E-mail address: [email protected] (S. Muzaferija).

https://doi.org/10.1016/j.ijthermalsci.2019.106201 Received 26 May 2019; Received in revised form 20 November 2019; Accepted 21 November 2019 Available online 5 December 2019 1290-0729/© 2019 Elsevier Masson SAS. All rights reserved.

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Superscripts

Nomenclature 𝑎 𝑎𝜙 , 𝑏𝜙 𝐴𝜙 , 𝐛𝜙 𝑐𝑣 𝐶𝑔 𝐶, 𝐷, 𝐸𝑉 , 𝐄𝑆 𝐝 𝐟𝑏 𝐢𝑘 𝐈 𝐺 𝑘 𝐿 𝑀 𝑚̇ 𝑛 𝑛𝑝ℎ 𝑝 𝐪ℎ 𝐪𝐺 𝐬 𝑆 𝑠𝛼𝑖 𝑡 𝑇 𝐓 𝐯 𝐯𝑆 𝑣𝑘 𝑉 𝑉̇ 𝑥, 𝑦 𝛼𝑖 𝛽𝑝 𝛾𝜙 𝜇 𝜌 𝜎 𝜎𝑠 𝜙

absorption coefficient coefficients of discretized equations coefficient matrix and source vector specific heat linear-anisotropic phase function coefficient coefficients in generic transport equation distance vector body force Cartesian basis vector unit tensor incident radiation thermal conductivity latent heat of fusion resistance coefficient mass flux refractive index number of phases pressure heat flux vector radiative energy flux surface vector surface source of phase 𝑖 time temperature Cauchy stress tensor velocity vector surface velocity vector Cartesian velocity component volume volume flux Cartesian coordinates phase volume fraction under-relaxation factor blending factor viscosity density Stefan–Boltzmann constant scattering coefficient generic variable

𝑚 𝑇 ′

interface between phases. Their results showed that radiation affects the dynamics of the phase-change process. Later Chan et al. [16] proposed a generalized phase-change model for melting and solidification with internal radiation which account for the existence of a two-phase zone in which partial phase-change can occur. Since the research of Chan et al. [16], two modeling approaches for simulating radiative heat transfer in flows with melting and solidification can be found in the literature. The first approach assumes a sharp interface between phases, and the second approach uses the mushy-zone model proposed by Chan et al. [16]. Leading researches of radiative heat transfer in multiphase flows with phase change which are based on a moving mesh method were performed by Brandon et al. [17–19] and Tsukada et al. [20]. When phase-change processes are solved using the moving-mesh method, the P–1 radiation model is most often used for the radiative heat transfer calculation in a combination with the finite-volume method [21–23]. On the other hand, in researches where the enthalpy method is used for melting and solidification and where the existence of mushy zone is admitted, the discrete ordinates [24,25] or the finite volume radiation models [26,27] are most commonly used for the radiative heat transfer calculations in a combination with the finite-volume discretization. Łapka and Furmański [28,29] proposed a new method for solidification processes based on the fixed grid front tracking method combined with the immersed boundary technique. In all the above researches where phase-change processes (solidification/melting) were encountered, it was assumed that phase change occurs at a constant temperature and the equilibrium solidification/melting model was used. Yao et al. [30,31] developed a nonequilibrium planar solidification model for emitting, absorbing, and isotropically scattering medium. A planar interface non-equilibrium solidification is assumed with crystalline phase nucleated on the surface at a given nucleation temperature, which may be significantly lower than the equilibrium melting temperature of the material. Yao et al. [32] then analyzed planar interface stability in non-equilibrium solidification of semitransparent materials. Their results suggested that, although the internal radiation leads to an undercooled melt in front of the interface, a planar interface can still be stable if there is a strong external heat transfer. The most recent research was done by Yang et al. [33] where they analyzed coupled conductive and radiative heat transfer during solidification and carried out experiments in order to obtain the boundary conditions and the radiative characteristics of the solid slags. For the phase-change processes the non-equilibrium model based on crystallization kinetics was used.

Subscripts 𝐵 𝑓 𝑖 𝑖𝑛 𝑗 𝑙, 𝑠 𝑃,𝑁 𝑤 𝜙

time-step counter transpose correction

boundary cell-face centers phase index initial phase change index liquid or solid phase cell centers wall solution variable

The modeling of radiative heat transfer in multiphase flows with phase change is mainly focused on the interface between solid and liquid phases, while the treatment of interfaces between liquid and gas or solid and gas phases is significantly simplified. Hasečić et al. [34] presented a mathematical model and a finite volume method which can be used to simulate the transport of mass, momentum and energy in multiphase flows at high temperatures without phase change. The interface-capturing method [4] was used to describe the motion of phases. In order to simplify the implementation of the radiation model, they focused on optically thick media where the radiative heat transfer could be described well with the P–1 model [35]. In this paper the method presented in [34,36] is extended to include melting/solidification processes.

and radiation. Abrams and Viskanta [15] studied the Stefan problem including combined radiative and conductive heat transfer (transient melting and solidification in translucent medium). They analyzed onedimensional energy transfer in a region of finite thickness with a planar 2

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

2. Mathematical model

A single set of governing equations (1)–(5) describes the transport process of a multiphase system, where each phase is associated with a set of physical properties (𝜌𝑖 , 𝜇𝑖 , 𝑘𝑖 , 𝑐𝑣𝑖 , 𝑎𝑖 , 𝑛𝑖 ). The local values of physical properties depend on the spatial distribution of phases represented by the volume fractions of phases 𝛼𝑖 : ∑ ∑ ∑ 𝜌= 𝜌𝑖 𝛼 𝑖 , 𝜇 = 𝜇𝑖 𝛼𝑖 , 𝑘 = 𝑘𝑖 𝛼𝑖 ,

In order to account for the phase change, a non-equilibrium rate driven melting/solidification model is incorporated into the mathematical model for solving the radiative transport in multiphase flows with free surfaces described in details in Ref. [34].

𝑖

2.1. Governing equations 𝑐𝑣 =

The governing equations describe the behavior of the mixture of phases. The momentum transport in the solid phases is simplified and the solid phase is modeled like a fluid immersed in the porous medium whose resistance is a function of the volume fraction of solid phases. The resulting equations describing mass, mass of phases, momentum, thermal energy, and incident radiation conservation in the volume 𝑉 , bounded by the surface 𝑆 with the outwards pointing surface vector 𝐬 are given below: d 𝜌 𝑑𝑉 + 𝜌 (𝐯 − 𝐯𝑆 ) ⋅ d𝐬 = 0 ∫𝑆 d𝑡 ∫𝑉 ) ( d 𝑠 d𝑉 𝜌 𝛼 d𝑉 + 𝜌 𝛼 𝐯 − 𝐯𝑆 ⋅ d𝐬 = ∫𝑉 𝛼𝑖 ∫𝑆 𝑖 𝑖 d𝑡 ∫𝑉 𝑖 𝑖 d 𝜌 𝐯 d𝑉 + 𝜌 𝐯(𝐯 − 𝐯𝑆 ) ⋅ d𝐬 = 𝐓 ⋅ d𝐬 + 𝐟 d𝑉 + 𝐟 d𝑉 ∫𝑆 ∫𝑆 ∫𝑉 𝑏 ∫𝑉 𝑠𝑟 d𝑡 ∫𝑉

𝑖

−

0=−

∫𝑆

𝐪𝐺 ⋅ d𝐬 +

∫𝑉 ∫𝑉

∫𝑉

(4𝑎𝑛2 𝜎𝑇 4 − 𝑎𝐺) d𝑉 +

(4𝑎𝑛2 𝜎𝑇 4 − 𝑎𝐺) d𝑉

∫𝑉

(1) (2) (3)

ℎ𝑚𝑠 d𝑉 = 0

(5)

where 𝑡 is the time, 𝜌 is the density, 𝐯 is the velocity vector, 𝐯𝑆 is the velocity of the surface 𝑆, 𝜌𝑖 and 𝛼𝑖 are the density and the volume fraction of the phase 𝑖, 𝑠𝛼𝑖 is the phase source or sink (e.g. when melting or solidification occurs), 𝐓 is the Cauchy stress tensor, 𝐟𝑏 is the body force, 𝐟𝑠𝑟 is the resistance force added when solid phases are present, 𝑐𝑣 is the specific heat, 𝑇 is the temperature, 𝐪ℎ is the heat flux vector, 𝑎 is the absorption coefficient, 𝑛 is the refractive index, 𝜎 is the Stefan– Boltzmann constant, 𝐺 is the incident radiation, ℎ𝑚𝑠 is the heat required for melting or solidification and 𝐪𝐺 is the radiative energy flux. If the problem to be solved involves volumes with moving boundaries (𝐯𝑆 ≠ 0), then, in addition to the above equations, the so-called space conservation law: d d𝑉 − 𝐯 ⋅ d𝐬 = 0 ∫𝑆 𝑆 d𝑡 ∫𝑉

(6)

which links the rate of change of volume 𝑉 and surface velocities 𝐯𝑆 of a deforming control volume, has to be satisfied [37]. System of Eqs. (1)–(5) is closed by employing the Stokes’ law, the Fourier’s law, and the constitutive relations for the P–1 radiation model: [ ] 𝐓 = 𝜇 ∇𝐯 + (∇𝐯)𝑇 − 𝑝 𝐈 (7)

𝐪ℎ = −𝑘 ∇𝑇

𝐪𝐺 = −

1 ∇𝐺 3(𝑎 + 𝜎𝑠𝑐 ) − 𝐶𝑔 𝜎𝑠𝑐

𝜌𝑖 , 𝑎= 𝜌

𝑖

√ 𝑎𝑖 𝛼 𝑖 , 𝑛 =

𝑖

∑

2 𝑖 𝑎𝑖 𝑛𝑖 𝛼𝑖

𝑎

(10)

There are two main numerical methods for solving phase-change problems, fixed and variable grid methods [38,39]. The most prominent fixed-grid method is the enthalpy porosity method which uses the field of enthalpy or temperature to estimate the volume fractions of liquid phases. A porosity model, where the local porosity depends on the volume fraction of liquid phases (melt fraction), is used to model the local porous resistance due to the presence of solid phases. The alternative approach uses moving grids to track interface between the liquid and solid phases. The methodology proposed in this work is somewhat similar to the enthalpy porosity method. The main difference is that the volume fraction of liquid, solid and gaseous phases are obtained solving their respective transport equations, and the phase-change rates appear as additional source terms in the phase transport equations. This approach can also be used to simulate phase change in liquid–gas systems, as it was done in [40] for example, and its extension to melting and solidification could allow us to simulate complex phase-change processes (e.g. modeling of evaporation and melting in a single simulation) using a unified modeling strategy. There are two possible modeling approaches when estimating phase-change rates: equilibrium models and non-equilibrium models, depending on whether solid overheating or melt undercooling is present or not. An equilibrium model assumes the local thermodynamic equilibrium condition at the interface. Such a model is presented for example in [7,41]. Depending on the treatment of the phase-change interface in the solidification and melting of semitransparent materials, planar interface model or mushy zone model can be used. Habib [13,14] and Abrams and Viskanta [15] who, among the first, solved the phasechange problem in semitransparent materials with the radiative heat transfer, used the planar interface model. The so-called mushy zone model was first proposed by Chan et al. [16] which argues that there should be a two-phase zone between the pure solid and liquid phase during the phase change of semitransparent materials. The formation of a mushy zone is derived from the thermodynamic equilibrium condition. Both of these models assume a local equilibrium condition at the solidification front and that assumption is valid only for the large solution domains with relatively small heat-transfer rate [31]. In advanced manufacturing processes, such as laser material processing, crystal growth or additive manufacturing, the treated material has a very small dimensions and high heat-transfer rate [42]. According to Yao [30], the combination of a small dimension and a fast rate of heat transfer leads to an extremely high velocity of the phase change interface and therefore a departure from the conventional local equilibrium condition. Because of that, the phase change will not occur at the equilibrium melting/solidification temperature and a large melt undercooling or solid overheating may exist at the interface. Also, local undercooling or overheating can occur due to an excessive internal radiative heat transfer in semitransparent materials experiencing solidification or melting [15] which cannot be described by an equilibrium model. In a non-equilibrium method phase-change kinetics is included which provides the relationship between the phase-change rate and the extent of local undercooling or overheating. Non-equilibrium phase

(4)

𝐓 ∶ ∇𝐯 d𝑉

𝑐𝑣𝑖 𝛼𝑖

𝑖

∑

2.2. Melting/solidification model

d 𝜌 𝑐 𝑇 d𝑉 + 𝜌 𝑐 𝑇 (𝐯 − 𝐯𝑠 ) ⋅ d𝐬 = − 𝐪 ⋅ d𝐬 ∫𝑆 𝑣 ∫𝑆 ℎ d𝑡 ∫𝑉 𝑣 +

∑

(8)

(9)

where 𝜇 is the viscosity, 𝑝 is the pressure, 𝐈 is the unit tensor, 𝑘 is the thermal conductivity, 𝜎𝑠𝑐 is the scattering coefficient, and 𝐶𝑔 is the linear-anisotropic phase function coefficient. 3

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

change processes have been researched by Grigoropoulos [43], Wang and Prasad [42], and recently Yao et al. [30–32] who were the first to propose the non-equilibrium model for semitransparent materials. The melting/solidification model used here is conceptually very similar to the homogeneous relaxation model [44] developed for cavitation and flash boiling. This modeling approach is attractive even if we intend to use it in situations where the equilibrium model is appropriate. The distribution of phases in the equilibrium model is completely determined by the energy transport and the distribution of the enthalpy. This complicates the solution of the phase transport equations for the liquid and solid phases. In the current approach, the phase transport equations for solid and liquid phases get additional sources and the phase transport equations are treated and solved in the way described in Ref. [34].

idea of the porous-resistance approach is to treat the whole solution domain as a porous medium, where porous resistance force is zero in cells where only liquid phase is present, and is non-zero where some solid phases are present. The additional resistance force is added to the momentum equation, and it stops the relative movement between solid and liquid phases. In this way the solid phase remains at rest or moves like a rigid body together with the neighboring wall boundaries to which it is attached. The resistance force 𝐟𝑠𝑟 is modeled as: √ √𝑚 √∑ 4 𝛼𝑗𝑠 (17) 𝐟𝑠𝑟 = −𝑀(𝐯 − 𝐯𝑠𝑡 ) √ 𝑗=1

where 𝑀 is a constant large enough to make the resistance force the most dominant term in the momentum transport where some solid phase is present, 𝐯 is the mixture velocity and 𝐯𝑠𝑡 is the target solid velocity, which is imposed by the boundaries to which the solid is attached. In order to increase the resistance force where the volume fraction of the solid phase is relatively small, the force is multiplied by the fourth root of the solid-phase volume fraction. In this way the force is zero where the volume fraction of solid phases is zero and grows quickly as soon as the volume fraction of solid phases starts increasing.

Contribution to phase transport equation. Every phase change 𝑗 is associated with a liquid phase 𝑖 = 𝑗𝑙 and a solid phase 𝑖 = 𝑗𝑠 , where 𝑗 = 1, 2, … , 𝑚 and 𝑚 is the total number of melting/solidification processes ongoing in the simulation. The phase sources 𝑠𝛼𝑖 , which are the result of a phase-change 𝑗, are approximated in the present work using the following expression: 𝑠𝛼𝑗 = 𝜌𝑗𝑙

𝛼 𝑗𝑙 − 𝛼𝑗∗

𝑙

𝛩

𝑙

𝑠𝛼𝑗 = −𝑠𝛼𝑗 𝑠

𝑓 (𝛼𝑗𝑙𝑠 )

(11)

Contribution to energy equation. The contribution to the energy equation due to the melting and solidification is defined as:

(12)

𝑙

ℎ𝑚𝑠 =

𝑗=1

where 𝛼 𝑗𝑙 and 𝛼𝑗∗ are the equilibrium and relative volume fraction of 𝑙 the liquid phase 𝑗𝑙 and 𝛩 is the relaxation time. Negative values of 𝑠𝛼𝑗 indicate solidification. The equilibrium volume fraction is typically 𝑙 a function of pressure, temperature, and composition of the material. For alloys, where the phase change takes place between the solidus 𝑇𝑠 and the liquidus temperature 𝑇𝑙 , 𝛼 𝑗𝑙 is often approximated as a linear function of temperature: 𝛼 𝑗𝑙 = max(0, min(1,

𝑇 − 𝑇𝑙 )) 𝑇𝑠 − 𝑇𝑙

where 𝛼𝑗𝑙𝑠 ⎧ ⎪𝛼𝑗 𝛼𝑗𝑙𝑠 = ⎨ 𝑙 ⎪𝛼𝑗𝑠 ⎩

2 is defined as if

𝑠𝛼𝑗 > 0

if

𝑠𝛼𝑗 < 0

𝑙

𝑠

3. Discretization and solution method In order to solve equations of the mathematical model the finitevolume method described in [45] is employed and will be, for completeness, briefly described here.

(13) 3.1. Discretization The discretization process consists of time, space and equation discretization. The time of interest is divided into a number of time steps 𝛿𝑡. The space is discretized by a finite number of contiguous non-overlapping polyhedral control volumes (CV) or cells with the computational points in their centers (see Fig. 1). The governing and constitutive equations form a closed set of 4+𝑛𝑝ℎ equations, where 𝑛𝑝ℎ is the number of phases present in a simulation. Instead of discretizing each of the governing equations, it is handy to discretize a generic transport equation d 𝜌𝐶𝜙 𝜙 d𝑉 + 𝜌𝐶 𝜙(𝐯 − 𝐯𝑆 ) ⋅ d𝐬 ∫𝑆 𝜙 d𝑡 ∫𝑉

Solidification can occur only in cells in which some liquid phase is present, and melting is possible only in cells which have some solid phase. Because of this, the melting rate 𝑠𝛼𝑗 should be multiplied with 𝑙 a function whose value is zero if the melting rate is positive and the volume fraction 𝛼𝑗𝑠 is zero or if the melting rate is negative and the volume fraction 𝛼𝑗𝑙 is zero and otherwise its value should be one. The use of the step function could cause some instabilities in the calculation, and in the present work the function 𝑓 is approximated as: 1 + tanh(30𝛼𝑗𝑙𝑠 − 3)

(18)

𝑠𝛼𝑗 𝐿𝑗

where 𝐿𝑗 is the latent heat of fusion of the phase-change 𝑗.

that could be a quite crude approximation in a general case. The same formula can also be used for pure materials where the phase change takes place at a constant temperature, only in this case it is assumed that the phase change takes place in a temperature interval which is much smaller than the global temperature scale of the problem. The relative volume fraction of the liquid phase represents the volume occupied by the liquid phase 𝑗𝑙 relative to the volume occupied by the liquid and the solid phases, i.e: 𝛼𝑗𝑙 𝛼𝑗∗ = (14) 𝑙 𝛼𝑗𝑙 + 𝛼𝑗𝑠

𝑓 (𝛼𝑗𝑙𝑠 ) =

𝑚 ∑

=

∫𝑆

𝐷𝜙 ∇𝜙 ⋅ d𝐬 +

∫𝑉

𝐸𝜙𝑉 d𝑉 +

∫𝑆

(19)

𝐄𝜙𝑆 ⋅ d𝐬

which represents the transport of phases, momentum, energy, and incident radiation. The continuity equation (1) is often combined with momentum equation (3) resulting in the pressure correction equation. 𝜙 in the generic transport equation (19) stands for the phase volume fraction 𝛼𝑖 (𝑖 = 1, … , 𝑛𝑝ℎ ), the velocity component 𝑣𝑘 (𝑘 = 1, 2, 3), the temperature 𝑇 , and the incident radiation 𝐺. The meaning of the terms 𝐶𝜙 , 𝐷𝜙 , 𝐸𝜙𝑉 , and 𝐄𝜙𝑆 are given in Table 1. Integrating generic transport equation (19) over a control volume results in the following (still exact) equation: ∑ ∑ d 𝜌 𝐶𝜙 𝜙(𝐯 − 𝐯𝑆 ) ⋅ d𝐬 = 𝐷 ∇𝜙 ⋅ d𝐬 + 𝜌𝐶𝜙 𝜙 d𝑉 + ∫𝑆𝑓 ∫𝑆𝑓 𝜙 dt ∫𝑉

(15)

(16)

𝑙

Rate of change

Contribution to momentum equation. The solid phase is modeled as a liquid phase immersed in a porous medium by introducing the additional porous resistance forces 𝐟𝑠𝑟 in the momentum equation. The basic

𝑓

∫𝑉 4

𝑓

Convection

𝐸𝜙𝑉 d𝑉 +

∑

∫𝑆𝑓 𝑓 Sources

Dif fusion

𝐄𝜙𝑆 ⋅ d𝐬

(20)

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Fig. 1. Control volume and geometric parameters. Table 1 The meaning of various terms in generic transport equation (19). 𝜙

𝑣𝑘 𝑇

𝐶𝜙 𝜌𝑖 𝜌 1 𝑐𝑣

𝐺

0

𝛼𝑖

𝐷𝜙

𝐸𝜙𝑉

𝐄𝜙𝑆

0

𝑠𝛼𝑖

0

𝜇 𝑘

𝑓 𝑏𝑘 𝐓 ∶ ∇𝐯 − (4𝑎𝑛2 𝜎𝑇 4 − 𝑎𝐺) + ℎ𝑚𝑠

(𝜇 ∇𝐯 − 𝑝 𝐈) ⋅ 𝐢𝑘 0

1 3(𝑎 + 𝜎𝑠 ) − 𝐶𝑔 𝜎𝑠

4𝑎𝑛2 𝜎𝑇 4 − 𝑎𝐺

0

∫𝑆𝑓

𝐷𝜙𝑃

𝜌(𝐯 − 𝐯𝑆 ) ⋅ d𝐬 ≈ 𝜌𝑓 (𝐯∗𝑓 ⋅ 𝐬𝑓 − 𝑉̇ 𝑓 )

∫𝑆𝑓

𝐯𝑆 ⋅ d𝐬 =

𝛿𝑉𝑓 𝛿𝑡

𝐝𝑓 ⋅ 𝐬𝑓

+ 𝐷𝜙𝑁

(27)

𝐝𝑃 ⋅ 𝐬𝑓 𝐝𝑓 ⋅ 𝐬𝑓

Source. The volume sources are approximated as: ∫𝑉

(28)

𝐸𝜙𝑉 d𝑉 ≈ (𝐸𝜙𝑉 )𝑃 𝑉𝑃

and the surface integrals involving vector 𝐄𝜙𝑆 are calculated as: ∫𝑆𝑓

(29)

𝐄𝜙𝑆 ⋅ d𝐬 ≈ (𝐄𝜙𝑆 )𝑓 ⋅ 𝐬𝑓

Resulting algebraic equations. After discretizing all the terms of Eq. (20), a non-linear algebraic equation is obtained for each control volume linking the value of 𝜙 at the CV center with 𝜙 values at the centers of the nearest neighbors: ∑ 𝑎𝜙 𝑃 𝜙 𝑃 − 𝑎𝜙𝑁 𝜙𝑁 = 𝑏𝜙 (30)

(22)

𝑁

where the summation is over all the neighboring CVs. The coefficients 𝑎𝜙𝑃 and 𝑎𝜙𝑁 and the source term 𝑏𝜙 result from the discretization and are summarized below: 𝐬𝑓 ⋅ 𝐬𝑓 ( ) − 𝐶𝜙𝑓 min 𝑚̇ 𝑓 , 0 𝑎𝜙𝑁 = 𝐷𝜙𝑓 𝐝𝑓 ⋅ 𝐬𝑓

(23)

𝑎𝜙𝑃 =

where 𝐯∗𝑓 is obtained from a special interpolation practice [46] which enables calculation of the pressure field (SIMPLE algorithm [47]) and ensures a strong coupling of velocity and pressure which leads to a stable solution procedure, as described later, and 𝑉̇ 𝑓 =

𝐝𝑁 ⋅ 𝐬𝑓

(21)

The cell-face values 𝜙𝑓 are calculated blending the first-order upwind (UD) and higher-order (HO) approximations (second-order central differencing, except for convective transport of phases where the high resolution interface capturing (HRIC) scheme [4] is used). The UD contribution is treated implicitly, while the difference between UD and HO approximations multiplied by a blending factor 𝛾𝜙 is treated explicitly. The mass flux is calculated as ∫𝑆𝑓

𝐷𝜙𝑃 𝐷𝜙𝑁

𝐷𝜙𝑓 =

Convection. The convection term is approximated as

𝑚̇ 𝑓 =

𝐬𝑓 ⋅ 𝐬𝑓 𝐝𝑓 ⋅ 𝐬𝑓

or by a harmonic interpolation in case of discontinuity of diffusion coefficient at contact surfaces between phases:

where the superscript 𝑚 is the time-step counter and the unsuperscripted terms refer to the current time 𝑡𝑚 .

𝜌 𝐶𝜙 𝜙 (𝐯 − 𝐯𝑆 ) ⋅ d𝐬 ≈ 𝑚̇ 𝑓 𝐶𝜙𝑓 𝜙𝑓

(25)

where 𝐷𝜙𝑓 is obtained by a linear interpolation between values calculated at nodes 𝑃 and 𝑁: 𝐝𝑁 ⋅ 𝐬𝑓 𝐝𝑃 ⋅ 𝐬𝑓 𝐷𝜙𝑓 = 𝐷𝜙𝑃 + 𝐷𝜙𝑁 (26) 𝐝𝑓 ⋅ 𝐬𝑓 𝐝𝑓 ⋅ 𝐬𝑓

Rate of change. The Euler implicit scheme is employed to approximate the rate of change term:

∫𝑆𝑓

𝐬𝑓 ⋅ 𝐬𝑓 𝐝𝑓 ⋅ 𝐬𝑓

+ 𝐷𝜙𝑓 (∇𝜙𝑁 ⋅ 𝛥𝐝𝑁 − ∇𝜙𝑃 ⋅ 𝛥𝐝𝑃 )

where the summation on 𝑓 is over all faces enclosing the CV. In what follows, the characteristic terms in Eq. (20) are approximated and a discrete form of the equation is obtained.

(𝜌𝐶𝜙 𝜙𝑉 )𝑃 − (𝜌𝐶𝜙 𝜙𝑉 )𝑚−1 d 𝑃 𝐶𝜙 𝜙 d𝑉 ≈ ∫ dt 𝑉 𝛿𝑡

𝐷𝜙 ∇𝜙 ⋅ d𝐬 ≈ 𝐷𝜙𝑓 (𝜙𝑁 − 𝜙𝑃 )

∑ 𝑁

𝑏𝜙 =

∑

𝑎𝜙𝑁 +

(𝜌𝐶𝜙 𝑉 )𝑚−1 𝑃

𝛿𝑡 ) ( UD 𝛾𝜙 𝐶𝜙𝑓 𝑚̇ 𝑓 𝜙HO 𝑓 − 𝜙𝑓

(31)

𝑓

+ (24)

∑

𝐷𝜙𝑓 (∇𝜙𝑁 ⋅ 𝛥𝐝𝑁 − ∇𝜙𝑃 ⋅ 𝛥𝐝𝑃 )

𝑓

+ (𝐸𝜙𝑉 )𝑃 𝑉𝑃 +

where 𝛿𝑉𝑓 is the volume swept by the cell-face 𝑓 during the time interval 𝛿𝑡.

𝐬𝑓 ⋅ 𝐬𝑓 𝐝𝑓 ⋅ 𝐬𝑓

∑ ∑ (𝜌𝐶𝜙 𝜙𝑉 )𝑚−1 𝑃 (𝐄𝜙𝑆 )𝑓 ⋅ 𝐬𝑓 + + 𝑎𝜙𝐵 𝜙𝐵 𝛿𝑡 𝑓 𝐵

where the summation on 𝐵 is over the boundary faces surrounding cell 𝑃.

Diffusion. The diffusion term is approximated as [48] 5

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Initial and boundary conditions. To start a transient analysis, all dependent variables have to be specified at the initial time 𝑡0 . The expressions (22), (25) and (29) for convection and diffusion fluxes and surface sources remain also valid for the boundary faces. In the case of Dirichlet boundary conditions, the value 𝜙𝑁 or 𝜙𝑓 are replaced with the known boundary values 𝜙𝐵 . In the case of Neumann boundary conditions, the boundary fluxes are known and as such are directly added to the source terms, while the boundary values 𝜙𝐵 are obtained by extrapolation from the interior of the solution domain. The boundary condition at wall boundaries for the P-1 radiation model has the following form [35]: 𝜖𝑤 (𝐺 − 4𝑛2 𝜎𝑇𝑤 4 ) |𝐬𝐵 | (32) 𝐪𝐺,𝑤 ⋅ 𝐬𝐵 = 2 (2 − 𝜖𝑤 ) 𝑤

Fig. 2. Solidification in a half-space.

Table 2 Values of the constant 𝐶(𝑆𝑡) for three 𝑆𝑡 numbers.

where 𝐬𝐵 is the boundary-face surface vector and 𝜖𝑤 and 𝐺𝑤 are the wall emissivity and incident radiation, respectively. 3.2. Solution algorithm

1.0

10.0

0.387

0.755

0.926

(33) 4.1. Solidification in a half-space, Stefan problem

for each solution variable. 𝐴𝜙 is a sparse, diagonally dominant matrix and vector 𝜱 contains values of dependent variable at CV centers. Eq. (33) is solved in turn for each solution variable using the conjugate gradient method with an incomplete Cholesky preconditioning. This procedure is repeated by updating the coefficient matrix and source terms until a converged solution is obtained, i.e. until the residual norm: ‖𝐫𝜙 ‖ = ‖𝐴𝜙 𝜱 − 𝐛𝜙 ‖

0.1

𝐶 (𝑆𝑡)

error irrespective of what dominates it, time or space discretization. In order to make this process more effective and avoid a situation where the discretization error in time dominates the overall discretization error, the time-step size used on the coarsest mesh should produce a result which is time-step size independent. The free-surface capturing and the HRIC discretization scheme require time-step sizes which respect the CFL number constrain (i.e. CFL < 1, so the free surface does not cross more than a cell in each time step) and this constrain often results in time-step sizes which already deliver a solution independent of the time-step size. In addition, the segregated solution algorithm and the highly nonlinear nature of radiative heat transfer and phase change could also require rather small time-step sizes in order to have stable and robust runs.

After discretizing the governing equations for all CVs and for all solution variables, a coupled system of non-linear algebraic equations is obtained. Those equations are solved by employing the following segregated iterative procedure. All dependent variables are given their initial values at the beginning of the first time step. The boundary conditions are applied and the sets of equations for each individual solution variable are linearized and temporarily decoupled by assuming that coefficients 𝑎𝜙𝑁 and source terms 𝑏𝜙 are known (calculated using variable values from the previous iteration) resulting in a system of linear algebraic equations of the form 𝐴𝜙 𝜱 = 𝐛𝜙

𝑆𝑡

Stefan problem is a standard test case for melting and solidification modeling. Computational domain is the half-space, 𝑥 > 0 (Fig. 2) in which the phase-change front moves through the solution domain as a function of the initial and boundary conditions and the physical properties of the material. The physical properties of solid and liquid phases are assumed to be the same (𝜌 = 1 kg∕m3 , 𝑘 = 1 W∕m K, and 𝑐𝑣 = 1 J∕kg K) as in the study by Teskeredžić et al. [41], who used an enthalpy melting/solidification model and a similar finite volume method. Initially, the whole domain is filled with the liquid at the temperature 𝑇𝑖𝑛 1 K larger than the temperature at which the solidification occurs. The wall on the left side of the solution domain is at temperature 𝑇𝑤 which is 1 K lower than the solidification temperature. As the liquid cools down and solidifies, the phase-change front travels in the positive 𝑥 direction. The exact position of the solidification front 𝑠 is a function of time and the Stefan number: 𝑐 𝛥𝑇 (35) 𝑆𝑡 = 𝑣 𝐿 where 𝛥𝑇 is the temperature difference between the solidification temperature 𝑇𝑠𝑜𝑙 and the cold wall temperature 𝑇𝑤 and 𝐿 is the latent heat. It is derived on the assumption that the interface is in the thermodynamic equilibrium [49]: √ 𝑠(𝑡) = 𝐶 (𝑆𝑡) 𝑡 (36)

(34)

is reduced by three to four orders of magnitude. In the next time step the whole procedure is repeated, except that the initial values are replaced by the values from the previous time step. Note that in order to promote stability of the solution method, an under-relaxation is often necessary. 4. Verification and application The verification work is focused on the implementation of the nonequilibrium phase-change model and its interaction with the previously developed multiphase model with radiation [34,36]. First, a Stefan problem, solidification in a square enclosure with radiative heat transfer, and solidification of a gas–liquid system inside a rotating tank are analyzed and then the method is applied to a blowing phase of a glass manufacturing process. The under-relaxation factors of 0.8 for momentum, 0.2 for pressure, 0.7 for energy and 1.0 for the transport of phases are used for all problems, expect for the solidification in a rotating tank, where the under-relaxation factors for the transport of momentum and phases are reduced to 0.7 and 0.9 respectively. All calculations are performed on a sequence of systematically refined meshes and time-step sizes in order to evaluate their influence on results. In the current study, the spatial and time coordinates are consistently refined (e.g. if the mesh spacing is halved in each direction, the time step is also halved). The procedure reduces the discretization

where 𝐶(𝑆𝑡) is the constant which dependents on the 𝑆𝑡 number and is given in Table 2. A number of tests are conducted with the same boundary and initial conditions and three different values of the latent heat (10, 1 and 0.1 J∕kg). The mesh and time-step size independent results are obtained on a uniform 1D numerical mesh made of 800 CVs using the time step of 0.01 s. All calculations are performed in double precision and the iteration error is reduced to the machine accuracy. Also, different 6

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Fig. 3. Solidification front position for the whole range of 𝛩 and Stefan numbers. Table 3 Properties of the solid and the liquid phase.

Solid Liquid

𝑎 (m−1 )

𝜎𝑠𝑐 (m−1 )

𝐶𝑔

𝑛

𝜌 (kg/m3 )

𝑘 (W/m K)

𝑐𝑣 (J/kg K)

3 3

0 0

0 0

1 1

1 1

7.56 4.536

2 2.4

values of the relaxation time are used (𝛩 = 100 , 10−2 , 10−4 , 10−6 , 10−8 , and 10−10 s) and the results are compared with the exact solution. Comparisons of the calculated and the analytical solutions are given in Fig. 3 for three different values of Stefan number and for the whole range of relaxation times 𝛩. As it can be seen from the figures, as 𝛩 gets lower the numerical results approach the analytical solution which assumes the thermo-dynamic equilibrium. This is an expected behavior since the non-equilibrium model should behave as the equilibrium model as 𝛩 → 0. The comparison shows a very good agreement with the exact solution for all three Stefan numbers.

Fig. 4. Solidification in a square enclosure. Geometry and boundary conditions.

4.2. Solidification in a square enclosure with radiative heat transfer

at temperature 𝑇1 = 500 K. Other three boundaries are kept at the initial temperature 𝑇2 = 1000 K. The phase change takes place in a temperature range, such that there is no solid above the liquidus temperature 𝑇𝑙 = 800 K and there is no liquid below the solidus temperature 𝑇𝑠 = 600 K. It is assumed that liquid volume-fraction varies linearly between 𝑇𝑠 and 𝑇𝑙 , i.e. 𝛼 𝑗𝑙 = (𝑇 − 𝑇𝑠 )∕(𝑇𝑙 − 𝑇𝑠 ). The relaxation time is set to 10−6 s.

The solidification of a material inside a two-dimensional black rectangular enclosure is considered. The geometry of the solution domain and boundary conditions are shown in Fig. 4, while the properties of phases are given in Table 3. The results are compared with the results obtained by Mishra et al. [26] who did the analysis of solidification of a two-dimensional semitransparent absorbing, emitting, and scattering medium. They solved the problem using the lattice Boltzmann method (LBM) and the finite volume method was used to compute the information about the radiative heat transfer required in the LBM formulation. At the beginning of the simulation, the material is at temperature 𝑇𝑖𝑛 = 1000 K. The upper wall of the rectangular enclosure is maintained

Conduction–radiation parameter 𝑁 = 𝑘𝑠 𝛽∕(4𝜎𝑇𝑖𝑛3 ) is 0.1, and the latent heat is 1000 J/kg. 𝛽 = 𝑎 + 𝜎𝑠𝑐 is the extinction coefficient, where 𝑎 is the absorption coefficient and 𝜎𝑠𝑐 is the scattering coefficient. The results presented here are obtained on a sequence of systematically refined uniform meshes and time-step sizes (17 × 17 CVs and 0.0002 s, 35 × 35 CVs and 0.0001 s, and 71 × 71 CVs and 0.00005 s). 7

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Fig. 6. Initial and boundary conditions of gas–liquid system in rotating tank.

Fig. 5. Variation of temperature and liquid volume fraction along the centerline at 𝑡 = 5 s.

Variation of the centerline normalized temperature 𝑇 ∕𝑇1 and liquid phase volume fraction at time 𝑡 = 5 s is shown in Fig. 5. The results are compared with results reported in [26]. It can be seen that solidification happens within the range of temperatures. Solidification begins at the 𝑇𝑙 = 800 K. Above the liquidus temperature the liquid volume fraction is unity. Solidification is completed at the 𝑇𝑠 = 600 K. Below that temperature the liquid volume fraction is zero. Between the solidus and liquids temperatures one can observe that both the solid and the liquid phase are present.

Fig. 7. Mass averaged system temperature as a function of mesh and time-step sizes.

Table 4 Material properties of gas, liquid and solid phases inside a rotating tank.

Gas Liquid Solid

𝑎 (m−1 )

𝜎𝑠𝑐 (m−1 )

𝐶𝑔

𝜌 (kg/m3 )

𝑘 (W/m K)

𝑐𝑣 (J/kg K)

10 200 5

0.1 0.1 0.1

0 0 0

1.188 4300 4300

0.0257 4 8

1007 800 800

4.3. Solidification of a gas–liquid system inside a rotating tank Hasečić et al. [34] analyzed behavior of a hot gas–liquid system which moves and cools down in a rotating tank. The study included the effects of radiation, but a solidification of the liquid phase was not analyzed. Adding the phase-change modeling to the test case results in some interesting results. A rectangular 0.2 m × 1 m tank is half filled with a liquid and half with a gas. Initial tank position, the axis and direction of rotation, and the initial and boundary conditions are shown in Fig. 6. Liquid phase solidifies at the temperature 1700 K and its latent heat is 𝐿 = 4.55 × 105 J/kg. The material properties of the three phases are given in Table 4. Since the walls of the tank rotate, the velocity of the solid phase 𝐯𝑠𝑡 from the Eq. (17) is equal to the rotational velocity of the tank: 𝐯𝑠𝑡 = 𝜔 × 𝐫

is then defined as: √ 𝐟𝑠𝑟 = −𝑀(𝐯 − 𝜔 × 𝐫) 4 𝛼𝑠

(38)

where 𝛼𝑠 is the volume fraction of the solid phase. The value of the resistance coefficient 𝑀 is 9 × 1012 . The total simulation time is 360 s, which corresponds to one complete tank revolution. A systematic mesh and time-step refinement, starting with 5x30 (150) CVs and 0.02 s and ending with 40x240 (9600) CVs and 0.0025 s, is used to evaluate the mesh and time-step size sensitivity of results. The system average temperature as a function of different mesh and timestep sizes is shown in Fig. 7. One cannot observe a gradual convergence of the results towards the grid and time-step independent solution. Fig. 8 shows the distribution of phases at 𝑡 = 45, 90, 180, 270, and 360 s after the beginning of rotation. One can see that at 𝑡 = 45 s the liquid phase is locally cooled to the solidification temperature at the

(37)

where 𝜔 is the angular speed of tank rotation and 𝐫 is the position vector relative to the center of rotation. The resistance force from Eq. (3) 8

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Fig. 8. Phases distribution at different times (liquid is colored yellow, gas is red and solid is blue).

cold walls. At 𝑡 = 90 s, the liquid phase occupies the bottom half of the tank and the newly formed solid phase remains attached to the tank wall. The continuous solidification which takes place along solid walls results in less and less of the liquid phase. At 𝑡 = 180 s the liquid phase occupies much less than a half of the tank. The liquid phase in contact with gas radiates the heat through the air and cools down faster than the liquid away from the interface and the walls. At some point in time, the liquid at the gas–liquid interface reaches the solidification temperature and solidifies. The consequence of this can be well observed at 𝑡 = 270 and 𝑡 = 360 s. The liquid phase is completely surrounded by the solid phase and it cannot freely move any more. At 𝑡 = 360 s the liquid phase occupies the top part of the tank, a result which probably does not match our initial expectations. An interesting question is how the results change if certain aspects of the original problem are neglected. The problem simplifications are necessary if one is not able to model all relevant physical processes. For example, one can neglect the rotation of the tank and the radiative heat transfer. We will refer to this simplest modeling as the reference case. The next problem description includes radiation but neglects the tank rotation. We will refer to this modeling as the case with radiation. We will also analyze the problem where the radiation is neglected but the tank rotation is modeled and we will refer to this case as the case with rotation. Fig. 9 shows the phases distributions at the end of calculation for different modeling approximation. It is not difficult to see that in this specific case neglecting some important aspects of the problem results in very large modeling errors. One can see that in the reference case the solid phase exists only in the lower half of the tank and in proximity of wall boundaries and the liquid phase occupies the bottom half of the tank. The liquid was in the contact with all boundaries in the case with rotation and consequently the solid phase also exist along all walls. The liquid phase occupies the bottom half of the tank at the end of simulation and the solid phase does not separate the liquid and gas phases. The solid and liquid phases also occupy the lower part of the tank in the case with radiation. The modeling of radiation is responsible that the interface with the gas phase cools down faster and as a consequence the solid phase is formed separating the liquid and gas phases. As one can see from Fig. 10, the most prominent qualitative

Fig. 9. Phases distributions at 𝑡 = 360 s for different modeling assumptions.

result, namely the final distribution of phases, does not changes with the systematic mesh refinement. The average system temperature and the percentage of the solid phase as a function of time and different modeling assumptions are presented in Fig. 11. The graphs show that the system cools down and consequently solidifies much faster if the radiation and rotation of the system are modeled. The influence of the radiative heat transfer is larger at the beginning when the system temperature is higher and the fluid motion is not so intensive. At this stage of simulation we see that the case with radiation cools down faster than the problem where radiation is neglected and rotation is modeled. As the liquid and solid 9

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Fig. 10. Final distribution of phases on three different meshes.

Fig. 11. Temperature and solid volume fraction as a function of time and different modeling assumptions. Table 5 Material properties of phases.

phases intensify their motion, the convective heat transfer starts making a significant contribution to the overall heat transfer. As expected, the problem where we model both rotation and radiation results in fastest cooling and solidification of the system.

𝑎 (m−1 ) 𝜎𝑠𝑐 (m−1 ) 𝐶𝑔 𝜌 (kg/m3 ) 𝑘 (W/m K) 𝑐𝑣 (J/kg K) 𝜇 (Pa s) Air 3 Liquid glass 10 Solid glass 200

0.1 0.1 0.1

0 0 0

1.188 2500 2500

0.0257 1.5 1.5

1007 1350 1350

1.824 × 10−5 1000 1000

4.4. Numerical modeling of a glass manufacturing process Dijkstra and Mattheij [50] have studied a blowing phase in the industrial manufacturing of glass containers. They used the boundary element method to analyze the glass blowing phase of several glass containers focusing on the glass shape and neglecting the heat transfer and phase change. The same problem is also analyzed in the present work, only the modeling now also includes the energy transport and solidification of glass. The shape of the mould and the initial and boundary conditions are adopted from Dijkstra and Mattheij [50] and are schematically shown in Fig. 12. Initial temperature of the air inside the mould is set to be equal to the mould temperature. The material properties of phases are given in Table 5. Since very often the mould is covered with a lubricating substance to improve the flow of glass, the walls of the mould are assumed to be slip. At the surface 𝑆 the glass is not allowed to move. The problem is axisymmetric, and the resulting 2D problem is solved using a number of meshes and time-step sizes (225 CVs and 0.0004 s, 900 CVs and 0.0002 s, 2025 CVs and 0.00015 s and 3600 CVs and 0.0001 s). The finest mesh is shown in Fig. 12. The mould is porous, and it is assumed that air can flow through it and glass cannot. The velocity normal to the mould is approximated as: 𝑣𝑛 =

𝛼𝑎𝑖𝑟 (𝑝 − 𝑝𝑎𝑚𝑏 ) 𝑅

and 𝑅 is a resistance coefficient which is chosen in a way that boundary conditions in this case match the boundary conditions described in [50] and its value is 10−4 . Fig. 13 shows the glass phase at several instances of time. The inlet pressure 𝑝𝑖𝑛 is first set equal to the ambient pressure. This lasts for 0.37 s and during this time the glass flows mainly downwards under the influence of gravity. One refers to this phase of the process as the sagging stage. Dijkstra and Mattheij [50] do not give any information about times when certain stages of the process start or end. Based on some suggestions found in literature [51] and observations from numerical experiments, it was decided to end the sagging stage at 𝑡 = 0.37 s. By that time, the glass sags towards and almost touches the bottom of the mould. As the sagging stage ends, the inlet pressure 𝑝𝑖𝑛 is increased to 1.38 bar and air starts to be blown into the mould. The glass blowing stage ends at 𝑡 = 0.44 s and at that time the parison follows the shape of the mould. It is observed that during the sagging stage the glass mainly moves in the vertical directions. During the blowing stage the glass moves in vertical and radial directions, and as a consequence of this, the glass is everywhere in contact with the mould at the end of the blowing stage.

(39)

The temperature fields at two different times and parison positions are shown on Fig. 14. At the 𝑡 = 0.37 s the parison is entering the blowing stage and the air at the ambient temperature starts to be blown

where 𝛼𝑎𝑖𝑟 is the volume fraction of air in the contact with the mould, 𝑝 is the local pressure, 𝑝𝑎𝑚𝑏 is the ambient pressure outside the mould 10

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

Fig. 12. Geometry, initial and boundary conditions (left) and numerical mesh (right)’.

Fig. 13. Shapes of glass at six instances of time.

Fig. 14. Temperature field at two instances of time.

into the mould through the inlet boundary. Since it is assumed that the top of the mould is opened and in a contact with the environment, the parison is at 𝑡 = 0.37 s already partially filled with the cold air. At 𝑡 = 0.42 s the air pocket between the parison and the mould, which was initially cold, is almost heated to the parison temperature.

air slowly enters the mould through the top inlet boundary, causing the mass averaged temperature to slowly decrease. The glass sagging stage ends at 𝑡 = 0.37 s when the cold air is blown into the mould. At this time, as one can see from Fig. 15, the average temperature starts dropping faster and this last until 𝑡 = 0.44 s when the blowing stage ends. From this time the cooling of the system slows down.

Fig. 15 shows the mass averaged temperature of the system as a function of time, mesh and time-step sizes. During the glass sagging stage, the parison sags towards the bottom of the mould, and the cold

The sagging and the blowing stages of the process are fast and last only 0.44 s. The liquid phase does not cool down much during this time 11

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

of phases could result in significant modeling errors. The developed model is also applied to a glass manufacturing process showing that the developed methodology is also fit for some real-life applications. The P-1 radiation model is relatively simple and its applicability is limited to optically thick media. It would be very interesting to formulate mathematical model and methodology based on more sophisticated radiation models, such as the discrete ordinates method (DOM) or the Monte Carlo method. This would expend the areas of applicability of the modeling and it would also allow one to model radiative phenomena taking place at the interfaces between phases (e.g. modeling of reflection and refraction). Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 15. Mass averaged temperature as a function of time, mesh and time-step size.

References [1] W.F. Noh, P. Woodward, SLIC (simple line interface calculation), Lecture Notes in Phys. 59 (1976) 330–340. [2] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [3] O. Ubbink, R.I. Issa, A method for capturing sharp fluid interfaces on arbitrary meshes, J. Comput. Phys. 153 (1999) 26–50. [4] S. Muzaferija, M. Perić, Computation of free surface flows using interface-tracking and interface-capturing methods, in: O. Mahrenholtz, M. Markiewicz (Eds.), Nonlinear Water Wave Interaction, Computational Mechanics Publications, Southampton, 1998, pp. 59–96 (Chapter 3). [5] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: Algorithms based on Hamilton–Jacobi formulation, J. Comput. Phys. 79 (1988) 12–49. [6] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1994) 146–159. [7] A. Teskeredžić, I. Demirdžić, S. Muzaferija, Numerical method for calculation of complete casting process – Part I: Theory, Numer. Heat Transfer B 68 (2015) 295–316. [8] A. Teskeredžić, I. Demirdžić, S. Muzaferija, Numerical method for calculation of complete casting process – Part II: Validation and application, Numer. Heat Transfer B 68 (2015) 317–335. [9] R. Gardon, Calculation of temperature distributions in glass plates undergoind heat treatment, J. Amer. Cer. Soc. 41 (1958) 200–209. [10] L.A. Tarshis, S. O’Hara, R. Viskanta, Heat transfer by simultaneous conduction and radiation for two absorbing media in intimate contact, Int. J. Heat Mass Transfer 12 (1969) 333–347. [11] H. Matshushima, R. Viskanta, Effects of internal radiative transfer on natural convection and heat transfer in a vertical crystal growth configuration, Int. J. Heat Mass Transfer 33 (1990) 1957–1968. [12] E.W. Larsen, G. Thommes, A. Klar, M. Seaid, T. Götz, Simplified P–N approximations to the equations of radiative heat transfer and applications, J. Comput. Phys. 183 (2002) 652–675. [13] I.S. Habib, Solidification of semitransparent materials by condcution and radiation, Int. J. Heat Mass Transfer 14 (1971) 2161–2164. [14] I.S. Habib, Solidification of semitransparent cylndrical medium by condcution and radiation, Int. J. Heat Mass Transfer 95 (1973) 37–41. [15] M. Abrams, R. Viskanta, The effect of radiative heat transfer upon the melting and solidification of semitransparent crystals, ASME J. Heat Transfer 96 (1974) 184–190. [16] S.H. Chan, D.H. Cho, G. Kocamustafaogullari, Melting and solidification with internal radiative transfer –A generalized phase change model, Int. J. Heat Mass Transfer 26 (1983) 621–633. [17] S. Brandon, J.J. Derby, Internal radiative transport in vertical bridgman growth of semitransparent crystals, J. Cryst. Growth 110 (1991) 481–500. [18] S. Brandon, J.J. Derby, Heat transfer in vertical bridgman growth of oxides: Effects of conduction, convection, and internal radiation, J. Cryst. Growth 121 (1992) 473–494. [19] S. Brandon, D. Gazit, A. Horowitz, Interface shapes and thermal fields during the gradient solidification method growth of sapphire single crystal, J. Cryst. Growth 167 (1996) 190–207. [20] T. Tsukada, K. Kakinoki, M. Hozawa, N. Imaishi, Effect of internal radiation within crystal and melt of czochralski crystal growth of oxide, Int. J. Heat Mass Transfer 38 (1995) 2707–2714. [21] C.W. Lan, C.Y. Tu, Y.F. Lee, Effects of internal radiation on heat flow and facet formation in bridgman growth of YAG crystals, Int. J. Heat Mass Transfer 46 (2003) 1629–1640.

Fig. 16. Mass fraction of the glass liquid phase as a function of time.

and the problem has to be analyzed much longer in order to see when the glass starts to solidify. Fig. 16 shows the mass fraction of the glass liquid phase as a function of time and, as one can see, the system needs almost 400 s to cool down to the point where the glass solidification stars. 5. Conclusions The multiphase flows at high temperatures are very often accompanied with phase-change processes. The research that treats radiative heat transfer in multiphase flows with phase change is mainly focused on the interface between solid and liquid phases, while the interfaces between liquid and gas or solid and gas phases are neglected or their treatment is significantly simplified. Here, a non-equilibrium rate-driven phase-change model for an arbitrary number of phases is developed and combined with the previously developed and validated method for simulating multiphase flows at high temperatures. A set of cases is used to verify the mathematical model, numerical discretization and solution algorithm. The phase-change model is verified on several examples and it is shown that it behaves like the equilibrium model if the relaxation times are small. A case where a liquid in contact with air solidifies in a rotating tank is designed to demonstrate the importance of a simultaneous modeling of all relevant transport processes. Unlike the cases that can be found in the literature, where the liquid–gas and solid–gas interfaces are considered to be fixed, here the interfaces between all phases could move. The designed case clearly indicates that neglecting radiative heat transfer and motion 12

International Journal of Thermal Sciences 149 (2020) 106201

A. Hasečić et al.

[36] A. Hasečić, Numerical Modeling of Radiative Heat Transfer in Multiphase Flows with Free Surface (Ph.D. thesis), Mechanical Engineering Faculty, University of Sarajevo, 2017. [37] I. Demirdžić, M. Perić, Space conservation law in finite volume calculations of fluid flow, Internat. J. Numer. Methods Fluids 8 (1988) 1037–1050. [38] V.R. Voller, Numerical methods for phase-change problems, in: W.J. Minkowycz, E.M. Sparrow, J.Y. Murthy (Eds.), Handbook on Numerical Heath Transfer, second ed., John Wiley & Sons, Inc., 2006, pp. 593–622 (Chapter 19). [39] K.R. Sultana, S.R. Dehghani, K. Pope, Y.S. Muzychka, Numerical techniques for solving solidification and melting phase change problems, Numer. Heat Transfer B 73 (2018) 129–145. [40] S. Muzaferija, D. Papoulias, M. Peric, VOF Simulations of hydrodynamic cavitation using the asymptotic and classical Rayleigh-Plesset models, in: Fifth International Symposium on Marine Propulsors, Espoo, Finland, June 12–15, 2017, pp. 50–57. [41] A. Teskeredžić, I. Demirdžić, S. Muzaferija, Numerical method for heat transfer, fluid flow, and stress analysis in phase-change problems, Numer. Heat Transfer B 42 (2002) 437–459. [42] G.-X. Wang, V. Prasad, Non-equilibrium phenomena in rapid solidification: Theoretical treatment for process modeling, Microscale Thermophys. Eng. 1 (1997) 143–157. [43] C.P. Grigoropoulos, Heat transfer in laser processing of thin films, Annu. Rev. Heat Transfer 7 (1994) 7–130. [44] Z. Bilicki, J. Kestin, Physical aspects of the relaxation model in two-phase flow, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 428 (1990) 379–397. [45] I. Demirdžić, S. Muzaferija, Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology, Comput. Methods Appl. Mech. Engrg. 125 (1995) 235–255. [46] C.M. Rhie, W.L. Chow, A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation, AIAA J. 21 (1983) 1525–1532. [47] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer 15 (1972) 1787–1806. [48] C.D. Perez-Segarra, C. Farre, J. Cadafalch, A. Oliva, Analysis of different numerical schemes for the resolution of convection-diffusion equations using finite-volume methods on three-dimensional unstructured grids. Part I: Discretization schemes, Numer. Heat Transfer B 49 (2006) 333–350. [49] M.N. Ozisik, Heat Conduction, Willey, New York, 1980. [50] W. Dijkstra, R.M.M. Mattheij, Numerical modeling of the blowing phase in the production of glass containers, Electron. J. Bound. Elem. 6 (2008) 1–23. [51] C.G. Giannopapa, J.A.W.M. Groot, A computer simulation model for the blowblow forming process of glass containers, in: Proceedings of PVP2007, 2007 ASME Pressure Vessels and Piping Division Conference, vol. 4, 2007, pp. 79–86.

[22] C.J. Jing, S. Ihara, K.I. Sugioka, T. Tsukada, M. Kobayashi, M. Mito, C. Yokoyama, Global analysis of heat transfer considering three-dimensional unsteady melt flow in CZ crystal growth of oxide, J. Cryst. Growth 307 (2007) 235–244. [23] C.J. Jing, S. Ihara, K.I. Sugioka, T. Tsukada, M. Kobayashi, Global analysis of heat transfer in CZ crystal growth of oxide taking into account three-dimensional unsteady melt convection: Effect of meniscus shape, J. Cryst. Growth 310 (2008) 204–213. [24] P.R. Parida, R. Raj, A. Prasad, S.C. Mishra, Solidification of a semitransparent planar layer subjected to radiative and convective cooling, J. Quant. Spectrosc. Radiat. Transfer 107 (2007) 226–235. [25] D. Makhanlall, L.H. Liu, Second law analysis of coupled conduction-radiation heat transfer with phase change, Int. J. Therm. Sci. 49 (2010) 1829–1836. [26] S.C. Mishra, N.C. Behara, A.K. Garg, A. Mishra, Solidification of a 2-d semitransparent medium using the lattice Boltzmann method and the finite volume method, Int. J. Heat Mass Transfer 51 (2008) 4447–4460. [27] B. Mondal, S.C. Mishra, Numerical analysis of solidification of a 3-d semitransparent medium in presence of volumetric radiation, Int. J. Therm. Sci. 48 (2009) 1116–1128. [28] P. Łapka, P. Furmański, Fixed cartesian grid based numerical model for solidification process of semi-transparent materials I: Modelling and verification, Int. J. Heat Mass Transfer 55 (2012) 4941–4952. [29] P. Łapka, P. Furmański, Fixed cartesian grid based numerical model for solidification process of semi-transparent materials II: Reflection and refraction or transmission of the thermal radiation at the solid-liquid interface, Int. J. Heat Mass Transfer 55 (2012) 4953–4964. [30] C. Yao, Combined Heat Transfer with Melting/Solidification in Semitransparent Materials (Ph.D. thesis), The Graduate Faculty of The University of Akron, 1999. [31] C. Yao, G.-X. Wang, B.T.F. Chung, Nonequilibrium planar interface model for solidification of semitransparent radiating model, J. Thermophys. Heat Transfer 14 (2000) 297–304. [32] C. Yao, B.T.F. Chung, G.-X. Wang, Mushy zone equilibrium solidification of a semitransparent layer subject to radiative and convective cooling, Int. J. Heat Mass Transfer 45 (2002) 2397–2405. [33] C. Yang, G. Wen, Q. Sun, P. Tang, Investigation of the coupled conductive and radiative heat transfer of molten slag in a cylindrical enclosure based on the zonal method, Int. J. Numer. Heat Transfer 110 (2017) 523–538. [34] A. Hasečić, S. Muzaferija, I. Demirdžić, Finite volume method for radiative transport in multiphase flows with free surfaces, Numer. Heat Transfer A 70 (2016) 347–365. [35] M.F. Modest, Radiative Heat Transfer, third ed., Academic Press, New York, 2013.

13