Modelling of multi-component droplet evaporation under cryogenic conditions

. The vaporization of drops of highly vaporizable liquids falling inside a cryogenic environment is far from being a trivial matter as it assumes harnessing specialized thermodynamics and physical equations. In this paper, a multi-component falling droplet evaporation model was developed for simulating the spray cooling process. The falling speed of the sprayed droplets was calculated with the momentum equations considering three forces (gravity, buoyancy and drag) applied to a droplet. To evaluate the mass and heat transfer between the sprayed droplet and the surrounding gas phase, a gaseous boundary ﬁ lm of suf ﬁ cient thinness was assumed to envelope the droplet, while the Peng-Robinson equation of state was used for estimating the phase equilibrium properties on the droplet ’ s surface. Based on the relevant conservation equations of mass and energy, the key properties (such as temperature, pressure and composition) of the liquid and gas phases in the tank during the spray process could be simulated. To conclude, the simulation algorithm is proposed.


Introduction
A wide body of studies exists on the vaporization of drops of multi-component liquids at moderate and high temperatures. By contrast, the same process has been the subject of relatively little research when it occurs at very low temperatures, especially under cryogenic conditions. This scarcity of literature is namely due to the need to resort to complex thermodynamics and access to specialized thermal and transport equations. The matter is all the more complex when the drops follow a falling motion that adds a dynamic dimension to the subject. The vaporization of drops of Liquefied Natural Gas (LNG) well exemplifies this statement.
Natural Gas (NG) is a flexible fuel that is used extensively in power generation, industrial and household consumption, as well as the production of advanced petrochemical derivatives. Compared with other fossil fuels, natural gas creates lower emissions of greenhouse gases and local pollutants, and is therefore expected to play a greater role in the future global energy mix. Natural gas can be delivered either by high pressure pipelines or, depending on the location of the gas field and the security of supply, it can be liquefied and then transported by Liquefied Natural Gas Carriers (LNGC) [1]. This part of the voyage is called the laden voyage. On arrival, the liquefied natural gas is unloaded and the LNGC travels to another loading terminal. This voyage is called the ballast voyage, due to the fact that the ship's ballast compartments are full of water while the LNG tanks are almost empty. Some LNG remains in the tanks, as the LNGC uses LNG as fuel. The Liquefied Natural Gas (LNG) is stored in highly insulated storage tanks at pressures slightly above atmospheric and temperatures close to boiling (% 111 K). During the laden and ballast voyages, some of the LNG will vaporize due to unavoidable heat ingress into the storage tank from its surroundings. Such a generated vapour is used as fuel gas in the main engines to drive the ship.
During the ballast voyage, the temperature of the tank significantly increases. In order to avoid excessive vapour generation during the next loading operation, the tanks are cooled down for a few days before loading. Spraying some LNG at the top of each tank cools down the vapour as well as part of the insulation system. Obviously, vapour and liquid temperatures, pressure and compositions have a significant impact on the spraying operation.
The aim of this paper is to propose a numerical method for simulating this spraying process. A multi-component falling droplet evaporation model is developed to estimate the mass and heat transfer between the sprayed droplets and the surrounding gas phase. The thermodynamic aspects are dealt with using the Peng-Robinson Equation 1. The set of connected storage tanks (as shown in Fig. 1) is considered to form an isolated system with no mass or heat exchange with its surroundings. The spraying time is normally quite short compared to the duration of LNG transportation, and thus the heat ingress from the surroundings is relatively small. 2. The NG phase has uniform properties (temperature, pressure and composition). 3. The LNG phase has also uniform properties (temperature, pressure and composition). The vaporization of the liquid phase is ignored, and the mass or heat exchange between the LNG phase and NG phase is also ignored. 4. Each droplet has uniform properties (temperature and composition), and its shape is assumed to be approximately spherical. As a simplification, the droplet distribution will be assumed to be monodisperse so that the droplets generated in a short interval of time (Dt layer ) can be grouped into a "droplet layer": all drops of a given layer will have the same vaporization history as well as the same falling kinetics and will come back to the LNG phase at the same time (in the case they are not totally vaporized before reaching the liquid phase). This notion of "droplet layer", in which all droplets have the same properties, will be useful to reduce computation time.

