Model based analysis of the boil-off gas management and control for LNG fuelled vessels

The immense pressure to decarbonise the maritime industry has led to the Lique ﬁ ed Natural Gas (LNG) uptake as a marine fuel and the LNG fuelled ships design. As LNG is stored in cryogenic conditions, the heat ingress from the ambient causes the boil-off gas (BOG) production, which, if not controlled, results in the tank overpressure with implications on the fuel storage system safe operation. This study aims at investigating an LNG storage tank behaviour for realistic operating conditions of an LNG fuelled ocean-going ship, targeting to identify the recommended control actions for avoiding tank overpressure. A dynamic model is developed by considering the mass and energy conservation in the liquid and vapour subsystems, the energy conservation in the tank walls, the vapour to liquid equilibrium (VLE), and real gas properties. Following the model validation for a holding test, simulation runs were performed for various operating conditions corresponding to typical long and short voyages of the investigated ship. The simulation results demonstrate that a boil-off gas (BOG) compressor capacity of 450 kg/h along with its on/off control setting the upper and lower limits of the tank absolute pressure at 5.5 and 4 bar leads to tank overpressure avoidance and the minimum number of BOG compressor activations. © 2022 The Authors. Published by


Background
Alternative fuels, such as Liquefied Natural Gas (LNG) [1], biofuels [2], hydrogen [3] and ammonia [4] can provide solutions in the short-term as well as to the medium to long-term, leading the maritime industry towards net-zero emissions operations.Among these alternative fuels, LNG is deemed a feasible solution contributing to the maritime industry transition towards decarbonisation in the short-term.The use of LNG in marine engines provides carbon dioxide (CO 2 ) emissions reduction by 25%, reduces the nitrogen oxides (NO x ) emissions by 75e90% depending on the engine type, whereas almost eliminates the sulphur oxides (SO x ) and particulate matter (PM) emissions [5,6].However, burning natural gas in marine gas or dual fuel engines increases the hydrocarbon emissions due to methane slip, although this issue has been partially addressed by improving the engine design and optimising the engine settings.
The commercial viability and the market acceptance of the LNG fuelled ships are expected to increase in the next years, as the shipboard LNG storage and process equipment further matures, yielding its end-use and cost advantages [5e8], compared to other alternative marine fuels.Nonetheless, several challenges for the LNG fuelled vessels are related to the lower fuel density (thus, the required greater storage volume), its storage in cryogenic conditions, as well as the boil-off gas (BOG) generation due to the inevitable heat ingress from the ambient [9].The latter can result in the tank overpressure, which, in turn, can have severe implications to the environmental and safety performance of the LNG fuelled ships.To address this challenge, the system behaviour and interactions must be well understood, whereas appropriate BOG management and pressure control systems are required.
LNG is a methane-rich mixture that also contains ethane, propane nitrogen and traces of heavier alkanes.Thus, due to the heat ingress, preferential vaporisation of the more volatile components takes place with the remaining liquid subsystem becoming richer in heavier components; this process is also called "weathering" (or "ageing") [9].The time variations of the LNG and the BOG compositions influence both its thermodynamic properties (boiling temperature, latent heat) and physical properties (heating value); the latter impacts the operation and performance of the ship machinery burning natural gas (main and auxiliary engines).For the LNG shore industry, the LNG weathering may impact the plant capability to provide LNG in the grid.The weathering process and its prediction are of particular significance for the LNG shore and shipping industries, as well as regasification plants, as accurate estimation of the stored LNG composition contributes to the effective management of the LNG cargo ships and terminals [9,10].As the LNG fuel weathering is expected to have varying effects in LNG fuelled ships [11], the sufficient understanding of the corresponding phenomena are essential for the effective management of the system operation and the control of the tank pressure.

Literature review
Early research studies focused on establishing accurate BOG rates by modelling the pure methane evaporation [12,13], and providing respective experimental data [14].More recent studies focused on developing models to represent various aspects of the LNG production, storage and transfer, including its weathering process considering both shore and ships LNG storage tanks [15e26].Studies pertaining to LNG cargo ships investigating the BOG generation and the transported LNG weathering were reported in Refs.[15e18].The BOG utilisation and processing in LNG terminals were investigated in Refs.[27,28], whereas the BOG system management and operation were studied in Refs.[29,30].The BOG re-liquefaction at LNG receiving terminals is investigated in Ref. [31].Several studies focused on modelling and investigating the rollover phenomena taking place in LNG storage tanks [32e35].The controllers design to adjust the LNG tanks pressure is reported in Refs.[36,37].Dimopoulos and Frangopoulos [16] developed a model to predict the LNG quantity and composition time variations, by coupling the non-linear Vapour-Liquid Equilibrium (VLE) with conservation equations, whilst making simplifying assumptions to keep the computational cost low.A typical LNG cargo vessel voyage of 25 days was investigated, and it was deduced that the BOG quantity and composition significantly vary during voyage.Miana et al. [17] presented a physical model and a data-driven model (based on Artificial Neural Networks) to evaluate the LNG cargo weathering process during typical ship voyages.The physical model was largely based on the mass and energy balances as well as several simplifying assumptions, such as equilibrium between the liquid and vapour subsystems during the voyage and constant evaporation rate, exhibiting considerable errors (up to 20%) compared to respective experimental data.Miana et al. [18] investigated the assumptions of constant evaporation rate and constant heat flow, whereas the derived simulation results were compared with pertinent literature data and historical data from 558 voyages.Errors up to 10% were reported for both modelling assumptions, concluding that the constant evaporation rate provided slightly more accurate results (compared to the constant heat flow).

