Directly-buried compensator residual life prediction method based on digital twinning

Abstract

This invention relates to the field of digital twin technology, specifically a method for predicting the remaining service life of directly buried compensators based on digital twins. The method includes collecting pipeline geometric data, compensator structural data, and soil physical parameters of a directly buried heating pipeline segment; constructing a three-dimensional mesh model and calculating the soil spring stiffness matrix based on the soil physical parameters to build a baseline finite element twin model; acquiring cyclic hardening curves and temperature-related yield strength data to generate a physical field coupling model; acquiring operating temperature time series and internal pressure fluctuation sequences to form a thermo-mechanical coupling simulation model; driving the thermo-mechanical coupling simulation model to perform nonlinear iterative calculations, outputting equivalent stress response time history and plastic strain range data; and calculating and outputting the predicted remaining service life value by modifying the fatigue damage assessment model. This invention constructs a strongly coupled analysis system of structure-thermal-soil multi-physics fields, improving the comprehensive reliability of remaining service life prediction results under complex coupling environments.

Description

Technical Field

[0001] This invention relates to the field of digital twin technology, specifically to a method for predicting the remaining life of directly buried compensators based on digital twins. Background Technology

[0002] In the safe operation and maintenance of urban centralized heating pipe network systems, conducting full life-cycle management and service status assessment of corrugated compensators in directly buried environments is an important technical application direction for ensuring heating safety. Current technologies generally use standard formulas from the American Expansion Joint Manufacturers Association (AEMA) or estimate service life based on the manufacturer's rated design cycle count. This involves statistically analyzing the rated displacement, working pressure, and geometric parameters of the bellows, combined with empirical fatigue curves, to derive its theoretical fatigue life. This is currently the mainstream technical method for assessing the service status of compensators in engineering design and operation and maintenance.

[0003] However, existing technologies have significant shortcomings in prediction accuracy and adaptability to operating conditions. Their computational models for assessing lifespan suffer from idealization and staticity. Existing methods are largely based on ideal boundary conditions under overhead or simple support structures, neglecting the nonlinear constraints of the soil on the pipeline system in directly buried environments (such as soil friction and clamping effects). This leads to a severe disconnect between the mechanical analysis boundary and actual operating conditions. The system can only perform one-time theoretical estimations based on parameters from the design phase, failing to reflect the dynamic evolution of materials over service time, such as stiffness decay, cyclic hardening, and corrosion aging. Furthermore, the output of assessment results suffers from lag, meaning it cannot update remaining lifespan data based on actual operating temperature and internal pressure fluctuations. This results in maintenance strategies lacking specificity and timeliness, easily leading to pipeline failures or over-maintenance.

[0004] To address this, a method for predicting the remaining life of directly buried compensators based on digital twins is proposed. Summary of the Invention

[0005] The purpose of this invention is to provide a method for predicting the remaining life of directly buried compensators based on digital twins, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for predicting the remaining life of directly buried compensators based on digital twins, comprising: Collect pipeline geometry data, compensator structure data, and soil physical parameters of the directly buried heating pipeline network; construct a three-dimensional mesh model based on the pipeline geometry data and compensator structure data, and calculate the soil spring stiffness matrix based on the soil physical parameters; map the soil spring stiffness matrix to the three-dimensional mesh model to construct a reference finite element twin model; Obtain cyclic hardening curves and temperature-related yield strength data, and map them to the baseline finite element twin model to generate a physical field coupling model; The operating temperature time series and internal pressure fluctuation series are obtained and converted into transient temperature field loads and internal wall pressure loads, respectively. The transient temperature field loads and internal wall pressure loads are applied to the physical field coupling model to form a thermo-mechanical coupling simulation model. The thermo-coupling simulation model is driven to perform nonlinear iterative calculations and outputs the equivalent stress response time history and plastic strain range data of the peak and trough nodes; Input the equivalent stress response time history and plastic strain range data into the modified fatigue damage assessment model to calculate the current cumulative damage degree; based on the current cumulative damage degree and the preset failure threshold, calculate and output the predicted value of the remaining service life.

[0007] Preferably, the pipeline geometric data, compensator structural data, and soil physical parameters specifically include: the pipeline geometric data includes the remaining wall thickness distribution matrix and pipe body out-of-roundness geometric deviation data obtained by scanning; the compensator structural data includes the initial waveform vector of the corrugated pipe extracted from the installation completion digital archive and the interlayer corrosion depth data based on eddy current nondestructive testing inversion; the soil physical parameters include lateral earth pressure response data collected by an earth pressure cell array pre-embedded in the backfill soil of the trench, soil moisture content collected by a time domain reflectometer, and soil dynamic shear modulus and soil dynamic unit weight parameters based on soil moisture content mapping.

[0008] Preferably, the specific construction process of the three-dimensional mesh model and the soil spring stiffness matrix includes: calling the pipe wall remaining thickness distribution matrix and pipe body out-of-roundness geometric deviation data from the pipe geometry data, superimposing the pipe body out-of-roundness geometric deviation data as radial displacement constraints onto the standard circular pipe model to generate a non-ideal pipe geometry topology; generating a corrugated pipe multi-layer solid structure using the corrugated pipe initial waveform vector and interlayer corrosion depth data from the compensator structure data; generating a three-dimensional solid model based on the non-ideal pipe geometry topology and the corrugated pipe multi-layer solid structure, and performing a hexahedral element discretization operation on the three-dimensional solid model to generate the three-dimensional mesh model; extracting soil dynamic shear modulus and lateral earth pressure response data from soil physical parameters, introducing a depth correction coefficient based on the burial depth difference, and calculating the axial, lateral, and vertical stiffness components in combination with the soil dynamic shear modulus, and assembling the axial, lateral, and vertical stiffness components to generate the soil spring stiffness matrix.

[0009] Preferably, the specific construction process of the benchmark finite element twin model includes: identifying the node coordinate sequence of the outer wall of the three-dimensional mesh model and defining the node coordinate sequence as the interactive interface; obtaining the spring element topology corresponding to the soil spring stiffness matrix and defining the spring element topology as the constraint boundary; constructing a nonlinear frictional contact pair connecting the interactive interface and the constraint boundary based on Coulomb friction theory; associating the soil spring stiffness matrix with the nonlinear frictional contact pair to construct the benchmark finite element twin model.

[0010] Preferably, the specific generation process of the physical field coupling model includes: calling the multi-temperature gradient uniaxial tensile test data of the compensator material, and introducing the material aging reduction coefficient in combination with the interlayer corrosion depth data to correct the nonlinear decay curve of the yield strength with temperature, and using the correction result as the temperature-related yield strength data; analyzing the stable hysteresis loop data of the strain-controlled cyclic loading experiment, extracting the parameters of the hybrid strengthening constitutive model including the initial yield surface size, back stress tensor and kinematic hardening modulus, and using them as the cyclic hardening curve; inputting the temperature-related yield strength data and the cyclic hardening curve into the mesh element integration points corresponding to the compensator structure data in the reference finite element twin model to generate the physical field coupling model.

[0011] Preferably, the specific generation process of the thermo-coupled simulation model includes: analyzing the operating temperature time series, dividing the pipe into an upper vapor phase region and a lower liquid phase region based on the condensate accumulation height, calculating the gas film heat transfer coefficient and liquid film heat transfer coefficient of the corresponding regions respectively, and mapping them to the inner surface nodes of the physical field coupling model to generate transient temperature field loads; converting the internal pressure fluctuation sequence into a set of distributed force vectors perpendicular to the inner wall unit surface of the physical field coupling model to generate the inner wall pressure load; and based on the unsteady time step criterion, synchronously assigning the transient temperature field load and the inner wall pressure load to the boundary condition set of the physical field coupling model to generate the thermo-coupled simulation model with time-varying thermodynamic boundaries.

[0012] Preferably, the specific process of the driving thermo-coupling simulation model performing nonlinear iterative calculation includes: performing nonlinear incremental iterative calculation based on the unsteady simulation load step size to calculate the transient trial deformation of the grid nodes under the current load increment step; extracting the radial displacement of the pipeline in the transient trial deformation of the grid nodes and updating the soil spring stiffness matrix according to the radial displacement of the pipeline; detecting the viscous state of the nonlinear friction contact pair and adjusting the tangential friction resistance of the nonlinear friction contact pair; iterative calculation until the transient trial deformation of the grid nodes meets the convergence equilibrium condition; in the convergence state, extracting the equivalent stress values ​​and equivalent plastic strain values ​​of the peak nodes and trough nodes corresponding to the compensator structure data and located in the lower liquid phase region, and generating equivalent stress response time history and plastic strain range data.