LNG spraying processthe modelling approach
The modelling approach uses mass and heat balances to predict the physical properties (temperature, pressure and density) and the amounts of the bulk gas phase (NG), liquid phase (LNG) and spraying droplets in the storage tanks. The approach will consist in estimating the mass and heat exchanges between the sprayed droplets and the bulk gas phase with the suitable models. In this section, the models and methods used to describe the sprayed droplets, the gas phase and the liquid phase are detailed.

Droplet modelling
In the spraying process, each droplet falls, and as it falls it exchanges mass and heat with the bulk gas. Therefore, a multi-component falling droplet evaporation model including motion equations and evaporation rate equations (for estimating heat and mass transfer between droplets and gas phase) is developed to predict the velocity, temperature, composition and mass of the sprayed droplets. While the general mass and heat transfer equations remain conventional ones, the main difficulties of this task reside in the requirement to dispose of an equation of state as well as properties correlations that fit the low temperature conditions.
(a) Motion of a droplet As shown in Figure 2, three forces, gravity (F G ), buoyancy or Archimedes' force (F A ), and drag force (F D ), act on a liquid droplet once it leaves the sprayer.
In Figure 2, u d is the velocity vector of the droplet, and h d is the angle between the droplet's motion and the direction of gravitational pull. According to Newton's second law, we obtain the following equations: where m d is the mass of a droplet. By assuming the droplets to be spherical, m d can be calculated from:  And F G , F A , and F D can be calculated using the following equations: where g is the standard gravitational acceleration, q g is the density of the gas mixture, and r d is the droplet's radius. The drag coefficient C d is classically correlated with Weber number "We" and Reynolds number "Re".
In this study, the empirical equations developed by Loth [2] and used: See equation (8) bottom of the page Re is obtained from: The Weber number (We), which indicates the degree of shape deformation of the falling droplet, is the ratio of continuous fluid stresses (which cause deformation) to surface tension stresses (which resist deformation): In the above equations, r d is the surface tension of the liquid (droplet) and the reduced dynamic viscosity l* is defined as: where l d and l g are the dynamic viscosities of the gas and liquid (droplet), respectively. In Loth's approach [2], C d,We ? 0 denotes the drag coefficient for a spherical solid particle, which is not deformable when falling, and C d,We ? 1 denotes the drag coefficient for a bubble with maximum deformation. As discussed in the paper by Loth [2], this method can accurately predict the drag coefficients of various airborne droplets under the conditions: We 12 and 400 Re 7000.

