System and method for visualizing state of muck in shield sealed cabin

Abstract

The invention relates to the technical field of particle state detection, in particular to a shield sealed cabin muck state visualization system and method.According to the shield sealed cabin muck state visualization system and method, through stress direction difference matching of an SIMPLE algorithm, inter-particle shear strength distribution is converted into a computable parameter set, quantitative extraction of shear characteristics can be achieved in a complex stress environment, and the detection accuracy is improved. Stable mechanical boundary conditions are provided for flow state calculation, a numerical coupling model of momentum exchange is established, judgment of turbulent fluctuation and slow flow partitions in a flow field does not depend on empirical thresholds any more, characteristic value changes of a mathematical matrix serve as partition bases, the adaptability and accuracy of flow field partition are enhanced, and the flow field partition efficiency is improved. A main load channel in a particle contact network can be automatically identified through combination of an intensity matrix and a centrality algorithm, the connectivity is used as a core calculation index for path sorting, the structural hierarchy of load transmission is clearly presented, and a stable matrix is formed based on weighted superposition of a path shearing and supporting force ratio. Extreme point extraction enables low-stability region identification to have quantification precision.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of particle state detection, in particular to a shield sealing cabin inner muck state visualization system and method. BACKGROUND

[0002] The technical field of particle state detection aims to establish a quantitative determination model of particle medium state by measuring and analyzing the physical characteristics of granular substances, to describe their morphology, distribution, density, water content and flowability under the conditions of stress, motion and environmental change, to visually characterize and numerically calculate the macroscopic morphology and microscopic state of the particle medium, and to realize real-time monitoring and evaluation of the running state of the particle medium.

[0003] The purpose of the shield sealing cabin inner muck state visualization system is to realize real-time monitoring and quantitative determination of the muck state inside the sealing cabin during shield construction, to obtain physical characteristic data such as the morphology, distribution, flowability and density of the muck inside the sealing cabin, to present the muck state in a visual form which cannot be directly observed, to provide accurate muck state feedback for shield construction, to assist the automatic control system in adjusting the tunneling parameters, to avoid water inrush, collapse and equipment damage problems caused by abnormal muck state, and to realize intelligent monitoring and efficient control of the construction process.

[0004] The existing particle state detection technology relies on single-point sensing signals and empirical parameter mapping in muck monitoring, and the data is distributed discretely, which is difficult to continuously reflect the overall response of the particle group under stress changes. The determination of muck morphology and density is mainly based on macroscopic indicators, and the stress chain structure characteristics between microscopic particles are not reflected, which leads to the inability to identify potential unstable trends in the early stage of abnormal state formation. The flow field and particle dynamics parameters are usually treated separately, lacking a unified calculation system, resulting in misjudgment of the muck state in the transition zone between turbulence and slow flow. When the tunneling rate changes and the cabin pressure fluctuates, the feedback lags behind. Since there is no spatial coupling determination mechanism for compaction and flow behavior, the evolution trend of the muck dense area can only be inferred from the late results, lacking real-time adjustment basis. When the load transmission path deviates, it is difficult to track the direction of contact force transmission between particles, and the identification accuracy of unstable areas is insufficient, which may lead to risks such as local collapse, water inrush and cutter damage. SUMMARY

[0005] The purpose of the present application is to solve the shortcomings in the prior art and to provide a shield sealing cabin inner muck state visualization system and method.

[0006] In order to achieve the above purpose, the present application adopts the following technical scheme: a shield sealing cabin inner muck state visualization system comprises: The particle dynamic response analysis module: based on the particle group in the shield sealing cabin, the particle velocity, displacement and stress are calculated in time sequence, the energy change is compared by step length, the abnormal stress area is identified, the SIMPLE algorithm is used, the stress direction difference between the particles in the abnormal area is matched and the shear strength distribution is extracted, and the particle stress distribution parameter set is formed; The cabin flow coupling calculation module: based on the particle stress distribution parameter set, the node system in the cabin is established, the flow velocity difference and pressure gradient are calculated, the momentum exchange value is merged to judge the turbulent and slow flow partition, the local shear and flux are matrix combined, and the cabin flow field response data is obtained; The muck density evolution determination module: based on the cabin flow field response data, the flow velocity gradient and density change rate are analyzed, the section standard deviation is calculated to judge the compaction state, the boundary gradient is interpolated and the spatial sequence is reorganized, and the muck dense structure parameter field is obtained; The contact transfer network identification module: based on the muck dense structure parameter field, the particle contact force and displacement deviation are extracted, the strength matrix is constructed and the node connectivity is analyzed, the centrality algorithm is used to identify the main load channel and sort the transmission path, and the particle load transfer structure set is established; The muck structure stability calculation module: based on the particle load transfer structure set, the path shear and support force ratio are calculated to judge the instability area, the low stability area range and direction are weighted and superimposed, the extreme value point is extracted to form the stability matrix, and the muck structure stability determination result is established.

[0007] As a further scheme of the application, the particle stress distribution parameter set includes particle velocity distribution parameters, displacement increment parameters and shear stress difference parameters, the cabin flow field response data includes flow velocity vector parameters, density gradient parameters and momentum flux parameters, the muck dense structure parameter field includes muck density distribution parameters, compaction gradient parameters and structure continuity parameters, the particle load transfer structure set includes node connectivity parameters, main force chain path parameters and stress center parameters, and the muck structure stability determination result includes stability coefficient parameters, instability area range parameters and direction strength parameters.

[0008] As a further scheme of the application, the particle dynamic response analysis module includes: The particle motion measurement sub-module: based on the particle group in the shield sealing cabin, the velocity, displacement and stress information of the particles are continuously recorded, the position change of each particle in the time sequence is mapped as a motion trajectory, the particle acceleration is calculated by comparing the velocity change amount at adjacent time points, and the stress direction and size are recorded synchronously, the three-dimensional motion data table is formed after being summarized, and the particle motion basic data set is generated; An energy difference determination sub-module: based on the particle motion basic data set, a kinetic energy change value of each particle in a continuous time period is calculated, an energy sudden increase and sudden decrease region is marked as an abnormal region, a force difference and a direction difference between particles in the region are extracted, an energy change rate and a force deviation amplitude are compared, a particle group with a shearing effect is identified, a stress distribution table is arranged, and a particle energy difference parameter set is obtained; A force distribution construction sub-module: based on the particle energy difference parameter set, a SIMPLE algorithm is used to superimpose a shearing strength, a normal force and a direction vector in a spatial coordinate, a stress transmission ratio between regions is calculated, a force direction distribution curve is drawn, a gradient difference between adjacent regions is spatially spliced, a complete stress network diagram is formed, and a particle force distribution parameter set is generated.