[0013] Preferably, the specific process for generating the remaining service life prediction value includes: calling the rainflow counting algorithm to analyze the equivalent stress response time history and identify the number of full-cycle stress cycles; combining the plastic strain range data and the number of full-cycle stress cycles to correct the low-cycle fatigue life equation and construct the corrected fatigue damage assessment model; driving the corrected fatigue damage assessment model to calculate the single-cycle damage increment, superimposing all the single-cycle damage increments to generate the current cumulative damage degree; retrieving the annual hardening trend of soil stiffness and the nonlinear growth trend of corrosion rate from the digital twin evolution database, and calculating the current damage rate in combination with the current cumulative damage degree; calculating the numerical difference between the preset failure threshold and the current cumulative damage degree as the remaining safety capacity; performing nonlinear trend evolution prediction on the remaining safety capacity based on the current damage rate to generate the remaining service life prediction value.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. By calling the pipe wall remaining thickness distribution matrix and pipe body out-of-roundness geometric deviation data from the pipe geometry data, a non-ideal pipe geometry topology is generated, so that the constructed physical field coupling model includes the geometric defect features obtained from actual scanning. When predicting the remaining life, calculations can be performed based on the local stress changes caused by pipe wall thinning and out-of-roundness, so that the final output evaluation data corresponds to the current actual geometric state of the pipe, rather than the ideal design state.

[0015] 2. By analyzing the operating temperature time series, the pipe is divided into an upper vapor phase region and a lower liquid phase region based on the condensate accumulation height. Finally, a transient temperature field load is generated, which transforms the condensate accumulation condition in the actual operating data into the model boundary conditions. This allows the thermo-mechanical coupling simulation model to take into account the influence of the difference in heat transfer coefficients between the gas and liquid phases on the thermal stress distribution of the pipe wall during calculation, thereby obtaining stress response time history data consistent with the actual operating temperature fluctuations and medium distribution.

[0016] 3. By extracting the radial displacement of the pipeline in the transient trial deformation of the grid nodes and updating the soil spring stiffness matrix based on the radial displacement of the pipeline, a coupling mechanism between deformation and constraint stiffness is introduced in the nonlinear iterative calculation process. Combined with the acquisition of the operating temperature time series and internal pressure fluctuation sequence, it is beneficial that the final output of the remaining service life prediction value is based on the actual load fluctuation and is calculated under the convergence state of dynamic updating of soil spring stiffness with pipeline deformation.

[0017] 4. By integrating the non-ideal geometric topology of the pipeline, transient temperature field loads, and dynamically updated soil spring stiffness matrix into a unified thermo-mechanical coupling simulation model and performing nonlinear iterative calculations, a strongly coupled analysis system of structure-thermal-soil multiphysics fields was constructed. This reduced the problem of the separation between geometric defects, thermal boundaries, and soil constraints, and facilitated the simultaneous reflection of the interaction between pipeline geometric distortion, gas-liquid temperature difference effects, and soil constraint reaction forces when calculating the equivalent stress response time history. This improved the overall credibility of the remaining life prediction results under complex coupling environments. Attached Figure Description

[0018] Figure 1 The flowchart shows the method for predicting the remaining life of a directly buried compensator based on digital twins, as proposed in an embodiment of this invention. Figure 2 This is a flowchart illustrating the construction process of the multi-source heterogeneous data fusion and physical field coupling model proposed in an embodiment of this invention. Figure 3 This is a flowchart of the thermodynamic coupling dynamic iteration and full-lifetime evolution prediction proposed in an embodiment of this invention application. Detailed Implementation

[0019] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] Please see Figures 1-3 The method for predicting the remaining life of directly buried compensators based on digital twins provided by this invention includes the following specific steps: Collect pipeline geometry data, compensator structure data, and soil physical parameters of the directly buried heating pipeline network; construct a three-dimensional mesh model based on the pipeline geometry data and compensator structure data, and calculate the soil spring stiffness matrix based on the soil physical parameters; map the soil spring stiffness matrix to the three-dimensional mesh model to construct a reference finite element twin model; Obtain cyclic hardening curves and temperature-related yield strength data, and map them to the baseline finite element twin model to generate a physical field coupling model; The operating temperature time series and internal pressure fluctuation series are obtained and converted into transient temperature field loads and internal wall pressure loads, respectively. The transient temperature field loads and internal wall pressure loads are applied to the physical field coupling model to form a thermo-mechanical coupling simulation model. The thermo-coupling simulation model is driven to perform nonlinear iterative calculations and outputs the equivalent stress response time history and plastic strain range data of the peak and trough nodes; Input the equivalent stress response time history and plastic strain range data into the modified fatigue damage assessment model to calculate the current cumulative damage degree; based on the current cumulative damage degree and the preset failure threshold, calculate and output the predicted value of the remaining service life.

[0021] The technical solution of the present invention will be further described in detail below with reference to specific embodiments.

[0022] Example 1

[0023] This application discloses a method for predicting the remaining life of directly buried compensators based on digital twins. (See attached document.) Figure 1 The specific steps proposed in this invention include: S1. Collecting pipeline geometric data, compensator structural data, and soil physical parameters of the buried heating pipeline network; constructing a three-dimensional mesh model based on the pipeline geometric data and compensator structural data, and calculating the soil spring stiffness matrix based on the soil physical parameters; mapping the soil spring stiffness matrix to the three-dimensional mesh model to construct a reference finite element twin model; S2. Obtaining cyclic hardening curves and temperature-related yield strength data, and mapping them to the reference finite element twin model to generate a physical field coupling model; S3. Obtaining the operating temperature time series and internal pressure fluctuation sequence, and converting them into transient temperature field loads and inner wall pressure loads, respectively; applying the transient temperature field loads and inner wall pressure loads to the physical field coupling model to form a thermo-mechanical coupling simulation model; S4. Driving the thermo-mechanical coupling simulation model to perform nonlinear iterative calculations, and outputting the equivalent stress response time history and plastic strain range data of the peak and trough nodes; S5. Inputting the equivalent stress response time history and plastic strain range data into the modified fatigue damage assessment model to calculate the current cumulative damage degree; calculating and outputting the predicted value of the remaining service life based on the current cumulative damage degree and the preset failure threshold.

[0024] Further, collect pipeline geometry data, compensator structural data, and soil physical parameters for the directly buried heating pipeline section; corresponding to step S1 above; see [link / reference]. Figure 2 The specific implementation process includes: The pipeline geometry data includes the remaining wall thickness distribution matrix and out-of-roundness geometric deviation data of the pipe body obtained by scanning; the compensator structure data includes the initial waveform vector of the bellows extracted from the installation completion digital archive and the interlayer corrosion depth data based on eddy current nondestructive testing; the soil physical parameters include the lateral earth pressure response data collected by the earth pressure cell array pre-embedded in the backfill soil of the pipe trench, the soil moisture content collected by the time domain reflectometer, and the soil dynamic shear modulus and soil dynamic unit weight parameters based on the soil moisture content mapping.

[0025] Specifically, the process of obtaining pipeline geometry data is as follows: In actual engineering implementation, the data acquisition process employs a multi-channel high-frequency ultrasonic scanning device. The detection cross-section is set at 50 mm increments along the axial direction of the directly buried heating pipeline, with 36 sampling points evenly distributed circumferentially in each cross-section to ensure full coverage. The acquired raw ultrasonic pulse echo signals are processed through data fusion to output a multi-dimensional matrix of the remaining pipe wall thickness distribution. For example, for a directly buried steel pipe with a nominal wall thickness of 8 mm, this matrix can record the precise location information of wall thickness reduction to 7.2 mm caused by localized electrochemical corrosion. The calculation process for the pipe body out-of-roundness geometric deviation data involves using a laser displacement sensor to obtain the maximum and minimum diameters of the actual cross-section, and quantifying the difference between the two diameters as a percentage of the nominal diameter, thus accurately reflecting the cross-sectional distortion state of the pipe body under soil load.

[0026] Specifically, the process of obtaining the compensator structure data is as follows: The initial waveform vector of the bellows is obtained by retrieving digital measurement files from the compensator manufacturing stage. Its data structure includes a two-dimensional coordinate sequence of corrugation height, pitch, crest radius, and trough radius. For multi-layered bellows, pulsed eddy current nondestructive testing (PDT) is used to acquire the induced current response signal generated by the probe on the bellows surface. The interlayer corrosion depth data is then retrieved based on the nonlinear attenuation law of the signal amplitude. For example, when the detection system detects a significant phase angle shift in the eddy current signal, the built-in inversion algorithm calculates the equivalent interlayer corrosion depth to be 0.12 mm; this data is defined as the geometric damage amount.

