Regular Article
Modelling of multicomponent droplet evaporation under cryogenic conditions
^{1}
Université de Lorraine, École Nationale Supérieure des Industries Chimiques, Laboratoire Réactions et Génie des Procédés (UMR CNRS 7274), 1 rue Grandville, 54000 Nancy, France
^{2} GTT, 1 Route de Versailles, 78470 SaintRémylèsChevreuse, France
^{*} Corresponding author: jeannoel.jaubert@univlorraine.fr
Received:
10
March
2020
Accepted:
9
September
2020
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 multicomponent 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 film of sufficient thinness was assumed to envelope the droplet, while the PengRobinson 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.
© X. Xu et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
A wide body of studies exists on the vaporization of drops of multicomponent 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 multicomponent 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 PengRobinson Equation of State (PREoS) in which the binary interaction parameters were relevant for low temperatures. A final calculation algorithm for simulation is then proposed.
2 General assumptions
The following hypotheses are assumed in modelling the spraying process:

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.
Fig. 1 Schematic of the LNG spraying process.
The NG phase has uniform properties (temperature, pressure and composition).
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.
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 (Δt_{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.
3 LNG spraying process – the 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.
3.1 Droplet modelling
In the spraying process, each droplet falls, and as it falls it exchanges mass and heat with the bulk gas. Therefore, a multicomponent 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.
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.
Fig. 2 The force field of a droplet. 
In Figure 2, u_{d} is the velocity vector of the droplet, and θ_{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:(1)(2)with(3)where m_{d} is the mass of a droplet. By assuming the droplets to be spherical, m_{d} can be calculated from:(4)
And F_{G}, F_{A}, and F_{D} can be calculated using the following equations:(5)(6)(7)where g is the standard gravitational acceleration, ρ_{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:(8)
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):(10)
In the above equations, σ_{d} is the surface tension of the liquid (droplet) and the reduced dynamic viscosity μ* is defined as:(11)where μ_{d} and μ_{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 → ∞} 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.
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.
Fig. 3 Schematic of a droplet enveloped by a gaseous film. 
If assuming the quasisteady species balance around the droplet, the instantaneous evaporation mass flow rate (ṁ_{d}) of the droplet can be estimated by equation [3]:(12)
In the above equation, N is the number of species in the droplet, is the evaporation mass flow rate of species i in the droplet. Note that would be negative in cases where the species i is condensed into the droplet. Furthermore: ρ_{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 ψ_{i} in the droplet, and it is defined as [3]:(13)ψ_{i} can be estimated with equation(14)where is the mole fraction of species i in the droplet, and is the molecular volume parameter of species i from UNIQUAC model [4].
is the modified Sherwood number of species i. accounts for the effect of a high mass transfer rate and it is estimated by equation [3, 5]:(15)
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]:(16)
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:(17)μ_{i,f} and ρ_{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:(18)Y_{i,S} and Y_{i,∞} 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,∞} is the mass fraction of species i in the bulk gas. Y_{i,S} can be estimated by PengRobinson Equation of State (PREoS) [9], which will be introduced in Section 3.4.
As illustrated in Figure 3, the accumulated thermal energy of the droplet is denoted . Its expression is basically:(19)where m_{d} is the mass of the droplet and is the specific heat capacity of the droplet at constant pressure.
In practice, equation (19) is used to express the derivative . To do so, it is necessary to dispose of another expression for . This heat flux can be deduced from the Spalding heat transfer number of species i, , the definition of which is [3]:(20)where is the contribution of species i to (this latter is obtained by summing the quantities over all the species); is the mean specific heat capacity of species i in the gas film. and are the droplet temperature and the bulkgas temperature, respectively. is the specific enthalpy of vaporization for species i.
Therefore, the knowledge of makes it possible to estimate :(21)
For estimating the Spalding heat transfer number of species i , we refer to Brenn et al. [3] and Abramzon and Sirignano [5] who postulate that is coupled with the mass transfer number through:(22)
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:(23)where λ_{i,f} is the thermal conduction coefficient of species i in the gas film.
In equation (22), 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:(24)
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 masstransfer rates:(25)Pr_{i} is the Prandtl number of species i defined by the equation:(26)
Note that, should be solved from equations (22) to (24) by iterations; its initial value can be assumed to be equal to B_{M,i}.
Finally, the rate of temperature change of the droplet can be calculated by equating equations (19) and (21):(27)
As illustrated in Figure 3, the heat flux exiting the film is simply given by:(28)
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:(29)
Such an equation makes it possible to calculate which is the heat flux from the bulk gas into the film. The average properties of the film (such as μ_{i,f}, D_{i,f}, C_{p,i,f}, λ_{i,f} and ρ_{i,f}) are estimated at the average temperature of the film T_{f}, given [3, 5] by the equation:(30)
It is worth noting that in the evaporation model can be positive (which indicates evaporation from the droplet) or negative (which indicates condensation to the droplet). Heat fluxes (, and ) can be positive or negative also. All computations described in this section are based upon one droplet.
3.2 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:(31)
The total heat flux from the bulk gas into the droplets can be computed with the following equation:(32)
The temperature change rate of the bulk gas in the tank can be calculated using the following equation:(33)where m_{g} is the total mass of the gas phase and is the specific heat capacity of the gas.
3.3 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:(34) and are the total molar amounts of the liquid phase before and after mixing, respectively. is the molar amount of the influx fluid. The new mole fraction of each species in the liquid can be computed with this equation:(35) and are the mole fractions of species i in the liquid phase before and after mixing, respectively. While 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:(36)
Assuming that the heat capacity of the liquid phase and the heat capacity of the influx liquid phase are constant (i.e. temperatureindependent), equation (36) can be converted to:(37)
This equation makes it possible to estimate the new liquid temperature .
3.4 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 PengRobinson Equation of State (PREoS) [7, 8] accounting for the specific behaviour of the various LNG molecules at low temperature and properly reflecting the corresponding molecular interactions.
At the dropletgas film interface (see Fig. 3), the condition of the vapourliquid phase equilibrium is given by the equations:(38)where x_{i,d} is the mole fraction of component i in the droplet, and y_{i,s} is the mole fraction of component i in the gas film adjacent to the droplet. The fugacity coefficient of each component in the liquid phase and in the vapour phase are calculated with the PREoS, the mathematical expression of which is:(39)
In this paper, the PREoS was not volumetranslated [9–11] and classical van der Waals mixing rules [12] were used:(40)
It is today acknowledged that the Soave [13] αfunction is nonconsistent [14–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 αfunction was chosen:(41)
The fugacity coefficient can be calculated [12] from the equation:(42)where(43)
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 covolume 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 ω_{i} is the experimental acentric factor of pure i. In this paper, the binary interaction parameters k_{ij} were predicted using the EPPR78 model [17–19].
4 Proposed algorithm
4.1 Modelling of one droplet in one simulation time step
The droplet’s properties are updated after each time step (Δt), which was preset for the simulation. The calculation steps for one droplet are listed below.
First we load the earlier properties of the droplet, i.e. velocity (u_{d}), angle between directions of motion and gravity (θ_{d}), height from the liquid phase surface (h_{d}), radius (r_{d}), temperature (T_{d}), composition (x_{i,d}) and amount (n_{d}). For the first step time, i.e., when the spraying starts, such data were calculated from the parameters and the initial values of the simulation (see details in Sect. 4.3), which need to be entered by user.
Then we compute the droplet’s new motion properties (u_{d} and θ_{d}) with the equations:
Next we compute the evaporation mass flow rates of each component () of the droplet with equation (12) (using Brenn’s method [3]).
If the droplet has vaporized completely, we compute the heat flux from the bulk gas and then skip the next steps.
Now we compute the heat fluxes (, and ) between the droplets, the gas film and the bulk gas with equations (21), (28) and (29). The quantity is obtained from equation (27) while is deduced from equation (33).
Then we compute the new mole amount (n_{d}) and composition (x_{i,d}) of the droplet with the equations:
Now we compute the new temperature (T_{d}) of the droplet using equation (27) and the following equation:
Now we compute the new molar volume (v_{m,d}) of the droplet using the PREoS 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 (Δt), a single droplet’s properties (i.e. velocity (u_{d}), angle between directions of motion and gravity (θ_{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 () and heat flux , 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.
Fig. 4 Flowchart detailing the modelling of a droplet layer and named LayerSim. 
4.2 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.
Compute the new spray droplet properties layer by layer, following the instructions in Section 4.1.
If the droplet layer reaches the liquid phase, update the liquid phase properties (T_{liq}, x_{i,liq}, and n_{liq}) with equations (34), (35), and (37). Once done, the volume of the liquid phase (V_{liq}) can be calculated with the help of the PREoS.
Compute the new mass (m_{g}) and compositions (y_{i,g}) of the gas phase using equation (31) and the following equations:
Compute the new temperature of the gas phase using equation (33) with and the following equation:
Compute the new volume of the gas phase as follows:
Calculate the gas pressure using the PREoS with the new temperature, volume and composition values.
If the gas phase is condensed, compute the properties of gas and condensed liquid with the PREoS. 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.
Fig. 5 Flowchart showing the modelling of one tank and named TankSim. 
4.3 Modelling overview
The algorithm for the master simulation program is shown in Figure 6. The sequence is as follows:
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.
Compute the properties of the gas, the liquid and the spray droplets tank by tank, following the instructions given in Section 4.2.
At each fixed time interval (Δt_{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).
Compute the equilibrations between tanks, if more than one tank is included in the spraying process.
Check the spraying time. If the spraying is not finished, go to the next simulation timestep.
Fig. 6 Flowchart of the master program. 
5 Results
5.1 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.
Fig. 7 An example case in which the droplet size is decreasing (due to evaporation). Temperature of gas = 173.15 K, temperature of droplet = 113.15 K, gas composition: 99% methane + 0.5% ethane + 0.5% nitrogen, droplet composition: 99.9% methane + 0.05% ethane + 0.05% nitrogen. 
Fig. 8 An example case in which the droplet size is increasing. Temperature of gas = 183.15 K, temperature of droplet = 103.15 K, gas composition: 9% methane + 80% ethane + 11% nitrogen, droplet composition: 19% methane + 80% ethane + 1% nitrogen. 
Fig. 9 An example case in which the droplet size is constant. Temperature of gas = 113.15 K, temperature of droplet = 113.15 K, gas composition: 99% methane + 0.5% ethane + 0.5% nitrogen, droplet composition: 94% methane + 5.9% ethane + 0.1% nitrogen. 
5.2 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 environment has been considered. As expected, the pressure and temperature of the gas phase decrease during the spraying process.
Fig. 10 Simulated gas phase temperature and pressure versus spraying time for an example case. 
This example highlights the usefulness of the spraying process for cooling down a liquidgas 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.
6 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 multicomponent 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.
References
 International Energy Agency (IEA). (2011) Are we entering a golden age of gas? Special report, World Energy Outlook 2011. https://www.iea.org/publications/freepublications/publication/WEO2011_GoldenAgeofGasReport.pdf [Google Scholar]
 Loth E. (2008) Quasisteady shape and drag of deformable bubbles and drops, Int. J. Multiph. Flow. 34, 523–546. https://doi.org/10.1016/j.ijmultiphaseflow.2007.08.010. [CrossRef] [Google Scholar]
 Brenn G., Deviprasath L.J., Durst F., Fink C. (2007) Evaporation of acoustically levitated multicomponent liquid droplets, Int. J. Heat Mass Transf. 50, 5073–5086. https://doi.org/10.1016/j.ijheatmasstransfer.2007.07.036. [Google Scholar]
 Fredenslund A., Jones R.L., Prausnitz J.M. (1975) Groupcontribution estimation of activity coefficients in nonideal liquid mixtures, AIChE J. 21, 1086–1099. https://doi.org/10.1002/aic.690210607. [Google Scholar]
 Abramzon B., Sirignano W.A. (1989) Droplet vaporization model for spray combustion calculations, Int. J. Heat Mass Transf. 32, 1605–1618. https://doi.org/10.1016/00179310(89)900434. [Google Scholar]
 Frössling N. (1938) Über die verdunstung fallender tropfen, Gerlands Beitr. Zur Geophys. 52, 170–216. [Google Scholar]
 Robinson D.B., Peng D.Y. (1978) The characterization of the heptanes and heavier fractions for the GPA PengRobinson programs (RR28), Res. Rep. GPA. 1–36. [Google Scholar]
 PinaMartinez A., Privat R., Jaubert J.N., Peng D.Y. (2019) Updated versions of the generalized Soave αfunction suitable for the RedlichKwong and PengRobinson equations of state, Fluid Phase Equilibria. 485, 264–269. https://doi.org/10.1016/j.fluid.2018.12.007. [Google Scholar]
 Jaubert J.N., Privat R., Le Guennec Y., Coniglio L. (2016) Note on the properties altered by application of a Pénelouxtype volume translation to an equation of state, Fluid Phase Equilibria. 419, 88–95. https://doi.org/10.1016/j.fluid.2016.03.012. [Google Scholar]
 Privat R., Jaubert J.N., Le Guennec Y. (2016) Incorporation of a volume translation in an equation of state for fluid mixtures: which combining rule? Which effect on properties of mixing? Fluid Phase Equilibria. 427, 414–420. https://doi.org/10.1016/j.fluid.2016.07.035. [Google Scholar]
 Péneloux A., Rauzy E., Freze R. (1982) A consistent correction for RedlichKwongSoave volumes, Fluid Phase Equilibria. 8, 7–23. https://doi.org/10.1016/03783812(82)800022. [Google Scholar]
 Kontogeorgis G.M., Privat R., Jaubert J.N. (2019) Taking another look at the van der waals equation of state–almost 150 years later, J. Chem. Eng. Data. 64, 11, 4619–4637. https://doi.org/10.1021/acs.jced.9b00264. [Google Scholar]
 Soave G. (1972) Equilibrium constants from a modified RedlichKwong equation of state, Chem. Eng. Sci. 27, 1197–1203. https://doi.org/10.1016/00092509(72)800964. [Google Scholar]
 Le Guennec Y., Lasala S., Privat R., Jaubert J.N. (2016) A consistency test for αfunctions of cubic equations of state, Fluid Phase Equilibria. 427, 513–538. https://doi.org/10.1016/j.fluid.2016.07.026. [Google Scholar]
 Le Guennec Y., Privat R., Jaubert J.N. (2016) Development of the translatedconsistent tcPR and tcRK cubic equations of state for a safe and accurate prediction of volumetric, energetic and saturation properties of pure compounds in the sub and supercritical domains, Fluid Phase Equilibria. 429, 301–312. https://doi.org/10.1016/j.fluid.2016.09.003. [Google Scholar]
 Le Guennec Y., Privat R., Lasala S., Jaubert J.N. (2017) On the imperative need to use a consistent αfunction for the prediction of purecompound supercritical properties with a cubic equation of state, Fluid Phase Equilibria. 445, 45–53. https://doi.org/10.1016/j.fluid.2017.04.015. [Google Scholar]
 Jaubert J.N., Mutelet F. (2004) VLE predictions with the PengRobinson equation of state and temperaturedependent k_{ij} calculated through a group contribution method, Fluid Phase Equilibria. 224, 285–304. https://doi.org/10.1016/j.fluid.2004.06.059. [Google Scholar]
 Qian J.W., Jaubert J.N., Privat R. (2013) Phase equilibria in hydrogencontaining binary systems modeled with the PengRobinson equation of state and temperaturedependent binary interaction parameters calculated through a groupcontribution method, J. Supercrit. Fluids. 75, 58–71. https://doi.org/10.1016/j.supflu.2012.12.014. [Google Scholar]
 Vitu S., Privat R., Jaubert J.N., Mutelet F. (2008) Predicting the phase equilibria of CO_{2} + hydrocarbon systems with the PPR78 model (PR EoS and k_{ij} calculated through a group contribution method), J. Supercrit. Fluids. 45, 1–26. https://doi.org/10.1016/j.supflu.2007.11.015. [Google Scholar]
All Figures
Fig. 1 Schematic of the LNG spraying process. 

In the text 
Fig. 2 The force field of a droplet. 

In the text 
Fig. 3 Schematic of a droplet enveloped by a gaseous film. 

In the text 
Fig. 4 Flowchart detailing the modelling of a droplet layer and named LayerSim. 

In the text 
Fig. 5 Flowchart showing the modelling of one tank and named TankSim. 

In the text 
Fig. 6 Flowchart of the master program. 

In the text 
Fig. 7 An example case in which the droplet size is decreasing (due to evaporation). Temperature of gas = 173.15 K, temperature of droplet = 113.15 K, gas composition: 99% methane + 0.5% ethane + 0.5% nitrogen, droplet composition: 99.9% methane + 0.05% ethane + 0.05% nitrogen. 

In the text 
Fig. 8 An example case in which the droplet size is increasing. Temperature of gas = 183.15 K, temperature of droplet = 103.15 K, gas composition: 9% methane + 80% ethane + 11% nitrogen, droplet composition: 19% methane + 80% ethane + 1% nitrogen. 

In the text 
Fig. 9 An example case in which the droplet size is constant. Temperature of gas = 113.15 K, temperature of droplet = 113.15 K, gas composition: 99% methane + 0.5% ethane + 0.5% nitrogen, droplet composition: 94% methane + 5.9% ethane + 0.1% nitrogen. 

In the text 
Fig. 10 Simulated gas phase temperature and pressure versus spraying time for an example case. 

In the text 