Metro tunnel grouting process simulation system based on digital twinning

Abstract

This invention relates to the field of digital twin technology, specifically a simulation system for the grouting process of a subway tunnel based on digital twins. The system includes: a grouting environment construction module, used to acquire geometric data and soil physical property parameters of the subway tunnel, define continuous phase field variables to characterize the soil integrity state, input a multiphysics solution environment, apply grouting hole pressure boundary conditions, and obtain the initial phase field boundary state. In this invention, a comprehensive instability risk factor is constructed by combining pressure load and virtual strain energy evolution, realizing the transformation from simple pressure numerical monitoring to quantifying the inherent stability of the structure. By identifying critical instability states exceeding safety thresholds and triggering targeted grouting stop commands, sudden structural instability caused by insufficient local stiffness is avoided.

Description

Technical Field

[0001] This invention relates to the field of digital twin technology, and in particular to a simulation system for the grouting process of subway tunnels based on digital twins. Background Technology

[0002] Digital twin technology utilizes data such as physical models, sensor updates, and operational history to integrate a simulation process involving multiple disciplines, physical quantities, scales, and probabilities, completing the mapping in virtual space to reflect the entire lifecycle of the corresponding physical equipment.

[0003] Existing technologies for simulating grouting processes in underground engineering often require high-frequency mesh reconstruction when dealing with complex geological fracturing phenomena, resulting in low computational efficiency and difficulty in guaranteeing accuracy. In the structural safety assessment phase, current monitoring methods mainly rely on pressure sensor readings or direct observation of surface settlement, which are passive, post-hoc, or critical-point monitoring methods. When the stiffness matrix of the tunnel lining structure becomes singular due to long-term loading or geological disturbance, and its resistance to deformation decreases significantly, the pressure reading may not yet have reached the design alarm value. This monitoring blind spot can easily mislead construction personnel to continue applying pressure, creating safety hazards. Therefore, improvements are needed. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a simulation system for the grouting process of subway tunnels based on digital twins.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a digital twin-based simulation system for subway tunnel grouting process includes:

[0006] The grouting environment construction module is used to acquire the geometric data and physical property parameters of the subway tunnel, define continuous phase field variables to characterize the soil integrity state, input the multiphysics solution environment, apply the grouting hole pressure boundary conditions, and obtain the initial phase field boundary state.

[0007] The fracture propagation simulation module is used to call the initial phase field boundary state, execute an interleaved solution sequence within a fixed time step to update the displacement field and phase field variables, assign the formation permeability coefficient in the spatial region where the phase field variable is higher than the threshold to the maximum limit value, generate a permeability evolution distribution field, and calculate the fluid velocity and pressure distribution in the updated fracture topology based on the permeability evolution distribution field to obtain the grouting fracture propagation topology;

[0008] The segment stiffness analysis module is used to map the fluid pressure distribution data in the grouting splitting extension topology to the nodes of the tunnel segment ring structure, extract the nonlinear tangential stiffness matrix of the segment ring structure under the current loading condition, apply a virtual disturbance vector to the existing load field, generate a virtual displacement response field, perform singular value decomposition on the stiffness components associated with the virtual displacement response field, and calculate the stiffness matrix condition number index.

[0009] The grouting control decision module is used to call the stiffness matrix condition number index to quantify the ability of the segment ring to resist grouting pressure fluctuations, calculate and generate a comprehensive instability risk factor, compare the comprehensive instability risk factor with a preset safety threshold, identify the critical instability state segment ring where the comprehensive instability risk factor exceeds the preset safety threshold, trigger a stop grouting command signal for the corresponding grouting hole, and obtain a grouting control command signal.

[0010] Preferably, the step of obtaining the initial phase field boundary state is as follows:

[0011] Read the geometric data and physical property parameters of the subway tunnel, extract the cross-sectional dimensions, centerline coordinates, grouting hole coordinates, density, elastic modulus, Poisson's ratio, permeability coefficient and fracture toughness, unify the units and complete the coordinate registration, establish the node number and unit relationship, and form a joint grid of the geometric data and physical property parameters of the subway tunnel.

[0012] Based on the combined grid of the subway tunnel geometric data and the physical property parameters of the strata, phase field variables are set to indicate the soil integrity state, a summation structure of elastic potential energy, fracture surface energy and fluid pressure work is constructed, the grouting hole pressure boundary and fixed support boundary are registered, and a set of phase field variable initialization rules and boundary conditions is formed.

[0013] According to the initialization rules and boundary conditions set of the phase field variables, the initial values ​​of the phase field variables are assigned to each node in the solution domain and the grouting hole pressure boundary is applied to the boundary nodes. The state assembly is performed to generate the initial field of displacement degree of freedom and phase field variable degree of freedom, thus forming the initial phase field boundary state.

[0014] Preferably, the step of obtaining the permeability evolution distribution field is as follows:

[0015] The initial phase field boundary state is invoked, and the displacement degrees of freedom and phase field variables are updated sequentially within a fixed time step. The spatial distribution of displacement field and spatial distribution of phase field variables are recorded according to the time step, and the one-to-one correspondence between the two types of distributions is maintained by the node index to obtain the update results of displacement field and phase field variables.

[0016] Based on the update results of the displacement field and phase field variables, the phase field variables are compared with the threshold cell by cell. For cells whose phase field variables are greater than the threshold, the formation permeability coefficient is assigned the maximum limit value. For cells whose phase field variables are not greater than the threshold, the original permeability coefficient is maintained, thus generating a permeability evolution distribution field.