[0009] As a further scheme of the present application, the SIMPLE algorithm first sets velocity components and pressure values of each node in an initial step, performs discretization calculation on a momentum term in a flow control equation, then uses a predicted velocity field to obtain a momentum residual, and corrects a pressure gradient to form a corrected pressure field, and then recalculates the node velocity components with the updated pressure values to obtain a corrected velocity field, and the difference between the new and old velocity and pressure is used for convergence judgment, and when the difference is lower than a set threshold, the iteration is ended, and the matching result of the velocity and the pressure is output, and the calculation of the particle force balance in the cabin is completed.

[0010] As a further scheme of the present application, the cabin flow coupling calculation module comprises: A node parameter construction sub-module: based on the particle force distribution parameter set, node regions are divided in the sealed cabin space, flow velocity, pressure and density parameters are assigned to each node, the velocity difference of adjacent nodes is compared and the pressure gradient is calculated, the node momentum information is arranged according to the spatial sequence, the numerical relationship of nodes in each layer is integrated, and a unified cabin node momentum data set is formed; A flow field difference calculation sub-module: based on the cabin node momentum data set, the momentum exchange amplitude of each node is compared to identify high-energy flow regions and slow flow regions, the shearing strength and flux value of different regions are fused, the result is mapped to the spatial coordinate, the flow direction change trend is extracted, a continuous flow distribution diagram is formed, and cabin flow field response data is generated.

[0011] As a further scheme of the present application, the slag densification evolution determination module comprises: A compaction feature recognition sub-module: based on the cabin flow field response data, the flow velocity gradient and the density change amplitude in the slag flow region are counted, the gradient difference of different sections is compared, the demarcation points of the compaction zone and the loose zone are located, the partition results are arranged according to the cabin position, a slag compaction distribution table is formed, and a slag compaction state parameter set is obtained. The compact structure generating submodule: based on the set of parameters of the compacted state of the slag, interpolation operation is performed on the boundary gradient value, the spatial sequence is rearranged to reflect the structural continuity, the node information of the compacted and loose areas is combined to generate a complete density matrix, and a three-dimensional structure mapping is formed after merging and arranging, so that a parameter field of the compacted structure of the slag is obtained.

[0012] As a further scheme of the present application, the contact transmission network identification module comprises: The contact feature extraction submodule: based on the compacted structure parameter field of the slag, the contact points between particles and their relative positions are identified, the normal force and tangential force components of each contact point are recorded, and the displacement offset is determined, the force data of adjacent particles is integrated through spatial aggregation, the high-strength contact area and the low-strength contact area are screened out, the corresponding data index is established, and a particle contact feature data set is generated; The connection path identification submodule: based on the particle contact feature data set, the contact points are gradually connected according to the continuity of the force direction and the contact strength, the particle node chain is formed, the main channel is determined by detecting the force consistency in the chain, the connectivity and hierarchical order between channels are identified, the structure path table is drawn, and a particle transmission path structure set is obtained; The load relationship modeling submodule: based on the particle transmission path structure set, the centrality algorithm is used to count the force value distribution of each path and the mutual support relationship between nodes, the load components are combined according to the path direction, the load concentration area and transmission strength between different channels are analyzed, the overall induction of the force direction relationship of each region is performed, a three-dimensional force topology diagram is established, and a particle load transmission structure set is formed.

[0013] As a further scheme of the present application, the centrality algorithm first establishes an adjacency matrix of the particle contact network, uses edge weight to represent the load strength between nodes, then traverses all paths, records the appearance frequency of nodes in each path, and uses it as the betweenness value, and then normalizes the node connectivity, betweenness value and adjacent node weight to calculate the centrality index, obtains the influence degree sorting of each node in the network, and draws a load relationship matrix according to the centrality value distribution.

[0014] As a further scheme of the present application, the slag structure stability calculation module comprises: The mechanical equilibrium calculation submodule: based on the particle load transmission structure set, the shear force and support force values of each path are extracted, the ratio of the two is calculated and the abnormal section is screened out, the low stability area is numbered and arranged into a force relationship table, which provides basic data for subsequent stability determination, and a path force ratio data set is obtained; The stable distribution evaluation sub-module: based on the path stress ratio data set, the range and direction information of the low stability area are summarized, the stability coefficient of each section is calculated and the extreme value node is determined, the extreme value area difference is processed to construct a stability matrix, the overall stress distribution structure in the cabin body is output, and the slag structure stability judgment result is established.

[0015] A shield sealing cabin internal slag state visualization method based on the above shield sealing cabin internal slag state visualization system is executed, comprising the following steps: S1: based on the shield sealing cabin internal slag particle group, record the particle velocity, displacement and stress, calculate the stress difference by comparing the time sequence change, extract the shear strength by using the SIMPLE algorithm balance iteration, and generate the particle stress distribution parameter set; S2: based on the particle stress distribution parameter set, establish the cabin node system, calculate the flow velocity difference and pressure gradient, correct the flow velocity and pressure by using the SIMPLE algorithm, judge the turbulent and slow flow partition and integrate the shear flux, and obtain the cabin flow field response data; S3: based on the cabin flow field response data, calculate the flow velocity gradient and density change rate, compare the section standard deviation to determine the compaction state, interpolate the boundary gradient to form a continuous sequence, and obtain the slag dense structure parameter field; S4: based on the slag dense structure parameter field, extract the particle contact force and displacement deviation, calculate the node weight by using the centrality algorithm to identify the main load channel, and establish the particle load transfer structure set; S5: based on the particle load transfer structure set, calculate the shear and support force ratio to identify the instability area, weight and integrate the range and direction of the low stability area, extract the extreme value node to form a stability matrix, and establish the slag structure stability judgment result.

[0016] Compared with the prior art, the advantages and positive effects of the present application are: In the present application, by introducing time sequence step comparison and energy change analysis in particle dynamics calculation, the dynamic tracking of particle group stress evolution is realized, the identification of abnormal stress area is converted from static threshold judgment to continuous time sequence judgment, and the sensitivity and response accuracy of stress abnormality detection are improved; In the present application, by matching the stress direction difference of the SIMPLE algorithm, the shear strength distribution between particles is converted into a calculable parameter set, which can realize quantitative extraction of shear characteristics in complex stress environment, provide stable mechanical boundary conditions for flow state calculation, and establish a numerical coupling model of momentum exchange through matrix operation of flow velocity difference and pressure gradient in the cabin, so that the judgment of turbulent and slow flow partition in the flow field is no longer dependent on empirical threshold, and the characteristic value change of mathematical matrix is used as the basis for partition, which enhances the adaptability and accuracy of flow field division; In the application, the combination of the strength matrix and the centrality algorithm enables the main load channel in the particle contact network to be automatically identified, the path ranking is based on the connectivity calculation index, the structural level of load transmission is clearly presented, the stable matrix is formed based on the weighted superposition of the path shear and support force ratio, and the low stability region identification has quantitative precision through extreme point extraction. BRIEF DESCRIPTION OF DRAWINGS