Nomenclature
Lin et al. [19] employed the equivalent thermal conductivity model to study the effect of the tank filling ratio to the boil-off rate, concluding that the heat transfer due to conduction dominates the BOG generation process.Pellegrini et al. [20] employed the assumption of vapour-liquid equilibrium, whereas the BOG rate was calculated directly from the heat flow to the tank by using the insulation layer characteristics.
In the studies discussed in the preceding paragraphs, it was assumed that the liquid zone is in thermal equilibrium with the vapour zone.Nonetheless, experimental studies supported that the vapour typically reaches the superheated state with its temperature being higher than the liquid temperature [9].This is accounted for in non-equilibrium models, which exhibit superior performance in predicting the liquid and vapour zones parameters variations [9].
More detailed models were developed in Refs.[9,38e45] for shore LNG tanks, by employing the VLE and calculating the heat ingress accounting for the daily and seasonal ambient temperature variations.These models considered thermodynamic equilibrium between the vapour and liquid subsystems, as well as different temperature and heat transfer between the two subsystems (from the hotter vapour to the colder liquid), concluding that both provide adequate accuracy; nonetheless, the thermodynamic equilibrium approach exhibits several advantages.Migliore et al. [9] was among the first studies that considered different temperature of the liquid and vapour, separately treating the heat ingress into each subsystem, whilst estimating the heat transferred from the hotter vapour to the colder liquid.The weathering effects of the long term LNG storage (52 weeks) for a 165,000 m 3 shore terminal tank were investigated, concluding that industrial evidence supported the hypothesis of the heat transfer between the vapour and liquid taking place by conduction; the inclusion of which provided more realistic results compared to the BOG rate over-estimation (in the case of equilibrium models).Huerta and Vesovic [38] developed a non-equilibrium model for large LNG storage tanks of constant pressure assuming the heat transfer by advection due to evaporated LNG upward flow, and conduction between vapour and liquid.It was concluded that the advective upward flow dominates the energy transfer within the vapour, whilst the natural convection within the vapour is negligible.
Using both equilibrium and non-equilibrium models, Hulsbosch-Dam et al. [39] studied the behaviour of LNG storage tank under fire sources, highlighting that the equilibrium model was able to describe the final heat ingress phenomena; however, the initial heat-up was not accurately captured, as the model did not include the tanks walls thermal inertia.The non-equilibrium model results, though, were in sufficient agreement with respective experimental data.Effendy et al. [40] studied an LNG storage tank at a regasification terminal by using a thermodynamic nonequilibrium model with the hypothesis of a thin film at the vapour-liquid interface in equilibrium with the liquid subsystem.This allowed for a formulation in which the LNG evaporation maintained the boiling point conditions.Although the model predictions were not validated against experimental data, the results demonstrated about 15 o C difference between the vapour and liquid temperatures during 15 days of operation.
In the studies discussed in the preceding paragraphs, the tank pressure was assumed constant.However, for LNG storage applications in both the automotive industry and LNG fuelled ships, the pressure time variation is critical for the tanks design and operation.Wang et al. [41] employed both the equilibrium and nonequilibrium models to evaluate the performance of vertical and horizontal LNG storage tanks in refuelling stations, based on pure methane considerations.By employing experimental data and respective results derived from Computational Fluid Dynamics (CFD) models, it was inferred that both models exhibited similar accuracy at steady state conditions.However, the non-equilibrium model provided predictions of adequate accuracy at dynamic conditions, where the equilibrium model failed.
Very few studies investigated the behaviour of storage tanks for LNG fuelled ships.Van Essen et al. [42] employed the model presented in Ref. [16] to predict the effect of natural and forced boil off gas composition on the knock resistance of dual fuel engines used in LNG carriers and LNG fuelled ships, providing insights on the fuel quality changes during voyages.However, this model does not consider the tank pressure variation.Thiaucourt et al. [43] employed a thermodynamic model considering three zones to predict the temperature gradients during bottom filling of LNG storage tanks.It was highlighted that the maximum tank filling ratio depends on the tank operating pressure and the opening of the BOG vent is required to achieve the complete tanks filling at the cost of strong thermal spatial and temporal gradients.Thiaucourt et al. [44] investigated the low pressure storage and feeding system of an LNG fuelled ship by employing a dynamic model, and estimated the methane number of the fuel supplied to the ship main engines during a typical voyage.A control approach is proposed to retain the methane number above its required lower limit; however, it was concluded that further studies are required for the system detailed design.
The investigation of the complex phenomena occurring in LNG storage tanks onboard ships require more detailed modelling approaches.Wu and Ju [45] reported the sloshing modelling of an LNG tank based on the volume of fluid method coupled with the mesh motion, concluding that the sloshing greatly affects the BOG generation.Such models can be employed in the tanks and its systems design phase, however, less computational intensive approaches are required for studies on the BOG management and pressure tank control.
Based on the preceding literature review, the following research gaps were identified: lack of studies focusing on storage tanks employed in LNG fuelled ships; sailing modes and operating conditions encountered in typical voyages are usually not accounted for; such considerations are needed to derive simulation results that realistically represent the tank actual operation; lack of simulation tools to support both the design process of the storage/feeding systems of such ships and the development of the BOG management and tank pressure control.

