Open Access
Numéro
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 75, 2020
Numéro d'article 81
Nombre de pages 10
DOI https://doi.org/10.2516/ogst/2020074
Publié en ligne 17 novembre 2020

© X. Xu et al., published by IFP Energies nouvelles, 2020

Licence Creative CommonsThis 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 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 of State (PR-EoS) 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:

  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.

    thumbnail Fig. 1

    Schematic of the LNG spraying process.

  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 (Δtlayer) 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 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.

  1. Motion of a droplet

As shown in Figure 2, three forces, gravity (FG), buoyancy or Archimedes’ force (FA), and drag force (FD), act on a liquid droplet once it leaves the sprayer.

thumbnail Fig. 2

The force field of a droplet.

In Figure 2, ud 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:mddud,xdt=-FDsinθd,$$ {m}_{\mathrm{d}}\frac{\mathrm{d}{u}_{\mathrm{d},x}}{\mathrm{d}t}=-{F}_D\mathrm{sin}{\theta }_{\mathrm{d}}, $$(1)mddud,ydt=-FDcosθd-FA+FG,$$ {m}_{\mathrm{d}}\frac{\mathrm{d}{u}_{\mathrm{d},y}}{\mathrm{d}t}=-{F}_{\mathrm{D}}\mathrm{cos}{\theta }_{\mathrm{d}}-{F}_{\mathrm{A}}+{F}_{\mathrm{G}}, $$(2)with{ud,x=udsinθdud,y=udcosθd,$$ \left\{\begin{array}{c}{u}_{\mathrm{d},x}={u}_{\mathrm{d}}\mathrm{sin}{\theta }_{\mathrm{d}}\\ {u}_{\mathrm{d},y}={u}_{\mathrm{d}}\mathrm{cos}{\theta }_{\mathrm{d}}\end{array}\right., $$(3)where md is the mass of a droplet. By assuming the droplets to be spherical, md can be calculated from:md=43πrd3ρd.$$ {m}_{\mathrm{d}}=\frac{4}{3}\pi {{r}_{\mathrm{d}}}^3{\rho }_{\mathrm{d}}. $$(4)

And FG, FA, and FD can be calculated using the following equations:FG=43πrd3ρdg,$$ {F}_{\mathrm{G}}=\frac{4}{3}\pi {r}_{\mathrm{d}}^3{\rho }_{\mathrm{d}}g, $$(5)FA=43πrd3ρgg,$$ {F}_{\mathrm{A}}=\frac{4}{3}\pi {r}_{\mathrm{d}}^3{\rho }_{\mathrm{g}}g, $$(6)FD=12πrd2ρgCdu2,$$ {F}_{\mathrm{D}}=\frac{1}{2}\pi {r}_{\mathrm{d}}^2{\rho }_{\mathrm{g}}{C}_{\mathrm{d}}{u}^2, $$(7)where g is the standard gravitational acceleration, ρg is the density of the gas mixture, and rd is the droplet’s radius. The drag coefficient Cd is classically correlated with Weber number “We” and Reynolds number “Re”. In this study, the empirical equations developed by Loth [2] and used:{Cd=Cd,We0+ΔCd*(Cd,We-Cd,We0)Cd,We0=24Re(1+0.15Re0.687)+0.421+42500Re1.16Cd,We=83+24Re(2+3μ*3+3μ*)ΔCd*=3.8×10-3(WeRe0.2)+3×10-5(WeRe0.2)2+9×10-7(WeRe0.2)3$$ \left\{\begin{array}{c}{C}_{\mathrm{d}}={C}_{\mathrm{d},\mathrm{We}\to 0}+\Delta {C}_{\mathrm{d}}^{\mathrm{*}}\left({C}_{\mathrm{d},\mathrm{We}\to \mathrm{\infty }}-{C}_{\mathrm{d},\mathrm{We}\to 0}\right)\\ {C}_{\mathrm{d},\mathrm{We}\to 0}=\frac{24}{\mathrm{Re}}\left(1+0.15{\mathrm{Re}}^{0.687}\right)+\frac{0.42}{1+\frac{42500}{{\mathrm{Re}}^{1.16}}}\\ \begin{array}{c}{C}_{\mathrm{d},\mathrm{We}\to \mathrm{\infty }}=\frac{8}{3}+\frac{24}{\mathrm{Re}}\left(\frac{2+3{\mu }^{\mathrm{*}}}{3+3{\mu }^{\mathrm{*}}}\right)\\ \Delta {C}_{\mathrm{d}}^{\mathrm{*}}=3.8\times 1{0}^{-3}\left(\mathrm{We}\cdot {\mathrm{Re}}^{0.2}\right)+3\times 1{0}^{-5}{\left(\mathrm{We}\cdot {\mathrm{Re}}^{0.2}\right)}^2+9\times 1{0}^{-7}{\left(\mathrm{We}\cdot {\mathrm{Re}}^{0.2}\right)}^3\end{array}\end{array}\right. $$(8)

Re is obtained from:Re=2rdρgudμg.$$ \mathrm{Re}=\frac{2{r}_{\mathrm{d}}{\rho }_{\mathrm{g}}{u}_{\mathrm{d}}}{{\mu }_{\mathrm{g}}}. $$(9)

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):We=2rdρgud2σd.$$ \mathrm{We}=\frac{2{r}_{\mathrm{d}}{\rho }_{\mathrm{g}}{u}_{\mathrm{d}}^2}{{\sigma }_{\mathrm{d}}}. $$(10)

In the above equations, σd is the surface tension of the liquid (droplet) and the reduced dynamic viscosity μ* is defined as:μ*=μdμg,$$ {\mu }^{\mathrm{*}}=\frac{{\mu }_{\mathrm{d}}}{{\mu }_{\mathrm{g}}}, $$(11)where μd and μg are the dynamic viscosities of the gas and liquid (droplet), respectively. In Loth’s approach [2], Cd,We → 0 denotes the drag coefficient for a spherical solid particle, which is not deformable when falling, and Cd,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.

  1. 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.

thumbnail Fig. 3

Schematic of a droplet enveloped by a gaseous film.

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]:ṁd=i=1Nṁi,d=2πi=1NrV,iρfDi,fShi*ln(1+BM,i).$$ {\dot{m}}_{\mathrm{d}}=\sum_{i=1}^N {\dot{m}}_{\mathrm{i},\mathrm{d}}=2\pi \sum_{i=1}^N {r}_{V,\mathrm{i}}{\rho }_{\mathrm{f}}{D}_{\mathrm{i},\mathrm{f}}{\mathrm{Sh}}_{\mathrm{i}}^{*}\mathrm{ln}\left(1+{B}_{\mathrm{M},\mathrm{i}}\right). $$(12)

In the above equation, N is the number of species in the droplet, ṁi,d$ {\dot{m}}_{i,\mathrm{d}}$ is the evaporation mass flow rate of species i in the droplet. Note that ṁi,d$ {\dot{m}}_{\mathrm{i},\mathrm{d}}$ would be negative in cases where the species i is condensed into the droplet. Furthermore: ρf is the density of the gas film, and Di,f is the mean diffusion coefficient of species i in the gas film. rV,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]:rV,i=rdψi1/3,$$ {r}_{V,\mathrm{i}}={r}_{\mathrm{d}}{\psi }_{\mathrm{i}}^{1/3}, $$(13)ψi can be estimated with equationψi=xi,dri,UNIQUACxi,dri,UNIQUAC,$$ {\psi }_{\mathrm{i}}=\frac{{x}_{\mathrm{i},\mathrm{d}}{r}_{\mathrm{i},\mathrm{UNIQUAC}}}{\sum {x}_{\mathrm{i},\mathrm{d}}{r}_{\mathrm{i},\mathrm{UNIQUAC}}}, $$(14)where xi,d$ {x}_{\mathrm{i},\mathrm{d}}$ is the mole fraction of species i in the droplet, and ri,UNIQUAC$ {r}_{\mathrm{i},\mathrm{UNIQUAC}}$ is the molecular volume parameter of species i from UNIQUAC model [4].

Shi*$ {\mathrm{Sh}}_{\mathrm{i}}^{\mathrm{*}}$ is the modified Sherwood number of species i. Shi*$ {\mathrm{Sh}}_{\mathrm{i}}^{\mathrm{*}}$ accounts for the effect of a high mass transfer rate and it is estimated by equation [3, 5]:{Shi*=2+Shi,0-2F(BM,i)F(BM,i)=(1+BM,i)0.7ln(1+BM,i)BM,i.$$ \left\{\begin{array}{c}{{\mathrm{Sh}}_{\mathrm{i}}}^{\mathrm{*}}=2+\frac{{\mathrm{Sh}}_{\mathrm{i},0}-2}{F({B}_{\mathrm{M},\mathrm{i}})}\\ F({B}_{\mathrm{M},\mathrm{i}})={\left(1+{B}_{\mathrm{M},\mathrm{i}}\right)}^{0.7}\frac{\mathrm{ln}\left(1+{B}_{\mathrm{M},\mathrm{i}}\right)}{{B}_{\mathrm{M},\mathrm{i}}}\end{array}\right.. $$(15)

Shi,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]:Shi,0=2+0.552Re1/2Sci1/3.$$ {\mathrm{Sh}}_{\mathrm{i},0}=2+0.552{\mathrm{Re}}^{1/2}{\mathrm{Sc}}_{\mathrm{i}}^{1/3}. $$(16)

Reynolds number (Re) can be computed with equation (9). Sci is the Schmidt number of species i in the gas film, and is defined by equation:Sci=μi,fρi,fDi,f,$$ {\mathrm{Sc}}_{\mathrm{i}}=\frac{{\mu }_{\mathrm{i},\mathrm{f}}}{{\rho }_{\mathrm{i},\mathrm{f}}{D}_{\mathrm{i},\mathrm{f}}}, $$(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), BM,i is the Spalding mass transfer number of species i, and is defined by equation:BM,i=Yi,S-Yi,1-Yi,S,$$ {B}_{\mathrm{M},\mathrm{i}}=\frac{{Y}_{\mathrm{i},\mathrm{S}}-{Y}_{\mathrm{i},\mathrm{\infty }}}{1-{Y}_{\mathrm{i},\mathrm{S}}}, $$(18)Yi,S and Yi,∞ 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. Yi,∞ is the mass fraction of species i in the bulk gas. Yi,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$ {\dot{Q}}_{\mathrm{d}}$. Its expression is basically:Q̇d=mdCp,ddTddt$$ {\dot{Q}}_{\mathrm{d}}={m}_{\mathrm{d}}{C}_{\mathrm{p},\mathrm{d}}\frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t} $$(19)where md is the mass of the droplet and Cp,d$ {C}_{\mathrm{p},\mathrm{d}}$ is the specific heat capacity of the droplet at constant pressure.

In practice, equation (19) is used to express the derivative dTddt$ \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}$. To do so, it is necessary to dispose of another expression for Q̇d$ {\dot{Q}}_{\mathrm{d}}$. This heat flux can be deduced from the Spalding heat transfer number of species i, BT,i$ {B}_{\mathrm{T},\mathrm{i}}$, the definition of which is [3]:BT,i=Cp,i,f(T-Td)Δvap,iH(Td)+Q̇d,iṁi,d$$ {B}_{\mathrm{T},\mathrm{i}}=\frac{{C}_{\mathrm{p},\mathrm{i},\mathrm{f}}\left({T}_{\mathrm{\infty }}-{T}_{\mathrm{d}}\right)}{{\Delta }_{\mathrm{vap},\mathrm{i}}H({T}_{\mathrm{d}})+\frac{{\dot{Q}}_{\mathrm{d},\mathrm{i}}}{{\dot{m}}_{\mathrm{i},\mathrm{d}}}} $$(20)where Q̇d,i$ {\dot{Q}}_{\mathrm{d},\mathrm{i}}$ is the contribution of species i to Q̇d$ {\dot{Q}}_{\mathrm{d}}$ (this latter is obtained by summing the Q̇d,i$ {\dot{Q}}_{\mathrm{d},\mathrm{i}}$ quantities over all the species); Cp,i,f$ {C}_{\mathrm{p},\mathrm{i},\mathrm{f}}$ is the mean specific heat capacity of species i in the gas film. Td$ {T}_{\mathrm{d}}$ and T$ {T}_{\mathrm{\infty }}$ are the droplet temperature and the bulk-gas temperature, respectively. Δvap,iH$ {\Delta }_{\mathrm{vap},\mathrm{i}}H$ is the specific enthalpy of vaporization for species i.