[0017] Preferably, the step of obtaining the grouting splitting propagation topology is as follows:

[0018] Based on the permeability evolution distribution field, candidate connected unit sequences of the fracture topology are extracted, the fluid velocity and pressure distribution of the candidate connected units are calculated, the energy functional values ​​are accumulated, and the crack propagation path is determined by the connected sequence corresponding to the minimum energy functional value. The crack propagation paths are merged into the fracture topology to form a grouting splitting propagation topology.

[0019] Preferably, the step of obtaining the virtual displacement response field is as follows:

[0020] The grouting splitting expansion topology is invoked, fluid pressure distribution data is read, the projection relationship between the node coordinates of the tunnel segment ring structure and the fracture surface element is established, the pressure values ​​are summarized according to the projection weight and converted to the node degrees of freedom to form the pressure load field of the tunnel segment ring structure node.

[0021] Based on the pressure load field of the tunnel segment ring structure nodes, the nonlinear tangent stiffness matrix under the current loading condition is extracted, a virtual disturbance vector with controlled amplitude is constructed and superimposed on the existing load field, the virtual displacement of each node is obtained by solving, and a virtual displacement response field is generated.

[0022] Preferably, the step of obtaining the stiffness matrix condition number index is as follows:

[0023] Based on the virtual displacement response field, the stiffness components associated with the virtual displacement response field are located and formed into an analysis submatrix. Singular value decomposition is performed to obtain the maximum and minimum singular values. The ratio of the maximum and minimum singular values ​​is calculated as an index value to obtain the stiffness matrix condition number index.

[0024] Preferably, the steps for obtaining the comprehensive instability risk factor are as follows:

[0025] The stiffness matrix condition number index is invoked to read the current pressure field data, aggregate the node pressure values ​​according to the segment ring number, and combine the stiffness matrix singular value information to establish a load-stiffness characteristic correspondence table for each segment ring, forming a comprehensive calculation input set;

[0026] Based on the comprehensive calculation input set, calculate the comprehensive instability risk factor for each segment ring.

[0027] Preferably, the step of acquiring the grouting control command signal is as follows:

[0028] Based on the comprehensive instability risk factor, threshold judgment is performed on each segment ring. Segment rings whose comprehensive instability risk factor exceeds the preset safety threshold are marked as critical instability states. The corresponding grouting hole numbers are extracted, and grouting control command signals are generated.

[0029] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0030] In this invention, by constructing continuous phase field variables to characterize the soil integrity state, and combining energy functionals that include elastic potential energy and fracture surface energy, the limitations of traditional mesh reconstruction on crack topological evolution can be overcome. This allows for the simulation of complex formation fracturing behavior caused by grouting under a fixed mesh, capturing crack propagation along arbitrary paths. An interleaved solution sequence is executed within a time step, and the formation permeability coefficient is dynamically updated based on the phase field variable threshold. Regions with high phase field values ​​are assigned limiting permeability, establishing a real-time mapping relationship between fracture aperture and fluid transport channels. This ensures that the fluid pressure distribution calculation accurately reflects the abrupt change in conductivity caused by grouting fracturing, solving the problem of seepage field lag in fluid-structure interaction analysis. To address the subsequent issues, fluid pressure data is mapped to the nodes of the tunnel segment ring structure, and the nonlinear tangential stiffness matrix is ​​extracted. Virtual disturbance vectors and virtual displacement response field analyses based on the principle of virtual work are introduced, which can detect the current stiffness degradation of the structure without damaging the actual structure. The ratio of the maximum to minimum singular values ​​of the stiffness matrix is ​​obtained using singular value decomposition technology, and a comprehensive instability risk factor is constructed by combining pressure load and virtual strain energy evolution. This realizes the transformation from simple pressure numerical monitoring to the quantification of the inherent stability of the structure. By identifying critical instability states that exceed the safety threshold and triggering targeted grouting stop commands, sudden structural instability caused by insufficient local stiffness is avoided. Attached Figure Description

[0031] Figure 1 This is a system flowchart of the present invention. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0033] Please see Figure 1 The present invention provides a technical solution: a simulation system for the grouting process of a subway tunnel based on digital twins, comprising:

[0034] The grouting environment construction module is used to acquire the geometric data and physical property parameters of the subway tunnel, define continuous phase field variables to characterize the soil integrity state, input the multiphysics solution environment, apply the grouting hole pressure boundary conditions, and obtain the initial phase field boundary state.

[0035] The splitting propagation simulation module is used to call the initial phase field boundary state, execute an interleaved solution sequence within a fixed time step to update the displacement field and phase field variables, assign the formation permeability coefficient in the spatial region where the phase field variable is higher than the threshold to the maximum limit value, generate the permeability evolution distribution field, and calculate the fluid velocity and pressure distribution in the updated fracture topology based on the permeability evolution distribution field to obtain the grouting splitting propagation topology.

[0036] The segment stiffness analysis module is used to map the fluid pressure distribution data in the grouting splitting propagation topology to the nodes of the tunnel segment ring structure, extract the nonlinear tangential stiffness matrix of the segment ring structure under the current loading conditions, apply a virtual disturbance vector to the existing load field, generate a virtual displacement response field, perform singular value decomposition on the stiffness components associated with the virtual displacement response field, and calculate the stiffness matrix condition number index.