[0027] Specifically, the process for obtaining soil physical parameters is as follows: During the trench backfilling stage, lateral earth pressure response data, in kilopascals, is collected and output in real time using a pre-embedded earth pressure cell array. Simultaneously, soil moisture content is output by measuring the soil dielectric constant using a time-domain reflectometer. The calculation process for the dynamic soil bulk density parameter is the product of the dry soil unit weight and the soil moisture content with the unit weight of water. The mapping calculation process for the dynamic soil shear modulus is the product of the reference shear modulus and the moisture content exponent. Specifically, the attenuation function is calculated by using the difference between the current volumetric moisture content and the optimum moisture content (e.g., 18%) as the exponent, calculating the negative exponent of the natural constant e, and introducing a structural damage factor (valued at 0.5) to reflect the softening effect of the skeleton after water immersion. For example, when the soil moisture content is monitored to increase from 15% to 22% due to external environmental influences, the mapped dynamic bulk density increases accordingly, while the dynamic shear modulus decreases from 15 MPa to 12 MPa due to soil softening.

[0028] By scanning and obtaining the remaining pipe wall thickness distribution matrix and pipe out-of-roundness geometric deviation data, the data includes local geometric features caused by manufacturing processes or service corrosion, thus enabling the location of specific stress concentration areas in the calculation. Simultaneously, by dynamically mapping the soil dynamic shear modulus based on soil moisture content, the model parameters can be dynamically adjusted with changes in the rainy season, dry season, or groundwater level, reflecting the actual fluctuations in the soil's constraint capacity on the pipeline at different times, which is beneficial to the realism and timeliness of the simulation calculation boundary conditions.

[0029] Furthermore, a three-dimensional mesh model is constructed based on the pipeline geometry data and compensator structural data, and the soil spring stiffness matrix is ​​calculated based on soil physical parameters; this corresponds to step S1 above; see [link / reference]. Figure 2 The specific implementation process includes: The remaining wall thickness distribution matrix and pipe body out-of-roundness geometric deviation data from the pipeline geometry data are called. The pipe body out-of-roundness geometric deviation data is superimposed on the standard circular pipe model as a radial displacement constraint to generate a non-ideal pipeline geometry topology. The initial waveform vector and interlayer corrosion depth data of the bellows in the compensator structure data are used to generate a multi-layer solid structure of the bellows. A three-dimensional solid model is generated based on the non-ideal pipeline geometry topology and the multi-layer solid structure of the bellows, and a hexahedral element discretization operation is performed on the three-dimensional solid model to generate the three-dimensional mesh model. The soil dynamic shear modulus and lateral earth pressure response data from the soil physical parameters are extracted. A depth correction coefficient is introduced based on the burial depth difference, and the axial, lateral and vertical stiffness components are calculated in combination with the soil dynamic shear modulus. The axial, lateral and vertical stiffness components are assembled to generate the soil spring stiffness matrix.

[0030] In this embodiment, the original matrix data structure and dimensions are kept unchanged. A Kelvin-Voyt viscoelastic model is introduced as a post-processing module to calculate the soil rheological relaxation coefficient and perform time-varying attenuation correction on each component of the original stiffness matrix. Specifically, after calculating the initial axial, lateral, and vertical stiffness components based on the soil dynamic shear modulus, the clay content index in the soil physical parameters is further read. If the clay content exceeds 15%, the viscoelastic correction procedure is initiated based on the original stiffness value. First, a parallel spring-damper mechanical unit, i.e., the Kelvin-Voyt model, is constructed, and the soil viscosity coefficient is set, for example, 250,000 Pascals-seconds for silty clay. Next, the complex stiffness modulus is calculated based on the loading frequency determined by the operating temperature time series. Subsequently, the Proni series is introduced to fit the stress relaxation function of the soil, and the relaxation time constant is set to 720 hours. This parameter should be determined according to soil consolidation experiments. 720 hours (approximately 30 days) corresponds to the typical characteristic period of pore water pressure dissipation and skeleton reconstruction of silty clay after groundwater level changes. When generating the final soil spring stiffness matrix, the diagonal structure of the original matrix is ​​not changed. Instead, the original static stiffness values ​​are multiplied by a time-dependent relaxation factor. The calculation logic for this relaxation factor is described as follows: 1 minus the product of the normalized shear modulus coefficient and the time decay term, where the time decay term is equal to the ratio of the current load duration to the relaxation time constant. The negative of this ratio is used as the exponent to calculate the natural constant e raised to this exponent. The normalized shear modulus coefficient is set to 0.15. This multiplication correction is performed on each stiffness component of the matrix, outputting a soil spring stiffness matrix with the same format as the original matrix but containing rheological properties in the values. This process, by introducing viscoelastic rheological correction, allows the soil spring stiffness to dynamically reflect the time-dependent changes in the soil's constraint force on the pipeline, improving the accuracy of remaining life prediction.

[0031] Specifically, the construction process of the 3D mesh model is as follows: The system uses a pre-stored 800x36 dimension pipe wall remaining thickness distribution matrix as initial input, where matrix elements represent the measured thickness values ​​of corresponding spatial nodes. First, a standard circular pipe geometric model with an outer diameter of 508 mm, a wall thickness of 8 mm, and a length of 10 m is established in the modeling space. Next, an interpolation algorithm maps the pipe's out-of-roundness geometric deviation data as radial displacement constraints to the circumferential nodes of the standard circular pipe. By superimposing displacement corrections on the radial coordinates of these nodes, the standard circular cross-section is transformed into a non-ideal ellipse or distorted cross-section. For example, if the out-of-roundness deviation at the right side of a cross-section is 1.5 mm, the radial coordinates of that node are shifted outward by the corresponding value. Subsequently, the values ​​in the wall thickness distribution matrix are assigned to the shell element thickness attributes at the corresponding coordinate positions, thus outputting a non-ideal pipe geometric topology that includes real geometric defects and uneven wall thickness characteristics.

[0032] Using the initial waveform vector of the bellows in the digital archive as the trajectory, a multi-layered concentric solid model is constructed in three-dimensional space through path sweeping technology. Taking a four-layer bellows structure as an example, the thickness of each layer is set to 0.8 mm, and the initial gap between layers is 0.05 mm. During the generation process, interlayer corrosion depth data based on eddy current inversion is read, and local Boolean subtraction operations are performed on the solid geometry of specific layers. For example, if the inversion data shows that the second layer of bellows has a pitting corrosion with a depth of 0.12 mm at the crest, then the corresponding volume of solid elements is removed within that local coordinate region. The final output multi-layered solid structure of the bellows has a complex sinusoidal or U-shaped corrugated shape and also integrates the geometric weakening characteristics caused by material damage.

[0033] The generated pipe topology and corrugated pipe entity are Boolean-assembled to form a unified three-dimensional solid model. The discretization process uses a mapped mesh generation technique. For the corrugated pipe region, high-order eight-node hexahedral linear reduced integral elements are used for refinement. At least 6 layers of elements are set in the corrugation height direction, and 24 elements are distributed per wave in the axial wave width direction to capture complex bending stress gradients and output a three-dimensional mesh model.

[0034] Specifically, the generation process of the soil spring stiffness matrix is ​​as follows: Using soil physical parameters as input, the core calculation logic is based on buried pipeline mechanics formulas. The process of calculating the maximum lateral soil resistance per unit length is described as follows: the maximum lateral soil resistance equals the product of the horizontal bearing capacity factor, the dynamic unit weight of the soil, the burial depth of the pipeline center, and the pipe diameter. The horizontal bearing capacity factor is defined as the natural constant raised to the power of the sum of the tangent of the soil's internal friction angle and pi. For silty clay with an internal friction angle of 20 degrees, this factor is set between 3.5 and 4.5. For the actual working condition where the burial depth fluctuates along the axial direction, a depth correction coefficient is introduced. The depth correction coefficient is defined as a power function of the ratio of the actual burial depth to the reference burial depth (e.g., 1.5 meters), with the exponent set as the soil pressure increment factor (0.5 to 0.8 for silty clay; 0.6 in this embodiment). Taking a sandy soil environment with a burial depth of 2.5 meters as an example, if the dynamic shear modulus of the soil is 20 MPa, the calculated maximum axial friction force is 35,000 N / m, and the lateral yield displacement is set to 4% of the pipe diameter. The calculated nonlinear stiffness values ​​in the four dimensions of axial, lateral, vertical upward pull and vertical downward compression are assembled into a diagonal matrix. Specifically, three mutually orthogonal linear spring elements are established between each foundation node and pipe wall node, and axial, lateral and vertical stiffness values ​​are assigned respectively. The off-diagonal elements of the stiffness matrix are set to zero to form a soil spring stiffness matrix, which is then used as a boundary condition to be associated with the outer surface nodes of the three-dimensional mesh model.