Aim and novelty
This study aims to investigate by simulation the behaviour of an IMO type C cryogenic tank in various realistic operating scenarios based on typical LNG fuelled vessel operating profiles.To accomplish this aim, a lumped parameter model is developed, based on the first principles of the mass and energy conservation as well as the VLE.The followed modelling approach exhibits relatively low computational cost, which renders it appropriate for studying the LNG tank boil-off gas management.The vapour and gas mixtures properties are calculated by employing the REFPROP database [46], whereas the material thermal conductivities and correlations for the heat transfer coefficients are based on the pertinent literature.The model results are verified utilising experimental data available in the pertinent literature [14], and is validated against respective data for a 30-day holding test.Subsequently, the tank operation in two typical scenarios (long and short voyages) of the considered LNG fuelled vessel is investigated.The derived results are analysed to conclude on the recommended BOG compressor capacity and control approach to retain the tank pressure within the normal operation range.
To the best of the authors' knowledge, this study is the first that investigates realistic scenarios for the considered LNG storage tank operation taking into account typical characteristics of the ship voyage and the power plant as well as the prevailing ambient conditions variation.Furthermore, the developed model extends previous approaches considering two-zones along with the respective wall segments and a massless interface layer at saturated conditions.This study provides decision support to the system designer and contributes to obtaining a better understanding for the LNG storage tank behaviour and interactions, by sufficiently describing the underlying phenomena and predicting the performance parameters with adequate accuracy.The results were employed to select the characteristics of the investigated system, in particular, the required BOG compressor capacity and the tank pressure control approach.
The remaining of this study is organised as follows.Section 2 underlines the problem statement and presents the developed modelling framework.Section 3 presents the case study of the investigated ocean-going vessel, the LNG storage tanks and considered operating scenarios.Section 4 describes the derived results and discusses the respective findings based on their analysis.Finally, Section 5 summarises the findings and conclusions drawn from this study.

System description
The schematic of the investigated LNG storage tank along with the modelled subsystems and wall layers are presented in Fig. 1.Two subsystems are employed to describe the liquid and vapour zones, respectively, whereas the tank wall consists of three layers: inner metal, insulation, and outer metal, as shown in Fig. 1(b).The typical tank operation includes the pumping of the LNG from the liquid subsystem to the ship fuel feeding system, which supplies fuel to the ship machinery.Furthermore, the BOG can also be removed from the vapour zone via the BOG compressor, to retain the tank pressure below its permissible upper limit.The heat ingress causes evaporation from the liquid zone and production of BOG, which flows to the vapour zone.The whole system receives heat from the surrounding environment through heat inflows into the vapour and liquid zones.

Model assumptions
Attention was paid to develop a model that sufficiently captures the occurring phenomena and adequately represent the LNG tank behaviour.To this end, the following assumptions were made: The tank is assumed to be at even keel, and completely quiescent.Therefore, any fluid movement is triggered by heat ingress to the tank, which is responsible for both free convention and boil-off production.The mixture in both the liquid and vapour subsystems consists of the following species: Methane (CH 4 ), Ethane (CH 3 ), Propane (C 3 H 8 ), n-Butane (C 4 H 10 ), and Nitrogen (N 2 ).The liquid subsystem state can be saturated or subcooled; the vapour subsystem state can by saturated or superheated.This study though considers that the liquid is at saturated conditions, whereas the vapour is at the superheated state.The vapour subsystem mixture behaves as a real gas.Both subsystems are considered homogeneous, and their constituents do not react with each other.Heat transfer between the tank and the surrounding environment occurs by conduction and convection only, whereas the effects of solar irradiation are neglected.The heat ingress from the surrounding environment into the liquid and vapour subsystems is accounted for each subsystem separately, assuming a uniform temperature at the corresponding inner tank walls; the conduction between the wall segments in contact with the liquid and vapour subsystems is neglected.Thermal expansion or contraction of the tank is ignored, and any volume occupied by any machinery (pumps, piping, etc.) within the tank is neglected.Due to the dominant effect of the tank insulation layer, the outer steel vessel (tank exterior wall) is neglected in the heat transfer calculations, whereas a constant convective heat transfer coefficient is considered to account for the heat transfer between the ambient and the tank outer surface.A massless film is considered at the vapoureliquid interface, which exchanges mass and heat with the liquid and vapour subsystems [40,43].This film is considered to be at the saturation temperature, whereas its liquid part has the same composition as the tank liquid subsystem.
Vapoureliquid equilibrium is considered in this film interface, which in combination with the energy and moles conservation provides a system of algebraic equations to calculate the film vapour subsystem composition and the evaporation molar rate (BOG molar rate).
The physical properties of the mixtures in both subsystems are calculated by utilising the REFPROP software package [46], which employs the GERG-2008 Equation Of State (EOS) for mixtures consisting of natural gas fluids [47].This EOS is valid for the homogeneous gas, the liquid, the supercritical regions, and the vapour-liquid equilibrium states.Despite its higher mathematical complexity, it was used as the reference EOS in the pertinent literature [15,44,48e50] to benchmark the results of other EOSs or to develop simplified models for calculating the natural gas mixtures properties.Thus, it is suitable for natural gas applications and particularly for the design of cryogenic LNG systems [51], where highly accurate calculation of the mixture thermodynamic properties is required [48,49].
The upper limit of the tank pressure was considered 5.4 barg (6.4 bara).