[0017] Figure 1 The system flowchart of the application is shown in the figure. Figure 2 The method step schematic diagram of the application is shown in the figure. DETAILED DESCRIPTION

[0018] In order to make the purpose, technical scheme and advantages of the application clearer, the application is further described in detail below in combination with the drawings and examples. It should be understood that the specific examples described herein are only used to explain the application and do not limit the application.

[0019] Example 1 Please refer to Figure 1 The application provides a technical scheme: a muck state visualization system in a shield sealing cabin includes: A particle dynamic response analysis module: based on the particle group in the muck in the shield sealing cabin, the particle velocity, displacement and stress are calculated in time sequence, the abnormal stress area is identified by comparing the energy change through the step length, the SIMPLE algorithm is used to match and extract the shear strength distribution of the stress direction difference between the abnormal area particles, and a particle stress distribution parameter set is formed; A cabin flow coupling calculation module: based on the particle stress distribution parameter set, a cabin node system is established, the flow velocity difference and pressure gradient are calculated, the momentum exchange value is merged to judge the turbulent and slow flow partition, the local shear and flux are matrix merged to obtain the cabin flow field response data; A muck densification evolution determination module: based on the cabin flow field response data, the flow velocity gradient and density change rate are analyzed, the section standard deviation is calculated to judge the compaction state, the boundary gradient is interpolated and the spatial sequence is reorganized to obtain the muck dense structure parameter field; A contact transmission network identification module: based on the muck dense structure parameter field, the particle contact force and displacement deviation are extracted, the strength matrix is constructed and the node connectivity is analyzed, the centrality algorithm is used to identify the main load channel and rank the transmission path, and a particle load transmission structure set is established; A muck structure stability calculation module: based on the particle load transmission structure set, the path shear and support force ratio are calculated to judge the instability area, the low stability area range and direction are weighted and superimposed, the extreme point is extracted to form a stable matrix, and a muck structure stability determination result is established.

[0020] The particle force distribution parameter set includes a particle velocity distribution parameter, a displacement increment parameter, and a shear stress difference parameter, the in-cabin flow field response data includes a flow velocity vector parameter, a density gradient parameter, and a momentum flux parameter, the muck dense structure parameter field includes a muck density distribution parameter, a compaction gradient parameter, and a structure continuity parameter, the particle load transmission structure set includes a node connectivity parameter, a main force chain path parameter, and a force center parameter, and the muck structure stability determination result includes a stability coefficient parameter, an instability region range parameter, and a directional strength parameter.

[0021] The particle dynamic response analysis module includes: The particle motion measurement submodule: based on the muck particle group in the shield sealed cabin, continuously records the velocity, displacement, and force information of the particles, maps the position changes of each particle in the time sequence into a motion trajectory, calculates the particle acceleration by comparing the velocity changes at adjacent time points, and synchronously records the force direction and size, and after summarizing, forms a three-dimensional motion data table to generate a particle motion basic data set; The energy difference determination submodule: based on the particle motion basic data set, calculates the kinetic energy change value of each particle in the continuous time period, marks the regions with sudden energy increases and drops as abnormal regions, extracts the force difference and direction difference between the particles in the regions, compares the energy change rate and the force deviation amplitude, identifies the particle group with shearing effect, organizes a stress distribution table, and obtains a particle energy difference parameter set; The force distribution construction submodule: based on the particle energy difference parameter set, uses the SIMPLE algorithm to superimpose the shear strength, normal force, and direction vector in the spatial coordinates, calculates the stress transmission ratio between each region, draws a force distribution curve, and spatially splices the gradient difference between adjacent regions to form a complete stress network diagram, and generates a particle force distribution parameter set; The particle motion measurement submodule: based on the muck particle group in the shield sealed cabin, uses the Euler-Lagrange hybrid solution algorithm to continuously record the particle velocity, displacement, and force information, sets the time step parameter to 0.002 seconds and the spatial resolution unit to 1.5×10 -3 The particle motion measurement submodule: based on the muck particle group in the shield sealed cabin, uses the Euler-Lagrange hybrid solution algorithm to continuously record the particle velocity, displacement, and force information, sets the time step parameter to 0.002 seconds and the spatial resolution unit to 1.5×10 Energy difference determination submodule: based on the particle motion basic data set, the kinetic energy change difference algorithm is used to calculate the kinetic energy change value of each particle in the continuous time period, the time period length is set to 0.01 second, the instantaneous kinetic energy of each particle is obtained by multiplying the mass by half of the square of the velocity, and then the energy change value is calculated by the kinetic energy difference of adjacent time periods. The absolute value of the energy change and the set threshold 5.0x10 -3 Joule compares and judges the energy mutation area, uses the logical judgment instruction to mark the energy sudden increase and sudden decrease particles and generates an abnormal area list, performs force difference and direction difference calculation between the particles in the abnormal area, the direction difference is obtained by using the unit direction vector angle function, the force difference is obtained by using the absolute value of the force value difference of the two particles, and the energy change rate and the force deviation amplitude are compared. After that, the strongest particle group is selected by the sorting command, the particle group is arranged in a matrix form to generate a stress distribution table, and the particle energy difference parameter set is output; Force distribution construction submodule: based on the particle energy difference parameter set, the SIMPLE algorithm is used to superimpose the shear strength, normal force and direction vector in the spatial coordinates, the density, viscosity, pressure and velocity four input variables are defined, the alternating update operation of the pressure field and the velocity field is performed, the pressure relaxation coefficient is set to 0.7, the velocity relaxation coefficient is set to 0.8, the loop iteration is used until the pressure residual is less than the set limit value, the stress transfer ratio between the calculation areas is calculated and recorded in the numerical matrix, the force distribution curve is drawn by using the drawing command, and then the gradient difference between the adjacent spatial units is performed. Three-dimensional splicing operation, complete spatial data field is generated by using coordinate grid command, gradient matrix after superposition is combined according to spatial direction, stress transfer matrix is constructed to obtain stress network diagram, and particle force distribution parameter set is generated.