(b) Evaporation of one droplet
The evaporation model is employed to estimate the heat and mass exchange between the droplet and the gas phase.
As shown in Figure 3, to evaluate mass and heat fluxes, a sufficiently thin gaseous boundary film is assumed to be around the droplet.
If assuming the quasi-steady species balance around the droplet, the instantaneous evaporation mass flow rate (ṁ d ) of the droplet can be estimated by equation [3]: In the above equation, N is the number of species in the droplet, _ m i;d is the evaporation mass flow rate of species i in the droplet. Note that _ m i;d would be negative in cases where the species i is condensed into the droplet. Furthermore: q f is the density of the gas film, and D i,f is the mean diffusion coefficient of species i in the gas film. r V,i is the volume equivalent partial radius of component i corresponding to its instantaneous volume fraction w i in the droplet, and it is defined as [3]: w i can be estimated with equation where x i;d is the mole fraction of species i in the droplet, and r i;UNIQUAC is the molecular volume parameter of species i from UNIQUAC model [4]. Sh Ã i is the modified Sherwood number of species i. Sh Ã i accounts for the effect of a high mass transfer rate and it is estimated by equation [3,5]: Sh i,0 is the Sherwood number of species i for the droplet with a small mass transfer rate, and it can be given by the Frössling correlation [6]: Reynolds number (Re) can be computed with equation (9). Sc i is the Schmidt number of species i in the gas film, and is defined by equation: l i,f and q i,f are the dynamic viscosity and the partial density of species i in the gas film, respectively.
In equations (12) and (15), B M,i is the Spalding mass transfer number of species i, and is defined by equation: Y i,S and Y i,1 are the mass fraction of species i at the inner sides (which is close to the droplet surface) and the outer side (which is at infinity) of the gaseous film (see Fig. 3), respectively. Y i,1 is the mass fraction of species i in the bulk gas. Y i,S can be estimated by Peng-Robinson Equation of State (PR-EoS) [9], which will be introduced in Section 3.4. As illustrated in Figure 3, the accumulated thermal energy of the droplet is denoted _ Q d . Its expression is basically: where m d is the mass of the droplet and C p;d is the specific heat capacity of the droplet at constant pressure. In practice, equation (19) is used to express the derivative dT d dt . To do so, it is necessary to dispose of another expression for _ Q d . This heat flux can be deduced from the Spalding heat transfer number of species i, B T;i , the definition of which is [3]: where _ Q d;i is the contribution of species i to _ Q d (this latter is obtained by summing the _ Q d;i quantities over all the species); C p;i;f is the mean specific heat capacity of species i in the gas film. T d and T 1 are the droplet temperature and the bulk-gas temperature, respectively. Á vap;i H is the specific enthalpy of vaporization for species i.
Therefore, the knowledge of B T ;i makes it possible to estimate _ Q d : For estimating the Spalding heat transfer number of species i B T;i , we refer to Brenn et al. [3] and Abramzon and Sirignano [5] who postulate that B T;i is coupled with the mass transfer number B M;i through: C p;i;g is the mean specific heat capacity of species i in the bulk gas. Le i is the Lewis number of species i defined [3,5] by: where k i,f is the thermal conduction coefficient of species i in the gas film.
In equation (22), Nu Ã i is the modified Nusselt number of species i estimated by equation (24). Following Abramzon and Sirignano [5], it was indeed assumed that such an equation could be used for the case of an evaporating droplet: Nu i;0 is the Nusselt number of species i for the nonevaporating droplets, that can be computed [3,6] using the following correlation especially developed for spheres from experimental data on mass-transfer rates: Pr i is the Prandtl number of species i defined by the equation: Note that, B T ;i should be solved from equations (22) (19) and (21): As illustrated in Figure 3, the heat flux exiting the film is simply given by: The film around the droplet can be assumed to be at steady state (no energy nor mass accumulation). Consequently, as shown in Figure 3, the energy balance equation for the film can be written as: Such an equation makes it possible to calculate _ Q g which is the heat flux from the bulk gas into the film. The average properties of the film (such as l i,f , D i,f , C p,i,f , k i,f and q i,f ) are estimated at the average temperature of the film T f , given [3,5] by the equation: It is worth noting that in the evaporation model _ m i;d can be positive (which indicates evaporation from the droplet) or negative (which indicates condensation to the droplet). Heat fluxes ( _ Q d , _ Q f and _ Q g ) can be positive or negative also. All computations described in this section are based upon one droplet.

Gas phase modelling
As noted in the above section, the mass and heat exchanges between a single droplet and the gas phase can be calculated with equations (12) and (29) respectively. Assuming N L layers of sprayed droplets in the tank, with each layer containing N d droplets, the total mass transferred into the bulk gas can be computed with the following equation: The total heat flux from the bulk gas into the droplets can be computed with the following equation: The temperature change rate of the bulk gas dTg dt in the tank can be calculated using the following equation: where m g is the total mass of the gas phase and C p;g is the specific heat capacity of the gas.