Tank geometry
The tank is modelled as a cylinder with hemispherical ends.The liquid subsystem volume (V L ) is calculated as a function of the liquid height (H) and the inner tank diameter (D) according to where V s and V c denote the liquid volume in the spherical and cylindrical parts of the tank, respectively.The total area of the inner steel vessel in contact with the liquid A L and vapour A V are calculated by with A s and A c being the areas of the spherical and cylindrical parts of the tank in contact with the liquid, respectively.Finally, the free surface area of the liquid A FS , and the corresponding circumference C FS are calculated by

Liquid subsystem
The conservation of the total liquid mixture moles provides the following equation where n L refers to the total moles in the liquid subsystem, _ n e corresponds to the evaporation molar flow rate, and _ n L;out refers to the molar LNG outflow from the tank.It must be noted that Equation ( 9) takes different forms depending on the LNG tank operating modes.For modelling the holding test, the first term of the righthand side is only employed.
The liquid mixture components molar balance is expressed by the following equation where y F,i refers to the molar fraction of the ith component in the thin vapour-liquid interface film, and x i corresponds to the molar fraction of ith component in the liquid subsystem.
For the liquid subsystem (considering a mixture of N components), the model employs the total moles conservation (Equation ( 9)), along with N À 1 equations representing the molar conservation of the respective components (Equation (10)).This approach conserves the total moles and reduces numerical errors due to the species with very small molar fractions [40].By integrating these equations, the total moles as well as the N À 1 components' moles are calculated, which also allows for calculating the last component moles amount.The nitrogen (N 2 ) was selected to be the species calculated implicitly, as this is the most volatile constituent and thus, it is expected to exhibit the lowest concentration in the liquid zone.Hence, the N 2 moles and composition of the liquid subsystem are calculated from the following equations The temperature of the liquid subsystem (T L ) is derived by applying the energy conservation, which after some manipulation provides the following equation where h L denotes the specific molar enthalpy of the liquid subsystem, and h V,F is the specific molar enthalpy of the vapour in the thin vapour-liquid interface that is in equilibrium with the liquid subsystem.
The model integrates Equation (12) to calculate the temperature of the liquid subsystem.Then, the liquid subsystem temperature and composition (from Equation (11)) are used in REFPROP to calculate the liquid subsystem properties.
By employing the definition of the liquid molar volume, i.e., moles divided by density, and considering the density as a function of the liquid subsystem temperature only, the following equation is derived for the liquid volume time derivative with r L being the molar density of the liquid subsystem.
The integration of Equation ( 13) provides the liquid subsystem volume.Considering that the tank volume remains constant, the vapour subsystem volume is calculated by The heat ingress into the liquid subsystem ð _ Q L Þ consists of the heat ingress through the respective wall segment (which is a combination of conduction and convection, and depends on the temperature difference between the surrounding air and the liquid subsystem) and the heat due to convection from the vapour subsystem (which is also a function of the respective temperature difference).The following equation applies with _ Q F being the heat transferred to the liquid from the vapour via the vapour-liquid interface, T LW being the temperature of the inner steel surface in contact with the liquid with a total area of A L , and U L being the overall heat transfer coefficient from the surrounding environment of the tank to the inner steel surface in contact with the liquid subsystem.
The overall heat transfer coefficient is calculated by the following equation where a amb is the ambient convective heat transfer coefficient, t s , t i refer to the thickness of the tank inner steel plate and insulation material, respectively, whilst k s , k i denote their corresponding thermal conductivities.It must be noted that both k s , k i are functions of temperature, as it is discussed in the next section.The application of the energy conservation in the inner wall segment in contact with the liquid provides the following differential equation for the calculation of this wall segment temperature (T LW ) where r s is the density of the steel vessel, and a L is the convective heat transfer coefficient from this wall segment to the liquid, which is evaluated by considering the following equation where k L denotes the thermal conductivity of the liquid calculated from REFPROP, L L ¼ A SF /C SF denotes the characteristic length that defines the convective domain, and Nu L is the Nusselt Number.The latter is evaluated according to the horizontal heated upwardfacing plate with uniform wall and liquid temperatures model [52], by employing the following equations Nu L;tur ¼ 0:14 C l ¼ 0:671 with Pr L , Ra L being the Prandtl and Rayleigh numbers, respectively.

