Numerical simulation of natural convection and boil-off in a small size pressurized LNG storage tank

Numerical simulation of natural convection and boil-off in a small size pressurized LNG storage tank J.L. Ferr´ına,b,∗ , L.J. P´erez-P´ereza,c aUniversidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain bInstituto Tecnol´ogico de Matem´atica Industrial, 15782 Santiago de Compostela, Spain cInstituto de Matem´aticas, 15782 Santiago de Compostela, Spain Abstract A numerical simulation of the flow of LNG stored in a small-sized cylindrical tank is presented. The main scope of this work is to characterize the heat ingress to the tank as well as the boil-off rate, depending on the filling level and the insulation layer thickness. The tank is assumed to be in a state of inactivity, such that the fluid phases are initially quiescent and are not released to the exterior. The proposed mathematical model consists of a conjugate heat transfer problem coupled with the SST K−ωturbulence model for the fluid phases. The free surface separating liquid and vapor in the tank is tracked using the Volume of Fluid method (VOF). The model is solved using the software ANSYS Fluent. It is shown that the filling level of the tank substantially influences the boiling rate and the degree of stratification, as well as the flow structures generated by free convection. The relation among the insulation thickness and total heat leak is established, showing that the obtained increased heat ingress due to lower insulation thickness values leads to a rise in pressurization and boil-off rates. Keywords: Numerical Simulation, LNG, boil-off, natural convection, conjugate heat transfer. 1. Introduction Natural gas is a fossil fuel mainly composed of methane, ethane and propane, with traces of heavier hydrocarbons and also a varying amount of nitrogen. Its composition is highly dependent on its origin. Usually, methane content is above 90%. Natural gas is currently the fastest-growing non-renewable source of energy [1] and is also expected5 to continue growing in the coming years, with a projected growth rate of 1.6% (see [2]). This is principally due to its cost being competitive and comparatively low compared to the historical one. Furthermore, it is cleaner than the remaining fossil alternatives, as its combustion produces substantially less carbon dioxide and also reduces the emissions of other highly contaminant gases, such as sulfur dioxide.10 ∗Corresponding author Email addresses: [email protected] (J.L. Ferr´ın), [email protected] (L.J. P´erez-P´erez) Preprint submitted to Computers and Chemical Engineering April 6, 2020 Nomenclature Greek symbols αVolume fraction βThermal expansion coefficient (K−1) ∆tTime step size (s) γEvaporation rate coefficient (s−1) µDynamic viscosity (Pa·s) µTTurbulent viscosity (Pa·s) ωSpecific dissipation rate (s−1) ρDensity (kg/m3) Latin symbols AlWall surface area in contact with liquid (m2) cpSpecific heat, at constant pressure (J/(kg·K)) ETotal specific energy (J/kg) Hvap Latent heat of vaporization (J/kg) hHeat transfer coefficient (W/(m2·K)) gGravitational acceleration (m/s2) MMolar mass (kg/kmol) mlLiquid mass (kg) ˙mlv Mass transfer source term (kg/(m3·s)) nOutward-pointing unit normal vector (m) KTurbulent kinetic energy (m2/s2) kThermal conductivity (W/(m·K)) kTTurbulent thermal conductivity (W/(m·K)) pPressure (Pa) qlAverage heat flux, liquid phase (W/m2) QTotal heat flux (W) TTemperature (K) rRadial coordinate (m) RIdeal gas constant (J/(mol·K)) TRReference temperature (K) TRb Reference temperature (Boussinesq model) (K) tTime (s) vSpecific volume (m3/kg) VmMolar volume (m3/mol) vVelocity (m/s) wInsulation layer thickness (mm) yVertical coordinate (m) Subscripts cCritical property fFluid lLiquid phase pPerlite sSolid sat Saturation value vVapor phase 2 In this context, Liquefied Natural Gas (LNG) becomes especially relevant, as natural gas may be stored by decreasing its temperature below its boiling point, approximately 111 K at standard atmospheric pressure. In this state, it is much denser and requires about 600 times less volume for storage. LNG enables profitable transportation by sea in LNG carriers or by road in trucks, eliminating the need for a costly gas network15 and, therefore, helping to diversify gas suppliers. Albeit its market is still immature, small scale storage of LNG has seen a rise and is projected to continue expanding in the coming years, as it may become a key commodity in different sectors, such as fuel for road transportation or to generate off-grid power (see [3]). Due to the cryogenic conditions required for its storage, heat ingress to LNG tanks is20 an unavoidable problem. This turns to be a major drawback, as it results in the heating of the stored LNG, which vaporizes as it reaches its boiling point following the process known as boil-off. The boil-off in the tank also triggers the auto-refrigeration process, which keeps LNG cold as it boils. However, as tank pressure rises (self-pressurization), the boiling point of LNG also increases. Thus, if tank pressure is not relieved, LNG25 starts to heat up. In small-sized LNG storage, the most commonly used tanks are type C pressure vessels, which can be designed to operate at high pressure values. This allows to hold LNG in a quiescent state for several days but, inevitably, if LNG inside the tank is not used or reliquefied, the generated boil-off gas (BOG) has to be vented out to maintain an acceptable tank pressure, with the consequent economic losses and30 environmental impact, as methane is one of the main greenhouse effect gases. To the best knowledge of the authors, even though there is a significant body of research performing CFD analyses of the boil-off process, there are not many studies specifically dealing with tank self-pressurization and temperature stratification in small and quiescent tanks. For LNG tankers, boil-off rate (BOR) can be estimated using35 empirical relations, as detailed in [4], where a mathematical expression relating tank LNG volume and heat ingress with the BOR was presented. This expression leads to reasonable estimations of BOG production during a ship voyage where pressure, tank volume and fluid properties are approximately constant. However, it fails to account for a number of phenomena that become much more relevant in smaller sized tanks, such as40 the effect that temperature and pressure build-up have on the boil-off rate. Other approaches rely on using state equations to relate thermodynamic variables governing the behavior of the fluid phases within the tank. In [5], models using Lee- Kesler-Plocker and Benedict-Webb-Rubin state equations were employed to predict BOG production in a tank assuming that the vapor and liquid phase were under thermody-45 namic equilibrium. Heat ingress was supplied as data, and the models were used in different tanks with varying LNG volume and operating pressures. In [6], a system of ODEs derived from global mass and enthalpy balances and the Soave-Redlich-Kwong (SRK) equation of state was proposed to model the boil-off phenomena, considering the different components of LNG. The model was able to predict50 the evolution of their concentrations separately. Consequently, the weathering process of LNG could be assessed, as its composition changes over time, with more volatile components like nitrogen boiling first. Moreover, in [7], a detailed study aiming to characterize LNG losses due to boil-off in all the steps of an LNG tanker voyage was performed, including the loading and unloading process of the tanks and the BOG generation during55 laden and ballast voyages, considering several parameters. One drawback with models relying on ODEs to model BOG generation in storage 3 tanks is that they may not completely account for the characteristics of the LNG flow inside the tank, driven by natural convection. In particular, the flow structures may depend on heat ingress, tank geometry or liquid level. Moreover, heat ingress is usually60 modeled in a simplified manner. Frequently, thermal equilibrium inside the tank is also assumed. Nevertheless, the actual conditions within small storage tanks under selfpressurization can substantially deviate from it. Thermal stratification in the liquid tends to appear as a warm layer at the free surface develops due to buoyancy (see [8]). In [9], the effect of the thermal aspect ratio on thermal stratification and self-pressurization was65 experimentally investigated in a liquid nitrogen storage tank. Their measurements show significant vapor superheating as well as stratification inside the liquid phase, which are shown to be dependent on the filling level of the tank. Comparable findings were obtained in [10] for a similar LN2tank. In addition, a thermal diffusion model was proposed to predict the self-pressurization curves, accounting for thermal stratification.70 To investigate the details of LNG flow inside tanks, PDE models that consider Navier- Stokes equations can be used, which are much more computationally expensive to solve numerically. Related to this kind of approach, in [11], the flow field inside a rectangular water tank subjected to lateral heating was visualized to study thermal stratification in cryogenic tanks. Then, the axisymmetric LH2flow in a cylindrical cavity considering75 only the liquid with no vaporization was solved. In [12], the flow inside a 2d cut of a membrane tank was solved, considering a two-phase flow. The boiling process was modeled using a mass transfer source term, assuming a constant saturation temperature. In [13], a detailed study of natural convection in an LNG tank was presented, aiming to investigate how tank gradual pressurization influences BOG generation. The proposed80 model solves the flow only for the liquid phase, given that the vapor is assumed to be at the same temperature as the liquid-vapor interface. It was concluded that the BOG generation rate is reduced as both pressure and tank size rises. In [14], a two-phase lumped vapor CFD model was used to investigate self-pressurization in a LH2storage tank for different filling levels, obtaining favorable agreement with the experimental85 observations from [15]. In [16], the flow due to free convection in a liquefied hydrogen tank was studied numerically using a 2d axisymmetric model. Boil-off was modeled using an expression based on the difference between saturation temperature and the temperature at the free surface, which was considered constant. The obtained stratification profiles for90 the temperature were compared for different tank aspect ratios. In [17], stratification was studied for a water tank under large incoming heat flux, replicating the conditions and observations made by [18], including a de-stratification process. A similar boiloff model was used, accounting for heat loss in the liquid and considering a pressuredependent saturation temperature. The proposed model was enhanced in [19], including95 the conjugate heat transfer problem among the solid wall of the tank and the fluid phases. CFD simulations were also used in [20] to investigate the diffusion and composition changes of leaked LNG through the insulation box of large storage tanks. In [21], a mixture model was used to study convection and self-pressurization in both a LN2and a LNG tanks, considering uniform heat flux values at the inner surface wall and different100 filling levels. Their results also show how thermal stratification is developed in the early stages of the considered cases. Moreover, in [22], a CFD simulation of a full-scale storage tank was presented, assessing BOG generation and flow dynamics. An estimation of the wall superheat needed for the inception of nucleate boiling was obtained, finding that it 4 (a) 3d view. Liquid Vapor Axis Solid wall 1.5 m 0.3315 m r y (b) Schematic of a cross section of the tank. Figure 1: Tank geometry. is around 2.5-2.8 K.105 At the present work, we aim to study the effect of pressurization in the boil-off process for a simplified prototype of a small-sized LNG storage tank using a CFD multiphase model. Different filling levels of the tank are considered, as well as the conjugate heat transfer problem involving the tank solid wall. Realistic data for material properties are used. The influence of the insulation layer thickness on heat ingress and flow structures is110 also investigated. In Section 2, we describe the physical problem, starting with the tank geometry. Then, the considered values for material properties and the complete mathematical model are presented. In Section 3, the numerical schemes, the computational grid and the procedure that are used to obtain the solution are described. In Section 4, the obtained results are discussed, both for varying filling levels and different insulation115 thickness values. Finally, in Section 5, some concluding remarks are addressed. 2. Mathematical model and problem statement 2.1. Problem description A 0.5 m3cylindrical LNG storage tank is considered, depicted in Figure 1a. The inner chamber, which holds the fluids, has a height of 1.5 m and a maximum interior radius of120 0.3315 m, as shown in Figure 1b. It is assumed to be fixed to the ground and completely quiescent. Thus, any fluid movement is triggered by heat ingress to the tank, which is responsible for both free convection and boil-off. The tank solid wall is composed of a wmm-thick evacuated perlite insulation layer and two 304 stainless steel vessels, each of them 3 mm-thick. In Figure 2, a schematic diagram with the different solid layers is125 displayed. For the filling level influence studies presented in Section 4.1, the considered value for the insulation thickness parameter is w= 11 mm. In Section 4.2, the influence of this parameter in the obtained results is investigated. 5 Inner steel vessel Outer steel vessel Evacuated perlite 3 mm 3 mm w mm Figure 2: Layers inside the solid wall. Given that the tank has a cylindrical structure, we assume that the flow also preserves this symmetry. Therefore, an axisymmetric problem is considered, such that the model130 is solved only in a 2d cut of the tank, as shown in Figure 1b. The tank is considered to be initially filled with LNG at constant temperature T= 113 K. Provided that the corresponding saturation pressure for this temperature is 1.13 bar, we choose this value as the initial pressure in the vapor phase. 2.2. Material properties135 The composition of the LNG is assumed to be pure methane, which is a frequent assumption in similar studies, such as [13] or [22]. Given that LNG has a content of methane that is typically over 90%, its overall physical properties are very similar to those of liquefied methane. Consequently, the assumption of pure methane instead of a more realistic composition does not significantly change the flow structures, temperature140 stratification profiles or heat ingress to the tank. Nonetheless, it is worth mentioning that there is a number of shortcomings derived from this simplifying assumption. Namely, rollover and LNG weathering phenomena cannot be predicted by the model, as they are related to the multicomponent nature of LNG. For more details on rollover, we refer to [23]. In addition, the boil-off rate can be substantially different for a more realistic145 LNG composition. This is particularly relevant when a significant amount of nitrogen is present in the mixture, as it is a more volatile fluid than methane (see [6] or [7]). At the liquid phase inside the tank, temperature variations are small. Thus, constant values for its properties are considered, as shown in Table 1. We use data corresponding to saturated liquid methane at 113 K, extracted from the National Institute of Standards150 and Technology (NIST) database, available in [24]. With ρ,k,cpand µ, we denote density, thermal conductivity, specific heat and dynamic viscosity, respectively. Note that although density may be considered constant, buoyancy forces must still be included in order to account for free convection in the liquid phase. For this purpose, we use the Boussinesq approximation. Hence, the body force at the liquid blis represented by:155 bl=ρl(1 −β(T−TRb))g,(1) where β= 3.54e−3 K−1denotes the thermal expansion coefficient and TRb is the reference temperature used for the approximation, which we consider equal to the initial temperature T0= 113 K. In addition, gis the gravitational acceleration vector. 6 Both latent heat of vaporization Hvap and saturation temperature Tsat are taken to be functions of pressure using the following expressions:160 Tsat(p) = 65.3180 + 2.3371p0.2592,(2) Hvap(p) = 5.4465e5−36.5741p0.5919,(3) which have been obtained fitting the data available in [24] for liquid methane. At the vapor phase, the specific heat is considered constant with value cp,v = 2230 J/(kg·K). Moreover, as a significant degree of superheating may be reached, large temperature variations can occur. Thus, thermal conductivity and dynamic viscosity are165 considered as linear functions of the temperature fitting the data available in [24]: µv(T)=3.610e−8T+ 5.639e−7,(4) kv(T) = 1.208e−4T−1.775e−3.(5) Provided that one of the main interests of this work is to study the impact of tank selfpressurization on boil-off rate, density at the vapor phase is considered as a dependent variable of temperature and pressure using the Peng-Robinson equation of state:170 p=RT Vm−b−aα V2 m−2bVm−b2,(6) where a=0.45724R2T2 c pc , b =0.7780RTc pc , α = 1 + κ 1−T Tc0.5!!2 ,(7) κ= 0.37464 + 1.54226ωac −0.26992ω2 ac.(8) Values with csubindices are the critical properties of the fluid, i.e., the pressure, temperature and density at the critical point. Moreover, ωac = 0.0114 is the accentric factor whereas Vmis the molar volume in m3/mol and Ris the ideal gas constant. The choice of175 the Peng-Robinson equation of state over a more traditional expression as the ideal gas law is motivated by the fact that the latter yields a worse approximation of the pressure values for the cryogenic temperatures inside the tank. Insulating properties of evacuated perlite are heavily dependent on temperature values. Thus, to model heat transfer inside the solid wall of the storage tank, perlite thermal180 conductivity is considered as a function of temperature using the following expression, gathered from [25]: kp(T) = 1.911e−4+3.476e−12T3.6783.(9) Both vessels of steel are built with the same material. However, the outer one, which is exposed to the environment, is approximately at ambient temperature (293 K). In contrast, the inner vessel is roughly at the temperature of the fluids in the tank, around185 113 K. Therefore, we consider constant properties for the 304 stainless steel corresponding to those temperatures for the different vessels, using the data available in [26], as shown in Table 1. In addition, the used values and functions for the fluids and perlite are summarized. Note that the vapor density is computed using the equation of state (6). 7 ρ[kg/m3]k[W/(m·K)] cp[J/(kg·K)] µ[Pa·s] Liquid phase 420.40 1.82e−1 3490 1.13e−4 Vapor phase - kv(T) 2230 µv(T) Perlite 50 kp(T) 837 - Steel, inner vessel 8030 9.8 304.7 - Steel, outer vessel 8030 15.1 470.5 - Table 1: Physical properties for the liquid, the vapor and the solid materials. Source: [24, 26, 27]. 2.3. Turbulence model190 To assess if turbulence is of relevance in the flow, a dimensional analysis may be performed. A characteristic length defined as the liquid volume divided by the wall surface in contact with the liquid is considered (L= 0.15 m). This choice is motivated by the work [28], where numerical studies that analyzed a cylindrical cavity subjected to natural convection for varying Rayleigh numbers were conducted. For an estimated heat195 flux of q= 20 W/m2, the Grashof number is equal to 2.67e10 and the Rayleigh number is 5.79e10. According to the criteria developed in [28], these numbers are indicative of a turbulent regime. Thus, in order to model the turbulent flow, we use the SST K−ω RANS turbulence model (see [29]), which uses a combination between the K−εin the free stream and a K−ωnear the boundary layers. This allows to retain the good behavior200 of the K−ωmodel close to the walls, which is reported to perform better than K−ε models in wall-bounded flows, as it can be extended to fully resolve the boundary layer without the use of damping functions (see [30]). This is a desirable feature to accurately model the heat transfer at the walls of the tank. The details concerning these models are extensively covered in the literature (see, for instance, [31] or [32]).205 2.4. Multiphase flow The flow in the tank is a two phase flow, with two distinct and immiscible phases separated by a free surface. To model these different phases, the Volume of Fluid method (VOF) is used, which is a front-capturing method proposed originally in [33] to solve free surface flows where a single fluid domain with fixed boundaries is considered. The210 computational domain is split in a fluid subdomain Ωfand the solid wall Ωs. Then, the fluid subdomain can be decomposed as Ωf= Ωl(t)∪Ωv(t)∪Γint(t), where Ωl(t) is the liquid phase subdomain, Ωv(t) is the vapor phase subdomain and Γint(t) is the free surface between both phases, which all depend on time t. A piecewise constant function ξlthat is discontinuous at Γint(t) is used to track the free surface position:215 ξl(x, t) = 1,x∈Ωl(t), 0,x6∈ Ωl(t).(10) This function is numerically approximated using the so-called volume fraction αl, which takes a constant value at each cell in the grid. The value of αlis interpreted as the fraction of it occupied by the liquid phase. Then, the vapor phase volume fraction is readily obtained, as it must satisfy that αv= 1 −αl. Values other than 0 or 1 for the volume fraction represent a mixture of the fluid phases and indicate that the free surface220 is located inside the corresponding cell. 8 With the VOF method, only one set of conservation equations is solved for both fluid phases, using mixture properties. For an arbitrary property field φit is considered that ˆ φ=αlφl+αvφv=αlφl+ (1 −αl)φv,(11) where φland φvdenote the values corresponding to the liquid and vapor phases, respectively.225 Therefore, all relevant fluid properties, such as density ρ, viscosity µ, specific heat cpor thermal conductivity k, are considered as mixture properties when solving the conservation equations. To ensure the accuracy of the method, the cells containing the vapor-liquid free surface should be sufficiently small. 2.5. Boil-off model230 To model the boiling process of the liquid phase, appropriate source terms are added to the conservation equations. We use an analogous model to the one described in [16, 17, 21, 22], known as Lee model. It uses the following source term to model mass transfer from the liquid to the vapor: ˙mlv =           γαlρl T−Tsat(p) Tsat(p),if T > Tsat(p), γαvρv T−Tsat(p) Tsat(p),if T≤Tsat(p), (12) where γis a coefficient related to the vaporization rate, with units s−1. The saturation235 temperature Tsat is obtained as a function of pressure, using (2). According to expression (12), temperature values below the saturation one result in condensation of the vapor whereas temperature values above saturation lead to the vaporization of the liquid. In the proposed case, heat ingress to the tank causes a temperature rise, which triggers the vaporization process. Condensation is negligible, unless the pressure rises rapidly, which240 can occur in real conditions when the available volume for vapor reduces significantly due to thermal expansion of the liquid. In our computations, the term due to condensation in (12) is not included in the model. The expected heat flux to the tank is around 20 W/m2. This heat flux value is well below the required amount to produce nucleate boiling in the liquid phase. Thus, the245 phase change should take place at the free surface separating the fluid phases. Consequently, the source term (12) is only applied in computational cells where the free surface is present, this is, those cells satisfying that αl6= 1 and αl6= 0. The latent heat of vaporization is absorbed as the liquid evaporates. To model this phenomenon, we introduce the following source term in the energy conservation equation:250 S=−Hvap(p) ˙mlv.(13) 2.6. Governing equations The flow inside the tank is modeled using the unsteady Reynolds-Averaged Navier- Stokes equations, with the VOF method to account for the multiphase flow. As described in Section 2.4, this implies that mixture properties are used and that only one set of255 9 during which this behavior is exhibited is short. 4. Results390 The mathematical model described in Section 2 is solved until reaching tend = 5 h using the software ANSYS Fluent 19.0. This end time value was considered enough to investigate the flow structures arising from natural convection inside the tank and to show the main trends of the behavior of the tank. The procedure detailed in Section 3 is followed for all the considered cases to solve the model. The simulations are run using395 16 cores at a computer cluster composed of 25 Dell PowerEdge machines, ranging from R815 ones, equipped with AMD Opteron 6174 processors, to C6220 II systems with Intel Xeon E5-2637 v4 processors. Convergence is achieved at each time step if the globally scaled residuals for all variables are below 1e−6. 4.1. Filling level influence study400 The filling levels reported in Table 2 are considered as initial condition for the model. In Figure 7, contours of the stream function are displayed for various time values for the 80% filling level. As shown in Figure 7a, the warm liquid rises along the side wall of the tank due to the buoyancy force and turns towards the tank axis after reaching the free surface. Then, it405 moves downward forming the primary circulation in the tank. A secondary circulation is developed near the free surface as time advances, as shown in Figure 7b. More complex and transient flow structures are generated in the bottom of the tank due to the heat leak coming from the bottom wall. Since the difference between stream function values at each time snapshot represents the mass flow rate between the streamlines that they410 define, it is observed that the primary circulation in the flow is more vigorous at the early stages and weakens over time. (a) t= 15 min. (b) t= 1 h. (c) t= 2 h. (d) t= 5 h. Figure 7: Stream function contours, 80% filling level. 16 Analogous results are obtained for the 60% and 40% filling levels, as depicted in Figure 8, where the stream function contours are shown for the 40% level, evidencing similar flow circulations. However, for the lower filling levels the results differ substantially,415 especially for the 10% level. At these levels, the circulations generated at the bottom of the tank now affect the whole bulk of the liquid, as shown in Figure 9, where the streamlines for the 10% level are displayed. In all the considered cases, the highest velocity values are obtained in the boundary layer next to the tank inner wall in both in the vapor and the liquid. In Figure 10, the420 contours of velocity magnitude are displayed at t= 5 h for various filling levels. The (a) t= 15 min. (b) t= 1 h. (c) t= 2 h. (d) t= 5 h. Figure 8: Stream function contours, 40% filling level. (a) t= 15 min. (b) t= 1 h. (c) t= 2 h. (d) t= 5 h. Figure 9: Stream function contours, 10% filling level. 17 scale hides the value of bulk fluid motion, as it is comparatively smaller than that in the boundary layer. Velocities are also significantly large at the generated circulations in the bottom of the tank and upon reaching the free surface. As time advances, heat ingress results in an increase of the internal energy in both425 phases and, consequently, tank pressure rises. Even though the initial temperature inside the tank is homogeneous, the temperature profiles along the vertical coordinate show early stratification in both phases. This is depicted for the different filling levels and the time values t= 2 h and t= 5 h in Figures 11 and 12 along the symmetry axis of the storage tank. To compare the results in the liquid phase, the vertical coordinate y430 is presented as a dimensionless height y/y0, where y0is the height of the free surface, which varies for each filling level. As the boil-off rates are low in all cases, the free surface barely moves from the initial position. Thus, the values gathered in Table 2 are used as y0. In the vapor phase, the dimensionless height is defined as (y−y0)/(L−y0), where (a) 10% filling level. (b) 20% filling level. (c) 60% filling level. (d) 80% filling level. Figure 10: Velocity magnitude contours, t= 5 h. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless height, liquid phase 114 114.5 115 115.5 116 116.5 117 117.5 118 Temperature (K) 10 % 20 % 40 % 60 % 80 % (a) t= 2 h. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless height, liquid phase 114 114.5 115 115.5 116 116.5 117 117.5 118 Temperature (K) 10 % 20 % 40 % 60 % 80 % (b) t= 5 h. Figure 11: Stratification at the liquid phase for different filling levels and time values. 18 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless height, vapor phase 115 120 125 130 135 140 145 150 155 160 Temperature (K) 10 % 20 % 40 % 60 % 80 % (a) t= 2 h. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless height, vapor phase 115 120 125 130 135 140 145 150 155 160 Temperature (K) 10 % 20 % 40 % 60 % 80 % (b) t= 5 h. Figure 12: Stratification at the vapor phase for the considered filling levels and time values. L= 1.5 m is the total height of the inner vessel of the tank.435 Both in the liquid and vapor phases the stratification degree grows as time advances and is especially relevant in the area closer to the free surface. In the liquid, stratification is highly dependent on the filling level of the tank. Significant stratification is observed in the 80%, 60% and 40% levels, with the bottom part of the tank being in a subcooled state. Nevertheless, for lower levels, stratification weakens substantially. In particular,440 at the 10% level, the liquid is essentially isothermal. The vapor phase displays significant superheating, following profiles with respect the dimensionless height that are very similar in all the considered filling levels. Slightly steeper temperature gradients close to the free surface are obtained in the lower filling levels. The numerically found temperature stratification profiles are in qualitative agreement445 with experimental results reported for different set-ups in small cryogenic tanks. In [9] and [10], the measured temperature profiles in a liquid nitrogen tank also showed vapor superheating. At higher filling levels, the temperature in the free surface was larger than the bulk temperature as a consequence of a bigger stratification degree, with stratification increasing over time. Similar temperature profiles were measured in a liquid hydrogen450 tank in [15]. Total heat ingress to the tank over time for the considered filling levels is depicted in Figure 13b. For higher filling levels, the obtained heat ingress is larger. Experimental studies in cryogenic vessels, such as [37], have reported similar findings. However, a more detailed analysis of the average heat flux qlper unit area, considering only the part455 of the wall in contact with the liquid phase, reveals larger values for the smaller filling levels, as displayed in Figure 13a. This is not in contradiction with the total heat ingress results, as the surface area in contact with the liquid is smaller for lower filling levels, as reported in Table 2. The larger heat flux values in these levels can be related to the flow structures, as liquid mixing becomes more efficient, especially at the 10% level.460 Computed boil-off rates are shown over time for the different filling levels in Figure 13c. A decreasing-increasing trend with decreasing filling level is observed. In particular, a larger boil-off rate is computed at the 80% level when compared to the 60% and 40% levels, which is mainly due to the larger amount of liquid mass. This decreasing trend is not observed when considering the boil-off rate in %/day (see Table 5). For465 19 0 1 2 3 4 5 time (hours) 0 5 10 15 20 25 30 35 40 45 Average heat flux (liquid phase)(W/m2) 10 % 20 % 40 % 60 % 80 % (a) Average heat flux to liquid phase. 0 1 2 3 4 5 time (hours) 0 10 20 30 40 50 60 70 80 90 Total heat ingress (W) 10 % 20 % 40 % 60 % 80 % (b) Total heat ingress to the tank. 0 1 2 3 4 5 time (hours) 0 1 2 3 4 5 6 7 Boil-off rate (kg/s) 10-6 10 % 20 % 40 % 60 % 80 % (c) Boil-off rate. 0 1 2 3 4 5 time (hours) 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Pressure (bar) 10 % 20 % 40 % 60 % 80 % (d) Vapor pressure. Figure 13: Evolution over time of several variables depending on tank filling level. the lower filling levels (10% and 20%), an increasing trend is observed, with the boiloff rate at the 10% becoming much higher than in the remaining levels. Due to the fast self-pressurization, most of the incoming energy is used to raise liquid temperature. Therefore, the high boil-off rate for the 10% may be partly explained by the especially worse heat ingress to liquid mass ratio (see Table 2), which is three times larger than in470 the 80% level. In addition, the enhanced liquid mixing may also play a role, resulting in larger average heat fluxes per unit area to the liquid phase, as discussed previously. In Figure 13d, the vapor pressure values are depicted over time. In the considered time interval, the vapor pressure profiles are almost linear. Similar pressurization rates are obtained for all the considered levels, showing that for lower filling levels self-pressurization475 is faster. Similar results were computationally obtained in [14]. Nevertheless, our numerical results contrast with the experimental findings in [9], which report higher pressurization rates for larger filling levels. This suggests that an accurate representation of the pressurization rate at the larger filling levels might need to include the thermal expansion of the liquid, as the available volume for the vapor may be significantly reduced. Further480 studies should aim to address this effect. In Table 5, we gather the boil-off rates, total heat flux and the average heat flux to the liquid phase for all the considered filling levels, computed at the end time t= 5 h. 20 BOR [%/day] BOR [kg/s] Q[W] ql[W/m2] 10% 2.98e+ 0 6.64e−6 29.52 40.99 20% 4.03e−1 1.91e−6 36.96 34.39 40% 9.62e−2 9.26e−7 51.55 30.87 60% 1.03e−1 1.50e−6 65.86 28.95 80% 1.01e−1 1.97e−6 80.20 27.97 Table 5: Boil-off rates, total heat flux, heat flux to the liquid phase and pressurization rate for the considered filling levels at t= 5 h. The values are consequent with the prior discussion, showing that for the 10% filling level, the average heat flux to the liquid phase qlis around 50% bigger. However, the485 total heat ingress to the tank Qis substantially lower. 4.2. Insulation layer thickness influence study The described mathematical model can be used to investigate the influence of the insulation layer thickness in the flow. In particular, heat leak and boil-off in the tank can be also characterized, as well as the pressurization rate. With this purpose, the490 CAD slice of the tank, shown in Figure 1b, is slightly changed increasing or decreasing the evacuated perlite layer thickness w, which was previously equal to 11 mm. Then, the mathematical model is solved, assuming an 80% filling level. We consider the values w= 5.5 mm, w= 8.25 mm, w= 16.5 mm and w= 22 mm. In Figure 14, the stream function contours for different insulation thicknesses are495 shown at t= 5 h. The flow structures are very similar to the previously described ones for w= 11 mm. However, the contour values indicate that the circulations in the tank become more vigorous as insulation thickness decreases. This is a consequence of the higher heat ingress to the tank, which grows for lower insulation thickness values, as (a) w= 5.5 mm. (b) w= 8.25 mm. (c) w= 16.5 mm. (d) w= 22 mm. Figure 14: Streamlines for 80% filling level at time t= 5 h and different insulation thickness values. 21 0 1 2 3 4 5 time (hours) 0 10 20 30 40 50 60 Average heat flux (liquid phase)(W/m2) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (a) Average heat flux to liquid phase. 0 1 2 3 4 5 time (hours) 0 20 40 60 80 100 120 140 160 Total heat ingress (W) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (b) Total heat ingress to the tank. 0 1 2 3 4 5 time (hours) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Boil-off rate (kg/s) 10-6 w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (c) Boil-off rate. 0 1 2 3 4 5 time (hours) 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 Pressure (bar) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (d) Vapor pressure. Figure 15: Evolution over time of several variables for different insulation layer thickness. 0 0.2 0.4 0.6 0.8 1 Height (m) 114 115 116 117 118 119 Temperature (K) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (a) t= 2 h. 0 0.2 0.4 0.6 0.8 1 Height (m) 114 115 116 117 118 119 Temperature (K) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (b) t= 5 h. Figure 16: Temperature stratification along the symmetry axis for 80% filling level at the liquid phase, considering different time and insulation layer thickness values. displayed in Figure 15b, where the total heat ingress to the tank is shown over time for500 all the considered values of w. The average heat flux to the liquid phase through the tank inner wall is shown over time in Figure 15a, which is also larger as thickness decreases. Larger heat ingress to the tank leads to increased boil-off rates, as depicted in Figure 22 1.2 1.25 1.3 1.35 1.4 1.45 1.5 Height (m) 120 130 140 150 160 170 180 Temperature (K) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (a) t= 2 h. 1.2 1.25 1.3 1.35 1.4 1.45 1.5 Height (m) 120 130 140 150 160 170 180 Temperature (K) w = 5.5 mm w = 8.25 mm w = 11 mm w = 16.5 mm w = 22 mm (b) t= 5 h. Figure 17: Temperature stratification along the symmetry axis for 80% filling level at the vapor phase, considering different time and insulation layer thickness values. 15c. Pressurization rates behave similarly, increasing with bigger heat leak into the tank, as displayed in Figure 15d. This behavior was observed experimentally in [15], where505 higher pressurization rates were measured for increased heat ingress to the tank. In Figures 16 and 17, the temperature values along the symmetry axis are compared for various insulation thicknesses and time values, both in the liquid and vapor phases. The results evidence that the stratification degree is increased as the heat ingress to the tank grows, resulting in larger temperature differences between the bottom of the tank510 and the free surface. In addition, larger heat ingress also results in faster superheating of the vapor. In Table 6, we gather the boil-off rates, total heat flux and average heat flux to the liquid for all the considered insulation thicknesses, computed at the final time t= 5 h. The relation among the heat ingress and the insulation thickness is displayed in Figure515 18a, where the obtained points are fitted using a function of the form f(Q) = AQB. The computed coefficients are A= 736.7 and B=−0.9249. Similar relations can be readily established among the insulation thickness and the boil-off rate, provided that the relation Insulation layer thickness (mm) 4 6 8 10 12 14 16 18 20 22 Total heat ingress (W) 50 100 150 200 Numerical result Fit (a) Insulation layer thickness. Total heat ingress (W) 50 60 70 80 90 100 110 120 130 140 150 Boil-off rate (kg/s) #10-6 1 1.5 2 2.5 3 3.5 4 (b) Boil-off rate. Figure 18: Obtained insulation layer thickness and boil-off rate with respect to the total heat ingress. 23 w[mm] BOR [%/day] BOR [kg/s] Q[W] ql[W/m2] 5.5 2.14e−1 4.19e−6 152.07 53.27 8.25 1.38e−1 2.69e−6 104.94 36.92 11 1.01e−1 1.97e−6 80.20 27.97 16.5 6.64e−2 1.30e−6 55.13 19.49 22 4.95e−2 9.66e−7 41.96 14.85 Table 6: Boil-off rates, total heat flux and heat flux to the liquid phase for the considered insulation thickness values at time t= 5 h. among the heat ingress and the boil-off rate is approximately linear, as shown in Figure 18b.520 5. Conclusions A systematic study of the flow and conjugate heat transfer inside a 0.5 m3LNG storage tank was presented using CFD. LNG was assumed to be composed of pure methane. Realistic values, extracted from the literature, were used for the properties of solid materials and fluid phases.525 Temperature stratification, boil-off and self-pressurization were studied for various filling levels and insulation layer thicknesses of the tank, as well as the dynamic behavior of both the liquid and vapor phases. The temperature profiles in the liquid phase were shown to be depend on the filling level and the insulation layer thickness. For low filling levels, almost isothermal profiles were obtained in the liquid. The degree of thermal530 stratification in the liquid phase grows with increasing filling level. The vapor shows significant superheating while the liquid remains in a subcooled state, which is consistent with the experimental findings in the literature. The obtained stratification profiles were found to be related to the flow structures due to natural convection. In particular, at low levels (10% and 20%) the circulations535 at the bottom of the tank result in better mixing of the liquid phase, which prevents thermal stratification from forming. In addition, increased average heat flux values to the liquid were obtained. However, the computed total heat ingress to the tank rises with increasing filling level. Using lower insulation thicknesses for the tank results in an increased heat leak. This540 also leads to a rise in pressurization and boil-off rates, as well as higher velocities inside the tank. The relation among insulation thickness and heat ingress was established fitting the obtained values. Further work should be conducted on addressing the effect of the thermal expansion of the liquid, which can influence self-pressurization, especially for large filling levels.545 Acknowledgements This work was partially supported by the Spanish Ministry of Science, Innovation and Universities through the Plan Nacional de I+D+i under grants MTM2015-68275-R and MTM2017-86459-R, the grant BES-2016-077228 and by FEDER and Xunta de Galicia funds under grant GRC GI-1563 ED431C 2017/60.550 24 The authors also thank the Energylab Technology Center, with whom they collaborated during the initial stages of this work. References [1] BP, BP Energy Outlook 2019 Edition, https://www.bp.com/en/global/corporate/ energy-economics/energy-outlook.html (2019).555 [2] International Energy Agency (IEA), World Energy Outlook 2018, https://www.iea.org/weo2018/ (2018). [3] G. Biscardini, R. Schmill, A. Del Maestro, Why small-scale LNG may be the next big wave, PWC, 2017. [4] D. Dobrota, B. Lali´c, I. Komar, Problem of boil-off in LNG supply chain, Trans. Marit. Sci. 2 (02)560 (2013) 91–100. [5] E. Adom, S. Z. Islam, X. Ji, Modelling of boil-off gas in LNG tanks: A case study, Int. J. Eng. Technol. 2 (4) (2010) 292–296. [6] L. A. Pellegrini, S. Moioli, F. Brignoli, C. Bellini, LNG technology: the weathering in above-ground storage tanks, Ind. Eng. Chem. Res. 53 (10) (2014) 3931–3937.565 [7] M. F. Hasan, A. M. Zheng, I. Karimi, Minimizing boil-off losses in liquefied natural gas transportation, Ind. Eng. Chem. Res. 48 (21) (2009) 9571–9580. [8] S. Gursu, S. Sherif, T. Veziroglu, J. Sheffield, Analysis and optimization of thermal stratification and self-pressurization effects in liquid hydrogen storage systems -— part 1: model development, J. Energy Res. Technol. 115 (1993) 221–227.570 [9] M. Kang, J. Kim, H. You, D. Chang, Experimental investigation of thermal stratification in cryogenic tanks, Exp. Therm. Fluid Sci. 96 (2018) 371–382. [10] M. Seo, S. Jeong, Analysis of self-pressurization phenomenon of cryogenic fluid storage tank with thermal diffusion model, Cryogenics 50 (9) (2010) 549–555. [11] S. Das, S. Chakraborty, P. Dutta, Studies on thermal stratification phenomenon in LH2 storage575 vessel, Heat Transf. Eng. 25 (4) (2004) 54–66. [12] M. S. Zakaria, K. Osman, M. N. A. Saadun, M. Z. A. Manaf, M. Hanafi, M. Hafidzal, Computational simulation of boil-off gas formation inside liquefied natural gas tank using evaporation model in ANSYS Fluent, in: Applied Mechanics and Materials, Vol. 393, Trans Tech Publ, 2013, pp. 839– 844.580 [13] S. Roh, G. Son, G. Song, J. Bae, Numerical study of transient natural convection in a pressurized LNG storage tank, Appl. Therm. Eng. 52 (1) (2013) 209–220. [14] S. Barsi, M. Kassemi, Numerical and experimental comparisons of the self-pressurization behavior of an LH2 tank in normal gravity, Cryogenics 48 (3) (2008) 122–129. [15] N. Van Dresar, C. Lin, M. Hasan, Self-pressurization of a flightweight liquid hydrogen tank: Effects585 of fill level at low wall heat flux, NASA TM 105411 (1992). [16] S. P. Kumar, B. Prasad, G. Venkatarathnam, K. Ramamurthi, S. S. Murthy, Influence of surface evaporation on stratification in liquid hydrogen tanks of different aspect ratios, Int. J. Hydrogen Energy 32 (12) (2007) 1954–1960. [17] J.-J. Ren, J.-Y. Shi, P. Liu, M.-S. Bi, K. Jia, Simulation on thermal stratification and de-590 stratification in liquefied gas tanks, Int. J. Hydrogen Energy 38 (10) (2013) 4017–4023. [18] J. Shi, M. Bi, X. Yang, Experimental research on thermal stratification of liquefied gas in tanks under external thermal attack, Exp. Therm. Fluid Sci. 41 (2012) 77–83. [19] J. Ren, H. Zhang, M. Bi, J. Yu, S. Sun, Numerical investigation of the coupled heat transfer of liquefied gas storage tanks, Int. J. Hydrogen Energy 42 (38) (2017) 24222–24228.595 [20] J. H. Lee, Y. J. Kim, S. Hwang, Computational study of LNG evaporation and heat diffusion through a LNG cargo tank membrane, Ocean Eng. 106 (2015) 77–86. [21] S. W. Choi, H. S. Kim, W. I. Lee, Analysis of leaked LNG flow and consequent thermal effect for safety in LNG cargo containment system, Ocean Eng. 113 (2016) 276–294. [22] A. Saleem, S. Farooq, I. A. Karimi, R. Banerjee, A CFD simulation study of boiling mechanism600 and BOG generation in a full-scale LNG storage tank, Comput. Chem. Eng. 115 (2018) 112–120. [23] M. Kulitsa, D. A. Wood, LNG rollover challenges and their mitigation on Floating Storage and Regasification Units: New perspectives in assessing rollover consequences, J. Loss Prevent. Proc. 54 (2018) 352–372. [24] E. W. Lemmon, M. O. McLinden, D. G. Friend, Thermophysical properties of fluid systems, in:605 P. Linstrom, W. Mallard (Eds.), NIST Chemistry WebBook, NIST Standard Reference Database 25