Liquid phase modelling
If the sprayed droplet reaches the liquid phase, or if the gas phase is condensed (due to a decrease in temperature), the liquid phase properties need to be updated. Assuming that the influx fluid (from either the falling droplet or the condensed gas) is instantaneously mixed with the liquid phase, the new total molar amount of the liquid phase is calculated as follows: N old liq and N new liq are the total molar amounts of the liquid phase before and after mixing, respectively. N in liq is the molar amount of the influx fluid. The new mole fraction of each species in the liquid can be computed with this equation: x old i;liq and x new i;liq are the mole fractions of species i in the liquid phase before and after mixing, respectively while x in i;liq is the mole fraction of species i in the influx fluid. If the mixing enthalpy is assumed to be negligible, the energy balance equation is expressed as follows: Assuming that the heat capacity of the liquid phase C p;liq and the heat capacity of the influx liquid phase C in p;liq are constant (i.e. temperature-independent), equation (36) can be converted to: This equation makes it possible to estimate the new liquid temperature T new liq .

Thermodynamic modelling
In the simulation, an adequate thermodynamic model is necessary to express the composition of the gas film around the droplet (i.e. Y i,S ) and the densities of the liquid phase (LNG), the bulk gas phase (NG), the droplets, and the gas film. This requires resorting to equation of states like the Peng-Robinson Equation of State (PR-EoS) [7,8] accounting for the specific behaviour of the various LNG molecules at low temperature and properly reflecting the corresponding molecular interactions. At the droplet-gas film interface (see Fig. 3), the condition of the vapour-liquid phase equilibrium is given by the equations: In this paper, the PR-EoS was not volume-translated [9][10][11] and classical van der Waals mixing rules [12] were used: It is today acknowledged that the Soave [13] a-function is non-consistent [14][15][16] since it diverges at very high temperatures. This study, carried out under cryogenic conditions, is however absolutely not concerned by such an inconsistency so that the Soave a-function was chosen:

See equation (41) bottom of the page
The fugacity coefficient can be calculated [12] from the equation: where In the above equations, P is the pressure, R is the gas constant, T is the temperature, a i and b i are the cohesive parameter and molar co-volume of pure component i, v is the molar volume, z i is the mole fraction of component i, T c,i is the experimental critical temperature, P c,i is the experimental critical pressure, and x i is the experimental acentric factor of pure i. In this paper, the binary interaction parameters k ij were predicted using the E-PPR78 model [17][18][19]. The droplet's properties are updated after each time step (Dt), which was preset for the simulation. The calculation steps for one droplet are listed below.  (1) and (2). 3. Next we compute the evaporation mass flow rates of each component ( _ m i;d ) of the droplet with equation (12) (using Brenn's method [3]). 4. If the droplet has vaporized completely, we compute the heat flux _ Q g from the bulk gas and then skip the next steps. 5. Now we compute the heat fluxes ( _ Q d , _ Q f and _ Q g ) between the droplets, the gas film and the bulk gas with equations (21), (28) and (29). The quantity dT d dt is obtained from equation (27) while dT d dt is deduced from equation (33). 6. Then we compute the new mole amount (n d ) and composition (x i,d ) of the droplet with the equations: 7. Now we compute the new temperature (T d ) of the droplet using equation (27) and the following equation: 8. Now we compute the new molar volume (v m,d ) of the droplet using the PR-EoS with the new temperature (T d ) and composition (x i,d ), then compute the new radius (r d ) of the droplet with the equation: In each simulation time step (Dt), a single droplet's properties (i.e. velocity (u d ), angle between directions of motion and gravity (h d ), height from the liquid phase surface (h d ), radius (r d ), temperature (T d ), composition (x i,d ) and amount (n d )) are updated, and would be used as earlier properties in the next time step. Meanwhile a single droplet's evaporation mass flow rates of each component ( _ m i;d ) and heat flux _ Q g , which would be used for computing the gas phase properties (according to Eqs. (31) and (33)), are outputted. The algorithm for computing the properties of one droplet is summarized in Figure 4. One shall remember that all droplets in the same layer are assumed to be identical.