Vapour subsystem
The total moles balance for the vapour subsystem provides the following equation where n V refers to the total moles in the vapour, and _ n V;out is the molar BOG outflow.
The individual components' moles balance provides the following equation where y i refers to the molar fraction of component i in the vapour subsystem.
Similarly to the liquid subsystem considerations, the composition of one of the components is implicitly calculated (the n-butane is the least volatile component that was selected for the vapour subsystem).Hence, the n-butane moles and composition of the vapour subsystem are calculated from the following equations As the mixture is considered a real gas, its properties depend on the composition, temperature, and pressure.To avoid a very complicated equation for the calculation of the vapour subsystem temperature time derivative, the model employs the total enthalpy of the vapour subsystem (n V h V ).In this case, the energy balance provides the following equation The model integrates Equation (26) to calculate the total enthalpy, which is subsequently divided by the total moles of the vapour subsystem to calculate the vapour molar specific enthalpy.This parameter along with the vapour composition (calculated by Equation ( 25)) is provided as input to REFPROP to calculate the tank pressure and the vapour subsystem temperature, as well as the other properties of the vapour subsystem.
The heat ingress into the vapour subsystem ð _ Q V Þ is calculated by the following equation, which is derived by considering the heat ingress from the wall segment in contact with the vapour and the heat transfer from the vapour subsystem to the liquid-vapour interface where T VW and A V denote the temperature and total area of the tank inner steel surface in conduct with the vapour zone, whereas U V is the overall heat transfer coefficient from the surrounding environment of the tank to this thanks inner steel surface.The latter is evaluated in a similar manner as in the liquid subsystem, according to The temperature of the inner vessel steel (inner wall segment) in contact with the vapour zone is evaluated by using Equation ( 29), which was derived by applying the energy conservation for the inner wall segment.
where a V is the convective heat transfer coefficient of the vapour, which is evaluated by considering the following equation with k V being the thermal conductivity of the vapour, L V ¼ V V /A V being the characteristic length that defines the convective domain, and Nu V being the Nusselt Number for the vapour subsystem.The latter is evaluated as a function of the Rayleigh number (Ra V ) as reported in Ref. [52], according to

Liquid-vapour film interface
The vapoureliquid equilibrium (VLE) applied at the film interface [16,53,54] provides the following equation where 4 is the fugacity coefficient (defined as the ratio of fugacity to pressure), x F,i is the film interface liquid composition (assumed equal to the liquid subsystem composition x i and provided as input to the VLE calculations), y F,i denotes the composition of the evaporated fluid at the film interface (which is mixed with the vapour subsystem), p is the tank pressure, whereas T sat denotes the saturated temperature of the liquid and vapour subsystems in the film interface.A system of N À 1 equations is formulated, from which the composition of methane in the evaporated fluid at the film interface is implicitly calculated, similarly to the vapour and liquid subsystems.
The film interface energy and moles balances govern the evaporation rate.It is considered that a liquid recirculation flow _ n r exists (that flows from the liquid subsystem to the interface and from the interface to the liquid subsystem) [55], whereas the evaporated fluid flows from the interface to the vapour subsystem.Heat flows from the vapour subsystem to the interface and causes the evaporation of the liquid subsystem at the interface.It is also assumed that the film interface does not accumulate mass (moles).After some manipulation of the respective equations, the following equation is derived for the calculation of the evaporated flow molar rate This approach reduces the number of equations and thus the required computational effort, without sacrificing accuracy.Equations (32) and (33) form a system of non-linear algebraic equations (similar to that described in Ref. [40]) that is solved to calculate y F,i and _ n e .The heat transfer rate from the vapour subsystem to the liquid interface is calculated by employing the following equation where T L,F denotes the temperature of the liquid interface, and a F refers to the convective heat transfer coefficient from the vapour subsystem to the interface.The latter is evaluated according to Equation ( 30), whereas the required properties are evaluated at a temperature equal to the average of T V and T sat .
Finally, the molar flow rate of the liquid subsystem recirculation is evaluated according to Ref. [55] by employing the following equation where r L , r L;F correspond to the molar densities of the liquid subsystem and the film interface liquid, respectively, whereas n L and k L refer to the kinematic viscosity and thermal conductivity of the liquid subsystem, respectively.

BOG compressor control
The tank pressure is controlled by employing an on/off controller that regulates the BOG outflow according to the following equation for deactivated BOG compressor at the previous timestep: for activated BOG compressor at the previous timestep: where p UL BOG , p LL BOG correspond to the upper and lower pressure limits for the activation and deactivation of the BOG compressor, whereas _ n nom V;out denotes the BOG compressor nominal molar flow rate.