[0022] SIMPLE algorithm, first, set the velocity component and pressure value of each node in the initial step, discrete calculation of the momentum term in the flow control equation is carried out, then the momentum residual is calculated by using the predicted velocity field, and the pressure gradient is corrected to form the corrected pressure field. The node velocity component is calculated again by using the updated pressure value to obtain the corrected velocity field. The difference between the new and old velocity and pressure is used for convergence judgment. When the difference is less than the set threshold, the iteration is ended, and the matching result of velocity and pressure is output, and the force balance calculation of particles in the cabin is completed. SIMPLE algorithm, according to the formula:

[0023] Among them: represents the predicted surface flux of the momentum equation of the slurry fluid in the sealed cabin, represents the Patankar momentum correction coefficient obtained by using the SIMPLE algorithm for discretization, represents the pressure correction term of the i-th spatial unit in the cabin, represents the pressure correction term of the i-th spatial unit in the cabin, Pressure correction term for space unit, This represents the numerical stability constant set to prevent iterative divergence. Represents the direction vector between the i-th and j-th units. The average shear strength of the projection, Indicates the first With the Spatial distance between unit centroids Indicates the first The normal force acting on the element, This represents the normal force acting on the j-th element. Indicates by the first Unit points to the first The unit direction vector of the element. This represents the anisotropy coefficient calculated based on the eigenvalue ratio of the tensor of the particle contact structure within the cabin. The shear dilatation geometric factor is used to describe the geometric coupling relationship between the particle volume fraction gradient and the contact orientation angle. The inelastic dissipation ratio of particles measures the degree of energy loss during a collision. This represents the micro-contact roughness coefficient of the particle surface to characterize the microscale friction enhancement effect. The local packing factor is used to reflect the density distribution of soil particles within the sealed compartment. This represents the shear strength weighting coefficient. Indicates the normal force weighting coefficient. This represents the anisotropy weighting coefficient. This represents the energy dissipation weighting coefficient. This represents the roughness weighting coefficient. This represents the stacking coefficient and weighting coefficient. This represents the calculated stress transfer ratio between the slag and soil regions within the sealed chamber. Execution process: First, a discrete control volume model of the slag fluid is established inside the sealed chamber, and the momentum equation is used for prediction and calculation. That is, the predicted mass flux values ​​on each control surface, and the pressure field of each spatial unit is corrected using the framework to obtain... and This is to reflect the pressure correction value of the slag inside the compartment under local stress difference, and at the same time, the Patankar coefficient is obtained based on the discrete scheme solution. To characterize the equilibrium relationship between momentum convection and diffusion, the positions of the centroids of adjacent units are then extracted in spatial coordinates for calculation. With direction vector The mean value is obtained by projecting the shear strength distribution along this direction. To characterize the flow tendency of the slag under shearing action, the normal force is then calculated based on the contact force sensing data of each unit node. With and the normal stress gradient is obtained by the difference of the two directions, then the eigenvalue ratio is calculated by the particle contact network structure tensor to get the anisotropy coefficient and the shear dilation geometric factor is constructed by combining the volume fraction gradient and the contact angle to reflect the local shear dilation coupling effect, and the inelastic dissipation ratio is calculated by the particle restitution coefficient and the energy loss ratio and its multiplication with corrects the flux attenuation characteristics, and the micro-contact roughness coefficient is obtained by the surface power spectrum density integration and are combined to reflect the micro-scale friction amplification effect, and then the local packing coefficient is calculated according to the inverse proportion of the Voronoi volume to correct the density influence, and then the design matrix is constructed by experimental and simulation data and the weight coefficient is determined by non-negative constrained ridge regression , the prediction error is minimized by cross-validation and the coefficient range is limited to 0 to 1, and all parameters are substituted into the formula to obtain , the stress transfer ratio matrix between regions is formed and spliced into a stress network diagram to realize the visual presentation of the stress state of the slurry in the sealed cabin.

[0024] The cabin flow coupling calculation module includes: Node parameter construction submodule: based on the particle stress distribution parameter set, the node area is divided in the sealed cabin space, the flow velocity, pressure and density parameters are assigned to each node, the velocity difference of adjacent nodes is compared and the pressure gradient is calculated, the node momentum information is arranged according to the spatial sequence, the numerical relationship of each layer of nodes is integrated, and a unified cabin node momentum data set is formed; Flow field difference calculation submodule: based on the cabin node momentum data set, the momentum exchange amplitude of each node is compared to identify high-energy flow areas and slow flow areas, the shear strength and flux value of different areas are fused, the results are mapped to the spatial coordinates, the flow direction change trend is extracted, and a continuous flow distribution map is formed to generate the cabin flow field response data; The node parameter construction submodule: based on the particle force distribution parameter set, the finite volume method solving algorithm is used to divide the node region in the sealed cabin space, the grid division step size is set to 0.05 meters for the space coordinate axes x, y and z respectively, a three-dimensional node matrix is generated and indexed, the flow velocity, pressure and density parameters of each node are assigned, the initial value of pressure is set to 101325 Pa, the value of density is 1850 kg / m3, the adjacent node parameters are read in sequence using the node index matrix, the velocity difference command is executed to obtain the velocity gradient by dividing the adjacent node velocity difference by the distance step, the pressure gradient matrix is obtained by dividing the adjacent node pressure difference by the space step using the pressure gradient calculation command, and the momentum value of each node is recorded as the density multiplied by the velocity vector using the loop command to traverse all nodes, the momentum information of all nodes is arranged in sequence according to the space, and the six-column matrix containing node number, coordinate, flow velocity, pressure, density and momentum is stored in the form of a three-dimensional array, the data splicing command is executed for each layer of nodes to integrate the numerical relationship of each layer into a unified data structure, and the cabin node momentum dataset is generated; The flow field difference calculation submodule: based on the cabin node momentum dataset, the turbulent flow vorticity difference algorithm is used to compare the momentum exchange amplitude of each node, first, the item-by-item momentum difference calculation command is executed on the node matrix to obtain the momentum change rate by dividing the adjacent node momentum difference by the time step 0.01 seconds, then the node energy level is judged according to the change rate and an identification table of high-energy flow area and slow flow area is established, the shear strength and flux value of different areas are fused, the two groups of data are superimposed in the form of weighted average using the matrix fusion command, and the flow field matrix is updated, the flux weight coefficient is set to 0.6, the shear strength weight coefficient is set to 0.4, the space mapping command is executed to map the fusion result to the three-dimensional space coordinate, the coordinate step is defined as 0.05 meters and arranged according to the node index number, the direction analysis command is executed on the mapped data to determine the flow direction change trend by the momentum vector difference of adjacent nodes, then the flow path is drawn in a continuous coordinate point sequence and a continuous flow distribution map is generated, and the cabin flow field response data is output.