[0037] The grouting control decision module is used to call the stiffness matrix condition number index to quantify the ability of the segment ring to resist grouting pressure fluctuations, calculate and generate a comprehensive instability risk factor, compare the comprehensive instability risk factor with a preset safety threshold, identify the critical instability state segment ring where the comprehensive instability risk factor exceeds the preset safety threshold, trigger a stop grouting command signal for the corresponding grouting hole, and obtain a grouting control command signal.

[0038] The steps for obtaining the initial phase field boundary state are as follows:

[0039] Read the geometric data and physical property parameters of the subway tunnel, extract the cross-sectional dimensions, centerline coordinates, grouting hole coordinates, density, elastic modulus, Poisson's ratio, permeability coefficient and fracture toughness, unify the units and complete the coordinate registration, establish the node number and unit relationship, and form a joint grid of the geometric data and physical property parameters of the subway tunnel.

[0040] Based on the joint grid of subway tunnel geometric data and stratum physical property parameters, phase field variables are set to indicate the soil integrity state, a summation structure of elastic potential energy, fracture surface energy and fluid pressure work is constructed, grouting hole pressure boundary and fixed support boundary are registered, and a set of phase field variable initialization rules and boundary conditions is formed.

[0041] Based on the initialization rules of phase field variables and the set of boundary conditions, initial values ​​of phase field variables are assigned to each node in the solution domain, and grouting hole pressure boundaries are applied to the boundary nodes. State assembly is performed to generate the initial field of displacement degrees of freedom and phase field variable degrees of freedom, thus forming the initial phase field boundary state.

[0042] Specifically, the process involves reading the geometric data and geological physical properties of the subway tunnel, and using computer-aided design software to analyze the tunnel design documents and geological survey reports. This allows for the extraction of the diameter of the circular profile, the three-dimensional coordinates of the center in the global coordinate system, the relative position angle of the grouting holes on the tunnel segments, and the density, elastic modulus, Poisson's ratio, permeability coefficient, and fracture toughness of the soil. All physical quantities are standardized, for example, length units are converted to meters, stress and modulus units to Pascals, and density units to kilograms per cubic meter. A coordinate transformation algorithm is then used to map the grouting hole coordinates from the tunnel's local coordinate system to the global geological coordinate system. The calculation formula is as follows: ,in, The transformed global x-coordinate. The transformed global ordinate. This represents the rotation angle of the tunnel cross-section relative to the global coordinate axes. Let x be the x-coordinate of the grouting hole in the local coordinate system. Let be the ordinate of the grouting hole in the local coordinate system. Let x be the x-coordinate of the tunnel center in the global coordinate system. The vertical coordinate of the tunnel center in the global coordinate system is used. After coordinate registration, the computational domain is discretized using the Delaunay triangulation algorithm. The upper limit of the grid size is set to 0.5 meters. For the dense area within 0.5 meters around the grouting holes, the upper limit of the grid size is set to 0.05 meters. A grid topology structure composed of multiple triangular elements is generated. All generated grid nodes are traversed and assigned unique integer index numbers. A topological connection matrix of element index and node index is established. Each row of this matrix stores the index numbers of the three nodes contained in a single element. The geological physical property parameters are mapped to the corresponding grid elements, such as the elastic modulus. Compared with Poisson Values ​​are assigned to each triangular element to form a joint mesh of subway tunnel geometric data and geological physical property parameters.

[0043] Based on the joint mesh of subway tunnel geometric data and geological physical property parameters, a scalar field with values ​​ranging from 0 to 1 is defined as the phase field variable. Where a value of 0 represents the soil in a completely intact state, and a value of 1 represents the soil in a completely fractured state, a total energy functional is constructed based on the variational principle, which includes elastic potential energy, fracture surface energy, and work done by fluid pressure. The expression for this total energy functional is: ,in, Represents the total energy functional. Represents the computational region. To prevent numerical singularities, the residual stiffness coefficient is set to a value of [value missing]. , The elastic strain energy density of the undamaged material. For strain tensor, This represents the critical energy release rate of the soil strata. The regularization length scale parameter, used to control crack spread width, is set to twice the average mesh size. Let be the gradient vector of the phase field variables. For the boundary region where fluid pressure is applied, The magnitude of the grouting fluid pressure, It is a displacement vector. Given the boundary outward normal vector, traverse all nodes in the mesh and calculate the Euclidean distance from each node to the center of the grouting hole. If the distance Smaller than the preset grouting hole radius If the value is 0.05 meters, then the node is marked as the grouting hole pressure boundary node. For nodes on the outer contour of the computational domain, they are marked as fixed support boundary nodes. These marking information are stored in the boundary condition list file to form the phase field variable initialization rules and boundary condition set.

[0044] Based on the initialization rules for phase field variables and the set of boundary conditions, an array space is allocated in computer memory to store the degree-of-freedom vectors of all nodes in the field. Each node is assigned three degrees of freedom: x-direction displacement, y-direction displacement, and phase field variable degrees of freedom. The initial values ​​of the phase field variables for all nodes are uniformly set to 0, indicating that the formation is undamaged at the initial moment. The boundary condition list file is read, and nodes marked as grouting hole pressure boundaries are identified. No mandatory assignment of values ​​is applied to the phase field variable degrees of freedom of these nodes, but nodal force loads are applied in subsequent load steps. Nodes marked as fixed support boundaries are identified, and their x-direction and y-direction displacement degrees of freedom are forcibly constrained to 0. The global stiffness matrix assembly process is executed, traversing each element and calculating the element stiffness matrix. It is accumulated into the global stiffness matrix. The corresponding position is calculated using the following formula: ,in, The global stiffness matrix, The total number of grid cells. For the first The degree-of-freedom extraction matrix of each unit. For the first The element stiffness matrix of each element. To extract the transpose of the matrix, and simultaneously initialize the global displacement vector. and phase field vector The zero vector represents the initial pressure value applied to the grouting hole boundary. Convert to equivalent nodal force vector The numerical model containing geometric information, material properties, boundary constraints and initial loads is constructed, and the initial field with displacement degree of freedom and phase field variable degree of freedom is generated, forming the initial phase field boundary state.