Model summary
In total, the model consists of 15 state variables ðy 2R 15 Þ and 1 control variable ðu 2RÞ, forming a system of Ordinary Differential Equations (ODEs) of the form y 0 ¼ F (t, y, u).The system of ODEs and the control variables are summarised in Table 1.The model was developed in the Matlab software [56].An explicit Runge-Kutta method based on the Dormand-Prince (4, 5) pair [57] was employed as the ODE solver.
The model employs as input the investigated LNG storage tank dimensions as well as the time variations of the ambient temperature and the liquid outflow rates.In addition, initial conditions must be provided for the state variables (operating pressure, liquid zone composition, filling ratio).At each time step, the heat transfer rates into the liquid and vapour subsystems are evaluated.The VLE calculations are performed leading to the estimation of the film interface compositions, which are subsequently employed in the time derivatives calculations taking place at the end of each timestep.
Prior to conducting the simulation runs described in the following sections, convergence tests were performed to estimate an appropriate maximum time step.It was deduced that the upper bound of the time-step for the typical voyage simulations is less than one day, under the assumption of constant environmental temperature.Similar findings were reported in the pertinent literature [9].

Case study
The investigated LNG storage tank is cylindrical and has a volume of 1750 m 3 .Its inner chamber has a length of 30.5 m and a maximum internal diameter of 9 m.The tank wall consists of two metallic layers (inner and outer) of stainless-steel each having a thickness of 25 mm, whereas the middle layer consists of perlite with 350 mm thickness.
The insulating properties of perlite are greatly dependent on temperature, according to Equation (37) as reported in Ref. [58].The developed model calculates the perlite's thermal conductivity as the average value of the conductivities at the insulation's hot and cold boundaries.
Both the tank steel vessels are built from the same material (9-Ni steel).The outer one is expected to be at approximately the ambient temperature, as it is exposed to the environment, whereas the inner vessel temperature at each segment (liquid and vapour) approaches the respective temperatures of the liquid and vapour subsystems, respectively.The steel thermal conductivity and specific heat are calculated as functions of temperature, as reported in Ref. [46], according to Equations 38 and 39 employing the coefficients values presented in Table 2. Due to the minor temperature effect on the steel vessels' density, the latter was assumed constant and equal to 7860 kg/m 3 .
The summary of all input parameters for the investigated tank, which are employed in the simulations presented in the following sections is provided in Table 3.The initial composition of the liquid subsystem is assumed as listed in Table 4.The temperature of the liquid subsystem is assumed to be the saturation temperature corresponding to the provided initial tank pressure, which is considered equal to 2 bar (abs) in most simulation runs.The initial composition of the vapour subsystem was calculated by considering liquid-vapour equilibrium conditions, which provided the molar fractions also listed in Table 4.The initial state of the vapour for deactivated BOG compressor at the previous timestep: for activated BOG compressor at the previous timestep: subsystem was assumed superheated with its temperature being 15 o C higher than the liquid subsystem temperature.This is a realistic assumption representing the tank conditions at the end of the bunkering phase.In most simulation runs, the tank filling ratio was taken 90% to represent the condition after the tank bunkering; however, several initial values for the filing ratio were also tested.

Model verification study
To verify the model, a 30-day holding test for the investigated LNG storage tank was simulated.The ambient temperature was considered constant at 20 o C, whereas the initial tank pressure was taken at 2 bar (abs).Data for a holding test was available by the tank manufacturer, providing adequate estimations of the time variations of the tank pressure, temperature, density and volume of both liquid and vapour subsystems, as well as the BOG rate and the total heat ingress.The model was set up for the same initial conditions and tank geometric characteristics, and the performed simulation run duration was taken as 30 days.The derived simulation results along with the manufacturer data are presented in Fig. 2.
The predictions of the model in general lie well within the 5% error with respect to the manufacturer estimates, for the entire time span of the simulation.The simulated tank pressure variation follows the manufacturer data; however, a steeper rise is observed after the 20th day.The liquid and vapour temperature variations are also sufficiently predicted, with the vapour temperature being around 10 K higher than the liquid temperature.It must be noted that the manufacturer data vapour temperature was calculated in REFPROP by using the provided tank pressure and vapour density.
The calculated densities of the liquid and vapour are slightly different than the manufacturer data.More specifically, the maximum error for the liquid density occurs at the 21st day, and is approximated to 2%.Likewise, the maximum error (of 3.9%) for the vapour density is also observed at the 21st day.The observed errors (especially for the liquid phase) are attributed to the different properties model employed by the manufacturer to estimate densities.However, REFPROP is the most validated properties database according to the pertinent literature, and therefore, the simulated predictions are considered trustful.The error in the heat ingress increases as the simulation time proceeds.This is attributed to the assumptions for the thermal conductivity of the insulation material, as this information was not available.
The model predicts condensation of vapour after the 20th day.This justifies the greater increase of the tank pressure in the simulated results, as the condensation phenomena were not modelled.The increase of the liquid volume due to the liquid thermal expansion (as the liquid temperature increases) as well as the resultant decrease of the vapour volume are sufficiently predicted.The model also predicts the BOG mass flow rate, although very high values are observed in the initial steps, which are attributed to the model parameters initialisation (and pertinent differences from the manufacturer data).The maximum deviation in the BOG mass flow rate of 5.3% is observed at the 6th day.
From the preceding discussion, it can be concluded that the model behaves reasonably and predicts the tank parameters adequately.A number of additional simulation runs were performed for different values for the filling ratios and insulation thickness, the results of which (not presented in this study for brevity) also verified that the model behaves reasonably.

