Tank leakage risk prediction method and system based on digital twinning
By using digital twin technology, combined with acoustic path analysis and nonlinear optimization, the problem of inaccurate location of corrosion defects in the tank bottom plate caused by sparse sensor arrangement was solved, enabling accurate assessment of tank bottom plate corrosion and dynamic prediction of leakage risk.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- HAOHE JINYANG (BEIJING) TECHNOLOGY CO LTD
- Filing Date
- 2026-04-03
- Publication Date
- 2026-07-03
Smart Images
Figure CN122333869A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of digital twin technology, specifically to a method and system for predicting tank leakage risks based on digital twins. Background Technology
[0002] Large vertical storage tanks are key equipment in the petrochemical industry for storing crude oil and refined oil products. Their bottom plates are in direct contact with the deposited media and the soil at the bottom of the tank, and are subjected to the combined effects of media corrosion and soil electrochemical corrosion over long periods. This makes them prone to hidden localized thinning and even perforation leaks, affecting the safe operation of the equipment. Because the bottom plate is covered by the media and can reach diameters of tens to hundreds of meters, traditional offline tank inspection requires cleaning and production shutdown, resulting in high inspection costs and disruption to production. Therefore, online bottom plate monitoring technology based on acoustic emission and ultrasonic guided waves has become the mainstream technology direction in the industry.
[0003] In practical engineering applications, online monitoring of the bottom plates of large open-air storage tanks is limited by installation costs and explosion-proof safety regulations. Sensor arrays can only be sparsely arranged along the outer edge of the tank, and the parameters of the dense internal mesh are inverted using limited path integral data from the edges, resulting in a serious underdetermined problem. Conventional algebraic reconstruction techniques (ART) and total variational (TV) regularization methods, when processing the collected sparse data, tend to evenly distribute single-path measurements along the path, producing strip-shaped artifacts radiating from the edges to the center. This makes it difficult to accurately distinguish between real corrosion areas and algorithm noise, and to quantify the true remaining load-bearing capacity of the bottom plate, thus hindering the engineering-scale application of online monitoring technology. Summary of the Invention
[0004] This application provides a method and system for predicting tank leakage risks based on digital twins, in order to solve the problem of insufficient accuracy in locating corrosion defects in the tank bottom plate caused by strip-shaped artifacts resulting from sparse sensor arrangement. The specific technical solution adopted is as follows: In a first aspect, one embodiment of this application provides a method for predicting tank leakage risk based on digital twins, the method comprising the following steps: Based on the liquid level height of the storage tank at the time of sampling, the sampling window is divided, and the liquid level height, sound wave path, sound time delay, and sound emission energy sequence of each grid of the tank bottom plate at each sampling time within the divided sampling window are collected. Based on the relationship between the time interval of different acquisition times, the difference in liquid level height and the acoustic time lag at the corresponding acquisition time, the acoustic transmission time difference-liquid level change rate of each acoustic path is calculated, and an acoustic transmission time difference-liquid level change rate vector is established. Based on the correlation between the acoustic emission energy sequence of the grid and the sequence composed of liquid level height, the acoustic emission activity index of the grid is assigned a value, and an acoustic emission feature vector is established. Establish the geometric mapping matrix and calculate the spatial smoothing constraint coefficient of the mesh. The spatial smoothing constraint coefficient is used to limit the degree of smoothing constraint on the mesh during the inversion process. Based on the theoretical response of the relative thickness of the bottom plate of the grid, the geometric mapping matrix, the sound transmission time difference-liquid level change rate vector, and the spatial smoothing constraint coefficient, a nonlinear optimization objective function is established, and the nonlinear optimization objective function is iteratively solved to obtain the relative optimal thickness vector of the bottom plate. Based on the relative optimal thickness of the base plate and the acoustic emission activity index, the actual physical thickness and strength parameters of the grid are calculated respectively. Different virtual liquid level heights are generated according to the liquid level height of the tank at the time of acquisition. Combined with the strength parameters, the structural yield risk of the grid at each virtual liquid level height is calculated. Based on all structural yield risk values, it is determined whether there is a leakage risk in the tank.
[0005] Furthermore, the method for dividing the acquisition window is as follows: When the rate of change of liquid level at the sampling moment is greater than the threshold of the rate of change of liquid level, the sampling moment is taken as the start moment of the sampling window, and the last sampling moment among the consecutive sampling moments with the rate of change of liquid level greater than the threshold of the rate of change of liquid level is taken as the end moment of the sampling window.
[0006] Furthermore, the specific calculation method for the acoustic transmission time difference-liquid level change rate of the sound wave path is as follows: For any acoustic path at any acquisition time within the acquisition window, the difference between the liquid level height at the acquisition time and the liquid level height at the start time of the acquisition window is denoted as the liquid level difference at the acquisition time. The liquid level difference at the acquisition time and the time interval between the acquisition time and the start time of the acquisition window are used as independent variables, and the acoustic time lag at the acquisition time is used as the dependent variable to establish a binary linear regression model. The coefficient values obtained by binary linear regression of the liquid level difference at the acquisition time are denoted as the acoustic transmission time difference - liquid level change rate of the corresponding acoustic wave path.
[0007] Furthermore, the specific method for determining the acoustic emission activity index of the grid is as follows: The Pearson correlation coefficient between the acoustic emission energy sequence of the grid and the liquid level sequence composed of the liquid level height is denoted as the acoustic emission activity correlation of the grid. When the acoustic emission activity correlation of the grid is less than the preset activity correlation threshold, the acoustic emission activity index of the grid is assigned the value of 0; otherwise, the normalized value of the cumulative energy value of the grid in the acquisition window is recorded as the acoustic emission activity index of the grid.
[0008] Furthermore, the specific method for establishing the geometric mapping matrix is as follows: Calculate the grid that each sound wave path passes through under the assumption of ideal straight-line propagation, and establish a geometric mapping matrix. The geometric mapping matrix contains the... line, number The elements of the column represent the first... The path passes through the first The length of the line segments of each grid, where, It means greater than or equal to 1 and less than or equal to 1. integers, It means greater than or equal to 1 and less than or equal to 1. integers, Indicates the number of sound wave paths. Indicates the number of grid cells.
[0009] Furthermore, the specific calculation method for the spatial smoothing constraint coefficient of the mesh is as follows: A geometric intersection iterative algorithm is used on the sound transmission time difference-liquid level change rate vector to obtain the score of each grid. The normalized grid score vector is denoted as the ray intersection confidence. The negative correlation processing result of the normalized grid score is denoted as the spatial smoothing constraint coefficient of the grid.
[0010] Furthermore, the actual physical thickness and strength parameters of the mesh are calculated as follows: The actual physical thickness of the mesh is the product of the original design thickness of the mesh and the relative optimal thickness of the base plate. The product of the acoustic emission activity index of the mesh and the preset strength reduction factor is recorded as the mesh correction product. The product of the difference between the number 1 and the mesh correction product and the standard yield strength of the mesh material is recorded as the mesh strength parameter.
[0011] Furthermore, the method for determining the structural yield risk of the grid at each virtual liquid level height is as follows: For each virtual liquid level height, the static pressure applied to the bottom plate corresponding to the virtual liquid level height is calculated according to the liquid static pressure formula. The equivalent stress distribution of each grid on the tank bottom plate is calculated. The ratio of the equivalent stress distribution of the grid to the strength parameter is recorded as the structural yield risk degree of the grid at the corresponding virtual liquid level height.
[0012] Furthermore, the specific steps for determining whether a storage tank has a leakage risk based on the yield risk of all structures are as follows: When the structural yield risk is greater than or equal to the preset safety warning threshold, the minimum value of all virtual liquid level heights corresponding to the structural yield risk that is greater than or equal to the preset safety warning threshold is recorded as the upper limit of the safe liquid level. When the liquid level height of the storage tank is greater than or equal to the upper limit of the safe liquid level, it is determined that the storage tank has a leakage risk. When the liquid level height of the storage tank is less than the upper limit of the safe liquid level, it is determined that the storage tank does not have a leakage risk. Conversely, it is determined that there is no risk of leakage from the storage tank.
[0013] Secondly, embodiments of this application also provide a tank leakage risk prediction system based on digital twins, including a memory, a processor, and a computer program stored in the memory and running on the processor, wherein the processor executes the computer program to implement the steps of any of the methods described above.
[0014] The beneficial effects of this application are: This application first divides the acquisition window based on the liquid level height of the storage tank at the acquisition time to ensure that the data acquired during the tank leakage risk prediction process mainly reflects changes in liquid level load rather than ambient temperature drift. It calculates the acoustic time lag caused by a unit change in liquid level, obtaining the acoustic transmission time lag minus the liquid level change rate. The larger the acoustic transmission time lag minus the liquid level change rate, the more severe the deformation of the bottom plate area traversed by the corresponding acoustic wave path under pressure, i.e., the stronger the acoustoelastic effect. The bottom plate area traversed by the acoustic wave path is more likely to correspond to a corrosion-thinned region. Based on the correlation between the acoustic emission energy sequence of the grid and the sequence composed of liquid level height, an acoustic emission feature vector is established. Utilizing the topological prior that real defects must be located at the geometric intersection of multiple high-response paths, and based on the acoustic wave path on an ideal straight line... Under the propagation assumption, a geometric mapping matrix is established for the mesh through which the sound wave passes. This matrix guides the subsequent inversion solution. Based on the possibility that the high response of the sound wave path is caused by real local defects, the spatial smoothness constraint coefficient of the mesh is calculated to limit the smoothness constraint on the mesh during the inversion process. Then, a nonlinear optimization objective function is established and iteratively solved to obtain the relatively optimal thickness vector of the bottom plate. Finally, using digital twin technology, geometric weak points and active damage points are mapped onto the finite element model. Through virtual working condition simulation, abstract monitoring data is transformed into quantitative indicators to guide production, obtaining the results of tank leakage risk prediction and solving the problem of insufficient accuracy in locating corrosion defects in the tank bottom plate caused by the strip-shaped artifacts generated by sparse sensor arrangement. Attached Figure Description
[0015] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0016] Figure 1 This is a flowchart illustrating a digital twin-based method for predicting tank leakage risks, provided as an embodiment of this application. Detailed Implementation
[0017] Next, the technical solutions in the embodiments of the present application will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present application without creative efforts shall fall within the protection scope of the present application.
[0018] Please refer to Figure 1 , which shows a flowchart of a method for predicting the leakage risk of a storage tank based on digital twin provided by an embodiment of the present application. The method includes the following steps: Step S001, according to the liquid level height of the storage tank at the acquisition moment, divide the acquisition window, and acquire the liquid level height, acoustic wave path, acoustic time lag amount, and acoustic emission energy sequence of each grid of the storage tank bottom plate within the divided acquisition window at each acquisition moment.
[0019] Real-time collect the data stream of the liquid level gauge of the storage tank, and continuously calculate the liquid level change rate at each acquisition moment.
[0020] Among them, the data stream of the liquid level gauge refers to a continuous data set related to the liquid level in the storage tank that is continuously output by the liquid level gauge supporting the storage tank and is accessed by the system in real time, including the liquid level height and the corresponding acquisition moment. The liquid level change rate is the ratio of the liquid level height to the time interval of adjacent acquisition moments.
[0021] In order to ensure that the data collected during the process of predicting the leakage risk of the storage tank mainly reflects the change of liquid level load rather than the environmental temperature drift, a preset liquid level change rate threshold is set. When the liquid level change rate at the acquisition moment is greater than the liquid level change rate threshold, it is determined that the storage tank enters the rapid feeding and discharging state, and the acquisition moment is used as the starting moment of the acquisition window. The determination of the acquisition window stops until the liquid level change rate at the acquisition moment is less than the liquid level change rate threshold, and the last acquisition moment with a liquid level change rate greater than the liquid level change rate threshold is used as the end moment of the acquisition window to determine an effective acquisition window.
[0022] Among them, in order to prevent the cumulative temperature drift from being too large due to too long acquisition time, the maximum duration of the acquisition window in this embodiment is set to 2 hours; the value of the liquid level change rate threshold in this embodiment is 0.5 m / h.
[0023] The determination of the acquisition window can utilize the working condition characteristics of short-term high flow rate to create a prerequisite for subsequent separation of temperature interference.
[0024] Collect the liquid level height at each acquisition moment within the acquisition window, and establish a liquid level sequence of the acquisition window.
[0025] Within the acquisition window, at each acquisition time, the acoustic wave path formed by each pair of transmitting and receiving sensors is acquired through the ultrasonic sensor array at the edge of the tank. The time difference between the arrival of the acoustic wave signal at the acquisition time and the arrival of the acoustic wave signal at the start time of the acquisition window is recorded as the acoustic time lag at the acquisition time. Simultaneously, the impact event data stream captured by all acoustic emission sensors within the acquisition window is intercepted. Using the TDOA time difference of arrival localization algorithm, the acquired acoustic emission events are mapped onto the discrete grid of the tank bottom plate to obtain the acoustic emission energy sequence of each grid on the tank bottom plate at each acquisition time within the acquisition window.
[0026] Among them, the ultrasonic sensor can be a PZT piezoelectric ceramic probe or an EMAT electromagnetic ultrasonic probe, as long as it can excite and receive plate waves; the acoustic emission sensor can be a resonant or broadband acoustic emission sensor, and the frequency band is recommended to cover 30kHz-150kHz to adapt to the crack signal characteristics of metal storage tanks.
[0027] The method for obtaining the discrete grid of the tank bottom plate is as follows: based on the geometric dimensions of the tank bottom plate, the tank bottom plate is divided into a uniform discrete grid. The number of discrete grids is determined by those skilled in the art; in this embodiment, the number of discrete grids is set to 3000. The time difference of arrival (TDOA) positioning algorithm is a well-known technology and will not be described in detail here.
[0028] At this point, the liquid level height, acoustic path, acoustic time delay, and acoustic emission energy sequence of each grid at each acquisition moment within the acquisition window are obtained.
[0029] Step S002: Based on the relationship between the time interval of different acquisition times, the difference in liquid level height and the acoustic time lag at the corresponding acquisition time, calculate the acoustic transmission time difference-liquid level change rate of each acoustic path, establish the acoustic transmission time difference-liquid level change rate vector, assign a value to the acoustic emission activity index of the grid based on the correlation between the acoustic emission energy sequence of the grid and the sequence composed of liquid level height, and establish the acoustic emission feature vector.
[0030] For any acoustic path at any acquisition time within the acquisition window, the difference between the liquid level height at the acquisition time and the liquid level height at the start time of the acquisition window is denoted as the liquid level difference at the acquisition time. The liquid level difference at the acquisition time and the time interval between the acquisition time and the start time of the acquisition window are used as independent variables, and the acoustic time lag at the acquisition time is used as the dependent variable. A binary linear regression model is established, and the least squares method is used to solve the binary linear regression model. The coefficient values obtained by the binary linear regression of the liquid level difference at the acquisition time are denoted as the acoustic transmission time difference-liquid level change rate of the corresponding acoustic path. The acoustic transmission time difference-liquid level change rate corresponds to the structural stress response characteristics of the acoustic path.
[0031] Since the propagation speed of sound waves is affected by both the stress changes in the bottom plate caused by the hydrostatic pressure of the liquid level and the thermal expansion and contraction of the material caused by changes in ambient temperature, directly ignoring the temperature term would lead to calculation errors. Furthermore, the stress changes in the bottom plate are related to the liquid level height, and the thermal expansion and contraction of the material drift slowly over time. Therefore, the binary linear regression model is as follows: In the formula, Indicates the sound wave path In the Each sampling time The acoustic time lag; Indicates the first Each sampling time The liquid level height; The liquid level height at the start of the acquisition window; Indicates the first sampling time; Indicates the sound wave path The sound transmission time difference minus the liquid level change rate, i.e., the sound wave path Obtained by binary linear regression The coefficient; Indicates the sound wave path The temperature drift coefficient, i.e., the sound wave path Obtained by binary linear regression The coefficient; This represents the constant intercept term in a binary linear regression model, which is used to absorb the fixed bias of the system. This represents the random observation error of the binary linear regression model.
[0032] The acoustic transmission time difference-liquid level change rate represents the acoustic time lag caused by a unit change in liquid level. The larger the acoustic transmission time difference-liquid level change rate, the more severe the deformation of the bottom plate area traversed by the corresponding acoustic wave path under pressure, that is, the stronger the acoustoelastic effect, and the more likely the bottom plate area traversed by the acoustic wave path is to correspond to a corrosion thinning area.
[0033] The temperature drift coefficient represents the amount of natural sound time drift caused by the passage of a unit of time. The temperature drift coefficient reflects the linear effect of changes in ambient temperature on the speed of sound.
[0034] It is understandable that the liquid level difference at the acquisition time is usually determined by the pumping power, exhibiting nonlinear or step characteristics. The time interval between the acquisition time and the start time of the acquisition window increases strictly linearly. The time interval and the liquid level difference have significant nonlinearity in the time domain. Therefore, by taking the liquid level difference at the acquisition time and the time interval between the acquisition time and the start time of the acquisition window as independent variables, and the acoustic time lag at the acquisition time as the dependent variable, the acoustic transmission time lag-liquid level change rate and temperature drift coefficient can be stably solved simultaneously using the least squares method.
[0035] The method of using least squares to solve a binary linear regression model is a well-known technique and will not be elaborated further.
[0036] Arrange the acoustic transmission time difference - liquid level change rate of all acoustic wave paths in sequence, and denote the obtained column vector as the acoustic transmission time difference - liquid level change rate vector.
[0037] Furthermore, the correlation between acoustic emission signals and liquid level load was analyzed.
[0038] The Pearson correlation coefficient between the cumulative acoustic emission energy sequence of the grid and the liquid level sequence composed of the liquid level height is denoted as the acoustic emission activity correlation of the grid. When the acoustic emission activity correlation of the grid is less than the preset activity correlation threshold, the acoustic emission activity index of the grid is assigned the value of 0. When the acoustic emission activity correlation of the grid is greater than or equal to the preset activity correlation threshold, the normalized value of the cumulative energy value of the grid in the acquisition window is denoted as the acoustic emission activity index of the grid.
[0039] In this embodiment, the threshold value for active correlation is set to 0.6; the calculation of the Pearson correlation coefficient is a well-known technique and will not be described in detail here. It should be noted that this embodiment uses the maximum-minimum normalization method to calculate the normalized value. In practical applications, implementers may use other methods of existing technology, such as the tanh function or the sigmoid function, to calculate the normalized value, and no limitation is made here.
[0040] Arrange the acoustic emission activity indices of the grid sequentially, and denote the resulting column vector as the acoustic emission feature vector.
[0041] Thus, the acoustic transmission time difference-liquid level change rate vector, acoustic emission feature vector, and the acoustic emission activity index of the grid are obtained.
[0042] Step S003: Establish the geometric mapping matrix and calculate the spatial smoothing constraint coefficient of the mesh. The spatial smoothing constraint coefficient is used to limit the degree of smoothing constraint on the mesh during the inversion process.
[0043] To address the inversion underdeterminacy caused by sparse sensor placement, it is necessary to utilize the topological prior that the real defect must be located at the geometric intersection of multiple high-response paths to construct a spatially adaptive regularized weight field to guide the subsequent inversion solution.
[0044] Based on the physical installation coordinates of the sensor array, the grid traversed by each sound wave path under the assumption of ideal straight-line propagation is calculated, and a grid of size is established. The geometric mapping matrix, where, Indicates the number of sound wave paths. Represents the number of meshes, the th in the geometry mapping matrix line, number Column elements Indicates the first The path passes through the first The length of each line segment in the grid is such that when the path does not cross a grid, the corresponding element in the geometric mapping matrix has a value of 0. It means greater than or equal to 1 and less than or equal to 1. integers, It means greater than or equal to 1 and less than or equal to 1. Integers.
[0045] To distinguish whether the high response of the acoustic path is caused by real local defects or by measurement noise, a geometric intersection iterative algorithm is used to process the acoustic transmission time difference-liquid level change rate vector. The confidence of each acoustic path and the score of each grid are obtained through iteration, and the normalized grid score vector is denoted as the ray intersection confidence.
[0046] The ray intersection confidence score is a column vector composed of the normalized scores of all grids. The larger the normalized score of a grid, the more likely the grid corresponds to the true defect center. The geometric intersection iterative algorithm is based on the ray acoustic approximation of high-frequency guided waves and can confirm the defect location through a multipath voting mechanism.
[0047] The negative correlation result of the normalized scores of the grid is denoted as the spatial smoothing constraint coefficient of the grid. The spatial smoothing constraint coefficients of all grids are arranged in order to obtain the spatial smoothing constraint coefficient vector.
[0048] It is understandable that the normalized score of the grid is negatively correlated, meaning that the normalized score of the grid is negatively correlated with the spatial smoothing constraint coefficient vector of the grid. It is also understood that the negative correlation in this application refers to the relationship between the independent and dependent variables, where the independent variable is the normalized score of the grid and the dependent variable is the spatial smoothing constraint coefficient vector of the grid. The negative correlation means that the dependent variable decreases (increases) as the independent variable increases (increases), and can be an inverse relationship, a subtraction relationship, etc.
[0049] Preferably, as an embodiment of this application, the negative of the product of the normalized score of the grid and a preset regularization adjustment factor is used as the exponent of an exponential function with the natural constant as the base, and the calculated value of the exponential function is recorded as the spatial smoothing constraint coefficient of the grid. In this embodiment, the regularization adjustment factor is set to 1.0.
[0050] The larger the normalized value of the mesh score, the greater the possibility of defects in the mesh. In this case, the smaller the spatial smoothing constraint coefficient of the mesh, the more relaxed the smoothing constraint should be in the subsequent inversion process, allowing more drastic changes in the mesh parameters, that is, allowing the existence of corrosion pits at the mesh location. Conversely, the smaller the normalized value of the mesh score, the smaller the possibility of defects in the mesh. In this case, the larger the spatial smoothing constraint coefficient of the mesh, the more strengthened the smoothing constraint should be in the subsequent inversion process, using the information of the surrounding mesh to smooth random noise, thereby suppressing artifacts.
[0051] At this point, the spatial smoothing constraint coefficients of the mesh are obtained.
[0052] Step S004: Based on the theoretical response of the relative thickness of the bottom plate of the mesh, the geometric mapping matrix, the sound transmission time difference-liquid level change rate vector, and the spatial smoothing constraint coefficient, establish a nonlinear optimization objective function, iteratively solve the nonlinear optimization objective function, and obtain the relative optimal thickness vector of the bottom plate.
[0053] Establish the following nonlinear optimization objective function The optimal relative thickness vector of the base plate is found by nonlinearly optimizing the objective function.
[0054] in, Indicates the first The relative thickness of the base plate of each grid, express -2; express The The theoretical response of acoustic transmission time difference - liquid level change rate for each grid; This represents a physical constant containing the material's acoustic elasticity and reference sound velocity, used to balance the relative thickness of the base plate. Its theoretical response The units between them are determined to ensure that the mapped predicted quantity has the same units as the observation vector; This represents the positive response operator, which is a column vector composed of the theoretical responses of the relative thicknesses of the base plates in each grid arranged sequentially. This represents the relative thickness vector of the base plates, which consists of the relative thicknesses of all the base plates. The initial value of the relative thickness vector of the base plates is a vector composed of the number 1. Represents the relative thickness vector of the base plate The corresponding nonlinear optimization objective function; Represents the geometric mapping matrix; This represents the sound transmission time difference minus the rate of change of liquid level vector; This represents the preset global regularization parameter; Represents a grid The thickness gradient at that point; Represents a grid The square of the magnitude of the thickness gradient at that point; Indicates the first Spatial smoothing constraint coefficients for each grid; Indicates the number of grid cells.
[0055] The global regularization parameter can be determined using the L-curve method based on the signal-to-noise ratio of the observed data. The global regularization parameter controls the overall smoothing strength, and its value should be greater than or equal to... and less than or equal to .
[0056] Spatial smoothing constraint coefficients are used for spatial weighting, providing stronger smoothing in non-defect areas and protecting edge features in defect areas.
[0057] Used to punish spatial mutations.
[0058] The nonlinear optimization objective function contains Therefore, the process of finding the optimal relative thickness vector of the base plate through the nonlinear optimization objective function is a nonconvex optimization problem. The nonlinear conjugate gradient method is selected to iteratively solve the nonlinear optimization objective function to obtain the optimal relative thickness vector of the base plate, that is, the relative optimal thickness vector of the base plate, which is composed of the relative optimal thicknesses of each base plate arranged in sequence.
[0059] During the iterative solution process, the value of the relative thickness vector of the base plate is updated along the conjugate direction. When the relative thickness of the base plate of the mesh is greater than the maximum relative thickness threshold, the relative thickness of the base plate of the mesh is assigned to the maximum relative thickness threshold. When the relative thickness of the base plate of the mesh is less than the preset minimum relative thickness threshold, the relative thickness of the base plate of the mesh is assigned to the minimum relative thickness threshold.
[0060] In this embodiment, the minimum relative thickness threshold and the maximum relative thickness threshold are set to 0.1 and 1, respectively; the function of the minimum relative thickness threshold and the maximum relative thickness threshold is to physically constrain the relative thickness of the base plate of the mesh.
[0061] The relative optimal thickness vector of the base plate directly reflects the degree of corrosion thinning of each grid.
[0062] At this point, the relative optimal thickness vector of the base plate is obtained.
[0063] Step S005: Based on the relative optimal thickness of the base plate and the acoustic emission activity index, calculate the actual physical thickness and strength parameters of the grid. Generate different virtual liquid level heights based on the liquid level height of the tank at the time of data acquisition. Combined with the strength parameters, calculate the structural yield risk of the grid at each virtual liquid level height. Based on all structural yield risk values, determine whether the tank has a leakage risk.
[0064] Finally, using digital twin technology, geometric weak points and active damage points are mapped onto the finite element model. Through virtual working condition simulation, abstract monitoring data is transformed into quantitative indicators that guide production—dynamic upper limit of safe liquid level.
[0065] The pre-established finite element model of the parameterized shell of the tank bottom plate is called. The mesh of the finite element model of the parameterized shell of the tank bottom plate is consistent with the inversion mesh. Then, the state of the finite element model of the parameterized shell of the tank bottom plate is updated using the relative optimal thickness vector of the bottom plate and the acoustic emission feature vector.
[0066] The specific steps for updating the state are as follows: 1. Geometric parameter update: The product of the original design thickness of the mesh and the relative optimal thickness of the base plate is denoted as the actual physical thickness of the mesh.
[0067] 2. Strength parameter correction: The product of the acoustic emission activity index of the mesh and the preset strength reduction factor is recorded as the mesh correction product. The product of the difference between the number 1 and the mesh correction product and the standard yield strength of the mesh material is recorded as the mesh strength parameter.
[0068] The calculation process of the actual physical thickness of the mesh is to map the corrosion thinning information obtained by inversion to the geometric properties of the model to realize the update of geometric parameters; the standard yield strength of the material is a known parameter, for example, the standard yield strength of Q235 steel is 235MPa; the value of the strength reduction factor should be greater than or equal to 0.1 and less than or equal to 0.3. The value range of the strength reduction factor corresponds to the 20%-30% safety factor margin commonly used in engineering design. In this embodiment, the value of the strength reduction factor is 0.2.
[0069] It is important to note that the calculation of the strength parameters of the mesh is not intended to accurately describe the physical degradation of the material's microscopic yield strength, but rather to construct an engineering assessment model based on risk aversion. In this model, the acoustic emission activity index of the mesh characterizes the crack activity within the mesh. By introducing a reduction term, the upper limit of the load-bearing capacity of highly active regions in the numerical simulation is artificially reduced, reserving a computational margin for potential failure risks. This forces the model to identify risks earlier in subsequent simulations, thereby providing more conservative and safer decision-making recommendations.
[0070] It is important to note that when the acoustic emission activity index of the mesh is 0, the strength parameter of the mesh is equal to the standard yield strength of the mesh material. In other words, the model degenerates into a conventional evaluation that only considers geometric thinning, which is logical.
[0071] After updating the state of the finite element model of the parameterized shell of the tank bottom plate, a virtual loading experiment is performed to predict the load-bearing capacity of the current structural state for future loads.
[0072] The specific steps of the virtual loading experiment are as follows: 1. Starting from the liquid level height of the storage tank at the time of data collection, the liquid level height is increased by a preset step size to generate a series of virtual liquid level heights until the maximum design liquid level of the storage tank is reached. In this embodiment, the preset step size is set to 0.1 meters.
[0073] 2. For each virtual liquid level height, calculate the static pressure applied to the bottom plate corresponding to the virtual liquid level height according to the liquid static pressure formula, call the nonlinear solver, and calculate the equivalent stress distribution of each grid of the tank bottom plate.
[0074] 3. For each virtual liquid level height, the ratio of the equivalent stress distribution of the grid to the strength parameter is denoted as the structural yield risk of the grid at the corresponding virtual liquid level height.
[0075] The structural yield risk directly reflects how close the mesh material is to the calculated failure point. Since the mesh strength parameters have been reduced according to the activity index, for meshes with active cracks, the risk may exceed the calculated failure point even if the equivalent stress distribution of the mesh has not yet reached the true physical limit of the material. This triggers an alarm.
[0076] Furthermore, the structural yield risk of all grids at all virtual liquid level heights is compared with a preset safety warning threshold. When a structural yield risk is greater than or equal to the preset safety warning threshold, the minimum value of all virtual liquid level heights corresponding to structural yield risk values greater than or equal to the preset safety warning threshold is recorded as the upper limit of the safe liquid level. When the liquid level of the storage tank is greater than or equal to the upper limit of the safe liquid level, it is determined that the storage tank has a leakage risk, triggering an audible and visual alarm and generating an operation command suggesting stopping feeding or suggesting lowering the liquid level. When the liquid level of the storage tank is less than the upper limit of the safe liquid level, it is determined that the storage tank has no leakage risk. When all structural yield risks are less than the preset safety warning thresholds, it is determined that the storage tank has no leakage risk.
[0077] This completes the prediction of storage tank leakage risks.
[0078] To further guide those skilled in the art in implementing this application, the following provides a detailed description of the key parameter configurations and initialization strategies involved. The selection of these parameters directly affects the convergence of the algorithm and the reliability of the prediction results.
[0079] In actual engineering deployments, the bottom plate dimensions, material properties, and environmental noise levels of different storage tanks may vary. This application fine-tunes the key algorithm parameters based on the actual situation. The key algorithm parameters include regularization adjustment factors, global regularization parameters, and intensity reduction factors. The recommended values and debugging methods for the key algorithm parameters are as follows.
[0080] Specifically, if there are still many artifacts in the inversion results, the value of the regularization adjustment factor can be appropriately increased. For example, the value of the regularization adjustment factor can be set to 1.5 to enhance the smoothing effect on non-intersection regions. If the inversion results lose small real defects, the value of the regularization adjustment factor can be appropriately decreased. For example, the value of the regularization adjustment factor can be set to 0.8.
[0081] Specifically, an excessively large global regularization parameter will result in overly smooth results, losing defect depth information; an excessively small global regularization parameter will result in unstable results, filled with noise. The global regularization parameter can be determined using the L-curve method based on the signal-to-noise ratio of the observed data. That is, a curve is plotted in a logarithmic coordinate system showing the variation of the data residual norm and the smoothness norm of the solution with the global regularization parameter, and the global regularization parameter value corresponding to the inflection point of the curve is selected as the optimal parameter.
[0082] Specifically, the strength reduction factor is used to control the conservatism of the safety strategy. The larger the strength reduction factor, the greater the strength reduction in the active defect area, the lower the output safe liquid level, and the more conservative the safety strategy. Therefore, the value of the strength reduction factor should refer to the safety factor in the tank design specification. If the tank stores high-risk media such as liquefied gas, it is recommended to take the upper limit of the strength reduction factor of 0.3. If the tank stores ordinary crude oil, it is recommended to take the lower limit of the strength reduction factor of 0.1.
[0083] Before conducting tank leakage risk prediction, initialization and cold start are required.
[0084] The initial value of the relative thickness vector of the bottom plate is a vector consisting entirely of the number 1, which means that it is assumed by default that the initial state of the tank bottom plate is intact, that is, the actual thickness of the tank bottom plate is equal to the design thickness. Then, the corrosion area is gradually discovered based on the observation data.
[0085] If no valid acoustic emission signal is detected within the acquisition window, the tank is in a low liquid level and static state, or the acoustic emission activity index of all grids after correlation screening is 0. The acoustic emission feature vector is assigned as a vector consisting entirely of 0s. In this case, the intensity parameters of the grids are not reduced to ensure that the prediction process can still operate normally when there is a lack of active source data and will not crash due to data loss.
[0086] In each ratio calculation process of this application, to avoid the denominator being zero, a preset value needs to be added to the denominator. The example of the preset value is [value to be filled in]. .
[0087] This completes the prediction of storage tank leakage risks.
[0088] Based on the same inventive concept as the above methods, this application also provides a tank leakage risk prediction system based on digital twins, including a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of any one of the above-described tank leakage risk prediction methods based on digital twins.
[0089] The above description is only a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the principles of this application should be included within the protection scope of this application.
Claims
1. A method for predicting tank leakage risk based on digital twins, characterized in that, The method includes the following steps: Based on the liquid level height of the storage tank at the time of sampling, the sampling window is divided, and the liquid level height, sound wave path, sound time delay, and sound emission energy sequence of each grid of the tank bottom plate at each sampling time within the divided sampling window are collected. Based on the relationship between the time interval of different acquisition times, the difference in liquid level height and the acoustic time lag at the corresponding acquisition time, the acoustic transmission time difference-liquid level change rate of each acoustic path is calculated, and an acoustic transmission time difference-liquid level change rate vector is established. Based on the correlation between the acoustic emission energy sequence of the grid and the sequence composed of liquid level height, the acoustic emission activity index of the grid is assigned a value, and an acoustic emission feature vector is established. Establish the geometric mapping matrix and calculate the spatial smoothing constraint coefficient of the mesh. The spatial smoothing constraint coefficient is used to limit the degree of smoothing constraint on the mesh during the inversion process. Based on the theoretical response of the relative thickness of the bottom plate of the grid, the geometric mapping matrix, the sound transmission time difference-liquid level change rate vector, and the spatial smoothing constraint coefficient, a nonlinear optimization objective function is established, and the nonlinear optimization objective function is iteratively solved to obtain the relative optimal thickness vector of the bottom plate. Based on the relative optimal thickness of the base plate and the acoustic emission activity index, the actual physical thickness and strength parameters of the grid are calculated respectively. Different virtual liquid level heights are generated according to the liquid level height of the tank at the time of acquisition. Combined with the strength parameters, the structural yield risk of the grid at each virtual liquid level height is calculated. Based on all structural yield risk values, it is determined whether there is a leakage risk in the tank.
2. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The method for dividing the acquisition window is as follows: When the rate of change of liquid level at the sampling moment is greater than the threshold of the rate of change of liquid level, the sampling moment is taken as the start moment of the sampling window, and the last sampling moment among the consecutive sampling moments with the rate of change of liquid level greater than the threshold of the rate of change of liquid level is taken as the end moment of the sampling window.
3. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The specific calculation method for the acoustic transmission time difference-liquid level change rate of the sound wave path is as follows: For any acoustic path at any acquisition time within the acquisition window, the difference between the liquid level height at the acquisition time and the liquid level height at the start time of the acquisition window is denoted as the liquid level difference at the acquisition time. The liquid level difference at the acquisition time and the time interval between the acquisition time and the start time of the acquisition window are used as independent variables, and the acoustic time lag at the acquisition time is used as the dependent variable to establish a binary linear regression model. The coefficient values obtained by binary linear regression of the liquid level difference at the acquisition time are denoted as the acoustic transmission time difference - liquid level change rate of the corresponding acoustic wave path.
4. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The specific method for determining the acoustic emission activity index of the grid is as follows: The Pearson correlation coefficient between the acoustic emission energy sequence of the grid and the liquid level sequence composed of the liquid level height is denoted as the acoustic emission activity correlation of the grid. When the acoustic emission activity correlation of the grid is less than the preset activity correlation threshold, the acoustic emission activity index of the grid is assigned the value 0. Conversely, the normalized value of the accumulated energy of the grid within the acquisition window is denoted as the acoustic emission activity index of the grid.
5. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The specific method for establishing the geometric mapping matrix is as follows: Calculate the grid that each sound wave path passes through under the assumption of ideal straight-line propagation, and establish a geometric mapping matrix. The geometric mapping matrix contains the... line, number The elements of the column represent the first... The path passes through the first The length of the line segments of each grid, where, It means greater than or equal to 1 and less than or equal to 1. integers, It means greater than or equal to 1 and less than or equal to 1. integers, Indicates the number of sound wave paths. Indicates the number of grid cells.
6. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The specific calculation method for the spatial smoothness constraint coefficient of the grid is as follows: A geometric intersection iterative algorithm is used on the sound transmission time difference-liquid level change rate vector to obtain the score of each grid. The normalized grid score vector is denoted as the ray intersection confidence. The negative correlation processing result of the normalized grid score is denoted as the spatial smoothing constraint coefficient of the grid.
7. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The actual physical thickness and strength parameters of the mesh are calculated as follows: The actual physical thickness of the mesh is the product of the original design thickness of the mesh and the relative optimal thickness of the base plate. The product of the acoustic emission activity index of the mesh and the preset strength reduction factor is recorded as the mesh correction product. The product of the difference between the number 1 and the mesh correction product and the standard yield strength of the mesh material is recorded as the mesh strength parameter.
8. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The method for determining the structural yield risk of the grid at each virtual liquid level height is as follows: For each virtual liquid level height, the static pressure applied to the bottom plate corresponding to the virtual liquid level height is calculated according to the liquid static pressure formula. The equivalent stress distribution of each grid on the tank bottom plate is calculated. The ratio of the equivalent stress distribution of the grid to the strength parameter is recorded as the structural yield risk degree of the grid at the corresponding virtual liquid level height.
9. The method for predicting tank leakage risk based on digital twins according to claim 1, characterized in that, The specific steps for determining whether a storage tank has a leakage risk based on the yield risk of all structures are as follows: When the structural yield risk is greater than or equal to the preset safety warning threshold, the minimum value of all virtual liquid level heights corresponding to the structural yield risk that is greater than or equal to the preset safety warning threshold is recorded as the upper limit of the safe liquid level. When the liquid level height of the storage tank is greater than or equal to the upper limit of the safe liquid level, it is determined that the storage tank has a leakage risk. When the liquid level height of the storage tank is less than the upper limit of the safe liquid level, it is determined that the storage tank does not have a leakage risk. Conversely, it is determined that there is no risk of leakage from the storage tank.
10. A digital twin-based tank leakage risk prediction system, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method as claimed in any one of claims 1-9.