[0025] The slag densification evolution judgment module includes: The compaction feature recognition submodule: based on the cabin flow field response data, the flow velocity gradient and density change amplitude in the slag flow area are counted, the gradient difference of different sections is compared, the demarcation point of the compaction area and the loose area is located, the partition result is arranged according to the cabin position, a slag compaction distribution table is formed, and a slag compaction state parameter set is obtained; The dense structure generation submodule: based on the slag compaction state parameter set, the boundary gradient value is subjected to interpolation operation, the space sequence is rearranged to reflect the structure continuity, the node information of the compaction and loose areas is combined to generate a complete density matrix, and a three-dimensional structure mapping is formed after merging and arranging, and a slag dense structure parameter field is obtained; Compaction feature recognition submodule: based on the flow field response data in the cabin, the multi-dimensional gradient difference recognition algorithm is used to identify the flow velocity gradient and the density change amplitude in the slag flow area. First, a three-dimensional node coordinate matrix is established with a spatial step of 0.05 meters, the flow velocity and density values of each node are read and arranged in node index order, the gradient calculation command is executed to obtain the flow velocity gradient matrix by dividing the flow velocity difference of adjacent nodes by the spatial step, the density change calculation command is executed to obtain the density change matrix by dividing the density difference of adjacent nodes by the node distance, the two matrices are paired and operated according to the node coordinate position, and the threshold judgment command is used to mark the compaction zone node when the gradient difference is greater than 0.8 and the density difference is greater than 120, and mark the loose zone node when the gradient difference is less than 0.2 and the density difference is less than 40. Execute loop traversal on all nodes and generate a boundary identification list, use the linear sequence splicing command to arrange the node index in the longitudinal direction of the cabin, compare the gradient difference of different sections, and determine the compaction and loose boundary points by the numerical comparison command. Record the cabin position, gradient difference value, density difference value and node number information in the form of a two-dimensional table for each partition result, and generate a slag compaction state parameter set; Dense structure generation submodule: based on the slag compaction state parameter set, a cubic spline interpolation algorithm is used to perform interpolation operation on the boundary gradient values. The interpolation step is set to 0.02 meters. The boundary nodes of the compaction zone and the loose zone are arranged in ascending order of spatial index. The interpolation command generates intermediate gradient values between adjacent nodes by using a spline function and writes them into an interpolation matrix. The spatial sequence is rearranged to reflect the continuity of the structure. The matrix rearrangement command recombines the node number, coordinates and gradient values into a new node sequence table. The merge command is executed on the compaction zone and loose zone node information to generate a complete density matrix. The matrix contains four columns of data: node number, coordinate position, interpolated gradient value and density value. The matrix is sorted by x, y and z axis indexes by the matrix merging command. The three-dimensional mapping command is executed on the merged matrix to generate a spatial body data structure. The three-dimensional structure mapping is output in the form of a voxel grid, and a slag dense structure parameter field is generated.

[0026] The contact transfer network identification module includes: Contact feature extraction submodule: based on the slag dense structure parameter field, the contact points between particles and their relative positions are identified. The normal force and tangential force components of each contact point are recorded, and the displacement offset is measured. The force data of adjacent particles is integrated through spatial aggregation. High and low intensity contact areas are selected and indexed to generate a particle contact feature data set. Connection path identification submodule: based on the particle contact feature data set, the contact points are connected step by step according to the continuity of force direction and contact strength to form a particle node chain. The main channel is determined by detecting the consistency of force direction in the chain. The connectivity and hierarchical order between channels are identified to draw a structure path table and obtain a particle transfer path structure set. Load relationship modeling submodule: based on the particle transmission path structure set, using the centrality algorithm, the force value distribution of each path and the mutual support relationship between nodes are counted, the load components are combined according to the path direction, the load concentration area and transmission strength between different channels are analyzed, the force direction relationship of each region is summarized as a whole, a three-dimensional stress topology diagram is established, and a particle load transmission structure set is formed; Contact feature extraction submodule: based on the dense structure parameter field of the slag, a three-dimensional contact point clustering recognition algorithm is used to identify the contact points between particles and their relative positions, the spatial resolution step is set to 0.01 meters, the coordinate index matrix of the input three-dimensional parameter field is established and traversed in the order of x, y and z three axes, the contact detection command is executed to determine whether the contact is through the distance between the centers of adjacent particles and the difference between the radii, the contact judgment threshold is set to 1.2 times the average particle size, the contact coordinates and contact numbers of the particles meeting the conditions are recorded, the normal force calculation command is executed to obtain the normal force component by multiplying the contact direction unit vector and the force scalar, the tangential force calculation command is executed to obtain the tangential force component by cross multiplying the force direction and the contact surface tangent direction vector, the displacement offset is determined by comparing the displacement value obtained by the difference between the positions of the centers of adjacent particles and the initial contact distance, a record matrix is established for each contact point information, including contact number, coordinate position, normal force value, tangential force value and displacement offset, the spatial aggregation command is executed to combine the particle groups with similar force directions in the three-dimensional coordinate into the same aggregation unit, the force data of the aggregation unit is summed and averaged to generate the force aggregation table, the loop filtering command is executed for all aggregation units to divide the high-intensity contact area with a force average greater than 500 Newton and the low-intensity contact area with a force average less than 100 Newton, and a high and low intensity partition data index table is established, and the output is a particle contact feature data set; Connection path identification submodule: based on the particle contact feature data set, a force direction continuous chain generation algorithm is used to connect the contact points according to the continuity of force direction and contact strength, the direction deviation threshold is set to 15 degrees, the direction vector and intensity value of each contact point are extracted, the connection condition is judged by the adjacent node direction vector angle calculation command, the nodes meeting the direction deviation threshold are connected in the form of adjacency matrix to establish a connection chain, the chain construction command is executed to extend the connection to the end node, the particle node chain is generated and each chain is assigned a unique number, the force direction consistency detection command is executed for each chain to calculate the standard deviation of the direction vector of each node in the chain, and the chain with a standard deviation less than 0.05 is marked as the main channel. The connectivity checking command is executed for the main channel and its adjacent channels to determine the hierarchical relationship by comparing the number of node overlaps and the connection order, the hierarchical structure is numbered and arranged to generate a channel hierarchy table, the table drawing command is executed to generate a structure path table, and the output is a particle transmission path structure set; The load relationship modeling submodule: based on the particle transmission path structure set, the weighted centrality modeling algorithm is used to count the force value distribution and the support relationship between nodes of each path, the centrality calculation weight parameter a is set to 0.65, the node force value of each path is extracted and sorted according to the path direction, the load merging command is executed to sum the adjacent node force values in the path direction by vector weighting, the path load vector sequence is generated, the aggregation command is executed on the load concentrated area between paths to aggregate the force values between nodes as the aggregation condition, the node coordinates, path numbers and concentrated force values of the aggregation results are recorded, the transmission strength analysis command is executed to calculate the transmission strength of the region through the ratio of the average force value and the length, the matrixing arrangement command is executed on the force direction relationship of all regions to establish the node force direction matrix, including path number, node coordinates, load vector and transmission ratio, the three-dimensional topological mapping operation is performed on the matrix to generate the node connection network and draw the stress direction diagram, and the particle load transmission structure set is generated.