[0045] The steps for obtaining the permeability evolution distribution field are as follows:

[0046] The initial phase field boundary state is invoked, and the displacement degrees of freedom and phase field variables are updated sequentially within a fixed time step. The spatial distribution of displacement field and spatial distribution of phase field variables are recorded according to the time step, and the one-to-one correspondence between the two types of distributions is maintained according to the node index to obtain the update results of displacement field and phase field variables.

[0047] Based on the update results of displacement field and phase field variables, the phase field variables are compared with the threshold cell by cell. For cells whose phase field variables are greater than the threshold, the formation permeability coefficient is assigned the maximum limit value. For cells whose phase field variables are not greater than the threshold, the original permeability coefficient is maintained, thus generating a permeability evolution distribution field.

[0048] Specifically, the initial phase field boundary state is invoked, and the total duration of time integration is set. The time step is 100 seconds. 0.5 seconds, total time steps The calculation consists of 200 steps. An interleaved solution algorithm is used to iteratively calculate the coupled physical fields. In each time step, the phase field variables are first fixed, the displacement field equilibrium equations are solved, and the Newton-Raphson iterative method is used to update the displacement degrees of freedom, thus constructing a linearized system of equations. ,in, Here is the tangent stiffness matrix under the current state. It is the displacement increment vector. The residual vector, formed by the difference between the external load vector and the internal force vector, is used to accumulate the solved displacement increments to obtain the displacement field of the current iteration step. until the norm of the residual vector is less than the preset convergence tolerance. The convergence tolerance Set as the initial residual norm The historical maximum tensile elastic energy density was then calculated based on the updated displacement field. The calculation formula is: ,in Let be the strain energy density function of the stretched portion. For a moment Given the strain tensor, fix the displacement field and solve the phase field evolution equation, then assemble the phase field stiffness matrix. With phase field load vector Solve the equation ,in Let be the vector of phase field variables at the nodes to be determined. Includes the regularization length parameter and crack surface energy The determined diffusion and reaction terms, Based on historical energy density After the source term of the driving function completes the convergence calculation of the current time step, it extracts the displacement vector components and phase field scalar values ​​of all nodes and stores them into the corresponding two-dimensional array matrix according to the global index number of the nodes. The row index of the matrix corresponds to the node number, and the column index corresponds to the current time step number, ensuring a strict mapping between spatial location and physical quantity values, and obtaining the update results of displacement field and phase field variables.

[0049] Based on the update results of displacement field and phase field variables, each triangular element in the joint mesh is traversed, the phase field variable values ​​of each vertex of the element are extracted, and the average phase field variable at the center of the element is calculated. Set a phase field threshold for determining complete soil fracture. The threshold is set based on the material stiffness degradation function. The critical point, where To prevent numerically singular values ​​of the residual stiffness coefficient A crack is considered to be fully penetrating when its stiffness degrades to less than 1% of its initial stiffness. Solve for the threshold Set to 0.9, if the average phase field variable of the current unit is... If the value is greater than 0.9, the element is determined to be a fracture channel element, and its corresponding maximum limiting permeability coefficient is calculated. The calculation formula uses the equivalent form of the cube law. ,in The equivalent hydraulic fracture width, estimated based on the phase field distribution characteristics, is taken as 0.1 times the average side length of the element. For example, when the element side length is 0.1 meters... Taking 0.01 meters, we obtain Approximately The maximum limit value is assigned to the permeability attribute of the current unit, calculated in square meters. If the average phase field variable of the current unit... If the value is not greater than 0.9, the unit is determined to be a matrix unit, and its initial permeability coefficient is maintained. Unchanged, for example, remaining unchanged After traversing all units and completing the assignment operation on the order of square meters, a scalar field data structure containing the permeability values ​​of all units in the field is constructed, and a permeability evolution distribution field is generated.

[0050] The steps to obtain the grouting splitting propagation topology are as follows:

[0051] Based on the permeability evolution distribution field, candidate connected unit sequences of fracture topology are extracted, fluid velocity and pressure distribution of candidate connected units are calculated, energy functional values ​​are accumulated, and crack propagation paths are determined based on the connected sequence corresponding to the minimum energy functional value. Crack propagation paths are merged into fracture topology to form grouting splitting propagation topology.