Investigated operating scenarios
This study considers an LNG fuelled product carrier of around 50,000 mt deadweight.Two identical LNG storage tanks are installed in the main deck of this ship and it is assumed that they are used interchangeably for 24 h each to provide natural gas to the following ship consumers: main engine, auxiliary engines, and boilers.A number of operating scenarios were investigated by simulation, in order to study the behaviour of the LNG storage tank in the following conditions: (a) typical operation in realistic voyage conditions; (b) holding test operation.The purpose of (a) was to provide recommendations for the BOG compressor capacity as well as the upper and lower limits of the tank pressure for the BOG compressor activation/deactivation.The purpose of (b) is to reveal the holding test outcome in realistic and somehow extreme conditions.For operating scenarios (a), the simulation runs consider two typical voyage scenarios (a long one and a short one) for the selected ocean-going LNG fuelled vessel taking into account various BOG compressors capacities.The investigated scenarios are listed in Table 5.
More specifically, the voyage scenarios cover various ship sailing modes including sea going with/without heating and washing, as well as cargo loading and unloading operations.Two different cases with duration of 40 days (long voyage; A1) and 30 days (short Fig. 3. LNG outflow rate utilised in voyages A 1 , A 2 .voyage; A2) were investigated corresponding to typical voyages of the selected ocean-going ship.The time variation of the total LNG consumption rate for one tank and each voyage (taking into account the ship machinery power demand) is presented in Fig. 3.A daily cyclical variation of the ambient temperature was considered as illustrated in Fig. 4 corresponding to typical conditions prevailing in these two voyages.Lastly, the holding test was simulated for a period of 60 days, with an initial filling ratio equal to 90% and a constant ambient temperature of 45 o C, to assess the time required for the tank to reach its critical pressure upper limit.

Simulation results of the investigated scenarios
This section provides an overview and discussion of the simulation results for the investigated operating scenarios.

Long voyage (A1) results and discussion
The results for the investigated scenarios of the long voyage are presented in Figs. 5 and 6.For scenarios A 1a and A 1b , the upper and lower limits (UL and LL) for the BOG compressor on/off control were set to 5.5 bar (abs) and 4 bar (abs), respectively.For scenario A 1c , the upper and lower limits for the BOG compressor on/off control were set to 5.5 bar (abs) and 3 bar (abs), respectively.The initial values for the tank pressure and filling ratio were taken as 2 bar (abs) and 90%, respectively.The following findings are drawn from the analysis of the results presented in Figs. 5 and 6.
In case of no BOG outflow, the upper pressure limit is exceeded at the 13th day.The tank filling ratio at the end of the voyage almost reaches its lower limit, which indicates that bunkering will be required at each leg of the long voyage.For scenario A 1a (BOG outflow of 100 kg/h); 5 activations of the BOG compressor are required to keep the tank pressure within the predetermined zone between the UL and LL; the first activation takes place at the 10th day, whereas the pressure reduction (from the UL to the LL) takes place within 10.5 h removing around 1.3 t of BOG.For scenario A 1b (BOG outflow of 250 kg/h); 5 activations of the BOG compressor are required to keep the tank pressure within the predetermined zone between the UL and LL; the first activation takes place at the 10th day, whereas the pressure reduction (from the UL to the LL) takes place within 9.6 h removing around 1.4 t of BOG.For scenario A 1c (BOG outflow of 450 kg/h); 4 activations of the BOG compressor are required to keep the tank pressure within the predetermined zone between the UL and LL (the LL was set to 3 bar (abs) in this scenario); the first activation takes place at the 10th day, whereas the pressure reduction (from the UL to the LL) takes place within 7.3 h removing around 3.3 t of BOG.