Modelling of one tank in one simulation time step
After calculating the new properties of all spray droplet layers, the properties of the gas phase and the liquid phase in the tank can be updated as well. The computation steps are listed below.
4. Compute the new temperature of the gas phase using equation (33) with m new g and the following equation: 5. Compute the new volume of the gas phase as follows: where V tank is the tank total volume. 6. Calculate the gas pressure using the PR-EoS with the new temperature, volume and composition values. 7. If the gas phase is condensed, compute the properties of gas and condensed liquid with the PR-EoS. We then update the liquid phase properties using equations (34), (35) and (37).  The algorithm for updating the properties of liquid and gas in the tank is summarized in Figure 5.

Modelling overview
The algorithm for the master simulation program is shown in Figure 6. The sequence is as follows: 1. Set up the parameters and the initial values of the simulation: time step, time interval of one droplet layer, volume of tank; height of tank; number of sprayers; sprayers orientation; diameter of sprayer's circular orifice; size of droplet; volume flow rate of droplets; volume of liquid phase in tank; volume of gas phase in tank; liquid level in tank; pressure in the tank; gas temperature; liquid temperature; number of components; gas and liquid compositions; number of droplets in each layer (computed from number of sprayers, volume flow rate and diameter of droplets, and time interval of one droplet layer); and initial velocity (determined by volume flow rate of droplets and diameter of sprayer's circular orifice), motion angle (depending on sprayers orientation), radius, height, temperature, mole amount and compositions of droplet in first layer. 2. Compute the properties of the gas, the liquid and the spray droplets tank by tank, following the instructions given in Section 4.2.   3. At each fixed time interval (Dt layer ), generate a new droplet layer (please note that the new droplet's initial velocity, motion angle, size, and height are determined by the parameters of sprayer and tank, and are therefore fixed during simulation, however its temperature, mole amount and compositions are changed with the updated properties of spraying liquid). 4. Compute the equilibrations between tanks, if more than one tank is included in the spraying process. 5. Check the spraying time. If the spraying is not finished, go to the next simulation time-step.

Size changes of the falling droplets
Depending on the initial conditions (such as bulk gas pressure, temperature and composition, droplet temperature and composition), the size of a sprayed droplet may decrease, remain constant or increase as it falls. A typical example showing the partial vaporization of a liquid droplet during spraying (falling inside a hotter gaseous environment) is shown in Figure 7. Conversely, the gas phase may condense onto a sprayed droplet, making it increase in size. An example is shown in Figure 8. As this figure shows, the droplet increases in size during the first 10 s of its fall before becoming almost constant as the droplet and the gas phase reach an equilibrium. Finally, if the droplet and the bulk gas phase are close to the equilibrium condition when the spraying starts, the size of this sprayed droplet remains almost constant during its fall, as shown in Figure 9.

Example case of spraying process
An example of simulated gas phase pressure and temperature variations over time is presented in Figure 10. For this simulation, only the effect on the vapour phase has been modelled: no heat ingress from the surrounding environ-ment has been considered. As expected, the pressure and temperature of the gas phase decrease during the spraying process. This example highlights the usefulness of the spraying process for cooling down a liquid-gas system made up of a cold liquid phase and a hot gas phase. The tool we propose is capable of predicting the kinetics of a cooling process and could be used to automatically control the temperature of a LNG tank during a laden voyage.

Conclusion
This paper was intended to address a relatively complex scenario that occurs in a cryogenic process where the cooling phase consists of drops of a highly vaporizable, multicomponent liquid that are sprayed into the gaseous environment they cool and evaporate during their fall. The proposed approach offers an overall simulation of all the thermal and kinetic aspects of this process, using a set of numerical methods combined with physics and thermodynamics matching the low temperature conditions. Such model describing the evaporation multi-component falling droplet enables the computation of the mass and heat exchanges between the droplet and the gas phase. This simulation allows the prediction of the properties of the gas and liquid phases in the tank during the spraying process.
The complete simulation algorithm that has been developed can be of interest to the Oil and Gas community facing similar heat and mass transfer within low temperature conditions.