[0052] Specifically, based on the permeability evolution distribution field, a breadth-first search algorithm is used to scan the entire grid cells to identify the permeability coefficient that equals the maximum limit value. The set of elements is used to establish an adjacency matrix based on the common node relationships of the elements. All connected element clusters are analyzed as candidate connected element sequences for the fracture topology. For each connected sequence, a local fluid flow finite element model is constructed, and the steady-state Darcy flow equation is solved. ,in For Hamiltonian operators, Per unit penetration rate The dynamic viscosity of the grout is set to 0.02 Pascals per second. For fluid pressure, the pressure distribution field and velocity vector distribution within the connected sequence are obtained by solving, and then the energy functional value of each candidate expansion path is calculated. The calculation formula is: ,in The elastic strain energy density stored within the cell. For the integration region, This represents the surface area of ​​the crack. For fracture surface energy, The fluid pressure within the fracture. This represents the jump in displacement as the crack surface opens. Given the crack surface normal vector, the energy functional values ​​corresponding to different candidate paths are compared. Based on the principle of minimum potential energy, the connected sequence corresponding to the minimum energy functional value is selected as the real physical crack propagation path at the current time step. The determined propagation path unit is marked and incorporated into the global crack topology list, and the geometric connection relationship of the crack network is updated to form the grouting splitting propagation topology.

[0053] The steps for obtaining the virtual displacement response field are as follows:

[0054] The grouting splitting expansion topology is invoked, fluid pressure distribution data is read, the projection relationship between the node coordinates of the tunnel segment ring structure and the fracture surface element is established, the pressure values ​​are summarized according to the projection weight and converted to the node degrees of freedom to form the pressure load field of the tunnel segment ring structure node.

[0055] Based on the pressure load field of the tunnel segment ring structure nodes, the nonlinear tangent stiffness matrix under the current loading condition is extracted, a virtual disturbance vector with controlled amplitude is constructed and superimposed on the existing load field, and the virtual displacement of each node is obtained by solving, thus generating a virtual displacement response field.

[0056] Specifically, the grouting and fracturing expansion topology is invoked to read fluid pressure distribution data. First, each cell storing fluid pressure in the fracture topology is traversed, and the three-dimensional spatial coordinates of the geometric center of the cell and the corresponding pressure scalar value are extracted to construct a source data point set containing spatial coordinates and pressure values. The KD-tree spatial indexing algorithm is then used to structure the source data point set, and a spatial search radius is set. The radius is set to 1.5 times the thickness of the tunnel segment, i.e., 0.525 meters, where the segment thickness is set to 0.35 meters. The model iterates through each finite element node in the tunnel segment ring structure model, using the current node coordinates as the center of the sphere. Within the specified range, a KD-tree search is performed to find the center point of all fracture elements falling within that range. If the search result is empty, the node is determined not to be directly affected by grouting pressure. If the search result is not empty, the equivalent fluid pressure of the node is calculated using the inverse distance weighted interpolation method. The calculation formula is as follows: ,in, The desired segment node pressure value is... This represents the total number of fracture elements found within the search range. For the first The fluid pressure value of each fracture element. For the segment node and the first To prevent the denominator from being zero, a local minimum is added to the distance term to calculate the Euclidean distance between the centers of each fracture element. After obtaining the nodal pressure scalar, calculate the associated force-bearing area of ​​each node on the outer surface of the segment. This area is obtained by calculating the area of ​​the dual element surrounding the node, converting the pressure scalar into a nodal force vector, using the following formula: ,in, For nodal force vectors, Let be the unit vector of the outward normal to the surface where the node is located. For the force-bearing area attached to the node, after traversing all the segment nodes to complete the calculation, all non-zero node force vectors are assembled into a global load vector. The projection relationship between the node coordinates of the tunnel segment ring structure and the crack surface element is established. The pressure values ​​are summarized according to the projection weight and converted to the node degrees of freedom to form the pressure load field of the tunnel segment ring structure node.

[0057] Based on the pressure load field at the nodes of the tunnel segment ring structure, the global stiffness matrix storage area after the convergence of the current nonlinear iteration step is accessed in the finite element solver. The sparse matrix structure updated by tangent stiffness is extracted. This matrix includes contributions from material nonlinearity and geometric nonlinearity. To evaluate the stability of the current equilibrium state, a virtual perturbation vector with the same degree of freedom dimension as the structure is constructed. The vector is not constructed entirely randomly, but is based on the current load field. The modulus is set with a small perturbation amplitude, and the perturbation coefficient is set. for That is, a perturbation ratio of one ten-thousandth, which generates a unit random vector whose components follow a standard normal distribution through a random generating function. The formula for calculating the virtual disturbance vector is: ,in, This is the constructed virtual perturbation vector. The disturbance coefficient is... Let Euclidean norm be the total pressure load vector currently applied to the segment. For the generated random vector, Given the magnitude of the random vector, this virtual perturbation vector is superimposed onto the existing load field. Essentially, without changing the tangent stiffness matrix, this involves solving a linearized incremental equation and constructing a system of linear equations. ,in, This is the currently extracted nonlinear tangent stiffness matrix. Let be the virtual displacement increment vector to be determined. As a virtual perturbation vector, the equation system is numerically solved using the preprocessed conjugate gradient method or a sparse direct solver. Since the stiffness matrix may be close to singular, the relative residual tolerance of the solver is set to [value missing]. The displacement increment response of all nodes under virtual disturbance is obtained by solving the problem. These displacement increments are arranged according to the node degree of freedom index. The nonlinear tangent stiffness matrix under the current loading condition is extracted. A virtual disturbance vector with controlled amplitude is constructed and superimposed on the existing load field. The virtual displacement of each node is obtained by solving the problem, and a virtual displacement response field is generated.

[0058] The steps for obtaining the condition number index of the stiffness matrix are as follows:

[0059] Based on the virtual displacement response field, the stiffness components associated with the virtual displacement response field are located and formed into an analysis submatrix. Singular value decomposition is performed to obtain the maximum and minimum singular values. The ratio of the maximum and minimum singular values ​​is calculated as an index value to obtain the condition number index of the stiffness matrix.

[0060] Specifically, based on the virtual displacement response field, a set of key nodes in the virtual displacement response whose displacement modulus exceeds a preset response threshold is identified. This response threshold is set to 5% of the maximum virtual displacement modulus. Sub-blocks corresponding to the rows and columns of the degrees of freedom of these key nodes are extracted from the global nonlinear tangent stiffness matrix to form a reduced stiffness matrix for stability analysis. Singular value decomposition (SVD) is performed on the reduced stiffness matrix, following the standard form in numerical algebra. ,in, For the stiffness submatrix used in the analysis, It is a left singular vector matrix. This is the transpose of the right singular vector matrix. For a diagonal matrix, its diagonal elements That is, it is a singular value, and satisfies From the diagonal matrix Extracting the maximum singular value directly from the middle and minimum singular value To avoid computational overflow caused by the minimum singular value approaching machine precision, a lower limit for numerical truncation is set. ,like Then a mandatory order The condition number index of the stiffness matrix is ​​calculated using the two extracted extreme values. The calculation formula is: ,in, The condition number index of the stiffness matrix. For the maximum singular value, The minimum singular value is the index that physically reflects the anisotropy of the segment structure's resistance to deformation and its proximity to instability. The stiffness components associated with the virtual displacement response field are located and formed into an analysis submatrix. Singular value decomposition is performed to obtain the maximum and minimum singular values. The ratio of the maximum and minimum singular values ​​is calculated as the index value, thus obtaining the stiffness matrix condition number index.

[0061] The steps for obtaining the comprehensive instability risk factors are as follows:

[0062] Call the stiffness matrix condition number index, read the current pressure field data, aggregate the node pressure values ​​according to the segment ring number, combine the stiffness matrix singular value information, establish a load-stiffness characteristic correspondence table for each segment ring, and form a comprehensive calculation input set;

[0063] Based on the comprehensive calculation input set, the comprehensive instability risk factor for each segment ring is calculated using the following formula:

[0064] ;

[0065] in, For the first The comprehensive instability risk factors of each tunnel segment ring, For the first Total pressure of each segment ring The equivalent radius of the segment ring, The current nonlinear tangent stiffness matrix The minimum singular value, This represents the global nonlinear tangent stiffness matrix under the current loading state. This is the reference nonlinear tangent stiffness matrix under initial loading conditions. ,in Let be the virtual displacement response field vector, representing the vector in the matrix. Virtual strain energy under action, , which serves as the baseline virtual strain energy.

[0066] Specifically, the process involves calling the stiffness matrix condition number index to read the current pressure field data. First, it accesses the real-time fluid pressure distribution array stored in memory. The index of this array strictly corresponds to the node number of the finite element mesh. Using a pre-built spatial hash mapping table, tens of thousands of nodes across the entire field are mapped to their corresponding tunnel segment ring numbers according to their spatial coordinates. For example, all nodes with longitudinal coordinates between 100 meters and 101.5 meters are marked as associated nodes of segment ring number 50. All identified segment ring objects are traversed. For each segment ring, the pressure scalar value (unit: Pascal) in its associated node set is extracted. A numerical integration algorithm is used to perform area-weighted summation of the pressures on these nodes. For each node, its associated area weight (unit: square meters) on the outer surface of the segment is obtained. This weight is determined by the area of ​​the dual element during the mesh generation stage. The node pressure value is multiplied by the associated area to obtain the nodal force (unit: Newton). Finally, the nodal forces of all nodes on the segment ring are calculated. The scalar values ​​are accumulated to aggregate the total pressure load data acting on the specific segment ring. Simultaneously, the singular value decomposition results of the stiffness matrix corresponding to the current segment ring are read from the output cache of the stiffness analysis module. In particular, the minimum singular value representing the weakest resistance direction of the structure is extracted, along with the baseline virtual strain energy and the current virtual strain energy data calculated in the previous steps. A dynamic structure array is constructed in memory, with the ring number of each segment ring as the primary key. The aggregated total pressure value, the corresponding minimum singular value, and the two types of virtual strain energy values ​​are filled in sequentially. During the filling process, double-precision floating-point verification of the data type is performed to ensure that there are no invalid values ​​(NaN) or infinite values ​​(Inf) caused by calculation overflow. For segment rings with missing data, interpolation is performed to complete the data or they are marked as skipped. Finally, the organized structured data rows are serialized into a standardized comprehensive calculation input set, and a load-stiffness characteristic correspondence table for each segment ring is established to form the comprehensive calculation input set.

[0067] In the formula for calculating the comprehensive instability risk factor, the first term measures the proportion of the current external load to the current remaining stiffness of the structure, directly reflecting the mechanical load rate of the structure. The second term uses the rate of change of virtual strain energy as an exponential correction term, which nonlinearly amplifies the impact of stiffness degradation on risk. That is, when the stiffness of the structure decreases due to damage, resulting in a change in strain energy under the same virtual disturbance, the risk factor will increase exponentially, thus sensitively capturing the precursors of brittle instability.