[0035] By superimposing the pipe's out-of-roundness geometric deviation data as radial displacement constraints onto the standard model, a non-ideal geometric topology of the pipeline with actual physical form is generated, allowing the finite element mesh to directly inherit the geometric defects of the solid. Furthermore, based on the burial depth difference, a depth correction coefficient is introduced, and axial, lateral, and vertical stiffness components are assembled to construct an anisotropic soil spring stiffness matrix. This enables high-fidelity simulation of the mechanical properties of the pipe-soil interaction under direct burial conditions, facilitating the differentiation and expression of the differences in soil constraint stiffness experienced by the pipeline in different directions.

[0036] Furthermore, the soil spring stiffness matrix is ​​mapped to a three-dimensional mesh model to construct a baseline finite element twin model; corresponding to step S1 above; see [link to relevant documentation]. Figure 2 The specific implementation process includes: Identify the node coordinate sequence on the outer wall of the 3D mesh model and define the node coordinate sequence as the interactive interface; obtain the spring element topology corresponding to the soil spring stiffness matrix and define the spring element topology as the constraint boundary; construct a nonlinear frictional contact pair connecting the interactive interface and the constraint boundary based on Coulomb friction theory; associate the soil spring stiffness matrix with the nonlinear frictional contact pair to construct a reference finite element twin model.

[0037] Specifically, the construction process of the benchmark finite element twin model is as follows:

[0038] Using the pipe's central axis as a reference, all nodes with a radial distance equal to the pipe's outer diameter radius are retrieved. For non-ideal geometric topologies with out-of-roundness deviations, a 0.5 mm extreme diameter tolerance range is set to capture the node coordinate sequence located on the outer surface of the pipe's insulation sheath. The output is an interactive node matrix. For example, for a typical pipe segment of 10 meters in length, if the axial step is set to 50 mm and 36 measuring points are evenly distributed around the circumference, approximately 7200 feature coordinate points will be identified and extracted.

[0039] Using the soil spring stiffness matrix as input data, and based on the spatial distribution of the interactive interface nodes, corresponding virtual foundation nodes are automatically generated outside the normal vector of each interface node (e.g., at an offset of 100 mm). These foundation nodes are assigned fully constrained properties with six degrees of freedom, i.e., displacement and rotation are both set to 0, thus constructing the spring element topology. The input soil spring stiffness matrix is ​​discretized and mapped to each pair of interactive nodes and foundation nodes, forming a one-to-one correspondence of mechanical support points, and outputting a set of constraint boundaries.

[0040] Then, a mechanical response criterion connecting the interactive interface and the constraint boundary is established. Based on Coulomb's friction theory, the key parameter, the tangential friction coefficient, is set to 0.4 (for the interface between silty clay and epoxy coating). Its mechanical logic is described as follows: when the tangential shear stress between the contact surfaces is less than the product of the friction coefficient and the normal effective pressure, the interactive interface is in a viscous state, and the deformation is controlled by the initial tangential stiffness of the soil spring; when the shear stress reaches this critical threshold, the nonlinear friction contact pair enters a slip state. In the nonlinear solution and incremental iteration process, the residual force norm is introduced as a convergence criterion for determining the equilibrium state. The criterion is set as: the ratio of the infinite norm of the unbalanced force to the current load increment is less than 0.1%. Then, the assembled axial, lateral, and vertical stiffness components are dynamically associated with the corresponding contact pair elements. Subsequently, self-weight and initial ground stress loads are applied to drive the solver to perform nonlinear incremental iteration, eliminating spurious unbalanced forces in the model assembly process. When the residual force satisfies the convergence criterion, the current deformation and contact state are locked and taken as the initial state at time zero, i.e., the baseline finite element twin model is constructed.

[0041] By identifying the node coordinate sequence on the outer wall of the 3D mesh model to define the interaction interface, and constructing nonlinear frictional contact pairs connecting the interaction interface and the constraint boundary based on Coulomb friction theory, a clear pipe-soil boundary interaction mechanism was established. Through this contact pair setup, the effective absorbed displacement of the corrugated compensator under actual working conditions was predicted, and the local buckling behavior of the pipe segment that may be induced by uneven friction force distribution was captured, ensuring that the boundary constraint conditions conform to the actual physical and mechanical processes.

[0042] Further, the cyclic hardening curve and temperature-dependent yield strength data are acquired and mapped to the reference finite element twin model to generate a physical field coupling model; corresponding to step S2 above; the specific implementation process includes: Multi-temperature gradient uniaxial tensile test data of the compensator material are used, and material aging reduction coefficients are introduced in combination with interlayer corrosion depth data to correct the nonlinear decay curve of yield strength with temperature. The correction result is used as the temperature-dependent yield strength data. The stable hysteresis loop data of strain-controlled cyclic loading test are analyzed, and parameters of the hybrid strengthening constitutive model including the initial yield surface size, back stress tensor, and kinematic hardening modulus are extracted and used as the cyclic hardening curve. The temperature-dependent yield strength data and the cyclic hardening curve are input into the mesh element integration points of the reference finite element twin model corresponding to the compensator structure data to generate a physical field coupling model.

[0043] Specifically, the construction process of the physical field coupling model is as follows: The dataset of uniaxial tensile tests on pipe steel within a multi-temperature gradient range of 20 to 650 degrees Celsius is used. The original curves show that the material's yield strength decreases non-linearly and exponentially with increasing temperature, for example, from 319 MPa at 20 degrees Celsius to 181 MPa at 300 degrees Celsius. Subsequently, interlayer corrosion depth data is used to quantify material damage by calculating a material aging reduction factor. This factor is calculated as follows: the material aging reduction factor equals 1 minus the ratio of the measured corrosion depth to the nominal wall thickness. If the nominal thickness of the bellows is 0.8 mm and the measured interlayer corrosion is 0.12 mm, then the material aging reduction factor is 0.85. Based on this factor, a double correction is performed: the original yield strength is multiplied by this factor to obtain the corrected yield strength.

[0044] Using strain-controlled cyclic loading experimental data at a constant temperature of 300 degrees Celsius as input, hysteresis loop features were extracted by identifying the stable hysteresis loop data of the stress-strain response in the 200th stable cycle. The experiment employed a nonlinear least squares fitting algorithm, the core of which was to optimize three sets of Chabosh kinematic hardening components to make the model's predictions approximate the experimental trajectory. During training, the sum of squared residuals of stress deviations was introduced as a loss function, and a convergence threshold of 10 was set. -4 The specific parameter settings employ a segmented coverage strategy: the first component describes the initial high hardening rate stage, with an initial yield surface size set to 220 MPa, a kinematic hardening modulus of 100350 MPa, and a coefficient of restitution of 2750; the second component describes the elastoplastic transition stage, with a kinematic hardening modulus of 15000 MPa and a coefficient of restitution of 350; the third component describes the large strain saturation stage, with a kinematic hardening modulus of 1500 MPa and a coefficient of restitution of 8. The output is a hybrid strengthening constitutive model parameter that includes the initial yield surface, the evolution law of the back stress tensor, and the kinematic hardening parameters. The kinematic hardening modulus is multiplied by the material aging reduction factor to generate a cyclic hardening curve reflecting the decline in plastic strengthening capacity after section damage.

[0045] Input the corrected yield strength data and cyclic hardening curve, and traverse the refined mesh region corresponding to the compensator corrugated structure in the baseline finite element twin model. Through the defined element attribute association logic, the temperature-dependent elastic modulus, Poisson's ratio, and corrected hardening parameters are assigned to the mesh element integration points of each solid element layer. Taking the four-layer corrugated pipe model as an example, the integration points of the second layer elements located in the most corroded area are assigned a lower strength limit, and the physical field coupling model is output.

[0046] By incorporating interlayer corrosion depth data and introducing a material aging reduction factor to correct the yield strength curve, the performance evolution of the material under long-term service conditions is integrated into the physical field coupling model, enabling the model parameters to objectively reflect the strength decay state of the metallic material. Simultaneously, by analyzing strain-controlled cyclic loading experimental data to extract the back stress tensor and kinematic hardening modulus, the output plastic strain range data is made more consistent with the true elastoplastic response of metallic materials under complex cyclic loading, thus improving the material mechanics accuracy of the fatigue assessment input data.