[0027] The centrality algorithm first establishes the adjacency matrix of the particle contact network, takes the edge weight as the load strength between nodes, then traverses all paths, records the appearance frequency of nodes in each path as the betweenness value, and then normalizes the node degree, betweenness value and adjacent node weight to calculate the centrality index, obtains the influence degree sorting of each node in the network, and draws the load relationship matrix according to the centrality value distribution; The centrality algorithm is according to the formula:

[0028] Among them: The betweenness centrality value of node v in the muck stress network in the shield sealing cabin, The shortest path count function from node to node under cost measurement, The load line integral quantity in the edge path direction, The support strength coefficient between node and , The path orientation coupling weight, The structure tightness factor, The path timing stability coefficient, The community bridging coefficient, The path curvature penalty term, The load fluctuation coefficient, The load line integral weight coefficient, The support strength weight coefficient, The orientation coupling weight coefficient, The structure tightness weight coefficient, The timing stability weight coefficient, This represents the bridging weight coefficient of the community. This represents the curvature penalty weighting coefficient. This represents the load fluctuation weighting coefficient. Represents the numerical stability constant. and This represents any starting and ending node in the force network of the waste soil. and Representing adjacent nodes in the network; Execution process: First, a weighted directed network graph is established based on the structure set of the waste soil transfer path. The contact points of particles inside the hopper are used as nodes, and edge relationships are constructed based on the load transfer direction. The load line integral on the edge is then extracted. Used to represent the bearing capacity of particles along the direction of force transmission, and the support strength coefficient is calculated based on the local particle support redundancy and compressive balance. This is used to measure the interparticle support stability, and then the path orientation coupling weight is calculated. The coupling degree is determined by the cosine of the angle between the principal stress direction and the path direction. At the same time, the local contact density and node connectivity of the network are extracted, and the structural compactness factor is calculated. To reflect the concentration of local force chains, the system tracks load changes over time and calculates the path temporal stability coefficient. This value represents the proportion of a path that is continuously carried over the time series. Then, by dividing the network community structure, the community bridging coefficient for cross-regional connections is calculated. This is used to measure the load propagation effect of the path in different soil accumulation zones, and the curvature penalty term is calculated based on the ratio of the total curvature to the average arc length of the three-dimensional path curve segment. This reflects the attenuation of load transfer efficiency due to path geometric tortuosity, and the load fluctuation coefficient is then calculated as the ratio of the load standard deviation to the mean within the sliding time window. This is used to evaluate load stability, then all parameters are normalized and a path cost function is constructed with each weight. To represent the importance of different physical features in path calculation, weight values ​​are determined by minimizing prediction error through elastic network regression with non-negative constraints combined with cross-validation, and then normalized to the range of 0 to 1. Path cost is then used as the weight calculation node. To node Number of shortest paths and the number of paths that pass through node v. Finally, the sum of the ratios yields the node betweenness centrality value. It is used to quantify the criticality of node V in the force chain transmission network of the slag and map it into a three-dimensional force topology diagram to realize the dynamic visualization of the stress concentration area, channel load intensity and stress network structure of the slag in the compartment.

[0029] The slag structure stability calculation module comprises: The mechanical balance calculation submodule: based on the particle load transfer structure set, the shear force and support force values of each path are extracted and calculated by using a numerical iterative balance algorithm. First, the path number, node coordinates, shear force value and support force value are read from the input structure set to establish a path index matrix. The shear force and support force of adjacent nodes are respectively subtracted and stored in a temporary vector table according to the node order. The force ratio parameter is obtained by dividing the shear force value by the support force value. All force ratio data are combined into a two-dimensional matrix according to the path number and node order. The abnormal screening threshold interval is set to 0.75 to 1.25. The node sequence with a ratio exceeding the interval is screened out by executing the condition judgment command. The low stability zone number list is generated and mapped with the path number as the primary key. The data arrangement command is used to write the low stability zone number, node coordinates, shear force value, support force value and ratio information into the stress relationship table according to the path number. The number sorting and index merging operation is performed on the data in the table to generate the path force ratio data set. The stability distribution evaluation submodule: based on the path force ratio data set, the range and direction information of the low stability zone are summarized and calculated by using a stability coefficient merging algorithm. The low stability zone number table is read and the node index is established according to the spatial coordinates. The range extraction command is executed for each low stability zone node to obtain the adjacent node coordinate difference and record the direction vector. The stability coefficient calculation command obtains the section stability coefficient value by multiplying the mean value of the adjacent node force ratio and the force direction angle deviation. The merging step is set to 0.05 meters. The stability coefficients of the continuous nodes in the same section are averaged and written into the stability coefficient matrix. The extreme value identification command extracts the maximum and minimum value nodes in each section stability coefficient sequence. The extreme value node number, coordinates and stability coefficient value are recorded into the extreme value table. The difference processing command takes the absolute value of the stability coefficient difference between adjacent extreme value sections to form a difference matrix. The difference matrix and node coordinates are indexed and matched to construct a stability matrix. The matrix includes node number, stability coefficient, direction vector and difference parameter. The three-dimensional drawing command outputs the overall stress distribution structure in the cabin body to generate the slag structure stability judgment result. The mechanical balance calculation submodule: based on the particle load transfer structure set, the shear force and support force values of each path are extracted and calculated by using a numerical iterative balance algorithm. First, the path number, node coordinates, shear force value and support force value are read from the input structure set to establish a path index matrix. The shear force and support force of adjacent nodes are respectively subtracted and stored in a temporary vector table according to the node order. The force ratio parameter is obtained by dividing the shear force value by the support force value. All force ratio data are combined into a two-dimensional matrix according to the path number and node order. The abnormal screening threshold interval is set to 0.75 to 1.25. The node sequence with a ratio exceeding the interval is screened out by executing the condition judgment command. The low stability zone number list is generated and mapped with the path number as the primary key. The data arrangement command is used to write the low stability zone number, node coordinates, shear force value, support force value and ratio information into the stress relationship table according to the path number. The number sorting and index merging operation is performed on the data in the table to generate the path force ratio data set. The stability distribution evaluation submodule: based on the path force ratio data set, the range and direction information of the low stability zone are summarized and calculated by using a stability coefficient merging algorithm. The low stability zone number table is read and the node index is established according to the spatial coordinates. The range extraction command is executed for each low stability zone node to obtain the adjacent node coordinate difference and record the direction vector. The stability coefficient calculation command obtains the section stability coefficient value by multiplying the mean value of the adjacent node force ratio and the force direction angle deviation. The merging step is set to 0.05 meters. The stability coefficients of the continuous nodes in the same section are averaged and written into the stability coefficient matrix. The extreme value identification command extracts the maximum and minimum value nodes in each section stability coefficient sequence. The extreme value node number, coordinates and stability coefficient value are recorded into the extreme value table. The difference processing command takes the absolute value of the stability coefficient difference between adjacent extreme value sections to form a difference matrix. The difference matrix and node coordinates are indexed and matched to construct a stability matrix. The matrix includes node number, stability coefficient, direction vector and difference parameter. The three-dimensional drawing command outputs the overall stress distribution structure in the cabin body to generate the slag structure stability judgment result.

