Liquid phase physical property field in storage tank and layering prediction method

Abstract

A storage tank internal liquid phase physical property field and layering prediction method comprises the steps that radial and axial gridding rules of an LNG storage tank model are formulated, historical data are collected to correct a preset model so as to generate a dynamic temperature initial value, temperature distribution models of a first layer and a last layer are constructed, the axial temperature of a storage tank is corrected, a boundary layer temperature distribution model is constructed, and then the storage tank internal liquid phase physical property field is obtained. The method comprises the following steps: correcting the radial temperature of a storage tank, collecting field actual measurement data, carrying out parameter estimation of a preset model through a least square method, calculating a temperature field, converting the temperature field into a physical property field, carrying out layering coefficient characterization, and finally constructing an interlayer dynamic interaction model to predict the dynamic temperature, composition and layering trend of each layer in the future. According to the method, through axial-radial grid modeling and a dynamic calibration mechanism, the layering prediction efficiency is remarkably improved, the physical property distribution change in the tank can be accurately captured, and safe operation of the storage tank is effectively guided.

Description

Technical Field

[0001] The present invention relates to the technical field of natural gas storage, and particularly to a method for predicting the liquid-phase physical property field and stratification in a storage tank. Background Art

[0002] Liquefied Natural Gas (LNG), as a clean and efficient energy source, has the characteristics of low temperature and high volatility, and its storage and transportation have extremely high technical requirements. An LNG storage tank is a special equipment for storing liquefied natural gas, usually a super-large equipment composed of an inner tank, an outer tank and an insulation layer, with a volume of more than 10000 m 3 Above, it can maintain the liquid state of natural gas at low temperature (about -162 °C). During the operation of an LNG storage tank, due to the evaporation and natural convection of the liquid, the differences in interlayer temperature, density and composition in the large space inside the tank often lead to the stratification of the liquid in the storage tank. This stratification phenomenon will not only affect the stability of the storage tank, but may also cause sudden changes in local temperature or pressure, and in severe cases, lead to the failure of the storage tank operation or safety accidents. Therefore, the accurate calculation and prediction of the physical property distribution and stratification degree inside the LNG storage tank have important engineering significance.

[0003] However, considering the configuration and maintenance costs of the instruments, it is very difficult for an LNG storage tank to be equipped with enough on-line instruments, and it is impossible to obtain the liquid-phase physical property field in the large space inside the tank in time, making it difficult to predict the liquid-phase dynamic interaction and stratification degree inside the storage tank. Therefore, there is an urgent need for a technical method that can combine on-site instruments, efficiently calculate the liquid-phase physical property distribution and accurately predict the liquid-phase stratification degree, in order to improve the safety of LNG storage tanks. Summary of the Invention

[0004] In order to overcome the deficiencies of the prior art, the present invention provides a method for predicting the liquid-phase physical property field and stratification in a storage tank, which is used to predict the dynamic physical property distribution of a super-large LNG storage tank and provide a stratification coefficient. Through axial-radial grid modeling and a dynamic calibration mechanism, the stratification prediction efficiency is significantly improved, the changes in the physical property distribution inside the tank can be accurately captured, and the safe operation of the storage tank can be effectively guided.

[0005] The purpose of the present invention is achieved through the following technical solutions:

[0006] A method for predicting the liquid-phase physical property field and stratification in a storage tank, the method comprising the following steps:

[0007] S1, complete the radial and axial meshing of the LNG storage tank according to the rules;

[0008] S2, collect the historical data of other storage tanks, generate a preset model, and generate the initial temperature value of the grid inside the storage tank;

[0009] S3. Construct the temperature distribution models for the first layer and the last layer, correct the axial temperature of the storage tank, construct the boundary layer temperature distribution model, and correct the radial temperature of the storage tank;

[0010] S4. Compare the historical data of the current storage tank with the model calculation values, and perform parameter estimation on the preset model;

[0011] S5. Convert the calculated temperature field into the corresponding physical property field and characterize it with a layering coefficient;

[0012] S6. Construct the dynamic model of the LNG storage tank to predict the future trends of the temperature field, physical property field, and layering coefficient.

[0013] Furthermore, the process of S1 is as follows:

[0014] (1.1) Axial grid division: Vertically divide the liquid phase area of the super-large vertical liquid storage tank into n layers evenly along the vertical direction, number them from top to bottom in sequence, and use k to represent the number, where k = 1, 2,..., n. The number of layers n is determined by the following formula, where H L is the total height of the liquid phase area in the LNG storage tank, with the unit of m:

[0015]

[0016] (1.2) Radial grid division: Define the boundary area of the storage tank as the area in each horizontal layer where the straight-line distance from the inner wall of the tank is not greater than 5% of the diameter D of the storage tank. The central area is the remaining part of each horizontal area except the boundary area, that is, the area where the distance from the inner wall of the tank is > 5% D, and no radial grid division is performed in the central area. Divide the boundary area into 5 layers non-uniformly along the horizontal direction, number them in sequence from the central axis to the tank wall direction, and use f to represent the number, where f = 0, 1,..., 5. When f = 0, it represents the central area of this layer. When f = 1, 2,..., 5, it represents the boundary area of this layer. The distance formula from the f-th layer of the boundary area to the tank wall is shown in formula (2);

[0017]

[0018] The thickness formula of the f-th layer of the boundary area is shown in formula (3);

[0019]

[0020] Define the radial grid label of the k-th layer as G k,f , and the schematic diagram of the grid division of the liquid phase area of the storage tank is as Figure 2 shown.

[0021] Furthermore, the process of S2 is as follows:

[0022] Collect historical data of other storage tanks. Each set of data includes the temperature T inside the storage tank, the corresponding liquid level height H, the feed temperature T in , the ambient temperature T a , and the heat resistance R of the tank wall b . Through regression analysis, a temperature generation preset model for the LNG storage tank is generated, as shown in Equation (4):

[0023]