Therefore, the knowledge of BT,i$ {B}_{T,\mathrm{i}}$ makes it possible to estimate Q̇d$ {\dot{Q}}_{\mathrm{d}}$:Q̇d=i=1Nṁi,d[Cp,i,f(T-Td)BT,i-Δvap,iH(Td)].$$ {\dot{Q}}_{\mathrm{d}}=\sum_{i=1}^N {\dot{m}}_{\mathrm{i},\mathrm{d}}\left[\frac{{C}_{\mathrm{p},\mathrm{i},\mathrm{f}}\left({T}_{\mathrm{\infty }}-{T}_{\mathrm{d}}\right)}{{B}_{T,\mathrm{i}}}-{\Delta }_{\mathrm{vap},\mathrm{i}}H({T}_{\mathrm{d}})\right]. $$(21)

For estimating the Spalding heat transfer number of species i BT,i$ {B}_{\mathrm{T},\mathrm{i}}$, we refer to Brenn et al. [3] and Abramzon and Sirignano [5] who postulate that BT,i$ {B}_{\mathrm{T},\mathrm{i}}$ is coupled with the mass transfer number BM,i$ {B}_{\mathrm{M},\mathrm{i}}$ through:BT,i=(1+BM,i)ϕi-1withϕi=Cp,i,fCp,i,gShi*Nui*1Lei.$$ {B}_{\mathrm{T},\mathrm{i}}={\left(1+{B}_{\mathrm{M},\mathrm{i}}\right)}^{{\phi }_{\mathrm{i}}}-1\hspace{1em}\mathrm{with}\hspace{1em}{\phi }_{\mathrm{i}}=\frac{{C}_{\mathrm{p},\mathrm{i},\mathrm{f}}}{{C}_{\mathrm{p},\mathrm{i},\mathrm{g}}}\frac{{\mathrm{Sh}}_{\mathrm{i}}^{\mathrm{*}}}{{\mathrm{Nu}}_{\mathrm{i}}^{\mathrm{*}}}\frac{1}{{\mathrm{Le}}_{\mathrm{i}}}. $$(22)