[0047] Further, the operating temperature time series and internal pressure fluctuation sequence are obtained and converted into transient temperature field loads and internal wall pressure loads, respectively; the transient temperature field loads and internal wall pressure loads are applied to the physical field coupling model to form a thermo-mechanical coupling simulation model; corresponding to step S3 above; see reference Figure 3 The specific implementation process includes: The operating temperature time series is analyzed, and the pipe is divided into an upper vapor phase region and a lower liquid phase region based on the condensate accumulation height. The gas film heat transfer coefficient and liquid film heat transfer coefficient of the corresponding regions are calculated and mapped to the inner surface nodes of the physical field coupling model to generate transient temperature field loads. The internal pressure fluctuation sequence is converted into a set of distributed force vectors perpendicular to the inner wall unit surface of the physical field coupling model to generate the inner wall pressure load. Based on the unsteady time step criterion, the transient temperature field load and the inner wall pressure load are synchronously assigned to the boundary condition set of the physical field coupling model to generate the thermo-mechanical coupling simulation model with time-varying thermodynamic boundaries.

[0048] Specifically, the construction process of the thermo-coupling simulation model is as follows: Using the operating temperature time-series data and condensate level monitoring data output by the acquisition system as input, the spatial coordinates of each grid node within the pipe cross-section are identified. Based on the condensate accumulation height (e.g., set to 150 mm), the inner wall of the pipe is divided into two physical boundary layers: the area above the liquid level is defined as the upper vapor phase region, and the area below the liquid level is defined as the lower liquid phase region. The calculation process is described as follows: the gas film heat transfer coefficient and liquid film heat transfer coefficient are calculated using the Nusselt formula based on the Dittosh-Belt correlation. Under typical operating conditions, if the steam temperature is 260 degrees Celsius, the gas film heat transfer coefficient is set between 6000 and 15000 Kelvin per square meter; for the lower liquid phase region, the liquid film heat transfer coefficient is set between 1000 and 6000 Kelvin per square meter. Due to the difference in heat transfer efficiency between the liquid and gas phases, these non-uniform convective heat transfer coefficients are mapped to the inner surface nodes of the physical field coupling model, generating transient temperature field loads that vary with spatial location.

[0049] In this embodiment, after mapping the overall basic heat transfer coefficient, the boundary between the upper vapor phase and the lower liquid phase is identified, a gas-liquid interface fluctuation thermal shock model is constructed, the local oscillation enhancement coefficient is calculated, and the existing heat transfer coefficients of the boundary region nodes are superimposed and corrected to generate a transient temperature field load containing local thermal shock characteristics. Specifically, after mapping the gas film heat transfer coefficient and liquid film heat transfer coefficient to the inner surface nodes of the physical field coupling model to form the basic heat transfer boundary, the rate of change of condensate level height is further tracked. The range of liquid level fluctuation, for example, within ±20 mm, is defined as the thermal shock sensitive zone. Within this region, the original heat transfer coefficient is not replaced, but a turbulent oscillation enhancement mechanism of the gas-liquid two-phase flow is introduced for numerical superposition. If the absolute value of the liquid level change rate exceeds five millimeters per second, thermal shock correction is triggered. The correction algorithm uses a gain function based on the Nusselt number to calculate an additional increment. The calculation of this increment is described as: the basic heat transfer coefficient multiplied by the turbulence intensity coefficient, and then multiplied by the natural logarithm of the liquid level fluctuation factor. The calculation process for the liquid level fluctuation factor is as follows: Obtain the absolute value of the current liquid level change rate, divide it by a preset reference change rate (set to 1 mm / s, using the noise limit value of the condensate level sensor under steady flow conditions as the reference), and obtain a dimensionless ratio. The turbulence intensity coefficient is set to 0.35, determined by pre-simulating the gas-liquid two-phase flow under the same pipe diameter and steam velocity using computational fluid dynamics software, extracting the average turbulent kinetic energy intensity at the gas-liquid interface. The value typically ranges from 0.2 to 0.5. For nodes located within the sensitive zone, the final heat transfer coefficient equals the original basic heat transfer coefficient plus the additional increment calculated above. For example, during a sudden and violent fluctuation in liquid level, the heat transfer coefficient at the interface node increases significantly from its original value, and an axially distributed gradient heat flux density vector is applied to this region. This process only updates local node attributes and does not change the data format of the transient temperature field load, which is beneficial for the subsequent thermo-mechanical coupling simulation model to load correctly. This process, by constructing a gas-liquid interface fluctuation thermal shock model, effectively captures the sudden change in local thermal stress caused by condensate sloshing inside the compensator and predicts the damage accumulation in high-incidence areas of corrosion fatigue, while retaining the original macroscopic partition calculation.

[0050] The input is a sequence of internal pressure fluctuations from the monitoring system. First, the normal vector of each element surface on the inner wall of the physics-field coupled model is extracted. The load transformation process is described as follows: for any inner wall element, the distributed force vector acting on it is equal to the current internal pressure value (e.g., fluctuating within the range of 0.8 to 1.2 MPa) multiplied by the surface area of ​​the element, and acts on the mesh nodes along the direction of the surface normal vector. The output is the dynamically changing internal wall pressure load. Taking a pipe with a nominal diameter of 500 mm as an example, if the internal pressure fluctuation is 0.1 MPa, the stress state of the inner wall nodes is automatically calculated and updated.

[0051] The input consists of transient temperature field loads and internal wall pressure loads. An automatic step-size algorithm is executed based on the severity of load fluctuations. The step-size setting criteria are as follows: Monitoring the inlet temperature change rate and internal pressure change rate, when the absolute value of the temperature change rate exceeds 2 degrees Celsius per minute, or the absolute value of the internal pressure change rate exceeds 0.05 MPa per minute, it is determined to be a period of severe fluctuation, and the step size is automatically increased to 2 to 5 minutes (e.g., set to 120 seconds) to capture transient thermal shock effects; when both of the above change rates are below 10% of the set threshold, it is determined to be a stable operating phase, and the step size is increased to 2 hours to improve computational efficiency. Through boundary condition synchronization logic, the time-varying temperature load and pressure load are synchronously applied to the model's boundary condition set, and a convergence criterion equivalent to a loss function is introduced, namely, the global plastic dissipation energy increment norm must be less than 1% of the work increment of the external load. The output is a thermo-mechanical coupling simulation model with time-varying thermo-mechanical boundaries.

[0052] By calculating the heat transfer coefficients of the upper vapor phase region and the lower liquid phase region respectively and mapping them to the inner surface nodes, a non-uniformly distributed transient temperature field load is constructed. Combined with the distributed force vector generated by the internal pressure fluctuation sequence, this thermo-hydraulic coupling simulation model can reproduce the stress state of the pipeline under complex thermal-hydraulic conditions, which is beneficial for load application and response analysis of the peak and trough nodes with the most severe stress concentration.

[0053] Furthermore, the thermo-coupling simulation model is driven to perform nonlinear iterative calculations, outputting the equivalent stress response time history and plastic strain range data for the peak and trough nodes; corresponding to step S4 above; the specific implementation process includes: Based on the unsteady simulation load step, nonlinear incremental iterative calculations are performed to calculate the transient trial deformation of the mesh nodes under the current load increment step; the radial displacement of the pipeline in the transient trial deformation of the mesh nodes is extracted, and the soil spring stiffness matrix is ​​updated according to the radial displacement of the pipeline; the viscous state of the nonlinear friction contact pair is detected, and the tangential friction resistance of the nonlinear friction contact pair is adjusted; iterative calculations are performed until the transient trial deformation of the mesh nodes meets the convergence equilibrium condition; in the convergence state, the equivalent stress values ​​and equivalent plastic strain values ​​of the peak nodes and trough nodes corresponding to the compensator structure data and located in the lower liquid phase region are extracted, and the equivalent stress response time history and plastic strain range data are generated.

[0054] Specifically, the detailed process of the nonlinear incremental iterative operation is as follows: The input time-varying load increment step is used, employing the Newton-Raphson iterative method. During the calculation, the unsteady simulation load step size is set to a hierarchical control logic: 3600 seconds for static operating conditions, and refined to 2 to 5 minutes (e.g., 120 seconds) for transient conditions with severe fluctuations. The radial displacement of the pipe at each mesh node is output by solving a set of nonlinear equations where the product of the tangent stiffness matrix and the displacement increment equals the residual force. A convergence control parameter, equivalent to a loss function, is introduced, with the convergence equilibrium condition setting the ratio of the residual force norm to the current load increment to be less than one-thousandth. For example, within an increment step where the internal pressure jumps from 0.8 MPa to 1.2 MPa, the radial displacement of the pipe at each node of the compensator corrugated structure is obtained through five sub-increment iterations.