[0024] where Z norm is the normalized independent variable set of the temperature calculation preset model, Z = [H norm , T in,norm , T a,norm , R b,norm . The normalization formula for each independent variable is shown in Equation (5). Define the parameter set as θ, θ = [0.15, 0.6, -0.3, 0.2, -0.1];

[0025]

[0026] Substitute the independent variable values collected on-site into the calculation model T model (Z norm ) to generate the initial temperature value of the central area of each layer

[0027] Furthermore, the process of S3 is as follows:

[0028] (3.1) Axial temperature correction model: Correct the temperature of the first layer , read the gas phase temperature T V on-site, and perform logarithmic mean calculation on T V and . Use the calculation result as the corrected temperature value of the first layer. Correct the temperature of the nth layer , and perform logarithmic mean calculation on the bottom temperature T B and . Use the calculation result as the corrected temperature value of the nth layer;

[0029] (3.2) Radial temperature correction: Input the tank wall material and the tank wall thickness δ b , calculate the inner wall surface temperature T k,b of each layer according to Fourier's law. The radial grid temperature in the boundary region is corrected using the smoothed particle hydrodynamics method. Use T k,0 and T k,bUsing the inner and outer boundary values respectively, calculate the grid temperature of the boundary region and update it to the temperature field. Here, the smoothed particle hydrodynamics method is a meshless numerical calculation method. By discretizing the continuous medium into a series of particles, each particle carries physical quantities such as mass, density, and temperature, and the physical quantities between particles are interpolated and approximately calculated through a kernel function. The role of the kernel function is to distribute the particle attributes within a certain spatial range, so that it affects the surrounding particles. Through the weighted summation of the kernel function, the physical quantity at any position can be estimated. In this step, the temperature field of the boundary region is discretized into particles, and a quadratic spline kernel function is used for interpolation calculation between particles. Its expression is shown in Equation (6);

[0030]

[0031] where W(r - r j , h) is the kernel function value, |r - r j | represents the distance between position r and particle j, and h is the kernel radius, which determines the range of particle action;

[0032] Calculate the temperature values of each particle in turn, and the calculation formula is shown in Equation (8);

[0033]

[0034] where m j is the mass of particle j, ρ j is the density of particle j, T j is the temperature of particle j. Repeat the above steps for all positions r in the boundary region to generate the temperature distribution at the current time.

[0035] Furthermore, the process of S4 is as follows:

[0036] Collect N groups of historical measurement data of the current storage tank. Each group of data includes the historical measurement temperature value of the current storage tank and the corresponding historical measurement liquid level H o , feed temperature T in,o , ambient temperature T a,o and tank wall thermal resistance R b , where R b is regarded as a constant value. Calculate the initial value through the grid division rule of step S1 and the preset model T model (Z norm ), and then combine the axial temperature correction model and radial temperature correction model of S3 to generate the comprehensive simulation temperature value T k,f,o of each grid. In this step, parameter estimation is only optimized for the linear parameters θ = [θ1, θ2, θ3, θ4, θ5] of T model (Z norm ). The optimization proposition is as follows:

[0037]

[0038]

[0039] ‖θ (t+1) -θ (t) ‖2≤10 -4

[0040] o = 1, 2, ···, N

[0041] The parameter θ optimized by the least squares method will be updated to T model (Z norm ), and the correction model in step S3 is a fixed algorithm and does not participate in parameter iteration.

[0042] The process of S5 is as follows:

[0043] (5.1) Input the feed composition x in , assuming that the molar compositions of each layer in the liquid phase region of the LNG storage tank are equal at the current time period and are the same as the molar composition of the feed, that is, x 1,f,t＝0 = x 2,f,t＝0 =... = x n,f,t＝0 = x in ;

[0044] (5.2) Input the Z value at the current time period, and calculate the current temperature field of the LNG storage tank through T model (Z norm ), the axial correction model and the radial correction model, that is, the temperature values T k,f of each grid;

[0045] (5.3) Substitute the temperature and molar composition of each layer into the empirical formula, and output the physical property field of the LNG storage tank, that is, the physical property parameters of each grid. The physical property parameters include density ρ, pressure P, bubble point BP, latent heat of vaporization ΔH vap , thermal conductivity λ and dynamic viscosity μ;

[0046] (5.4) Quantify the overall stratification degree and local stratification instability of the storage tank according to the physical property field output in (5.3), and calculate the current stratification coefficient of the storage tank with a comprehensive index. The formula for quantifying the overall stratification degree of the storage tank is as follows. β = 1 indicates that the liquid is in a saturated state, and the larger β is, the higher the overall stratification tendency of the liquid phase region of the storage tank:

[0047]

[0048] where P is the arithmetic mean of the pressures in the liquid phase region of the storage tank, and the calculation formula is as follows:

[0049]

[0050] Where P sat is the mass-average liquid temperature Calculated with the average mole fraction The saturation pressure, The calculation formula is as follows:

[0051]

[0052] The calculation formula is as follows:

[0053]

[0054] The formula for quantifying the local stratification instability of a tank is as follows, R k is the quantitative value of the k-th layer’s stratification instability, β T is the volume expansion coefficient caused by the temperature of the Kth layer, β S is the volume expansion coefficient caused by the concentration of the kth layer:

[0055]

[0056] β T (k), β S The calculation formula of (k) is as follows:

[0057]

[0058] The calculation formula of ΔT(k) is as follows:

[0059]

[0060] ΔT(k)=T k,0 -T k+1,0

[0061] The formula for the stratification coefficient I used to assess the stratification degree of LNG storage tanks and provide stratification warning is as follows:

[0062]

[0063] ω(R k ) is the local stability ratio R k The weight function, when (R k <0, β T (k)ΔT(k)>0) or (|R k |<1, R k >0) or (|R k |>1. R k >0), ω(R k )=10, in other cases ω(R k )=1;

[0064] Different ranges of I value correspond to different degrees of stratification severity. When I > 5, on-site personnel should promptly adjust the operating parameters to suppress stratification.

[0065] The process of S6 is as follows:

[0066] (6.1) Using the temperature and molar composition of each region output by S5 as the initial values, and the physical properties of each region output by S5 as the basis for calculating the heat terms in the dynamic model, construct a dynamic model of the LNG storage tank, and input the static parameters and dynamic parameters on-site. The static parameters include the storage tank diameter D, the height H of the liquid phase region of the storage tank L , the tank wall material, the tank wall thickness δ b , the tank bottom material, the tank bottom thickness δ B , the circulating pipeline diameter D loop , the circulating pipeline length L loop , the installation height H of the feed pipeline in , the installation height H of the discharge pipeline out , and the dynamic parameters include the feed rate the discharge rate the circulation rate the gas phase temperature T V , the ambient temperature T a , the tank bottom temperature T B , * represents the molar rate. The MESH equation is a mathematical model of the equilibrium stage separation process, which is established based on mass conservation, phase equilibrium relationship, heat balance, and summation of mole fractions. Use the MESH equation to establish an interlayer dynamic interaction model of the LNG storage tank. Assume that the temperature, molar composition, and physical properties of each layer of the dynamic model are calculated based on the value of the central region G of each layer k,0 , set the total time step and the unit time step, which are used to predict the dynamic temperature and dynamic molar composition of each layer in the LNG storage tank, and evaluate the stratification trend according to the prediction results;

[0067] In the above (6.1), the specific dynamic model of each layer is as follows:

[0068] (6.1.1) In the calculation of the dynamic model of a general single horizontal layer, heat transfer is the convective heat transfer between liquid layers and heat transfer through the tank wall, and mass transfer is the liquid-liquid interface mass transfer based on molecular diffusion. The formula is as follows, where the subscripts A and B represent methane and nitrogen respectively, and the subscript TRANS represents mass transfer based on molecular diffusion represents the specific enthalpy of methane in the central region of the k + 1 layer, Q LL,kRepresents the convective heat transfer rate from the k-th layer region to the (k - 1)-th layer region, Q b,k Is the heat transfer rate of the tank wall in the k-th layer region, and C is the number of components of LNG;

[0069] Material balance equation:

[0070]

[0071] Heat balance equation:

[0072]

[0073] Mole fraction summation formula:

[0074]

[0075] (6.1.2) In the dynamic model calculation of the first layer, heat transfer is the convective heat transfer between the gas-liquid interface and the liquid layer, and mass transfer is the gas-liquid interface mass transfer based on the double-film theory and the liquid-liquid interface mass transfer based on molecular diffusion, where the subscript ABS represents the mass transfer based on the double-film theory, Q VL Is the convective heat transfer rate from the gas phase region of the LNG storage tank to the liquid phase region of the first layer;

[0076] Material balance equation:

[0077]

[0078] Heat balance equation:

[0079]

[0080] Mole fraction summation formula:

[0081]

[0082] (6.1.3) In the dynamic model calculation of the first layer, it is necessary to judge whether evaporation occurs in the liquid phase of the first layer. If T 1,0 < BP 1,0 , then evaporation does not occur, and the dynamic model of the first layer is consistent with (4.2). If T 1,0 ≥BP 1,0 , then the evaporation rate needs to be calculated The formula is Q is the net heat receiving rate of the first layer. Update the composition and enthalpy value of the evaporated gas into the material balance and heat balance equations of the first layer dynamic model respectively. The formulas are as follows:

[0083] Material balance equation:

[0084]

[0085] Heat balance equation:

[0086]

[0087] Mole fraction summation formula:

[0088]

[0089] Phase equilibrium equation:

[0090] y i,E = K i x i,1,0 (26)

[0091] (6.1.4) In the dynamic model calculation of the first layer, if there is a circulation pipeline set above the storage tank, it is necessary to use the force balance formula to judge whether the circulating liquid can be mixed with the liquid phase of the first layer during the filling of the circulating liquid. The formula of the force balance formula is as follows. The first item on the left side is the inertial force promoting mixing, the first item on the right side is the buoyancy force resisting mixing, and the second item on the right side is the viscous stress resisting mixing. Where F is a constant, and its value range is 6 - 50, ρ h is the original LNG density, v is the internal convection velocity of LNG, and the calculation formula of the internal convection velocity of LNG in the first layer is shown in formula (29), Δρ c is the density difference between the circulating liquid density ρ RE and the original LNG density. The circulating liquid density takes the liquid phase density of the horizontal layer where the discharge pipeline is located. g is the acceleration due to gravity, D′ is the diameter of the impact area where the circulating liquid contacts the original LNG, μ is the viscosity of the circulating liquid, l s Characteristic length of the eddy of the first layer of LNG:

[0092]

[0093] If the force balance formula is not satisfied, it is regarded that the circulating liquid forms an independent layer and does not mix with the liquid phase of the first layer; if the force balance formula is satisfied, it is regarded that the circulating liquid mixes with the liquid phase of the first layer, and the composition and enthalpy value of the circulating liquid are updated to the mass transfer and heat transfer calculations of the first layer dynamic model. The formula is as follows:

[0094] Material balance equation:

[0095]

[0096] Heat balance equation:

[0097]

[0098] Mole fraction summation formula:

[0099]

[0100] Among them, the lag caused by the time for the recirculating liquid to flow back to the first layer through the pipeline needs to be considered. The lag equation is as follows, where R represents the corresponding layer that provides the recirculating liquid:

[0101]

[0102] The calculation formula for the recirculation time is as follows:

[0103]

[0104] (6.1.5) In the dynamic model calculation of the horizontal layer where the inlet and outlet pipelines are located, the composition and heat of the inlet and outlet liquids are respectively updated to the mass transfer and heat transfer calculations of the dynamic model of this layer. The layers where the inlet pipeline and the outlet pipeline are located are numbered as k in 、k out respectively. Taking the dynamic model of the horizontal layer where the inlet pipeline is located as an example, the formula is as follows:

[0105] Material balance equation:

[0106]

[0107] Heat balance equation:

[0108]

[0109] Mole fraction summation formula:

[0110]

[0111] (6.2) A calibration trigger mechanism is set up for the predicted temperature values of each region, and it is calibrated once every 24 hours. The calibration trigger mechanism is as follows:

[0112] Assume that there are p unit time steps in this period, and the temperature corresponding to each unit time step of the grid G k,f is T q,k,f , then the average temperature T k,f of the grid G avg,k,f in this period is:

[0113]

[0114] Collect the latest historical data of the current storage tank When , trigger temperature calibration, and re-estimate the parameters of the preset model until calibration is no longer triggered;

[0115] (6.3) Repeat step S5, and substitute the temperature T p,k,0 and the molar composition x p,k,0 corresponding to each time node within the total time step into the empirical formula, and output ρ at each time nodep,k,0 , P p,k,0 , BP p,k,0 , ΔH vap,p,k,0 , λ p,k,0 , μ p,k,0 and the layering coefficient I p,k,0 , if the value of the layering coefficient at each time node continuously shows I > 5, measures need to be taken in advance to prevent layering.

[0116] The beneficial effects of the present invention are as follows:

[0117] 1. The method for predicting the liquid-phase physical property field and layering in a storage tank of the present invention avoids the complex numerical simulation process and the data dependence of the simple interpolation method, realizes the efficient calculation of the dynamic physical property distribution in the storage tank, and significantly improves the calculation accuracy and efficiency;

[0118] 2. The method for predicting the liquid-phase physical property field and layering in a storage tank of the present invention considers the multiple effects of the temperature influence factor in the on-site data on the dynamic temperature distribution of the LNG storage tank, combines key parameters such as the liquid level height, and corrects the temperature calculation model through screening and regression analysis, and can utilize multiple sets of on-site operation data simultaneously, improving the adaptability and calculation accuracy of the model under different working conditions;

[0119] 3. The method for predicting the liquid-phase physical property field and layering in a storage tank of the present invention combines the local layering stability ratio R value, can evaluate the overall layering tendency and local layering instability in the storage tank in real time, and provides a layering warning mechanism through dynamic trend prediction, providing clear operation guidance for on-site operators and avoiding potential safety risks during the operation of the storage tank; Description of the Drawings

[0120] Figure 1 is a calculation flow chart of a method for predicting the liquid-phase physical property field and layering in a storage tank;

[0121] Figure 2 is a schematic diagram of the grid division of the liquid phase region of the storage tank;

[0122] Figure 3 is the dynamic temperature calculation result (first 24h) of the fifth horizontal layer of the LNG storage tank provided by the embodiment of the present invention;

[0123] Figure 4 is the dynamic methane content calculation result (first 24h) of the fifth horizontal layer of the LNG storage tank provided by the embodiment of the present invention;

[0124] Figure 5 is the dynamic trend (first 24h) of the density of the central region G 1,0 of the first horizontal layer provided by the embodiment of the present invention;

[0125] Figure 6The horizontal fifth - layer boundary region G provided by the embodiments of the present invention 5,f (when f≠0) the dynamic trend of density (for the first 24 h);

[0126] Figure 7 The calculation result of the dynamic stratification coefficient of the LNG storage tank provided by the embodiments of the present invention (for the first 24 h). Specific embodiments

[0127] The present invention will be further described below with reference to the accompanying drawings.

[0128] Refer to Figures 1 to 7 , a method for predicting the liquid - phase physical property field and stratification in a storage tank, the method comprising the following steps:

[0129] S1, complete the radial and axial meshing of the LNG storage tank according to the rules;

[0130] S2, collect the historical data of other storage tanks, generate a preset model, and generate the initial temperature values of the grids in the storage tank;

[0131] S3, construct the temperature distribution models of the first layer and the last layer, correct the axial temperature of the storage tank, construct the boundary - layer temperature distribution model, and correct the radial temperature of the storage tank;

[0132] S4, compare the historical data of the current storage tank with the model calculation values, and perform parameter estimation on the preset model;

[0133] S5, convert the calculated temperature field into the corresponding physical property field and characterize it with a stratification coefficient;

[0134] S6, construct a dynamic model of the LNG storage tank to predict the future trends of the temperature field, physical property field, and stratification coefficient.

[0135] Example: Calculate the dynamic physical property distribution and dynamic stratification degree of the storage tank of a certain LNG receiving terminal. The on - site parameters of the storage tank are shown in Table 1:

[0136] Table 1 On - site parameters of the storage tank of a certain LNG receiving terminal

[0137]

[0138] The method for predicting the liquid - phase physical property field and stratification in the storage tank includes the following steps:

[0139] S1, formulate the radial and axial meshing rules of the LNG storage tank model, and the process is as follows:

[0140] (1.1) Axial grid division: For a vertical super - large liquid - phase storage tank, the liquid - phase area is evenly divided into n layers along the vertical direction, numbered from top to bottom in sequence, and the numbering is represented by k, where k = 1, 2, …, n. The number of layers n is determined by the following formula, where H Lis the total height of the liquid phase region in the LNG storage tank, with the unit of m:

[0141]

[0142] (1.2) Radial grid division: Define the boundary region of the storage tank as the region where the straight-line distance from the inner wall of the tank in each horizontal layer is not greater than 5% of the diameter D of the storage tank. The central region is the remaining part of each horizontal region except the boundary region, that is, the region where the distance from the inner wall of the tank > 5% D, and no radial grid division is performed in the central region. The boundary region is non-uniformly divided into 5 layers along the horizontal direction and numbered in sequence from the central axis to the tank wall direction. The number is represented by f, f = 0, 1, …, 5. When f = 0, it represents the central region of this layer. When f = 1, 2 …, 5, it represents the boundary region of this layer. The distance formula from the f-th layer of the boundary region to the tank wall is as follows:

[0143]

[0144] The thickness formula of the f-th layer of the boundary region is as follows:

[0145]

[0146] Define the radial grid label of the k-th layer as G k,f , input the diameter D of the storage tank and the total height H of the liquid phase region L , perform axial and radial grid division. The storage tank is axially divided into layers, and the boundary region of the LNG storage tank is non-uniformly divided into 5 layers along the radial direction. The boundary region is the region where the straight-line distance from the inner wall of the tank in each horizontal layer is not greater than , and numbering is performed in sequence from the central axis to the tank wall direction. The number is represented by f (f = 1, 2, …, 5). The distance formula from the f-th layer to the tank wall is as follows:

[0147]

[0148] The thicknesses of each layer of the boundary region are: Δr1 = 1.43m, Δr2 = 1.26m, Δr3 = 1.06m, Δr4 = 0.83m, Δr5 = 0.50m;

[0149] S2, collect the historical data on site, correct the preset model, and generate the dynamic initial temperature value. The process is as follows:

[0150] Collect historical data. Each group of data includes the temperature T inside the storage tank, the liquid level height H corresponding to this temperature, the feed temperature T in , the ambient temperature T a , the heat resistance R of the tank wall b , and generate the preset temperature calculation model of the LNG storage tank through regression analysis:

[0151]

[0152] where Z norm is the set of normalized independent variables of the temperature calculation preset model, Z = [H norm , T in,norm , T a,norm , R b,norm . The normalization formulas for each independent variable are as follows. Define the parameter set as θ, θ = [0.15, 0.6, -0.3, 0.2, -0.1];

[0153]

[0154] Substitute the independent variable values collected on-site into the calculation model T model (Z norm ) to generate the initial values of the central region temperature of each layer where the R b term is calculated according to the following formula:

[0155]

[0156] S3. Construct the temperature distribution models of the first layer and the last layer, correct the axial temperature of the storage tank, construct the boundary layer temperature distribution model, and correct the radial temperature of the storage tank; The process is as follows:

[0157] (3.1) Axial temperature correction model: Correct the temperature of the first layer , read the on-site gas phase temperature value T V , perform a logarithmic mean calculation on T V and , and use the calculation result as the corrected temperature value of the first layer. Correct the temperature of the nth layer , perform a logarithmic mean calculation on the bottom temperature of the tank T B and , and use the calculation result as the corrected temperature value of the nth layer;

[0158] (3.2) Radial temperature correction: Input the tank wall material and the tank wall thickness δ b , calculate the inner wall surface temperature T of each layer according to Fourier's law k,b , and use the smoothed particle hydrodynamics method to correct the radial grid temperature in the boundary region. Use T k,0 and T k,bTaking them as the inner and outer boundary values respectively, calculate the grid temperature of the boundary region and update it to the temperature field. Among them, the smoothed particle hydrodynamics method is a meshless numerical calculation method. By discretizing the continuous medium into a series of particles, each particle carries physical quantities such as mass, density, and temperature, and interpolates and approximately calculates the physical quantities between particles through a kernel function. The role of the kernel function is to distribute the particle attributes within a certain spatial range so that it affects the surrounding particles. Through the weighted summation of the kernel function, the physical quantities at any position can be estimated. In this step, the temperature field of the boundary region is discretized into particles, and a quadratic spline kernel function is used for interpolation calculation between particles. Its expression is as follows:

[0159]

[0160] Among them, W(r - r j , h) is the value of the kernel function, |r - r j | represents the distance between position r and particle j, and h is the kernel radius, which determines the range of particle action;

[0161] Calculate the temperature values of each particle in turn. The calculation formula is as follows:

[0162]

[0163] Among them, m j is the mass of particle j, ρ j is the density of particle j, T j is the temperature of particle j. Repeat the above steps for all positions r in the boundary region to generate the temperature distribution at the current time;

[0164] S4. Compare the current measurement points with the model calculation values to estimate the parameters of the preset model;

[0165] Collect 50 groups of on-site measured data at different time points. Each group of data includes the measured temperature value and the corresponding measured liquid level H o , feed temperature T in,o , ambient temperature T a,o and tank wall thermal resistance R b , where R b is regarded as a constant value. Through the grid division rule in step S1, the preset model T model (Z norm ) calculates the initial value, and then combines the axial temperature correction model and the radial temperature correction model in S3 to generate the comprehensive simulated temperature value T k,f,o of each grid. In this step, the parameter estimation only optimizes the linear parameters θ = [θ1, θ2, θ3, θ4, θ5] of T model (Z norm ). The optimization proposition is as follows:

[0166]

[0167] ‖θ (t+1) -θ (t) ‖2≤10 -4

[0168] o = 1, 2, ···, 50

[0169] The parameter θ optimized by the least squares method will be updated to T model (Z norm ) where the correction model in step S3 is a fixed algorithm and does not participate in parameter iteration.

[0170] The corrected model parameters are θ = [0.14, 0.55, -0.28, 0.2, -0.1];

[0171] S5. After the model is corrected, the corrected temperature field is converted into the corresponding physical property field and characterized by the layering coefficient; the process is as follows:

[0172] (5.1) Input the feed composition x in , assuming that the molar composition of each layer in the liquid phase region of the LNG storage tank is equal at the current time period and is the same as the molar composition of the feed, that is, x 1,f,t＝0 = x 2,f,t＝0 =... = x n,f,t＝0 = x in ;

[0173] The feed molar composition x in of this storage tank is shown in Table 2:

[0174] Table 2 Feed Molar Composition of this Storage Tank (unit: mol%)

[0175] methane ethane nitrogen propane n-butane 88.25 5.62 3.87 1.56 0.70

[0176] (5.2) Input the Z value at the current time period, and calculate the corrected temperature field through T model (Z norm ), the axial correction model and the radial correction model. The temperature calculation results of the grid G k,0 in the central region of each layer are shown in Table 3 as follows:

[0177] Table 3 Temperature Calculation Results of the Grid G k,0 in the Central Region of Each Layer

[0178] position <![CDATA[G 1,0 > <![CDATA[G 2,0 > <![CDATA[G 3,0 > <![CDATA[G 4,0 > <![CDATA[G 5,0 > <![CDATA[G 6,0 > <![CDATA[G 7,0 > <![CDATA[G 8,0 > <![CDATA[G 9,0 > <![CDATA[G 10,0 > temperature 110.48 110.78 111.09 111.39 111.69 111.99 112.30 112.60 112.90 113.20

[0179] Taking the boundary region of the fifth layer as an example, the temperature at the center of the storage tank is 111.69 K calculated by the temperature calculation model of the LNG storage tank, and the temperature of the inner wall of the storage tank is 112.363 K calculated according to Fourier's law. The temperature calculation results of the grid in the boundary region of the fifth layer are shown in Table 4:

[0180] Table 4 Grid G in the boundary region of the fifth layer 5,f (f≠0) Temperature calculation results

[0181]

[0182]

[0183] (5.3) Substitute the corrected temperature and material composition into the empirical formula to output the physical property field of the LNG storage tank, that is, the physical property parameters of each grid. The physical property parameters include density ρ, pressure P, bubble point BP, latent heat of vaporization ΔH vap , thermal conductivity λ and dynamic viscosity μ. Among them, the density ρ, thermal conductivity λ and dynamic viscosity μ are estimated by the empirical formula in the physical property handbook. The pressure P of the kth layer k is calculated using P hydro,i is the static pressure of the ith layer, i = 1, 2,... k - 1. The bubble point BP and latent heat of vaporization ΔH vap are calculated using the Peng-Robinson equation;

[0184] Calculate the physical property parameters of each region in turn. Taking the grid G in the central region of the fifth layer as an example 5,0 , the density ρ, thermal conductivity λ and dynamic viscosity μ are estimated by the empirical formula in the physical property handbook, and the pressure P 5,0 is calculated using P hydro,i is the static pressure of the ith layer (i = 1, 2, 3, 4). The bubble point BP and latent heat of vaporization ΔH vap are calculated using the Peng-Robinson equation. Read P V = 116 kPa from the on-site pressure sensor. The calculation results are shown in Table 5:

[0185] Table 5 Calculation results of physical property parameters of grid G in the central region of the fifth layer 5,0

[0186]

[0187] (5.4) Quantify the overall stratification degree and local stratification instability of the storage tank according to the physical property field output in (5.3), and calculate the current stratification degree of the storage tank with a comprehensive index;

[0188] ​According to the calculated physical property distribution, the arithmetic mean pressure P of the liquid phase region in the storage tank is 162.52 kPa, and the saturation pressure P sat = 167.24 kPa. The overall stratification tendency is quantified as follows:

[0189]

[0190] Calculate the stratification stability ratio R value between each horizontal layer to quantify the local stratification tendency in the liquid phase region of the storage tank. Taking the local stratification tendency between the first layer and the second layer at the current time point as an example, the derivative is obtained according to the density empirical formula. The density empirical formula is as follows:

[0191]

[0192] The result of taking the partial derivative of the mixture density with respect to temperature is as follows:

[0193]

[0194] The result of taking the partial derivative of the mixture density with respect to the mole fraction of methane component is as follows:

[0195]

[0196] The calculation result of the R value between the first layer and the second layer is as follows:

[0197]

[0198] R1 satisfies |R1| < 1 and R1 > 0, and stratification occurs between the first layer and the second layer. Calculate the R k value in turn. When (R k < 0, β T (k)ΔT(k) > 0) or (|R k | < 1, R k > 0) or (|R k | > 1, R k > 0), stratification occurs at this liquid level; otherwise, stratification does not occur at this liquid level;

[0199] Taking the comprehensive stratification coefficient I at the current time point as an example, the calculation result is as follows:

[0200]

[0201] The comprehensive stratification coefficient at the current time point is within the safe range. When I > 5, on-site personnel should promptly adjust the operation parameters to suppress stratification;

[0202] S6. Construct a dynamic model of the LNG storage tank to predict the future temperature field, physical property field, and stratification coefficient trend; the process is as follows:

[0203] (6.1) Taking the temperature and molar composition of each region output by S5 as the initial values, and using the physical properties of each region output by S5 as the basis for calculating the heat terms in the dynamic model, a dynamic model of the LNG storage tank is constructed, and the static parameters and dynamic parameters of the site are input. The static parameters include the storage tank diameter D, the height H of the liquid phase region of the storage tank L , the tank wall material, the tank wall thickness δ b , the tank bottom material, the tank bottom thickness δ B , the circulation pipeline diameter D loop , the circulation pipeline length L loop , the installation height H of the feed pipeline in , the installation height H of the discharge pipeline out , and the dynamic parameters include the feed rate the discharge rate the circulation rate the gas phase temperature T V , the ambient temperature T a , the tank bottom temperature T B , * represents the molar rate. The MESH equation is used to establish an interlayer dynamic interaction model of the LNG storage tank. It is assumed that the temperature, molar composition and physical properties of each layer of the dynamic model are calculated based on the value of the central region G of each layer k,0 . The total time step is set to 720 h, and the unit time step is 10 min. The specific dynamic models of each layer are as follows:

[0204] The dynamic model of the general layer is as follows:

[0205] Material balance equation:

[0206]

[0207] Heat balance equation:

[0208]

[0209] Mole fraction summation formula:

[0210]

[0211] According to the calculation results at the current time point, the bubble point of the first layer of the liquid phase region is 111.7 K, satisfying T 1,0 ≥BP 1,0 . It is necessary to consider the evaporation term in (6.1.3). The installation height of the discharge pipeline is 15 m, corresponding to the sixth layer of the liquid phase region. Therefore, the temperature and composition of the circulating liquid corresponding to the sixth layer are taken. The constant F = 25 and the gravitational acceleration g = 9.81 m / s 2 , the original LNG density ρ h = ρ1 = 445.71 kg / m 3 , the density difference Δρ c= ρ6 - ρ1 = 0.25 kg / m 3 , the diameter D of the impact area where the circulating liquid contacts the original liquid is 0.3 m. The calculated LNG internal convection velocity is v = 0.6 m / s. The viscous stress of LNG is small, and the viscous stress term can be ignored, satisfying Fρ h v 2 ≥ Δρ c gD. The circulating liquid mixes with the first-layer liquid phase region, and the circulating liquid term in (6.1.4) needs to be considered;

[0212] The dynamic model of the first layer is as follows:

[0213] Material balance equation:

[0214]

[0215] Heat balance equation:

[0216]

[0217] Mole fraction summation formula:

[0218]

[0219] Phase equilibrium equation:

[0220] y i,E = K i x i,1,0

[0221] Considering the lag caused by the time when the recirculating liquid returns to the first layer through the pipeline, the lag equation is as follows:

[0222]

[0223] The recirculation time calculation formula is as follows:

[0224]

[0225] Taking the current time point as an example, the recirculation time is:

[0226]

[0227] In the dynamic model calculation of the horizontal layer where the inlet and outlet pipelines are located, the composition and heat of the inlet and outlet liquids are updated to the mass transfer and heat transfer calculations of the dynamic model of this layer respectively. The corresponding liquid level heights of the inlet and outlet pipelines under this condition are 15 m, both in the sixth layer. Then the dynamic model of the sixth layer is as follows:

[0228] Material balance equation:

[0229]

[0230] Heat balance equation:

[0231]

[0232] Mole fraction summation formula:

[0233]

[0234] Based on the dynamic models of each layer, use programming software for simulation to predict the temperature and composition change trends at each time point. According to the temperature and composition of each layer corresponding to each time point, taking the dynamic calculation results of the first 24 hours of the fifth layer as an example, the calculation results are respectively as Figure 3 and Figure 4 shown;

[0235] (6.2) Set up a calibration trigger mechanism for the temperature prediction values of each region, and calibrate once every 24 hours. The calibration trigger mechanism is as follows:

[0236] Assume that there are p unit time steps in this period, and the temperature corresponding to each unit time step of grid G k,f is T q,k,f , then the average temperature T k,f of grid G avg,k,f in this period is:

[0237]

[0238] Taking the dynamic calculation results of the first 24 hours of the horizontal fifth layer as an example, there are a total of unit time steps in this period, and the average temperature T avg,5,0 of this period is:

[0239]

[0240] Collect the latest historical data of the current storage tank Compare the measured values and calculated values in turn. Taking the central area of the fifth layer as an example, it is calculated that:

[0241]

[0242] Trigger calibration, repeat step S4 for parameter estimation, and the parameter estimation results are as follows:

[0243] θ = [0.13, 0.54, -0.28, 0.2, -0.1]

[0244] (6.3) Repeat step S5, substitute the temperatures and molar fractions in each layer corresponding to each time node within the total time step into the empirical formula, and output the physical property distributions and stratification coefficients at each time node. Taking the density distribution at each time point in the first 24 hours as an example, for the central region G of the first horizontal layer 1,0 the dynamic trend of the density is as Figure 5 shown. For the boundary region G of the fifth horizontal layer 5,f (f ≠ 0), the dynamic trend of the density is as Figure 6 shown. Taking the calculation results of the dynamic stratification coefficients in the first 24 hours corresponding to this working condition as an example, the dynamic trend of the stratification coefficient is as Figure 7 shown. The predicted I values within the total time step are all lower than 5, indicating that the current operating parameters are reasonable and the stratification risk is low. If the predicted index value continuously shows I > 5, measures need to be taken in advance to prevent stratification;

[0245] The content described in the embodiments of this specification is only a list of the implementation forms of the inventive concept and is only for illustrative purposes. The protection scope of the present invention should not be regarded as limited to the specific forms stated in this embodiment. The protection scope of the present invention also extends to equivalent technical means that can be conceived by those of ordinary skill in the art based on the inventive concept of the present invention.

Claims

1. A method for predicting the liquid-phase physical property field and stratification in a storage tank, characterized in that, The method includes the following steps: S1. Complete the radial and axial meshing of the LNG storage tank according to the rules; S2. Collect the historical data of other storage tanks, generate a preset model, and generate the initial temperature values of the grids inside the storage tank; S3. Construct the temperature distribution models of the first layer and the last layer, correct the axial temperature of the storage tank, construct the boundary layer temperature distribution model, and correct the radial temperature of the storage tank; S4. Compare the historical data of the current storage tank with the calculated values of the model, and perform parameter estimation on the preset model; S5. Convert the calculated temperature field into the corresponding physical property field and characterize it with a layering coefficient; S6. Construct a dynamic model of the LNG storage tank to predict the future trends of the temperature field, physical property field, and layering coefficient; 2. The method for predicting the liquid-phase physical property field and stratification in a storage tank according to claim 1, characterized in that The process of S1 is as follows: (1.1) Axial grid division: The liquid phase region of a vertical super-large liquid storage tank is evenly divided into n layers along the vertical direction, numbered sequentially from top to bottom, with the number represented by k, where k = 1, 2, …, n. The number of layers n is determined by the following formula, where H L is the total height of the liquid phase region in the LNG storage tank, with the unit of m: (1.2) Radial grid division: Define the boundary region of the storage tank as the region where the straight-line distance from the inner wall of the tank in each horizontal layer is not greater than 5% of the diameter D of the storage tank. The central region is the remaining part of each horizontal region except the boundary region, that is, the region where the distance from the inner wall of the tank is > 5% D, and no radial grid division is performed in the central region. The boundary region is non-uniformly divided into 5 layers in the horizontal direction and numbered in sequence from the central axis to the tank wall direction. The numbering is represented by f, where f = 0, 1, …, 5. When f = 0, it represents the central region of this layer. When f = 1, 2, …, 5, it represents the boundary region of this layer. The distance formula of the f-th layer of the boundary region from the tank wall is shown in Equation (2): The thickness formula of the f-th layer of the boundary region is shown in Equation (3): Define the radial grid label of the k-th layer as G k,f 。 3. A method for predicting the liquid-phase physical property field and stratification in a storage tank according to claim 1 or 2, characterized in that The process of S2 is as follows: Collect historical data of other storage tanks. Each group of data includes the temperature T inside the storage tank, the liquid level height H corresponding to this temperature, the feed temperature T in , the ambient temperature T a , the heat resistance R of the tank wall b . Through regression analysis, a temperature generation preset model for the LNG storage tank is generated, as shown in Equation (4): Among which Z norm is the set of normalized independent variables of the temperature calculation preset model, Z = [H norm , T in,norm , T a,norm , R b,norm , and the normalization formula for each independent variable is shown in Equation (5). Define the parameter set as θ, θ = [0.15, 0.6, -0.3, 0.2, -0.1]; Substitute the independent variable values collected on-site into the calculation model T model (Z norm ) to generate the initial temperature values of the central regions of each layer 4. A method for predicting the liquid-phase physical property field and stratification in a storage tank according to claim 1 or 2, characterized in that, The process of S3 is as follows: (3.1) Axial temperature correction model: Correct the temperature of the first layer and read the gas phase temperature value T at the site V . Take T V and to perform logarithmic mean calculation, and use the calculation result as the temperature value of the corrected first layer. Correct the temperature of the nth layer and take the bottom temperature of the tank T B and to perform logarithmic mean calculation, and use the calculation result as the temperature value of the corrected nth layer; (3.2) Radial temperature correction: Input the tank wall material and the tank wall thickness δ b , and calculate the temperature T of each inner wall surface according to Fourier's law k,b . For the radial grid temperature in the boundary region, the smoothed particle hydrodynamics method is used for temperature correction. Using T k,0 and T k,b as the inner and outer boundary values respectively, calculate the grid temperature in the boundary region and update it to the temperature field. Among them, the smoothed particle hydrodynamics method is a meshless numerical calculation method. By discretizing the continuous medium into a series of particles, each particle carries physical quantities such as mass, density, and temperature, and a kernel function is used to interpolate and approximately calculate the physical quantities between particles.

5. A method for predicting the liquid-phase physical property field and stratification in a storage tank according to claim 1 or 2, characterized in that, The process of S4 is as follows: Collect N groups of historical measurement data of the current storage tank, and each group of data includes the historical measured temperature value of the current storage tank and the corresponding historical measured liquid level H o 、feed temperature T in,o 、ambient temperature T a,o and the tank wall thermal resistance R b ,where R b is regarded as a constant value. Through the grid division rule of step S1 and the preset model T model (Z norm ) calculate the initial value, and then combine the axial temperature correction model and the radial temperature correction model of S3 to generate the comprehensive simulated temperature value T k,f,o of each grid. In this step, the parameter estimation is only for the linear parameter θ = [θ1, θ2, θ3, θ4, θ5] of T model (Z norm ) for optimization, and the optimization proposition is as follows: The parameter θ optimized by the least squares method will be updated to T model (Z norm ) In this case, the correction model in step S3 is a fixed algorithm and does not participate in parameter iteration.

6. The method for predicting the liquid-phase physical property field and stratification in a storage tank according to claim 1 or 2, characterized in that, The specific process of S5 is as follows: (5.1) Input feed composition x in , assuming that the molar composition of each layer in the liquid phase region of the LNG storage tank is equal at the current time period and is the same as the molar composition of the feed, i.e., x 1,f,t＝0 = x 2,f,t＝0 =... = x n,f,t＝0 = x in ; (5.2) Input the Z value of the current time period, and calculate the current temperature field of the LNG storage tank, i.e., the temperature value T of each grid, through T model (Z norm ), the axial correction model and the radial correction model k,f ; (5.3) Substitute the temperature and molar composition of each layer into the empirical formula to output the physical property field of the LNG storage tank, that is, the physical property parameters of each grid. The physical property parameters include density ρ, pressure P, bubble point BP, latent heat of vaporization ΔH vap , thermal conductivity λ, and dynamic viscosity μ; (5.4) According to the physical property field output in (5.3), quantify the overall layering degree and local layering instability of the storage tank, and calculate the current coefficient of the storage tank with a comprehensive index.

7. A method for predicting the liquid-phase physical property field and stratification in a storage tank according to claim 1 or 2, characterized in that, The process of S6 is as follows: (6.1) Using the temperature and molar composition of each region output by S5 as the initial values, and the physical properties of each region output by S5 as the basis for calculating the heat terms in the dynamic model, a dynamic model of the LNG storage tank is constructed. Input the static parameters and dynamic parameters of the site. The static parameters include the storage tank diameter D, the height H of the liquid phase region of the storage tank L , the tank wall material, the tank wall thickness δ b , the tank bottom material, the tank bottom thickness δ B , the circulation pipeline diameter D loop , the circulation pipeline length L loop , the installation height H of the feed pipeline in , the installation height H of the discharge pipeline out , and the dynamic parameters include the feed rate the discharge rate the circulation rate the gas phase temperature T V , the ambient temperature T a , the tank bottom temperature T B , * represents the molar rate. The MESH equation is a mathematical model for the equilibrium stage separation process, established based on mass conservation, phase equilibrium relationships, heat conservation, and the summation relationship of molar fractions. The MESH equation is used to establish the interlayer dynamic interaction model of the LNG storage tank. It is assumed that the temperature, molar composition, and physical properties of each layer of the dynamic model are calculated based on the values of the central region G k,0 of each layer. Set the total time step and the unit time step, which are used to predict the dynamic temperature and dynamic molar composition of each layer in the LNG storage tank, and evaluate the stratification trend based on the prediction results; (6.2) Set up a calibration trigger mechanism for the predicted values of temperature and molar composition in each region, calibrate once every 24 hours. The calibration method is to collect the latest historical data of the current storage tank, re-estimate the parameters of the preset model, and repeat the above calculations until the calibration is no longer triggered; (6.3) Repeat step S5 to substitute the temperatures T within each layer corresponding to each time node within the total time step p,k,0 and the molar composition x p,k,0 into the empirical formula, and output ρ p,k,0 , P p,k,0 , BP p,k,0 , ΔH vap,p,k,0 , λ p,k,0 , μ p,k,0 and the layering coefficient I p,k,0 . If the values of the layering coefficient at each time node continuously show I > 5, measures need to be taken in advance to prevent layering.

Images