Cp,i,g$ {C}_{\mathrm{p},\mathrm{i},\mathrm{g}}$ is the mean specific heat capacity of species i in the bulk gas. Lei is the Lewis number of species i defined [3, 5] by:Lei=λi,fρfDi,fCp,i,g,$$ {\mathrm{Le}}_{\mathrm{i}}=\frac{{\lambda }_{\mathrm{i},\mathrm{f}}}{{\rho }_{\mathrm{f}}{D}_{\mathrm{i},\mathrm{f}}{C}_{\mathrm{p},\mathrm{i},\mathrm{g}}}, $$(23)where λi,f is the thermal conduction coefficient of species i in the gas film.

In equation (22), Nui*$ {\mathrm{Nu}}_{\mathrm{i}}^{\mathrm{*}}$ 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:{Nui*=2+Nui,0-2F(BT,i)F(BT,i)=(1+BT,i)0.7ln(1+BT,i)BT,i.$$ \left\{\begin{array}{c}{{\mathrm{Nu}}_{\mathrm{i}}}^{\mathrm{*}}=2+\frac{{\mathrm{Nu}}_{\mathrm{i},0}-2}{F({B}_{T,\mathrm{i}})}\\ F({B}_{T,\mathrm{i}})={\left(1+{B}_{T,\mathrm{i}}\right)}^{0.7}\frac{\mathrm{ln}\left(1+{B}_{T,\mathrm{i}}\right)}{{B}_{T,\mathrm{i}}}\end{array}\right.. $$(24)

Nui,0$ \mathrm{N}{\mathrm{u}}_{\mathrm{i},0}$ is the Nusselt number of species i for the non-evaporating droplets, that can be computed [3, 6] using the following correlation especially developed for spheres from experimental data on mass-transfer rates:Nui,0=2+0.552Re1/2Pri1/3.$$ {\mathrm{Nu}}_{\mathrm{i},0}=2+0.552{\mathrm{Re}}^{1/2}{\mathrm{Pr}}_{\mathrm{i}}^{1/3}. $$(25)Pri is the Prandtl number of species i defined by the equation:Pri=Cp,i,fμi,fλi,f.$$ {\mathrm{Pr}}_{\mathrm{i}}=\frac{{C}_{\mathrm{p},\mathrm{i},\mathrm{f}}{\mu }_{\mathrm{i},\mathrm{f}}}{{\lambda }_{\mathrm{i},\mathrm{f}}}. $$(26)

Note that, BT,i$ {B}_{T,\mathrm{i}}$ should be solved from equations (22) to (24) by iterations; its initial value can be assumed to be equal to BM,i.

Finally, the rate of temperature change of the droplet dTddt$ \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}$ can be calculated by equating equations (19) and (21):dTddt=1mdCp,di=1Nṁi,d[Cp,i,f(T-Td)BT,i-Δvap,iH(Td)].$$ \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}=\frac{1}{{m}_{\mathrm{d}}{C}_{\mathrm{p},\mathrm{d}}}\sum_{i=1}^N {\dot{m}}_{\mathrm{i},\mathrm{d}}\left[\frac{{C}_{\mathrm{p},\mathrm{i},\mathrm{f}}\left({T}_{\mathrm{\infty }}-{T}_{\mathrm{d}}\right)}{{B}_{T,\mathrm{i}}}-{\Delta }_{\mathrm{vap},\mathrm{i}}H({T}_{\mathrm{d}})\right]. $$(27)