[0055] The radial displacement of the pipeline at each node on the outer wall is extracted and substituted into the nonlinear earth pressure-displacement curve equation. The update process is described as follows: when the radial displacement of the pipeline is within the elastic limit (e.g., 3 mm), the initial stiffness value remains unchanged; if the displacement exceeds this threshold and enters the soil yielding stage, the tangent stiffness slope is automatically calculated, and the corresponding component in the soil spring stiffness matrix is ​​softened. The specific softening calculation formula is: the updated stiffness equals the initial stiffness multiplied by the softening factor. The softening factor is calculated as: the yield displacement threshold divided by the current actual radial displacement. This calculation is beneficial because as the displacement increases, the provided constraint force tends towards the asymptote of the ultimate soil resistance. For example, if a 5 mm radial bulge is detected at the outlet end of the compensator, the lateral stiffness of the soil spring will be automatically reduced from 20 MN / m to 15 MN / m.

[0056] Based on Coulomb's law of friction, the input is normal pressure distribution data. The system detects the tangential shear stress values ​​of each nonlinear friction contact pair on the contact surface and applies a criterion for determining viscous and slip states. The determination process is as follows: if the tangential shear stress is less than the product of the friction coefficient (set to 0.4) and the normal pressure, the node is determined to be in a viscous state, and a displacement constraint force is applied. If the shear stress reaches this critical value, the system switches to a slip state. In this state, the tangential frictional resistance is constantly equal to the product of the normal contact pressure and the dynamic friction coefficient, but its direction is updated to be opposite to the current relative slip velocity vector. This nonlinear state switching can capture the micro-slippage generated between the bellows and the soil during thermal expansion and contraction. In a practical case, when the pipe wall temperature rises to 260 degrees Celsius, if the axial slippage reaches 12 mm, the frictional resistance vector direction will be automatically adjusted, thus outputting highly accurate axial force time history data of the compensator.

[0057] Specifically, the process for generating the equivalent stress response time history and plastic strain range data is as follows: The input is the full-field response matrix after iterative convergence. The data is located in the defined lower liquid phase region, and the geometrically significant peak and trough nodes are identified based on the initial waveform vector of the bellows. The data extraction process is described as follows: the von Mises equivalent stress algorithm is called to calculate the combined values ​​of the six stress components at the specific nodes, and the evolution of the equivalent plastic strain over time is recorded. Taking a compensator that has been in operation for 5 years as an example, the extraction results show that the equivalent stress at the trough node reaches 420 MPa under the maximum thermal cycle step, and the equivalent plastic strain range generated by a single thermal load cycle is 0.45%. These data are formatted and output as the equivalent stress response time history and signed plastic strain range data, with the sign consistent with the sign of the first principal stress at the current moment.

[0058] During the nonlinear iteration process, the radial displacement of the grid nodes is detected, and the stiffness parameters of the soil spring and the tangential friction resistance of the nonlinear friction contact pair are adjusted accordingly. This achieves a strong two-way coupling between the structural deformation field and the environmental constraint field, which is beneficial for each iteration calculation to tend towards the dynamic mechanical equilibrium state of the structure and soil, and ensures the mechanical reliability of the fatigue assessment basic data.

[0059] Furthermore, the equivalent stress response time history and plastic strain range data are input into the modified fatigue damage assessment model to calculate the current cumulative damage level; based on the current cumulative damage level and the preset failure threshold, the predicted value of remaining service life is calculated and output; corresponding to step S5 above; the specific implementation process includes: The equivalent stress response time history is analyzed using a rainflow counting algorithm to identify the number of full-cycle stress cycles. Combining the plastic strain range data and the number of full-cycle stress cycles, the low-cycle fatigue life equation is corrected, and the corrected fatigue damage assessment model is constructed. The corrected fatigue damage assessment model is driven to calculate the single-cycle damage increment, and all single-cycle damage increments are superimposed to generate the current cumulative damage degree. The annual hardening trend of soil stiffness and the nonlinear growth trend of corrosion rate are retrieved from a digital twin evolution database, and the current damage rate is calculated based on the current cumulative damage degree. The difference between the preset failure threshold and the current cumulative damage degree is calculated as the remaining safety capacity. Based on the current damage rate, a nonlinear trend evolution prediction is performed on the remaining safety capacity to generate a predicted value for the remaining service life.

[0060] In this embodiment, the modified fatigue damage assessment model is constructed using a multi-level logical architecture, specifically including: (1) Data preprocessing layer: as the model input, it receives the equivalent stress response time history of the peak and trough nodes, calls the rainflow counting algorithm to identify the number of full-cycle stress cycles, and extracts the stress amplitude and average stress of each cycle; (2) Basic damage calculation layer: it receives the plastic strain range data and the cycle parameters output by the data preprocessing layer, and calculates the mechanical damage increment of a single cycle based on the low-cycle fatigue life equation with the introduction of the Moro correction term; (3) Environmental coupling correction layer: it calculates the environmental acceleration coefficient based on the Gutmann mechanical-chemical coupling theory, performs a doubling correction on the mechanical damage increment, and generates the total damage increment of a single cycle; (4) Cumulative damage assessment layer: as the model output, it superimposes the total damage increment of all single cycles based on the Mainner linear damage criterion and outputs the current cumulative damage degree.

[0061] Specifically, the process of constructing the modified fatigue damage assessment model is as follows: The input is equivalent stress response time history data located at the trough node in the lower liquid phase region, following the international standard ASTM E1049-85. The processing logic is as follows: First, hysteresis filtering is used to remove non-damaging fluctuations with stress amplitudes below 5% of the material's fatigue limit, discretizing the complex continuous stress time history into a sequence of extreme points. A dynamic double-ended queue buffer is set up, and new stress extreme points are pushed in. Only when the extreme points in the buffer satisfy the three-point counting criterion to form a closed hysteresis loop is the three-point loop counting logic of the rainflow counting algorithm executed. By detecting continuous stress peaks and valleys, closed stress-strain hysteresis loops representing material energy dissipation are identified and extracted. Taking a DN500 compensator at the outlet of a heating station as an example, in the stress history of one heating cycle (3600 hours), 12 large-amplitude full-cycle stress cycles caused by violent start-stop operations and 150 secondary cycles caused by internal pressure regulation were identified. The output is a fatigue load spectrum containing stress amplitude, average stress, and corresponding timestamps.

[0062] The inputs are fatigue load spectrum and equivalent plastic strain range data. The core of the modified fatigue damage assessment model is the low-cycle fatigue life equation based on the Manson-Coffin equation. The calculation formula is described as follows: the total strain amplitude equals the fatigue strength coefficient divided by the calculated elastic modulus, multiplied by the fatigue strength exponent raised to twice the number of failure cycles, plus the fatigue ductility coefficient multiplied by the fatigue ductility exponent raised to twice the number of failure cycles. For 304L stainless steel at 300 degrees Celsius, the specific parameter settings are as follows: fatigue strength coefficient is set to 660.17 MPa, fatigue ductility coefficient is set to 2.43, fatigue strength exponent is set to -0.11, and fatigue ductility exponent is set to -0.16. To account for the influence of mean stress evolution on the ratcheting effect, a Moro correction term is introduced to compensate for the elastic component. The specific process is as follows: the fatigue strength coefficient is subtracted from the mean stress value of the current cycle, and the difference is used as the corrected fatigue strength coefficient, replacing the original fatigue strength coefficient in the numerator of the elastic strain component in the low-cycle fatigue life equation. The output is the damage increment value corresponding to each identified cycle.

[0063] In this embodiment, after calculating the basic mechanical fatigue damage using the modified low-cycle fatigue life equation, a stress-corrosion coupling factor is introduced based on the Gutmann mechanochemical effect theory to calculate the environmental acceleration factor. This factor is then used to multiply and correct the basic damage value, constructing an evaluation model that incorporates the mechanochemical coupling effect. Specifically, after calculating the pure mechanical fatigue damage increment for a single cycle using the modified low-cycle fatigue equation, the basic calculation logic is not changed; instead, a fully coupled correction algorithm is introduced. First, the instantaneous equivalent stress at the trough node is extracted, and the mechanochemical activation energy increment is calculated. According to the Gutmann equation, the specific calculation logic is as follows: multiply the instantaneous equivalent stress at the trough node by the molar volume of the metal (7.1 cubic centimeters per mole for stainless steel), and then divide by the product of the gas constant and the absolute temperature to obtain the dimensionless mechanochemical activation energy increment exponent. Subsequently, the natural constant is calculated to the power of this exponent, and the result is the multiplication factor characterizing the stress-induced anodic dissolution rate. Based on this, an environmental acceleration factor greater than 1 is calculated. This factor is calculated as 1 plus the product of the coupling sensitivity coefficient and the stress-assisted factor. The stress-assisted factor is equal to the ratio of instantaneous equivalent stress to the material's yield strength. This ratio is used as an exponent to calculate the natural constant raised to that power. The coupling sensitivity coefficient (e.g., 0.05) is obtained by fitting the slope of corrosion current density as stress increases based on slow strain rate tensile tests under specific electrochemical environments (such as chloride ion concentration). The final modified fatigue damage assessment model defines the total damage per cycle as the original mechanical fatigue damage increment multiplied by the environmental acceleration factor. For example, this factor significantly increases the damage assessment value when the stress at the trough exceeds a threshold. This process quantifies the accelerated anodic dissolution phenomenon on metal surfaces under high stress conditions, improving the predictive ability for sudden equipment failures in corrosive environments.