[0030] Please refer to Figure 2 A method for visualizing the state of muck in a shield sealing cabin, the method for visualizing the state of muck in a shield sealing cabin is executed based on the above-mentioned system for visualizing the state of muck in a shield sealing cabin, and comprises the following steps: S1: Based on the muck particle group in the shield sealing cabin, record the particle velocity, displacement and force, calculate the force difference by comparing the time sequence change, extract the shear strength by using the SIMPLE algorithm balance iteration, and generate the particle force distribution parameter set; S2: Based on the particle force distribution parameter set, establish the cabin node system, calculate the flow velocity difference and pressure gradient, correct the flow velocity and pressure by using the SIMPLE algorithm, judge the turbulent and slow flow partition and integrate the shear flux, and obtain the cabin flow field response data; S3: Based on the cabin flow field response data, calculate the flow velocity gradient and density change rate, compare the section standard deviation to determine the compaction state, interpolate the boundary gradient to form a continuous sequence, and obtain the muck dense structure parameter field; S4: Based on the muck dense structure parameter field, extract the particle contact force and displacement deviation, calculate the node weight to identify the main load channel by using the centrality algorithm, and establish the particle load transfer structure set; S5: Based on the particle load transfer structure set, calculate the shear and support force ratio to identify the instability area, weight and integrate the low stability area range and direction, extract the extreme value node to form the stability matrix, and establish the muck structure stability judgment result.

[0031] The above is only a preferred embodiment of the present application, and does not limit the form of the present application, any skilled person in the art can use the disclosed technical content to make changes or modifications as equivalent embodiments applied to other fields, but any simple modification, equivalent change and modification made according to the technical essence of the present application to the above embodiments without departing from the technical solution content of the present application, still belongs to the protection scope of the technical solution of the present application.

Claims

1. A system for visualizing the state of muck in a shield seal capsule, characterized in that, The system comprises: A granular dynamic response analysis module: based on the granular group in the shield sealing cabin, the particle velocity, displacement and stress are calculated in time sequence, the abnormal area of stress is identified by comparing the energy changes, the SIMPLE algorithm is used to match the stress direction difference between the abnormal area particles and extract the shear strength distribution, and the particle stress distribution parameter set is formed; A cabin flow coupling calculation module: based on the particle stress distribution parameter set, a cabin node system is established, the flow velocity difference and pressure gradient are calculated, the momentum exchange value is merged to judge the turbulent and slow flow partition, the local shear and flux are matrix combined to obtain the cabin flow field response data; A muck density evolution judgment module: based on the cabin flow field response data, the flow velocity gradient and density change rate are analyzed, the section standard deviation is calculated to judge the compaction state, the boundary gradient is interpolated and the spatial sequence is reorganized to obtain the muck dense structure parameter field; A contact transfer network identification module: based on the muck dense structure parameter field, the particle contact force and displacement deviation are extracted, the strength matrix is constructed and the node connectivity is analyzed, the centrality algorithm is used to identify the main load channel and sort the transfer path, and the particle load transfer structure set is established; A muck structure stability calculation module: based on the particle load transfer structure set, the path shear and support force ratio are calculated to judge the instability area, the low stability area range and direction are weighted and superimposed to extract the extreme value point to form the stability matrix, and the muck structure stability judgment result is established.

2. The slurry state visualization system within a shield tunnel according to claim 1, characterized in that, The particle stress distribution parameter set includes particle velocity distribution parameters, displacement increment parameters and shear stress difference parameters, the cabin flow field response data includes flow velocity vector parameters, density gradient parameters and momentum flux parameters, the muck dense structure parameter field includes muck density distribution parameters, compaction gradient parameters and structure continuity parameters, the particle load transfer structure set includes node connectivity parameters, main force chain path parameters and stress center parameters, and the muck structure stability judgment result includes stability coefficient parameters, instability area range parameters and direction strength parameters.

3. The slurry state visualization system within a shield tunnel according to claim 1, wherein, The granular dynamic response analysis module comprises: A granular motion measurement submodule: based on the granular group in the shield sealing cabin, the velocity, displacement and stress information of the particles are continuously recorded, the position changes of each particle in the time sequence are mapped as motion trajectories, the particle acceleration is calculated by comparing the velocity changes at adjacent time points, and the stress direction and size are recorded synchronously, and after being summarized, a three-dimensional motion data table is formed to generate a particle motion basic data set; An energy difference judgment submodule: based on the particle motion basic data set, the kinetic energy change value of each particle in the continuous period is calculated, the regions with sudden energy increase and sudden decrease are marked as abnormal areas, the stress difference and direction difference between the particles in the region are extracted, the energy change rate and stress deviation amplitude are compared, the particle group with shear effect is identified, and a stress distribution table is arranged to obtain a particle energy difference parameter set; The stress distribution construction submodule: based on the particle energy difference parameter set, using the SIMPLE algorithm, superimposing the shear strength, normal force and direction vector in the spatial coordinates, calculating the stress transmission ratio between each region, drawing the force distribution curve, and spatially splicing the gradient difference between adjacent regions to form a complete stress network diagram, and generating a particle stress distribution parameter set.