[0068] Representing the The total radial grouting pressure currently borne by each segment ring is expressed in Newtons (N). The steps for obtaining this parameter are as follows: First, collect the data in real-time or extract the data acting on the segment ring at the moment of impact. Fluid pressure scalar at all finite element nodes on the outer surface of each segment ring Simultaneously read each node Corresponding effective force-bearing area This area is determined by the element geometry and topology during mesh generation. For example, for the 50th segment ring, which contains 200 stress nodes, the average nodal pressure is... Pa (i.e., 0.5 MPa), the average attached area of ​​each node is If m², then the result of this parameter calculation is: N.

[0069] The equivalent geometric radius of the tunnel segment ring, expressed in meters (m), is obtained as follows: Based on the design and construction drawings or Building Information Modeling (BIM) data of the subway tunnel, the cylindrical fitting radius of the tunnel segment structure is extracted. For a standard circular tunnel cross-section, its outer diameter is directly taken. For non-circular cross-sections or those undergoing convergence deformation, the equivalent radius is obtained by performing a circle fitting on the point cloud of the segment's outer contour using the least squares method. For example, for a standard subway shield tunnel section with a designed outer diameter of 6.2 meters, this parameter... It is set to 3.1 m.

[0070] Represents the global nonlinear tangent stiffness matrix under the current loading state. The minimum singular value, expressed in Newtons per meter (N / m), is obtained by: calling the singular value decomposition result generated in the aforementioned stiffness matrix condition number index calculation step, and directly reading the diagonal matrix. The smallest non-zero diagonal element in the mean square root represents the tangent stiffness modulus of the structure under the most easily deformable mode, reflecting the structure's minimum resistance to instability. For example, considering the actual stiffness characteristics of concrete segment structures and the existing local microcrack damage, the smallest singular value obtained through SVD decomposition is... N / m.

[0071] This represents the baseline virtual strain energy under initial loading conditions, expressed in joules (J). The steps to obtain this parameter are: call the baseline stiffness matrix for the initial simulation stage before damage occurs. and the unit virtual displacement response vector generated in the current step. Using the quadratic form formula Matrix multiplication yields this parameter, which serves as a benchmark for measuring the original energy storage capacity of a structure. For example, by substituting the initial stiffness matrix and the current virtual displacement vector into the calculation, the benchmark virtual strain energy value is obtained as follows: J.

[0072] This parameter represents the virtual strain energy under the current loading state, expressed in joules (J). The steps to obtain this parameter are: using the currently updated nonlinear tangent stiffness matrix... This matrix already incorporates the effects of material damage and geometric nonlinearity, and is the same virtual displacement response vector used to calculate the baseline energy. Perform quadratic operations, that is This quantifies the energy storage characteristics of the current structural state. For example, since the structural stiffness decreases due to damage, the current virtual strain energy value obtained by substituting this value into the calculation is... J.

[0073] Calculations based on parameters:

[0074] First, calculate the load stiffness ratio:

[0075] ;

[0076] Next, calculate the stiffness degradation exponent:

[0077] ;

[0078] ;

[0079] Finally, calculate the overall instability risk factor:

[0080] ;

[0081] The result indicates that the comprehensive instability risk factor of the current No. 50 segment ring is 0.0098. Although the absolute value is small, for a dimensionless risk index, this value is close to the preset critical level, reflecting that under the current grouting pressure, the stiffness reserve of the segment ring is being consumed, and the nonlinear effect caused by damage is beginning to appear. It is necessary to enter the subsequent threshold determination stage to confirm whether the shutdown mechanism is triggered.

[0082] The steps for obtaining the grouting control command signal are as follows:

[0083] Based on the comprehensive instability risk factor, threshold judgment is performed on each segment ring. Segment rings whose comprehensive instability risk factor exceeds the preset safety threshold are marked as critical instability states. The corresponding grouting hole numbers are extracted, and grouting control command signals are generated.

[0084] Specifically, based on the comprehensive instability risk factor, the system reads the safety threshold parameter pre-stored in the system configuration file. The setup method is as follows: Collect data from 100 historical grouting operations in this tunnel section, screen out "critical failure" samples where minor segment misalignment or crack development has occurred, calculate the comprehensive instability risk factor corresponding to these samples, obtain a set of critical factor sequences, and calculate the average value of the sequence. and standard deviation To ensure sufficient safety redundancy, the lower confidence interval is used as the safety limit, and the following settings are made: The currently calculated comprehensive instability risk factor (For example, the aforementioned calculation result of 0.0098) is compared with the preset safety threshold of 0.009. The judgment logic is as follows: if If this is the case, the current segment ring is determined to be in an unstable state. In the logical judgment, because... If the conditions are met, the segment ring is immediately marked as "critical instability" and the alarm flag is activated. Then, the geometric mapping database of the segment ring and grouting holes is queried. Based on the segment ring number, the grouting hole numbers that are having a pressure impact on the ring and the adjacent area (e.g., grouting holes No. 3 and No. 4) are retrieved in reverse index. Stop grouting command signals are generated for these specific grouting holes. The signal contains the target device ID and action type code (STOP). The signal is sent to the grouting machine controller via the industrial Ethernet bus. The segment ring with the comprehensive instability risk factor exceeding the preset safety threshold is marked as critical instability state. The corresponding grouting hole number is extracted and a grouting control command signal is generated.

Claims