[0064] Specifically, the process for generating the remaining service life prediction is as follows: Based on the Mainner linear damage criterion, all detected single-cycle damage increments are algebraically superimposed. If the superimposed value of the currently identified historical operational data is 0.42, this value is defined as the current cumulative damage degree. Subsequently, environmental degradation characteristics are dynamically retrieved from the digital twin evolutionary database. The input includes the annual hardening trend of soil stiffness with service time (set as a 2% annual increase in stiffness due to soil compaction and water loss). The actual soil spring stiffness value obtained by inversion at the current time step is compared with the hardening trend. If the deviation exceeds 5%, the least squares method is used to correct the exponential parameter of the hardening trend curve online, and the corrected evolution law is written back to the digital twin evolutionary database. It also includes the nonlinear growth trend of corrosion rate based on the exponential power law model. The power law exponent is set to 1.5, that is, the corrosion depth is proportional to the 1.5th power of the service time. This is based on the soil corrosivity level, referring to standards such as pipeline corrosion protection engineering inspection or historical plate test data. Usually, the power exponent of pitting corrosion is greater than 1. The calculation logic of the current damage rate is described as follows: the damage rate is equal to the product of the fatigue damage increment in the current time step and the environmental coupling acceleration factor.

[0065] The inputs are a preset failure threshold (set to 1.0, representing fatigue life exhaustion) and the current cumulative damage level. The calculation process is described as follows: the difference between the failure threshold and the current cumulative damage level is defined as the remaining safety capacity. Taking a current damage level of 0.42 as an example, the remaining safety capacity is 0.58. Then, a nonlinear extrapolation algorithm based on second-order Taylor series expansion is driven, combining the current damage rate and its second-order evolution acceleration due to environmental degradation, to perform time-axis projection. The backward difference method is used to calculate the second-order evolution acceleration, specifically: the damage rate at the current time step is subtracted by twice the damage rate at the previous time step, and then the damage rates from the previous two time steps are added. The result is then divided by the square of the time step. An error correction based on historical prediction bias is introduced during the prediction process. Specifically, the mean of the accuracy ratio between the predicted and measured values ​​over the past 5 time steps (i.e., measured damage increment divided by predicted damage increment) is calculated, and the reciprocal of this mean ratio is used as a correction coefficient multiplied by the remaining life prediction value obtained through nonlinear extrapolation. For example, if the prediction curve shows that damage accumulation is non-linearly accelerating due to the combined effects of corrosion and hardening, the critical moment when the remaining capacity decreases to zero can be obtained by integrating the Taylor formula. The final output is a predicted remaining service life, such as 5.2 years. This data is fed back to the operation and maintenance decision system in real time through the twin platform, providing maintenance cycle recommendations based on the actual damage status of the physical structure.

[0066] The digital twin evolutionary database adopts a time-series database architecture. During initialization, the baseline corrosion rate model parameters (such as power law exponent and initial coefficients) are entered based on the soil corrosivity survey report, and the hardening theoretical curve of soil stiffness changing over time is entered based on soil consolidation experimental data. During operation, the database stores the pipe-soil contact pressure, local corrosion depth inversion value, and cumulative damage degree calculated in each iteration. The sliding window algorithm is used to periodically correct the curve parameters of the hardening trend and corrosion growth trend. Specifically, the sliding window algorithm uses a fixed-length time window, corresponding to several recent operating condition sampling periods (e.g., the operating time of the most recent heating season). Within each time window, the measured soil stiffness value is fitted with the theoretical hardening curve using least squares, and the exponent parameters of the hardening trend curve are updated. Similarly, the corrosion depth data over time obtained based on eddy current detection is fitted with the preset power law corrosion model using least squares, and the power law coefficient and exponent parameters are corrected. The corrected parameters replace the original model parameters, thereby forming a dynamically updated evolutionary rule base.

[0067] By extracting the annual hardening trend of soil stiffness and the nonlinear growth trend of corrosion rate, and combining them with the current cumulative damage, the current damage rate is calculated. Based on this rate, a nonlinear trend evolution prediction is performed on the remaining safe capacity. The long-term evolution of environmental parameters and the nonlinear degradation of material properties are incorporated into the life prediction model, which helps operation and maintenance personnel identify the inflection point of accelerated life decay, formulate scientific maintenance and replacement plans, and extend the safe service life of the pipeline network.

[0068] This invention provides a digital twin-based method for predicting the remaining service life of buried compensators. By collecting pipeline geometric data, compensator structural data, and soil physical parameters, a baseline finite element twin model is constructed, incorporating actual physical characteristics and environmental constraint parameters. Based on this, a thermo-mechanical coupling simulation model is further generated by integrating cyclic hardening curves and transient temperature field loads. This model can directly calculate the stress concentration effects caused by the combined effects of pipe wall thinning, geometric non-roundness, soil constraint changes, and nonlinear material constitutive properties. The thermo-mechanical coupling simulation model is driven to perform nonlinear iterative calculations and output the equivalent stress response time histories of peak and trough nodes, transforming macroscopic operating conditions into microscopic material damage response data. Finally, by combining a modified fatigue damage assessment model with the current cumulative damage level, the remaining service life is calculated. This modeling and calculation method, driven by comprehensive data, enables dynamic quantitative assessment of the service status of buried pipeline networks, ensuring the reliability of pipeline operation and maintenance.

[0069] Example 2

[0070] This embodiment uses a digital twin-based method for predicting the remaining life of buried compensators applied to a buried steam pipeline compensator located 2000 meters from the outlet of a heating trunk line in a northern city. The pipe section has a nominal diameter of 508 mm, a nominal wall thickness of 8 mm, and is made of 304 stainless steel; it has been in operation for a cumulative period of five years.

[0071] During the data acquisition phase, an internal inspection robot carrying 36 ultrasonic probes performed a full scan at a speed of 0.2 meters per second. The inspection results, including the remaining pipe wall thickness distribution matrix, revealed localized erosion at the six o'clock position at the bottom of the pipe, with the thinnest value being 7.2 mm. A radial collapse of 1.5 mm was also measured at the inlet section of the compensator, indicating a geometric deviation in pipe out-of-roundness. Digital as-built documentation showed that the corrugated pipe has a four-layer structure, with a nominal thickness of 0.8 mm per layer and an initial wave pitch of 40 mm. Using a pulsed eddy current non-destructive testing system, inversion was performed without removing the insulation layer, determining the interlayer corrosion depth of the outer corrugated pipe at the third wave trough to be 0.12 mm. Simultaneously, the pre-embedded earth pressure cell array provided real-time feedback of lateral earth pressure at 105 kPa, and the time-domain reflectometry measured the volumetric moisture content of the backfill silty clay to be 22%.

[0072] The aforementioned moisture content maps to a dynamic soil bulk density of 19,500 N / m³ and a dynamic shear modulus of 12 MPa. A 1.5 mm radial displacement constraint is superimposed onto the baseline circular pipe to generate a non-ideal topological model incorporating real geometric defects. Subsequently, the 3D solid is discretized using hexahedral elements, with the mesh size for the corrugated pipe region set to 0.5 mm, ensuring six layers of elements are distributed along the thickness direction, resulting in a total of 7200 mesh nodes. The node coordinate sequence of the outer wall of the insulation jacket is identified and defined as the interaction interface. A six-degree-of-freedom fully constrained foundation node is established 100 mm outward in the normal direction, constructing a nonlinear contact pair with a tangential friction coefficient of 0.4, thereby building the baseline finite element twin model.