Short voyage (A2) results and discussion
The results for the investigated scenarios of the short voyage are presented in Figs.7 and 8.For scenarios A 2a and A 2b , the upper and lower limits (UL and LL) for the BOG compressor on/off control were set to 5.5 bar (abs) and 4 bar (abs), respectively.For scenario 2c , the upper and lower limits for the BOG compressor on/off control were set to 5.5 bar (abs) and 3 bar (abs), respectively.The initial values for the tank pressure and filling ratio were taken as 2 bar (abs) and 90%, respectively.
The following findings are drawn from the analysis of the results of Figs. 7 and 8.
For the case of no BOG outflow, the upper pressure limit is approached at the 13th day.However, due to the considerable LNG consumption required at this day, the pressure UL is exceeded at the 18th day.The tank filling ratio at the end of the voyage reduces to around 30%; slightly greater filling ratio is expected in actual cases considering that the produced BOG will be fed to the ship machinery.For scenario A 2a (BOG outflow of 100 kg/h); 4 activations of the BOG compressor are required to keep the tank pressure within the predetermined zone between the UL and LL; the first activation takes place at the 10th day, whereas the pressure reduction (from the UL to the LL) takes place within 14.5 h removing around 1.45 t of BOG.For scenario A 2b (BOG outflow of 250 kg/h); 4 activations of the BOG compressor are required to keep the tank pressure within the predetermined zone between the UL and LL; the first activation takes place at the 10th day, whereas the pressure reduction (from the UL to the LL) takes place within 4.9 h removing around 1.3 t of BOG.For scenario A 2c (BOG outflow of 450 kg/h); 3 activations of the BOG compressor are required to keep the tank pressure within the predetermined zone between the UL and LL (the LL was set to 3 bar (abs) in this scenario); the first activation takes place at the 10th day, whereas the pressure reduction (from the UL to the LL) takes place within 5 h removing around 2.0 t of BOG.
The preceding discussion supports the findings of the previous section (scenario A 1 ) with regard to the selection of the BOG compressor nominal capacity of 450 kg/h as well as the tank absolute pressure upper and lower limits (5.5 and 3 bar), as only 3 activations were required for the BOG compressor and the tank pressure was kept at lower levels for longer periods.
The computational time of each simulation run was approximately 4 min for scenario A 1 (simulated duration of 40 days), and 3 min for scenario A 2 (simulated duration of 30 days), on a Windows PC equipped with a Xeon Silver 4208 processor and RAM of 32 GB.This demonstrates that the developed model computational effort is relatively low, and therefore, the model can be employed in applications, such as BOG management.

Holding test (B)
Finally, the simulation results from the holding test, taking into account ambient temperature 45 o C, initial pressure of 2 bar (abs) and initial filling ratio of 90%, are presented in Fig. 9.The simulation results demonstrate that the tank pressure limit is reached at the 42nd day, whereas the upper limit for the filling ratio (95%) is reached 6 days later, at the 48th day.Condensation conditions are predicted after the 30th day of operation, thus the pressure rise may increase in a lower rate (providing some additional days to the pressure limit reach).However, the model does not consider the condensation conditions.From the simulation runs at several values of filling ratio, it was confirmed that the model provided reasonable results in a filling ratio range, however, the model results uncertainty increases when the filling ratio lies outside the limits zones (less than 10%, more than 90%).Similar findings have been reported in the pertinent literature [19,22].

Conclusions
The analysis of the boil-off gas (BOG) formation and its effects on the tank operating parameters was presented considering an LNG fuelled ship actual operation.A dynamic model was developed that adequately represents the tank behaviour and the involved phenomena taking place during the ship voyages as well as the tank holding.The model governing equations were derived by considering the liquid and vapour subsystems along with the corresponding wall segments and by applying the following first principles: (a) the moles and energy conservation in both the liquid and vapour subsystems; (b) the energy conservation in the respective tank walls segments; (c) the vapour to liquid equilibrium (VLE) in the interface between the vapour and liquid; and (d) equation of state for real gases for the vapour mixture.The required heat transfer coefficients were estimated based on relevant correlations in the pertinent literature.
Following the model validation against available data for a tank holding test, simulation runs were performed for various operating scenarios and the derived results were employed to provide recommendations for the tank design and operation.
The main findings of this study are summarised as follows.
Despite the made assumptions and simplifications, the developed model was proved an effective tool and provided useful insights, thus contributing to the better understanding of the involved phenomena and interactions in LNG storage tanks.The developed model is of relatively low computational cost, hence it can be employed in applications, such as the BOG management.
A BOG compressor capacity around 450 kg/h along with its on/ off control is effective for avoiding the tank overpressure, as it resulted in 3 or 4 activations of the BOG compressor for the typical short and long ship voyages, respectively.
The recommended upper and lower limits for the activation/ deactivation of the BOG compressor were 5.5 bar (abs) and 4 bar (abs), respectively.Considering the holding test of a 90% filled tank in relatively extreme ambient temperature of 45 o C, it was found that the upper pressure limit will be exceeded after the 42nd day.Based on the expected ship operation, this outcome is characterised as acceptable.
The results derived by this study can be used to provide decision support for the LNG storage tanks design and operation.The developed model can be integrated in wider simulation platforms and employed to study the operations of LNG powered vessels.Lastly, the developed model and the derived results can be further expanded in future studies focusing on the development of shorebased intelligent LNG management systems for providing decision support on actual operations of LNG fuelled ships.

Table 1
Summary of the state variables and differential equations of the mathematical model.

Table 2
Values of the coefficients utilised inEquations 38 and 39

Table 3
Investigated LNG tank particulars and parameters employed in this study.

Table 4
Initial Compositions of the liquid and vapour subsystems.

Table 5
Investigated operating scenarios.