1. A simulation system for grouting process in subway tunnels based on digital twins, characterized in that, The system includes: The grouting environment construction module is used to acquire the geometric data and physical property parameters of the subway tunnel, define continuous phase field variables to characterize the soil integrity state, input the multiphysics solution environment, apply the grouting hole pressure boundary conditions, and obtain the initial phase field boundary state. The fracture propagation simulation module is used to call the initial phase field boundary state, execute an interleaved solution sequence within a fixed time step to update the displacement field and phase field variables, assign the formation permeability coefficient in the spatial region where the phase field variable is higher than the threshold to the maximum limit value, generate a permeability evolution distribution field, and calculate the fluid velocity and pressure distribution in the updated fracture topology based on the permeability evolution distribution field to obtain the grouting fracture propagation topology; The segment stiffness analysis module is used to map the fluid pressure distribution data in the grouting splitting extension topology to the nodes of the tunnel segment ring structure, extract the nonlinear tangential stiffness matrix of the segment ring structure under the current loading condition, apply a virtual disturbance vector to the existing load field, generate a virtual displacement response field, perform singular value decomposition on the stiffness components associated with the virtual displacement response field, and calculate the stiffness matrix condition number index. The grouting control decision module is used to call the stiffness matrix condition number index to quantify the ability of the segment ring to resist grouting pressure fluctuations, calculate and generate a comprehensive instability risk factor, compare the comprehensive instability risk factor with a preset safety threshold, identify the critical instability state segment ring where the comprehensive instability risk factor exceeds the preset safety threshold, trigger a stop grouting command signal for the corresponding grouting hole, and obtain a grouting control command signal.

2. The simulation system for grouting process in subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the initial phase field boundary state are as follows: Read the geometric data and physical property parameters of the subway tunnel, extract the cross-sectional dimensions, centerline coordinates, grouting hole coordinates, density, elastic modulus, Poisson's ratio, permeability coefficient and fracture toughness, unify the units and complete the coordinate registration, establish the node number and unit relationship, and form a joint grid of the geometric data and physical property parameters of the subway tunnel. Based on the combined grid of the subway tunnel geometric data and the physical property parameters of the strata, phase field variables are set to indicate the soil integrity state, a summation structure of elastic potential energy, fracture surface energy and fluid pressure work is constructed, the grouting hole pressure boundary and fixed support boundary are registered, and a set of phase field variable initialization rules and boundary conditions is formed. According to the initialization rules and boundary conditions set of the phase field variables, the initial values ​​of the phase field variables are assigned to each node in the solution domain and the grouting hole pressure boundary is applied to the boundary nodes. The state assembly is performed to generate the initial field of displacement degree of freedom and phase field variable degree of freedom, thus forming the initial phase field boundary state.

3. The simulation system for grouting process in subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the permeability evolution distribution field are as follows: The initial phase field boundary state is invoked, and the displacement degrees of freedom and phase field variables are updated sequentially within a fixed time step. The spatial distribution of displacement field and spatial distribution of phase field variables are recorded according to the time step, and the one-to-one correspondence between the two types of distributions is maintained by the node index to obtain the update results of displacement field and phase field variables. Based on the update results of the displacement field and phase field variables, the phase field variables are compared with the threshold cell by cell. For cells whose phase field variables are greater than the threshold, the formation permeability coefficient is assigned the maximum limit value. For cells whose phase field variables are not greater than the threshold, the original permeability coefficient is maintained, thus generating a permeability evolution distribution field.

4. The simulation system for grouting process in subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the grouting splitting propagation topology are as follows: Based on the permeability evolution distribution field, candidate connected unit sequences of the fracture topology are extracted, the fluid velocity and pressure distribution of the candidate connected units are calculated, the energy functional values ​​are accumulated, and the crack propagation path is determined by the connected sequence corresponding to the minimum energy functional value. The crack propagation paths are merged into the fracture topology to form a grouting splitting propagation topology.

5. The simulation system for grouting process in subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the virtual displacement response field are as follows: The grouting splitting expansion topology is invoked, fluid pressure distribution data is read, the projection relationship between the node coordinates of the tunnel segment ring structure and the fracture surface element is established, the pressure values ​​are summarized according to the projection weight and converted to the node degrees of freedom to form the pressure load field of the tunnel segment ring structure node. Based on the pressure load field of the tunnel segment ring structure nodes, the nonlinear tangent stiffness matrix under the current loading condition is extracted, a virtual disturbance vector with controlled amplitude is constructed and superimposed on the existing load field, the virtual displacement of each node is obtained by solving, and a virtual displacement response field is generated.

6. The simulation system for grouting process in subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the condition number index of the stiffness matrix are as follows: Based on the virtual displacement response field, the stiffness components associated with the virtual displacement response field are located and formed into an analysis submatrix. Singular value decomposition is performed to obtain the maximum and minimum singular values. The ratio of the maximum and minimum singular values ​​is calculated as an index value to obtain the stiffness matrix condition number index.

7. The simulation system for grouting process of subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the comprehensive instability risk factor are as follows: The stiffness matrix condition number index is invoked to read the current pressure field data, aggregate the node pressure values ​​according to the segment ring number, and combine the stiffness matrix singular value information to establish a load-stiffness characteristic correspondence table for each segment ring, forming a comprehensive calculation input set; Based on the comprehensive calculation input set, calculate the comprehensive instability risk factor for each segment ring.

8. The simulation system for grouting process in subway tunnels based on digital twins according to claim 1, characterized in that, The steps for obtaining the grouting control command signal are as follows: Based on the comprehensive instability risk factor, threshold judgment is performed on each segment ring. Segment rings whose comprehensive instability risk factor exceeds the preset safety threshold are marked as critical instability states. The corresponding grouting hole numbers are extracted, and grouting control command signals are generated.

Images