[0073] The yield strength benchmark of 304 stainless steel at 300°C is 181 MPa. Combined with a corrosion depth of 0.12 mm, the aging reduction factor is calculated to be 0.85, and the corrected temperature-dependent yield strength is determined to be 153.85 MPa. Parameters of the hybrid strengthening constitutive model are extracted, setting the initial yield surface size to 220 MPa, the first kinematic hardening modulus to 100350 MPa, and the kinematic hardening recovery coefficient to 2750. The 120-day operating sequence from the previous heating season is analyzed, with steam temperature fluctuating between 180°C and 260°C, and internal pressure cycling between 0.8 MPa and 1.2 MPa. Based on a condensate accumulation height of 150 mm, the film heat transfer coefficient of the upper vapor phase region is calculated to be 12000 W / m² Kelvin, and the liquid film heat transfer coefficient of the lower liquid phase region is calculated to be 3000 W / m² Kelvin.

[0074] During the iterative calculation phase, the steady-state step size was set to 3600 seconds, while the transient fluctuation step size was refined to 3 minutes. The driving model performed nonlinear incremental calculations, and the softened soil spring stiffness matrix was automatically adjusted when the radial displacement at the bottom of the pipe reached 5 mm. After the iteration converged to less than 0.1% of the unbalanced force residual, the peak value of the von Mises equivalent stress at the trough node in the liquid phase region was extracted to be 420 MPa, and the plastic strain range data for a single temperature difference cycle was output as 0.45%.

[0075] The rainflow counting algorithm identified 12 large-amplitude start-stop cycles and 150 secondary pressure fluctuation cycles in the annual operating history. The modified fatigue damage assessment model calculated the annual damage increment to be 0.085, and the current cumulative damage level, generated by combining five years of historical operating conditions, is 0.42. The power-law growth trend of the soil hardening coefficient and corrosion rate (2% annually) was retrieved from the evolution database, and the current damage rate was calculated to be 0.092 per year. Using 1 as the failure threshold, the remaining safe capacity was calculated to be 0.58. Based on nonlinear trend extrapolation and after excluding interference terms, the final prediction result shows that the remaining service life of the directly buried compensator is predicted to be 5.2 years.

[0076] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for predicting the remaining life of directly buried compensators based on digital twins, characterized in that, include: Collect pipeline geometry data, compensator structure data, and soil physical parameters for directly buried heating pipeline sections; A three-dimensional mesh model is constructed based on pipeline geometry data and compensator structural data, and the soil spring stiffness matrix is ​​calculated based on soil physical parameters; the soil spring stiffness matrix is ​​mapped to the three-dimensional mesh model to construct a benchmark finite element twin model; Obtain cyclic hardening curves and temperature-related yield strength data, and map them to the baseline finite element twin model to generate a physical field coupling model; The operating temperature time series and internal pressure fluctuation series are obtained and converted into transient temperature field loads and internal wall pressure loads, respectively. The transient temperature field loads and internal wall pressure loads are applied to the physical field coupling model to form a thermo-mechanical coupling simulation model. The thermo-coupling simulation model is driven to perform nonlinear iterative calculations and outputs the equivalent stress response time history and plastic strain range data of the peak and trough nodes; Input the equivalent stress response time history and plastic strain range data into the modified fatigue damage assessment model to calculate the current cumulative damage degree; Based on the current cumulative damage level and the preset failure threshold, calculate and output the predicted value of the remaining service life.

2. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The pipeline geometric data, compensator structural data, and soil physical parameters specifically include: the pipeline geometric data includes the remaining wall thickness distribution matrix and pipe body out-of-roundness geometric deviation data obtained through scanning; the compensator structural data includes the initial waveform vector of the corrugated pipe extracted from the installation completion digital archive and the interlayer corrosion depth data based on eddy current nondestructive testing; the soil physical parameters include the lateral earth pressure response data collected by the earth pressure cell array pre-embedded in the backfill soil of the trench, the soil moisture content collected by the time domain reflectometer, and the soil dynamic shear modulus and soil dynamic unit weight parameters based on the soil moisture content mapping.

3. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The specific construction process of the three-dimensional mesh model and the soil spring stiffness matrix includes: calling the pipe wall remaining thickness distribution matrix and pipe body out-of-roundness geometric deviation data from the pipe geometry data, superimposing the pipe body out-of-roundness geometric deviation data as radial displacement constraints onto the standard circular pipe model to generate a non-ideal pipe geometry topology; generating a corrugated pipe multi-layer solid structure using the corrugated pipe initial waveform vector and interlayer corrosion depth data from the compensator structure data; generating a three-dimensional solid model based on the non-ideal pipe geometry topology and the corrugated pipe multi-layer solid structure, and performing a hexahedral element discretization operation on the three-dimensional solid model to generate the three-dimensional mesh model; extracting soil dynamic shear modulus and lateral earth pressure response data from soil physical parameters, introducing a depth correction coefficient based on the burial depth difference, and calculating the axial, lateral, and vertical stiffness components in combination with the soil dynamic shear modulus, and assembling the axial, lateral, and vertical stiffness components to generate the soil spring stiffness matrix.

4. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The specific construction process of the reference finite element twin model includes: identifying the node coordinate sequence of the outer wall of the three-dimensional mesh model and defining the node coordinate sequence as the interactive interface; obtaining the spring element topology corresponding to the soil spring stiffness matrix and defining the spring element topology as the constraint boundary; constructing a nonlinear frictional contact pair connecting the interactive interface and the constraint boundary based on Coulomb friction theory; associating the soil spring stiffness matrix with the nonlinear frictional contact pair to construct the reference finite element twin model.

5. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The specific generation process of the physical field coupling model includes: calling the multi-temperature gradient uniaxial tensile test data of the compensator material, and introducing the material aging reduction coefficient in combination with the interlayer corrosion depth data to correct the nonlinear decay curve of yield strength with temperature, and using the correction result as the temperature-related yield strength data; analyzing the stable hysteresis loop data of the strain-controlled cyclic loading experiment, extracting the parameters of the hybrid strengthening constitutive model including the initial yield surface size, back stress tensor and kinematic hardening modulus, and using them as the cyclic hardening curve; inputting the temperature-related yield strength data and the cyclic hardening curve into the mesh element integration points of the reference finite element twin model corresponding to the compensator structure data to generate the physical field coupling model.

6. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The specific generation process of the thermo-coupled simulation model includes: analyzing the operating temperature time series, dividing the pipe into an upper vapor phase region and a lower liquid phase region based on the condensate accumulation height, calculating the gas film heat transfer coefficient and liquid film heat transfer coefficient of the corresponding regions respectively, and mapping them to the inner surface nodes of the physical field coupling model to generate transient temperature field loads; converting the internal pressure fluctuation sequence into a set of distributed force vectors perpendicular to the inner wall unit surface of the physical field coupling model to generate the inner wall pressure load; and based on the non-steady time step criterion, synchronously assigning the transient temperature field load and the inner wall pressure load to the boundary condition set of the physical field coupling model to generate the thermo-coupled simulation model with time-varying thermodynamic boundaries.

7. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The specific process of the nonlinear iterative calculation performed by the driving thermo-coupling simulation model includes: performing nonlinear incremental iterative calculation based on the unsteady simulation load step size to calculate the transient trial deformation of the grid nodes under the current load increment step; extracting the radial displacement of the pipeline in the transient trial deformation of the grid nodes and updating the soil spring stiffness matrix according to the radial displacement of the pipeline; detecting the viscous state of the nonlinear friction contact pair and adjusting the tangential friction resistance of the nonlinear friction contact pair; iterative calculation until the transient trial deformation of the grid nodes meets the convergence equilibrium condition; in the convergence state, extracting the equivalent stress and equivalent plastic strain values ​​of the peak nodes and trough nodes corresponding to the compensator structure data and located in the lower liquid phase region, and generating equivalent stress response time history and plastic strain range data.

8. The method for predicting the remaining life of a directly buried compensator based on digital twins according to claim 1, characterized in that, The specific process for generating the remaining service life prediction value includes: calling the rainflow counting algorithm to analyze the equivalent stress response time history and identify the number of full-cycle stress cycles; combining the plastic strain range data and the number of full-cycle stress cycles to correct the low-cycle fatigue life equation and construct the corrected fatigue damage assessment model; driving the corrected fatigue damage assessment model to calculate the single-cycle damage increment, superimposing all the single-cycle damage increments to generate the current cumulative damage degree; retrieving the annual hardening trend of soil stiffness and the nonlinear growth trend of corrosion rate from the digital twin evolution database, and calculating the current damage rate in combination with the current cumulative damage degree; calculating the difference between the preset failure threshold and the current cumulative damage degree as the remaining safety capacity; performing nonlinear trend evolution prediction on the remaining safety capacity based on the current damage rate to generate the remaining service life prediction value.

Images