As illustrated in Figure 3, the heat flux exiting the film is simply given by:Q̇f=Q̇d+i=1Nṁi,dΔvap,iH.$$ {\dot{Q}}_{\mathrm{f}}={\dot{Q}}_{\mathrm{d}}+\sum_{i=1}^N {\dot{m}}_{\mathrm{i},\mathrm{d}}{\Delta }_{\mathrm{vap},\mathrm{i}}H. $$(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:Q̇g=Q̇f+(net advectiveflux)=Q̇f+i=1Nṁi,dCp,i,f(T-Td).$$ {\dot{Q}}_{\mathrm{g}}={\dot{Q}}_{\mathrm{f}}+\left(\mathrm{net}\enspace \mathrm{advectiveflux}\right)={\dot{Q}}_{\mathrm{f}}+\sum_{i=1}^N {\dot{m}}_{\mathrm{i},\mathrm{d}}{C}_{\mathrm{p},\mathrm{i},\mathrm{f}}\left({T}_{\mathrm{\infty }}-{T}_{\mathrm{d}}\right). $$(29)

Such an equation makes it possible to calculate Q̇g$ {\dot{Q}}_{\mathrm{g}}$ which is the heat flux from the bulk gas into the film. The average properties of the film (such as μi,f, Di,f, Cp,i,f, λi,f and ρi,f) are estimated at the average temperature of the film Tf, given [3, 5] by the equation:Tf=Td+T-Td3.$$ {T}_{\mathrm{f}}={T}_{\mathrm{d}}+\frac{{T}_{\mathrm{\infty }}-{T}_{\mathrm{d}}}{3}. $$(30)

It is worth noting that in the evaporation model ṁi,d$ {\dot{m}}_{\mathrm{i},\mathrm{d}}$ can be positive (which indicates evaporation from the droplet) or negative (which indicates condensation to the droplet). Heat fluxes (Q̇d$ {\dot{Q}}_{\mathrm{d}}$, Q̇f$ {\dot{Q}}_{\mathrm{f}}$ and Q̇g$ {\dot{Q}}_{\mathrm{g}}$) 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 NL layers of sprayed droplets in the tank, with each layer containing Nd droplets, the total mass transferred into the bulk gas can be computed with the following equation:ṁg=i=1Nṁi,g=i=1NL=1NL(Ndṁi,d)L.$$ {\dot{m}}_{\mathrm{g}}=\sum_{i=1}^N {\dot{m}}_{i,\mathrm{g}}=\sum_{i=1}^N \sum_{L=1}^{{N}_L} {\left({N}_{\mathrm{d}}{\dot{m}}_{\mathrm{i},\mathrm{d}}\right)}_{\mathrm{L}}. $$(31)

The total heat flux from the bulk gas into the droplets can be computed with the following equation:Q̇gtot=L=1NL(NdQ̇g)L.$$ {\dot{Q}}_{\mathrm{g}}^{\mathrm{tot}}=\sum_{L=1}^{{N}_L} {\left({N}_{\mathrm{d}}{\dot{Q}}_{\mathrm{g}}\right)}_{\mathrm{L}}. $$(32)

The temperature change rate of the bulk gas dTgdt$ \frac{\mathrm{d}{T}_{\mathrm{g}}}{\mathrm{d}t}$ in the tank can be calculated using the following equation:Q̇gtot=L=1NL(NdQ̇g)L=-mgCp,gdTgdt,$$ {\dot{Q}}_{\mathrm{g}}^{\mathrm{tot}}=\sum_{L=1}^{{N}_L} {\left({N}_{\mathrm{d}}{\dot{Q}}_{\mathrm{g}}\right)}_L=-{m}_{\mathrm{g}}{C}_{\mathrm{p},\mathrm{g}}\frac{\mathrm{d}{T}_{\mathrm{g}}}{\mathrm{d}t}, $$(33)where mg is the total mass of the gas phase and Cp,g$ {C}_{\mathrm{p},\mathrm{g}}$ 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:Nliqnew=Nliqold+Nliqin.$$ {N}_{\mathrm{liq}}^{\mathrm{new}}={N}_{\mathrm{liq}}^{\mathrm{old}}+{N}_{\mathrm{liq}}^{\mathrm{in}}. $$(34)Nliqold$ {N}_{\mathrm{liq}}^{\mathrm{old}}$ and Nliqnew$ {N}_{\mathrm{liq}}^{\mathrm{new}}$ are the total molar amounts of the liquid phase before and after mixing, respectively. Nliqin$ {N}_{\mathrm{liq}}^{\mathrm{in}}$ is the molar amount of the influx fluid. The new mole fraction of each species in the liquid can be computed with this equation:xi,liqnew=xi,liqoldNliqold+xi,liqinNliqin(xi,liqoldNliqold+xi,liqinNliqin),$$ {x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{new}}=\frac{{x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{old}}{N}_{\mathrm{liq}}^{\mathrm{old}}+{x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{in}}{N}_{\mathrm{liq}}^{\mathrm{in}}}{\sum \left({x}_{i,\mathrm{liq}}^{\mathrm{old}}{N}_{\mathrm{liq}}^{\mathrm{old}}+{x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{in}}{N}_{\mathrm{liq}}^{\mathrm{in}}\right)}, $$(35)xi,liqold$ {x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{old}}$ and xi,liqnew$ {x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{new}}$ are the mole fractions of species i in the liquid phase before and after mixing, respectively. While xi,liqin$ {x}_{\mathrm{i},\mathrm{liq}}^{\mathrm{in}}$ 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:NliqnewTliqoldTliqnewCp,liqdT+NliqinTliqinTliqoldCp,liqindT=0.$$ {N}_{\mathrm{liq}}^{\mathrm{new}}\underset{{T}_{\mathrm{liq}}^{\mathrm{old}}}{\overset{{T}_{\mathrm{liq}}^{\mathrm{new}}}{\int }} {C}_{\mathrm{p},\mathrm{liq}}\mathrm{d}T+{N}_{\mathrm{liq}}^{\mathrm{in}}\underset{{T}_{\mathrm{liq}}^{\mathrm{in}}}{\overset{{T}_{\mathrm{liq}}^{\mathrm{old}}}{\int }} {C}_{\mathrm{p},\mathrm{liq}}^{\mathrm{in}}\mathrm{d}T=0. $$(36)

Assuming that the heat capacity of the liquid phase Cp,liq$ {C}_{\mathrm{p},\mathrm{liq}}$ and the heat capacity of the influx liquid phase Cp,liqin$ {C}_{\mathrm{p},\mathrm{liq}}^{\mathrm{in}}$ are constant (i.e. temperature-independent), equation (36) can be converted to:NliqnewCp,liq(Tliqnew-Tliqold)+NliqinCp,liqin(Tliqold-Tliqin)=0.$$ {N}_{\mathrm{liq}}^{\mathrm{new}}{C}_{\mathrm{p},\mathrm{liq}}\left({T}_{\mathrm{liq}}^{\mathrm{new}}-{T}_{\mathrm{liq}}^{\mathrm{old}}\right)+{N}_{\mathrm{liq}}^{\mathrm{in}}{C}_{\mathrm{p},\mathrm{liq}}^{\mathrm{in}}\left({T}_{\mathrm{liq}}^{\mathrm{old}}-{T}_{\mathrm{liq}}^{\mathrm{in}}\right)=0. $$(37)

This equation makes it possible to estimate the new liquid temperature Tliqnew$ {T}_{\mathrm{liq}}^{\mathrm{new}}$.

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. Yi,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:xi,dφiliq(T,p,2x)=yi,sφivap(T,p,2y)(i=1,,N),$$ {x}_{\mathrm{i},\mathrm{d}}\cdot {\phi }_{\mathrm{i}}^{\mathrm{liq}}\left(T,p,2x\right)={y}_{\mathrm{i},s}\cdot {\phi }_{\mathrm{i}}^{\mathrm{vap}}\left(T,p,2y\right)\hspace{1em}\left(\mathrm{i}=1,\cdots,N\right), $$(38)where xi,d is the mole fraction of component i in the droplet, and yi,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 φiliq$ {\phi }_{\mathrm{i}}^{\mathrm{liq}}$ and in the vapour phase φivap$ {\phi }_{\mathrm{i}}^{\mathrm{vap}}$ are calculated with the PR-EoS, the mathematical expression of which is:P=RTv-bm-amv(v+bm)+bm(v-bm).$$ P=\frac{{RT}}{v-{b}_{\mathrm{m}}}-\frac{{a}_{\mathrm{m}}}{v(v+{b}_{\mathrm{m}})+{b}_{\mathrm{m}}(v-{b}_{\mathrm{m}})}. $$(39)

In this paper, the PR-EoS was not volume-translated [911] and classical van der Waals mixing rules [12] were used:{am(T,z)=i=1Nj=1Nzizjai(T)aj(T)[1-kij(T)]bm(z)=i=1Nzibi.$$ \left\{\begin{array}{c}{a}_{\mathrm{m}}(T,z)=\sum_{i=1}^N \sum_{j=1}^N {z}_{\mathrm{i}}{z}_j\sqrt{{a}_i(T)\cdot {a}_j(T)}\left[1-{k}_{{ij}}(T)\right]\\ {b}_{\mathrm{m}}(z)=\sum_{i=1}^N {z}_i{b}_i\end{array}\right.. $$(40)

It is today acknowledged that the Soave [13] α-function is non-consistent [1416] 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:{R=8.314472 mol-1 K-1X=[1+34-22+34+22]-10.253076587bi=ΩbRTc,iPc,iwith Ωb=XX+30.0777960739ai(T)=ac,iαi(T)with{ac,i=ΩaR2Tc,i2Pc,i, and Ωa=8(5X+1)49-37X0.457235529αi(T)=[1+mi(1-TTc,i)]2if ωi0.491mi=0.37464+1.54226ωi-0.26992ωi2if ωi>0.491mi=0.379642+1.48503ωi-0.164423ωi2+0.016666ωi3.$$ \left\{\begin{array}{c}R=8.314472\mathrm{\enspace J\enspace mo}{\mathrm{l}}^{-1}\mathrm{\enspace }{\mathrm{K}}^{-1}\\ X={\left[1+3\sqrt{4-2\sqrt{2}}+3\sqrt{4+2\sqrt{2}}\right]}^{-1}\approx 0.253076587\\ \begin{array}{c}{b}_{\mathrm{i}}={\mathrm{\Omega }}_b\frac{R{T}_{\mathrm{c},\mathrm{i}}}{{P}_{\mathrm{c},\mathrm{i}}}\hspace{1em}\mathrm{with}\enspace {\mathrm{\Omega }}_{\mathrm{b}}=\frac{X}{X+3}\approx 0.0777960739\\ {a}_{\mathrm{i}}(T)={a}_{\mathrm{c},\mathrm{i}}{\alpha }_{\mathrm{i}}(T)\hspace{1em}\mathrm{with}\left\{\begin{array}{c}{a}_{\mathrm{c},\mathrm{i}}={\mathrm{\Omega }}_{\mathrm{a}}\frac{{R}^2{T}_{\mathrm{c},\mathrm{i}}^2}{{P}_{\mathrm{c},\mathrm{i}}},\enspace \mathrm{and}\enspace {\mathrm{\Omega }}_{\mathrm{a}}=\frac{8\left(5X+1\right)}{49-37X}\approx 0.457235529\\ {\alpha }_{\mathrm{i}}(T)={\left[1+{m}_{\mathrm{i}}\left(1-\sqrt{\frac{T}{{T}_{\mathrm{c},\mathrm{i}}}}\right)\right]}^2\end{array}\right.\\ \begin{array}{c}\mathrm{if}\enspace {\omega }_{\mathrm{i}}\le 0.491\hspace{1em}{m}_{\mathrm{i}}=0.37464+1.54226{\omega }_{\mathrm{i}}-0.26992{\omega }_{\mathrm{i}}^2\\ \mathrm{if}{\enspace {\omega }}_{\mathrm{i}}>0.491\hspace{1em}{m}_{\mathrm{i}}=0.379642+1.48503{\omega }_{\mathrm{i}}-0.164423{\omega }_{\mathrm{i}}^2+0.016666{\omega }_{\mathrm{i}}^3\end{array}\end{array}\end{array}\right.. $$(41)

The fugacity coefficient can be calculated [12] from the equation:lnφi=bibm(PvRT-1)-ln[PRT(v-bm)]+am22 RTbm(δi-bibm)ln[v+(1-2)bmv+(1+2)bm],$$ \mathrm{ln}{\phi }_{\mathrm{i}}=\frac{{b}_{\mathrm{i}}}{{b}_{\mathrm{m}}}\left(\frac{{Pv}}{{RT}}-1\right)-\mathrm{ln}\left[\frac{P}{{RT}}\left(v-{b}_{\mathrm{m}}\right)\right]+\frac{{a}_{\mathrm{m}}}{2\sqrt{2}\enspace {RT}{b}_{\mathrm{m}}}\left({\delta }_{\mathrm{i}}-\frac{{b}_{\mathrm{i}}}{{b}_{\mathrm{m}}}\right)\cdot \mathrm{ln}\left[\frac{v+(1-\sqrt{2}){b}_{\mathrm{m}}}{v+(1+\sqrt{2}){b}_{\mathrm{m}}}\right], $$(42)whereδi=2j=1Nczjaiaj(1-kij)am.$$ {\delta }_{\mathrm{i}}=\frac{2\sum_{j=1}^{{N}_{\mathrm{c}}} {z}_j\sqrt{{a}_i{a}_j}\left(1-{k}_{{ij}}\right)}{{a}_{\mathrm{m}}}. $$(43)

In the above equations, P is the pressure, R is the gas constant, T is the temperature, ai and bi are the cohesive parameter and molar co-volume of pure component i, v is the molar volume, zi is the mole fraction of component i, Tc,i is the experimental critical temperature, Pc,i is the experimental critical pressure, and ωi is the experimental acentric factor of pure i. In this paper, the binary interaction parameters kij were predicted using the E-PPR78 model [1719].

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.

  1. First we load the earlier properties of the droplet, i.e. velocity (ud), angle between directions of motion and gravity (θd), height from the liquid phase surface (hd), radius (rd), temperature (Td), composition (xi,d) and amount (nd). 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.

  2. Then we compute the droplet’s new motion properties (ud and θd) with the equations:

{udnew=(ud,xnew)2+(ud,ynew)2θdnew=arctan(ud,xnewud,ynew)ud,xnew=ud,x+dud,xdtΔtud,ynew=ud,y+dud,ydtΔt,$$ \left\{\begin{array}{c}{u}_{\mathrm{d}}^{\mathrm{new}}=\sqrt{{\left({u}_{\mathrm{d},x}^{\mathrm{new}}\right)}^2+{\left({u}_{\mathrm{d},y}^{\mathrm{new}}\right)}^2}\\ {\theta }_{\mathrm{d}}^{\mathrm{new}}=\mathrm{arctan}\left(\frac{{u}_{\mathrm{d},x}^{\mathrm{new}}}{{u}_{\mathrm{d},y}^{\mathrm{new}}}\right)\\ \begin{array}{c}{u}_{\mathrm{d},x}^{\mathrm{new}}={u}_{\mathrm{d},x}+\frac{\mathrm{d}{u}_{\mathrm{d},x}}{\mathrm{d}t}\Delta t\\ {u}_{\mathrm{d},y}^{\mathrm{new}}={u}_{\mathrm{d},y}+\frac{\mathrm{d}{u}_{\mathrm{d},y}}{\mathrm{d}t}\Delta t\end{array}\end{array}\right., $$(44)and a new height of droplet (hd) from liquid surface with the equation:hdnew=hdold-ud,yold+ud,ynew2Δt,$$ {h}_{\mathrm{d}}^{\mathrm{new}}={h}_{\mathrm{d}}^{\mathrm{old}}-\frac{{u}_{\mathrm{d},y}^{\mathrm{old}}+{u}_{\mathrm{d},y}^{\mathrm{new}}}{2}\Delta t, $$(45)where dud,xdt$ \frac{\mathrm{d}{u}_{\mathrm{d},x}}{\mathrm{d}t}$ and dud,ydt$ \frac{\mathrm{d}{u}_{\mathrm{d},y}}{\mathrm{d}t}$ are calculated from equations (1) and (2).
  1. Next we compute the evaporation mass flow rates of each component (ṁi,d$ {\dot{m}}_{\mathrm{i},\mathrm{d}}$) of the droplet with equation (12) (using Brenn’s method [3]).

  2. If the droplet has vaporized completely, we compute the heat flux Q̇g$ {\dot{Q}}_{\mathrm{g}}$ from the bulk gas and then skip the next steps.

  3. Now we compute the heat fluxes (Q̇d$ {\dot{Q}}_{\mathrm{d}}$, Q̇f$ {\dot{Q}}_{\mathrm{f}}$ and Q̇g$ {\dot{Q}}_{\mathrm{g}}$) between the droplets, the gas film and the bulk gas with equations (21), (28) and (29). The quantity dTddt$ \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}$ is obtained from equation (27) while dTddt$ \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}$ is deduced from equation (33).

  4. Then we compute the new mole amount (nd) and composition (xi,d) of the droplet with the equations:

{ni,dnew=ndoldxi,dold-ṁi,dMWiΔtndnew=i=1Nni,dnewxi,dnew=ni,dnewndnew.$$ \left\{\begin{array}{c}{n}_{\mathrm{i},\mathrm{d}}^{\mathrm{new}}={n}_{\mathrm{d}}^{\mathrm{old}}{x}_{\mathrm{i},\mathrm{d}}^{\mathrm{old}}-\frac{{\dot{m}}_{\mathrm{i},\mathrm{d}}}{M{W}_{\mathrm{i}}}\Delta t\\ {n}_{\mathrm{d}}^{\mathrm{new}}=\sum_{\mathrm{i}=1}^N {n}_{\mathrm{i},\mathrm{d}}^{\mathrm{new}}\\ {x}_{\mathrm{i},\mathrm{d}}^{\mathrm{new}}=\frac{{n}_{\mathrm{i},\mathrm{d}}^{\mathrm{new}}}{{n}_{\mathrm{d}}^{\mathrm{new}}}\end{array}.\right. $$(46)
  1. Now we compute the new temperature (Td) of the droplet using equation (27) and the following equation:

Tdnew=Tdold+dTddtΔt$$ {T}_{\mathrm{d}}^{\mathrm{new}}={T}_{\mathrm{d}}^{\mathrm{old}}+\frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}\Delta t $$(47)
  1. Now we compute the new molar volume (vm,d) of the droplet using the PR-EoS with the new temperature (Td) and composition (xi,d), then compute the new radius (rd) of the droplet with the equation:

rdnew=12(6ndnewvm,dnewπ)1/3.$$ {r}_{\mathrm{d}}^{\mathrm{new}}=\frac{1}{2}{\left(\frac{6{n}_{\mathrm{d}}^{\mathrm{new}}{v}_{\mathrm{m},\mathrm{d}}^{\mathrm{new}}}{\pi }\right)}^{1/3}. $$(48)

In each simulation time step (Δt), a single droplet’s properties (i.e. velocity (ud), angle between directions of motion and gravity (θd), height from the liquid phase surface (hd), radius (rd), temperature (Td), composition (xi,d) and amount (nd)) 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 (ṁi,d$ {\dot{m}}_{\mathrm{i},\mathrm{d}}$) and heat flux Q̇g$ {\dot{Q}}_{\mathrm{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.

thumbnailthumbnail 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.

  1. Compute the new spray droplet properties layer by layer, following the instructions in Section 4.1.

  2. If the droplet layer reaches the liquid phase, update the liquid phase properties (Tliq, xi,liq, and nliq) with equations (34), (35), and (37). Once done, the volume of the liquid phase (Vliq) can be calculated with the help of the PR-EoS.

  3. Compute the new mass (mg) and compositions (yi,g) of the gas phase using equation (31) and the following equations:

{mi,gnew=mi,g+ṁi,gΔtmgnew=iNmi,gnew,$$ \left\{\begin{array}{c}{m}_{\mathrm{i},\mathrm{g}}^{\mathrm{new}}={m}_{\mathrm{i},\mathrm{g}}+{\dot{m}}_{\mathrm{i},\mathrm{g}}\Delta t\\ {m}_{\mathrm{g}}^{\mathrm{new}}=\sum_{\mathrm{i}}^N {m}_{\mathrm{i},\mathrm{g}}^{\mathrm{new}}\end{array}\right., $$(49)yi,gnew=mi,gnew/MWik=1N(mk,gnew/MWk).$$ {y}_{\mathrm{i},\mathrm{g}}^{\mathrm{new}}=\frac{{m}_{\mathrm{i},\mathrm{g}}^{\mathrm{new}}/M{W}_{\mathrm{i}}}{\sum_{k=1}^N \left({m}_{k,\mathrm{g}}^{\mathrm{new}}/M{W}_k\right)}. $$(50)
  1. Compute the new temperature of the gas phase using equation (33) with mgnew$ {m}_{\mathrm{g}}^{\mathrm{new}}$ and the following equation:

Tgnew=Tg+dTgdtΔt.$$ {T}_{\mathrm{g}}^{\mathrm{new}}={T}_{\mathrm{g}}+\frac{\mathrm{d}{T}_{\mathrm{g}}}{\mathrm{d}t}\Delta {t}. $$(51)
  1. Compute the new volume of the gas phase as follows:

Vg=Vtank-Vliq-L=1NL(Ndndvm,d)L,$$ {V}_{\mathrm{g}}={V}_{\mathrm{tank}}-{V}_{\mathrm{liq}}-\sum_{L=1}^{{N}_L} {\left({N}_{\mathrm{d}}{n}_{\mathrm{d}}{v}_{\mathrm{m},\mathrm{d}}\right)}_L, $$(52)where Vtank is the tank total volume.
  1. Calculate the gas pressure using the PR-EoS with the new temperature, volume and composition values.

  2. 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.

thumbnail 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:

  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 (Δtlayer), 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.

thumbnail 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.

thumbnail 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.

thumbnail 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.

thumbnail 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.

thumbnail 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 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.

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, multi-component 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.

References

All Figures

thumbnail Fig. 1

Schematic of the LNG spraying process.

In the text
thumbnail Fig. 2

The force field of a droplet.

In the text
thumbnail Fig. 3

Schematic of a droplet enveloped by a gaseous film.

In the text
thumbnailthumbnail Fig. 4

Flowchart detailing the modelling of a droplet layer and named LayerSim.

In the text
thumbnail Fig. 5

Flowchart showing the modelling of one tank and named TankSim.

In the text
thumbnail Fig. 6

Flowchart of the master program.

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 10

Simulated gas phase temperature and pressure versus spraying time for an example case.

In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.