4. The slurry state visualization system within a shield tunnel according to claim 3, characterized in that, The SIMPLE algorithm first sets the velocity component and pressure value of each node in the initial step, discretizes and calculates the momentum term in the flow control equation, then uses the predicted velocity field to calculate the momentum residual, and corrects the pressure gradient to form a corrected pressure field, and then recalculates the node velocity component with the updated pressure value to obtain the corrected velocity field. The difference between the new and old velocity and pressure is judged for convergence, and when the difference is less than the set threshold, the iteration is ended, and the matching result of velocity and pressure is output, completing the calculation of the particle stress balance in the cabin.

5. The slurry condition visualization system within a shield tunnel according to claim 1, wherein, The cabin flow coupling calculation module comprises: The node parameter construction submodule: based on the particle stress distribution parameter set, dividing the node region in the sealed cabin space, assigning flow velocity, pressure and density parameters to each node, comparing the velocity difference of adjacent nodes and calculating the pressure gradient, arranging the node momentum information according to the spatial sequence, integrating the numerical relationship of nodes in each layer, and forming a unified cabin node momentum data set; The flow field difference calculation submodule: based on the cabin node momentum data set, comparing the momentum exchange amplitude of each node to identify high-energy flow areas and slow flow areas, fusing the shear strength and flux values of different regions, mapping the results to the spatial coordinates, extracting the flow direction change trend, forming a continuous flow distribution map, and generating cabin flow field response data.

6. The slurry state visualization system within a shield tunnel according to claim 1, wherein, The slag densification evolution determination module comprises: The compaction feature recognition submodule: based on the cabin flow field response data, the flow velocity gradient and density change amplitude in the slag flow region are counted, the gradient difference of different sections is compared, the demarcation point of the compaction zone and the loose zone is located, the partition result is arranged according to the cabin position, a slag compaction distribution table is formed, and a slag compaction state parameter set is obtained. The dense structure generation submodule: based on the slag compaction state parameter set, the boundary gradient numerical value is subjected to interpolation operation, the spatial sequence is rearranged to reflect the continuity of the structure, the node information of the compaction and loose zones is combined to generate a complete density matrix, and a three-dimensional structure mapping is formed after merging and arranging to obtain a slag dense structure parameter field.

7. The slurry condition visualization system within a shield tunnel according to claim 1, wherein, The contact transmission network identification module comprises: The contact feature extraction submodule: based on the slag dense structure parameter field, the contact points and their relative positions between particles are identified, the normal force and tangential force components of each contact point are recorded, and the displacement offset is measured. The stress data of adjacent particles are integrated through spatial aggregation, high and low intensity contact areas are selected, corresponding data indexes are established, and a particle contact feature data set is generated. The contact feature extraction submodule: based on the slag dense structure parameter field, the contact points and their relative positions between particles are identified, the normal force and tangential force components of each contact point are recorded, and the displacement offset is measured. The stress data of adjacent particles are integrated through spatial aggregation, high and low intensity contact areas are selected, corresponding data indexes are established, and a particle contact feature data set is generated. The connection path identification sub-module: based on the particle contact feature data set, the contact points are gradually connected according to the continuity of the stress direction and the contact strength, the particle node chain is formed, the main channel is determined by detecting the force direction consistency in the chain, the connectivity and hierarchical order between channels are identified, the structure path table is drawn, and the particle transfer path structure set is obtained; The load relationship modeling sub-module: based on the particle transfer path structure set, the centrality algorithm is used to calculate the force value distribution of each path and the mutual support relationship between nodes, combine the load components according to the path direction, analyze the load concentration area and transfer strength between different channels, and overall induce the force direction relationship of each region to establish a three-dimensional stress topology diagram and form a particle load transfer structure set.

8. The slurry state visualization system within a shield tunnel according to claim 7, characterized in that, The centrality algorithm first establishes the adjacency matrix of the particle contact network, uses edge weight to represent the load strength between nodes, then traverses all paths, records the frequency of node appearance in each path as the betweenness value, and then normalizes the node degree, betweenness value and adjacent node weight to calculate the centrality index and obtain the influence degree sorting of each node in the network. The load relationship matrix is drawn according to the centrality value distribution.

9. The slurry condition visualization system within a shield tunnel according to claim 1, wherein, The slag structure stability calculation module includes: The mechanical equilibrium calculation sub-module: based on the particle load transfer structure set, the shear force and support force values of each path are extracted, the ratio of the two is calculated and abnormal sections are screened, low stability regions are numbered and sorted into a stress relationship table to provide basic data for subsequent stability judgment, and a path stress ratio data set is obtained; The stable distribution evaluation sub-module: based on the path stress ratio data set, the range and direction information of the low stability region are summarized, the stability coefficients of each section are calculated and the extreme value node is determined, the extreme value region difference is processed to construct a stability matrix, the overall stress distribution structure in the cabin is output, and the slag structure stability judgment result is established.

10. A method of visualizing the state of the spoil in a shielded capsule, characterized by, The shield sealing cabin slag state visualization system according to any one of claims 1-9 is executed, including the following steps: S1: Based on the particle group in the shield sealing cabin, record the particle velocity, displacement and stress, calculate the stress difference by comparing the time sequence changes, extract the shear strength by using the SIMPLE algorithm balance iteration, and generate the particle stress distribution parameter set; S2: Based on the particle stress distribution parameter set, establish the node system in the cabin, calculate the flow velocity difference and pressure gradient, correct the flow velocity and pressure by using the SIMPLE algorithm, judge the turbulent and slow flow partition and integrate the shear flux to obtain the cabin flow field response data; S3: Based on the cabin flow field response data, calculate the flow velocity gradient and density change rate, compare the section standard deviation to determine the compaction state, interpolate the boundary gradient to form a continuous sequence, and obtain the slag compaction structure parameter field; S4: Based on the slag compaction structure parameter field, extract the particle contact force and displacement deviation, identify the main load channel by using the centrality algorithm to calculate the node weight, and establish the particle load transfer structure set; S5: Based on the particle load transfer structure set, calculate the shear and support force ratio to identify the instability region, weight and integrate the low stability region range and direction, extract the extreme value node to form a stability matrix, and establish the slag structure stability judgment result.

Images