Continuous-discontinuous numerical model coupling method

Abstract

The application discloses a continuous-discontinuous numerical model connecting method, comprising the following steps: step S1: importing a built continuous numerical model into a FLAC3D6.0 platform, finding out a model range by traversing nodes, and setting left, right, front, back and bottom boundaries of the model as viscous boundaries to absorb energy; step S2: setting a corresponding group of the discontinuous numerical model, and enlarging a group boundary to ensure that an enlarged range exceeds an actual area by at least 2 average particle sizes; step S3: generating a circular sphere with an average particle size r in the model area, uniformly balancing the particle system by using a particle servo method, and making the particle system tend to be balanced; and step S4: calculating the displacement, velocity and acceleration of the overlapping circular sphere by using an element interpolation method for the overlapping circular sphere overlapping with a solid element in the model, and transmitting the load to the discrete model range through the force transmission between the overlapping circular sphere particles. The method can effectively realize smooth transmission of stress waves in the continuous numerical model to the discontinuous model area.

Description

Technical Field

[0001] The present application relates to the field of geotechnical engineering technology, and in particular to a continuous-discontinuous numerical model connection method. Background Art

[0002] In geotechnical engineering, numerical simulation methods are primarily used to study the propagation of seismic waves and the damage process under dynamic influence. Continuous numerical simulation methods, such as the finite element method and the finite difference method, have clear mechanical significance and are the most widely used. Their wave fields are continuous and effective, but they are less effective in simulating geotechnical failure. Discontinuous numerical simulation methods, such as the particle discrete element method, can simulate damage caused by large deformations, but their displacement and stress fields are discontinuous.

[0003] Since continuous numerical simulation methods and discontinuous numerical simulation methods have their own advantages and disadvantages, treating the non-large deformation area as a continuous numerical model and the large deformation failure area as a discontinuous model such as particle discrete element has become an important means of studying dynamic failure such as rock and soil earthquake.

[0004] When using the continuous-discontinuous numerical coupling method to carry out numerical simulation, the most common method is to set up walls between the continuous part and the discontinuous part, let the discrete balls act on the walls, and the vertices that constitute the walls match the unit nodes of the continuous part and move in unison, thereby realizing the force-displacement transmission of the continuous-discontinuous numerical model.

[0005] However, the wall acting as an intermediary between the continuous and discontinuous numerical models does not satisfy Newton's second law. Therefore, stress waves at the interface between the continuous and discontinuous numerical models cannot reflect normal stress wave transmission conditions. As a result, the results obtained by using the existing model are inconsistent with the earthquake dynamic load and seismic mechanics mechanism. Summary of the Invention

[0006] In response to the above technical problems, the present application provides a continuous-discontinuous numerical model connection method to improve the simulation effect of stress waves under earthquake action and enhance the authenticity and reference value of the simulation.

[0007] The present application provides a method for connecting a continuous-discontinuous numerical model, comprising the following steps:

[0008] Step S1: Import the built continuous numerical model into the FLAC3D6.0 platform, traverse the nodes to find the model range, and set the left, right, front, back, and bottom boundaries of the model as viscous boundaries to absorb energy;

[0009] Step S2: setting the group corresponding to the discontinuous numerical model, and enlarging the group boundary to ensure that the enlarged range exceeds the actual area by at least 2 average particle sizes;

[0010] Step S3: generating spheres of a specific particle size in the model region, and homogenizing the particle system using a particle servo method to bring the particle system into equilibrium;

[0011] Step S3 includes the following steps:

[0012] Step S31: The boundary geometry of the discontinuous area defined in step S2 is defined as a geometry group. The name can be specified by you. Here, it is called "exterior". At the same time, it is converted into a wall and the group name is also "exterior".

[0013] Step S32: Using the random ball placement function in the FLAC3D platform, a series of balls (radius r) are generated under the control of the geometry group "exterior" generated in step S2. The r value is between the minimum radius rmin and the maximum radius rmax of the ball, showing a uniform random distribution. The average value of all particle radii is approximately equal to the average value of the minimum and maximum radii, denoted as rave. The ratio of the maximum particle radius to the minimum particle radius is preferably between 1.0 and 3.0.

[0014] Step S33: Set the contact between all spheres in the model area to a linear contact model, set the normal stiffness parameter kn = 1e8 N / m, and the tangential stiffness parameter ks = 1e8 N / m, and run the particle system with a time step of 1.0 for a certain number of iterations, such as 20,000 steps, until the particle system is basically balanced;

[0015] Step S34: traverse the vertices of the obtained geometry group "exterior" and calculate the minimum cuboid that accommodates the geometry range, which is defined by xmin0, xmax0, ymin0, ymax0, zmin0, and zmax0. Within the cuboid, set the measurement circle radius to 5 times the average particle radius to generate a series of regular measurement circles;

[0016] Step S35: Control the stress magnitude by measuring the circle. By reducing and enlarging the radius of all particles, the average stress value of the measuring circle is set to 100 kPa so that the particle system tends to equilibrium. The equilibrium model is the discontinuous particle system within the control domain of the discontinuous model.

[0017] Step S4: For overlapping spheres that overlap with solid elements in the model, the displacement, velocity, and acceleration of the overlapping spheres are calculated using the unit interpolation method, and the load is transferred to the discrete model range through the force transmission between the overlapping sphere particles.

[0018] Step S5: Calibrate the micromechanical parameters of the overlap area, transfer the states in the discrete model and the continuous numerical model to each other according to Newton's second law, gradually calculate the seismic wave propagation process, record the displacement, velocity, and acceleration peak of the nodes and discrete spheres, and finally obtain the response spectrum data.

[0019] Step S5 includes the following steps:

[0020] Step S51: When calibrating the microscopic parameters of the particle system, the following conditions should be met: the macroscopic properties of the particle parameters in the model overlap region are consistent or nearly consistent with the macroscopic properties of the continuous model, so as to ensure the stability of the stress wave in the model overlap region;

[0021] Step S52: Set the calculation time step of the particle system to a small value dt, such as 1.0e-7 (small enough to ensure the stability of the stress wave). When the velocity of each node of the model at any time t is known, the velocity of the particle position in the overlapping area is calculated by interpolation, and the obtained velocity is assigned to the ball i at the corresponding position, so that the fluctuation generated by the stress wave is transmitted to the particle system in the discontinuous area.

[0022] Step S53: After the force exerted on the particle system at time t+dt is transmitted to the ball i, the force is transmitted as a concentrated force to the continuous numerical model for calculation. Steps S52 to S53 are repeated until the stress wave is continuously and smoothly transmitted on the boundary, thereby achieving stress wave transmission;

[0023] Step S54: Compare the displacement, velocity, and acceleration peaks of the node and the discrete sphere at each time step, and record the maximum value to obtain response spectrum data.

[0024] Preferably, step S1 includes the following steps:

[0025] Step S11: Treat the discontinuous domain as a continuous model and import the established continuous model into the FLAC3D 6.0 or higher version platform. The model consists of node coordinates and units, where the node coordinates are (X, Y, Z) and the units are (each unit contains 8 nodes and unit group names);

[0026] Step S12: In principle, the model is projected as a regular rectangle on the horizontal plane, so the minimum value of the model range in the x direction is set to x min (initial value 100000), the maximum value is x max (Initial value is -10000); Model range: minimum value in y direction is y min (initial value 100000), the maximum value is ymax (initial value -100000); the model has a minimum value z in the z direction min (initial value 100000), traverse all nodes of the model and update x min 、x max 、y min 、y max 、z min .

[0027] Preferably, the updating step is: if the coordinate of node x is less than x min , then replace x with xmin ; If the node coordinate x value is greater than x max , then replace x with x max ; If the node coordinate y value is less than y min , then replace y with y min ; If the node coordinate y is greater than y max , then replace y with y max ; If the node coordinate z is less than z min , then replace z with z min .

[0028] Preferably, step S12 further includes: setting a boundary tolerance error=0.05, setting the x coordinate of the left boundary of the model to [xmin-error, xmin+error], the x coordinate of the right boundary to [xmax-error, xmax+error], the y coordinate of the front boundary to [ymin-error, ymin+error], the y coordinate to [ymax-error, ymax+error], and the z coordinate of the bottom boundary to [zmin-error, zmin+error]; and setting independent dampers in the normal direction and horizontal direction of the model to absorb incident waves generated by internal vibration of the model;

[0029] Normal viscous force t provided by the damper n and tangential viscosity t s As shown in formula (1):

[0030] t n =-ρC p v n

[0031] t s =-ρC s v s Formula (1)

[0032] Where: v n ,v s are the normal and tangential velocity components on the model boundary, ρ is the medium density, C p ,C s are the p-wave and s-wave velocities of the model, respectively.

[0033] Preferably, the group boundary enlargement operation in step S2 includes the following steps:

[0034] Step S21: group the continuous models that are at risk of damage, output their outer boundary geometry as a discontinuous model boundary, which is composed of a series of triangles, and estimate the value of the average particle radius rave0 of the continuous model;

[0035] Step S22: Normal judgment is performed on the non-continuous model boundary geometry, and the judgment condition and corresponding operation are as follows:

[0036] 1) If the normal vector of the model boundary is positive, the position of the lower half of the boundary is lowered by at least 2 times the average particle radius;

[0037] 2) If the upper boundary of the model is a slope surface, the upper boundary cannot be changed;

[0038] 3) If the upper boundary of the model is not the outer surface of the model but an intermediate layer, the upper boundary can also be translated upward by at least 2 times the average particle radius.

[0039] Preferably, step S4 comprises the following steps:

[0040] Step S41: All entity elements zone and particles ball are traversed, and the model is divided into a continuous model area, a discrete model area, and a model overlap area according to the geometric positions of the entity elements and the particles, the particles in the rock bridge area of the model are boundary control particles, and the area is an overlapping area of the entity elements and the particles;

[0041] Step S42: Since one entity grid element can overlap with one or more particles, the transmission of displacement of the overlapping area is calculated as follows:

[0042]

[0043] In the formula, α is equal to 1 and β is equal to 0 in the discrete area; in the continuous area, α is equal to 0 and β is equal to 1; for the ball-zone coupling method, α and β linearly change from 0 to 1 in the overlapping area, α j and β i are respectively the coefficients of the discrete particles and the continuous elements, and the value range is 0.01 to 0.99; m j is the mass of the discrete particle, d j is the displacement of the discrete particle, F j tot is the external force acting on the discrete particle, λ j is the Lagrange multiplier of the discrete particle, n j is the number of discrete particles; m i is the mass of the continuous element, u i is the displacement of the continuous element, F i tot is the external force acting on the continuous element, λ k is the Lagrange multiplier at the node of the continuous element, n i is the number of continuous elements; K is the motion matrix, k jk is the element motion matrix defined by the classical interpolation function, u kis the displacement function of the eight nodes (k) of the continuous element i surrounding the discrete particle j.

[0044] The beneficial effects of this application include:

[0045] 1) The continuous-discontinuous numerical model connection method provided in this application can effectively realize the smooth transfer of stress waves in the continuous numerical model to the discontinuous model area, making the transfer of seismic dynamic loads more realistic and the results obtained more consistent with mechanical mechanisms.

[0046] 2) The continuous-discontinuous numerical model connection method provided in this application repeatedly uses the boundary control particles at the connection between the continuous numerical model and the discrete numerical model in different time steps, so that the movement and force of different regions are transmitted to each other, and the stress waves between the continuous numerical model and the discontinuous numerical model are propagated according to their propagation characteristics, thereby obtaining more realistic and reliable simulation results and improving the authenticity and effectiveness of the simulation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 A schematic flow chart of the continuous-discontinuous numerical model connection method provided in this application;

[0048] Figure 2 It is a continuous-discontinuous coupled slope numerical model in the embodiment of the present application;

[0049] FIG3 is a diagram of a single unit node and its degenerate form in an embodiment of the present application;

[0050] Figure 4 is the boundary of the continuous-discontinuous numerical model in the embodiment of the present application;

[0051] Figure 5 is the discontinuous numerical model boundary expanded in the embodiment of the present application;

[0052] FIG6 is a measuring ball of a discontinuous particle state monitoring arrangement according to an embodiment of the present application;

[0053] Figure 7 is a monitoring curve of the equilibrium state of the discontinuous numerical model in the embodiment of the present application;

[0054] FIG8 is a diagram illustrating the interpolation processing of boundary overlap sphere displacement in an embodiment of the present application; a) is a schematic diagram of the overlap of each zone; b) is a schematic diagram of the position of any discrete particle j within a unit;

[0055] Figure 9 is the calculation result of the coupled boundary seismic response in the embodiment of the present application; DETAILED DESCRIPTION

[0056] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions of the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings herein can be arranged and designed in various different configurations.

[0057] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention as claimed, but rather merely represents selected embodiments of the present invention. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of the present invention without creative effort are also within the scope of protection of the present invention.

[0058] The controller used in this embodiment is an existing structure, and the control circuit can be implemented through simple programming by technicians in this field. It is common knowledge in this field and is only used without modification. Therefore, the control method and circuit connection will not be described in detail.

[0059] The technical means that are not described in detail in this application and are not used to solve the technical problems of this application are all set according to the common knowledge in this field, and can be implemented in a variety of common knowledge settings.

[0060] See also Figure 1 The continuous-discontinuous numerical model connection method provided in this application includes the following steps:

[0061] Step S1: Import the built continuous numerical model into the FLAC3D6.0 platform, traverse the nodes to find the model range, and set the left, right, front, back, and bottom boundaries of the model as viscous boundaries to absorb energy;

[0062] Step S2: setting the group corresponding to the discontinuous numerical model, and enlarging the group boundary to ensure that the enlarged range exceeds the actual area by at least 2 average particle sizes;

[0063] Step S3: generating spheres of a specific particle size in the model region, and homogenizing the particle system using a particle servo method to bring the particle system into equilibrium;

[0064] Step S3 includes the following steps:

[0065] Step S31: The boundary geometry of the discontinuous area defined in step S2 is defined as a geometry group. The name can be specified by you. Here, it is called "exterior". At the same time, it is converted into a wall and the group name is also "exterior".

[0066] Step S32: Using the random ball placement function in the FLAC3D platform, a series of balls (radius r) are generated under the control of the geometry group "exterior" generated in step S2. The r value is between the minimum radius rmin and the maximum radius rmax of the ball, showing a uniform random distribution. The average value of all particle radii is approximately equal to the average value of the minimum and maximum radii, denoted as rave. The ratio of the maximum particle radius to the minimum particle radius is preferably between 1.0 and 3.0.

[0067] Step S33: Set the contact between all spheres in the model area to a linear contact model, set the normal stiffness parameter kn = 1e8 N / m, and the tangential stiffness parameter ks = 1e8 N / m, and run the particle system with a time step of 1.0 for a certain number of iterations, such as 20,000 steps, until the particle system is basically balanced;

[0068] Step S34: traverse the vertices of the obtained geometry group "exterior" and calculate the minimum cuboid that accommodates the geometry range, which is defined by xmin0, xmax0, ymin0, ymax0, zmin0, and zmax0. Within the cuboid, set the measurement circle radius to 5 times the average particle radius to generate a series of regular measurement circles;

[0069] Step S35: Control the stress magnitude by measuring the circle. By reducing and enlarging the radius of all particles, the average stress value of the measuring circle is set to 100 kPa so that the particle system tends to equilibrium. The equilibrium model is the discontinuous particle system within the control domain of the discontinuous model.

[0070] Step S4: For overlapping spheres that overlap with solid elements in the model, the displacement, velocity, and acceleration of the overlapping spheres are calculated using the unit interpolation method, and the load is transferred to the discrete model range through the force transmission between the overlapping sphere particles.

[0071] Step S5: Calibrate the micromechanical parameters of the overlap area, transfer the states in the discrete model and the continuous numerical model to each other according to Newton's second law, gradually calculate the seismic wave propagation process, record the displacement, velocity, and acceleration peak of the nodes and discrete spheres, and finally obtain the response spectrum data.

[0072] Step S5 includes the following steps:

[0073] Step S51: When calibrating the microscopic parameters of the particle system, the following conditions should be met: the macroscopic properties of the particle parameters in the model overlap region are consistent or nearly consistent with the macroscopic properties of the continuous model, so as to ensure the stability of the stress wave in the model overlap region;

[0074] Step S52: When the velocity of each node of the model at any time t is known, the velocity of the particle in the overlapping region is calculated by interpolation, and the obtained velocity is assigned to the ball i at the corresponding position, thereby transmitting the fluctuation generated by the stress wave to the particle system in the discontinuous region;

[0075] Step S53: After the force exerted on the particle system at time t+dt is transmitted to the ball i, the force is transmitted as a concentrated force to the continuous numerical model for calculation. Steps S52 to S53 are repeated until the stress wave is continuously and smoothly transmitted on the boundary, thereby achieving stress wave transmission;

[0076] Step S54: Compare the displacement, velocity, and acceleration peaks of the node and the discrete sphere at each time step, and record the maximum value to obtain response spectrum data.

[0077] The method provided in this application can generate spheres of particle size r in the model area, and perform the above operations on the obtained spheres, so as to transfer the load to the discrete model range in the form of particles through the force transmission between overlapping spherical particles, thereby effectively improving the force conduction simulation accuracy of the model.

[0078] Preferably, step S1 includes the following steps:

[0079] Step S11: Treat the discontinuous domain as a continuous model and import the established continuous model into the FLAC3D 6.0 or higher version platform. The model consists of node coordinates and units, where the node coordinates are (X, Y, Z) and the units are (each unit contains 8 nodes and unit group names);

[0080] Step S12: In principle, the model is projected as a regular rectangle on the horizontal plane, so the minimum value of the model range in the x direction is set to x min (initial value 100000), the maximum value is x max (Initial value is -10000); Model range: minimum value in y direction is y min (initial value 100000), the maximum value is ymax (initial value -100000); the model has a minimum value z in the z direction min (initial value 100000), traverse all nodes of the model and update x min 、x max 、y min ,ymax,z min .

[0081] Preferably, the updating step is: if the coordinate of node x is less than x min , then replace x with x min ; If the node coordinate x value is greater than x max , then replace x with x max ; If the node coordinate y value is less than ymin , then replace y with y min ; If the node coordinate y is greater than y max , then replace y with y max ; If the node coordinate z is less than z min , then replace z with z min .

[0082] Preferably, step S12 further includes: setting a boundary tolerance error=0.05, setting the x coordinate of the left boundary of the model to [xmin-error, xmin+error], the x coordinate of the right boundary to [xmax-error, xmax+error], the y coordinate of the front boundary to [ymin-error, ymin+error], the y coordinate to [ymax-error, ymax+error], and the z coordinate of the bottom boundary to [zmin-error, zmin+error]; and setting independent dampers in the normal direction and horizontal direction of the model to absorb incident waves generated by internal vibration of the model;

[0083] Normal viscous force t provided by the damper n and tangential viscosity t s As shown in formula (1):

[0084] t n =-ρC p v n

[0085] t s =-ρC s v s Formula (1)

[0086] Where: v n ,v s are the normal and tangential velocity components on the model boundary, ρ is the medium density, C p ,C s are the p-wave and s-wave velocities of the model, respectively.

[0087] Preferably, the group boundary enlargement operation in step S2 includes the following steps:

[0088] Step S21: group the continuous models that are at risk of damage, output their outer boundary geometry as a discontinuous model boundary, which is composed of a series of triangles, and estimate the value of the average particle radius rave0 of the continuous model;

[0089] Step S22: Perform normal judgment on the boundary geometry of the discontinuous model, judgment conditions and corresponding operations:

[0090] 1) If the normal vector of the model boundary is positive, move the lower half of the boundary down by at least 2 times the average particle radius;

[0091] 2) If the upper boundary of the model is a slope surface, the upper boundary cannot be changed;

[0092] 3) If the upper boundary of the model is not the outer surface of the model but the middle layer, the upper boundary can also be translated upward by at least 2 times the average particle radius.

[0093] Preferably, step S4 includes the following steps:

[0094] Step S41: Traverse all entity unit zones and particle balls, and divide the model into the following areas according to the geometric positions of the entity units and particles: continuous model area, discrete model area, and model overlap area. The particles in the rock bridge area of ​​the model are boundary control particles, and this area is the overlapping area where the entity units and particles overlap;

[0095] Step S42: Since a solid grid cell may overlap with one or more particles, the displacement transfer in the overlapping area is calculated as follows:

[0096]

[0097] Where: α is equal to 1 and β is equal to 0 in the discrete region; α is equal to 0 and β is equal to 1 in the continuous region; for the ball-zone coupling method, α and β vary linearly from 0 to 1 in the overlapping region, and α j and β i are the coefficients of discrete particles and continuous elements, ranging from 0.01 to 0.99; m j is the mass of discrete particles, d j is the displacement of discrete particles, F j to t is the external force acting on the discrete particles, λ j is the Lagrange multiplier of the discrete particle, n j is the number of discrete particles; m i is the mass of the continuous unit, u i is the displacement of the continuous element, F i to t is the external force acting on the continuous element, λ k is the Lagrange multiplier on the continuous element nodes, n i is the number of continuous units; K is the motion matrix, k jk is the unit motion matrix defined by the classical interpolation function, u k is the displacement function of the eight nodes (k) of the continuous element i surrounding the discrete particle j.

[0098] Example

[0099] See also Figures 2 to 9 The slope of a certain hydropower project is close to 50 million cubic meters in size and nearly a kilometer in horizontal dimension. It is a huge scale. The propagation of stress waves under earthquake is complex, and it is necessary to focus on the dynamic response under earthquake and the deformation and failure mechanism of the accumulation body. However, due to the large size of the model, a discontinuous numerical method (PFC3D) is used to simulate the accumulation body, and a continuous numerical simulation method (FLAC3D6.0) is used to simulate other rock and soil bodies. However, how to deal with the stress wave transmission between the continuous model and the discontinuous model is very important to the calculation results. A transmission boundary processing method at the connection of a continuous-discontinuous numerical model of the present invention is used to implement it, and the steps are as follows:

[0100] (1) Based on the FLAC3D6.0 platform, which can be coupled with PFC3D, it can simultaneously realize the calculation of continuous numerical model and discontinuous numerical model. First, the built continuous numerical model (including the part of the landslide body that needs to be considered as discontinuous) is imported, such as Figure 2 The model is divided into four groups. From the bottom of the model upward, group "1" is slightly weathered rock mass, group "2" is weathered rock mass, group "3" is the sliding zone soil between the accumulation body and the weathered rock mass, and group "4" is the landslide body.

[0101] The model consists of 276,827 nodes and 531,815 units. Each unit can be Figure 3a ) or the degenerate triangular prism element ( Figure 3b )), tetrahedral element ( Figure 3c )).

[0102] Traverse all node coordinates to find the model range: the left interface xmin = -1400; the right interface xmax = 960; the front boundary ymin = -500; the rear boundary ymax = 2450; the bottom boundary zmax = 1600. Set the left, right, front, rear, and bottom boundaries as sticky boundaries. The sticky force is applied using the following formula and updated at each time step:

[0103] t n =-ρC p v n

[0104] t s =-ρC s v s

[0105] (2) For the continuous model landslide group “4” with the risk of destruction, the geometric operation function of the software platform is used to output the outer boundary geometry of group “4” as a discontinuous model boundary, which is composed of a series of triangles, such as Figure 4 shown.

[0106] According to experience, the landslide body is set as a sphere with a radius of 1.5 to 3.0 m, and the estimated value of the average particle radius rave is 2.25 m.

[0107] The normal vector of the discontinuous model boundary geometry is determined. If the normal vector is positive, it indicates the lower half of the boundary. Its position is shifted downward by 2 to 10 times the average particle radius, with a practical value of 20 m. If the upper boundary is a slope surface, the upper boundary cannot be changed, so the upper boundary surface translation distance is 0.0 m.

[0108] (2) Traverse the geometric nodes within the landslide model range and find the minimum x-axis value of -862.0m, the maximum x-axis value of 452.0m, the minimum y-axis value of -262.0m, the maximum y-axis value of 153.0m, the minimum z-axis value of 2132m, and the maximum z-axis value of 3248m. Using the rectangular region constrained along the three coordinate axes (x-axis, y-axis, and z-axis) and the discontinuous model boundary geometry, a total of 125,378 Gaussian-distributed spheres with radii of 1.5m to 3m were generated within the discrete model boundary control area.

[0109] The numerical model consists of the geometry after the discontinuous domain boundary is enlarged, the continuous unit (zone), and the discontinuous ball (ball). Figure 5 shown.

[0110] Generate a measurement sphere with a radius of 10m in the above rectangular area, such as Figure 6a ), delete the measuring balls with porosity greater than 0.5 or whose center is no longer within the range of the discontinuous model, and the remaining measuring balls are used for particle system monitoring, such as Figure 6b ).

[0111] Based on the discrete element calculation platform PFC3D, the effective modulus of the spherical particle system is set to 5e8MPa, the stiffness ratio is 2.0, the particle system is balanced, and the particle servo method is used to make the particle system meet the homogenization and the average stress of all measurement circles is 0.5MPa. The state of the particle system gradually approaches the servo stress. Figure 7 shown.

[0112] (3) Traverse all the entity unit zones and particle balls, and divide the model into three parts according to their geometric positions: continuous model area, discrete model area, and model overlap area, such as Figure 8a The overlap region of the model is the boundary control particle. This region is an overlapping area, and a solid grid unit may overlap with one or more particles.

[0113] Specific operations such as Figure 8b ), without loss of generality, assume that a unit consists of 1, 2, 3, 4, 5, 6, 7, and 8 nodes, and the displacements corresponding to these 8 nodes are The speed corresponding to these 8 nodes The displacement of the sphere (numbered j) located in the unit can be obtained by interpolation using the following formula.

[0114] The transfer calculation of the displacement in the overlapping area is shown as follows:

[0115]

[0116] (5) When calibrating the microscopic parameters of the particle system, it should be noted that the macroscopic properties reflected by the particle parameters in the model overlap area should be consistent or nearly consistent with the macroscopic properties of the continuous model to ensure the stability of the stress wave in the model overlap area.

[0117] When the velocity of the node t at any time is known, the velocity of the particle in the overlapping area is calculated by interpolation, and the velocity is assigned to the ball j at the corresponding position, and then transferred to the particle system in the discontinuous area.

[0118] After the force on the particle system at time t+dt is transmitted to i, this force is transferred as a concentrated force to the continuous numerical model for calculation. This process is repeated, and the stress wave can be transmitted continuously and smoothly on this boundary.

[0119] After assigning the parameters, a typical seismic wave is applied, and the stress wave propagation at the typical moment obtained by the above principle is as follows: Figure 9 As shown, it can be seen that the stress on the boundary of the continuous-discontinuous numerical model can be transferred smoothly, indicating that the data obtained by the model obtained by the method provided in this application can effectively improve the simulation effect of stress waves under earthquake action. The data obtained by fitting the model can be more consistent with the actual transmission of earthquake dynamic loads, and the results obtained are more consistent with the mechanical mechanism.

[0120] At each moment, the displacement, velocity, and acceleration of the continuous grid are compared with the historical maximum displacement, velocity, and acceleration of the node. If the current displacement, velocity, or acceleration is greater than the historical maximum, the current value replaces the historical maximum. The maximum displacement, velocity, and acceleration of each node at the end of the calculation are the displacement, velocity, and acceleration response spectrum values.

[0121] Although the present invention has been described in detail with reference to the aforementioned embodiments, it is still possible for those skilled in the art to modify the technical solutions described in the aforementioned embodiments, or to make equivalent substitutions for some of the technical features therein. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for connecting continuous and discontinuous numerical models, characterized in that: The following steps are involved: Step S1: Import the built continuous numerical model into the FLAC3D6.0 platform, traverse the nodes to find the model range, and set the left, right, front, back, and bottom boundaries of the model as viscous boundaries to absorb energy; Step S2: setting the group corresponding to the discontinuous numerical model, and enlarging the group boundary to ensure that the enlarged range exceeds the actual area by at least 2 average particle sizes; Step S3: generating spheres of a specific particle size in the model region, and homogenizing the particle system using a particle servo method to bring the particle system into equilibrium; Step S3 includes the following steps: Step S31: The boundary geometry of the discontinuous area defined in step S2 is used to form a geometry group, here called "exterior", and it is converted into a wall, and the group name is also "exterior"; Step S32: Using the random ball placement function in the FLAC3D platform, a series of balls are generated under the control of the geometry group "exterior" generated in step S2. The radius r is between the minimum radius rmin and the maximum radius rmax of the sphere, showing a uniform random distribution. The average value of all particle radii is approximately equal to the average value of the minimum and maximum radii, denoted as rave, where the ratio of the maximum particle radius to the minimum particle radius is between 1.0 and 3.0; Step S33: Set the contact between all spheres in the model area to a linear contact model, set the normal stiffness parameter kn = 1e8 N / m, and the tangential stiffness parameter ks = 1e8 N / m, and run the particle system for 20,000 iterations with a time step of 1.0 until the particle system is basically balanced; Step S34: Traverse the vertices of the obtained geometry group "exterior" and calculate the minimum cuboid that contains the geometry range, which is defined by xmin0, xmax0, ymin0, ymax0, zmin0, and zmax0. Within the cuboid, set the measurement circle radius to 5 times the average particle radius to generate a series of regular measurement circles; Step S35: Control the stress magnitude by measuring the circle. By reducing and enlarging the radius of all particles, the average stress value of the measuring circle is set to 100 kPa so that the particle system tends to equilibrium. The equilibrium model is the discontinuous particle system within the control domain of the discontinuous model. Step S4: For the overlapping spheres that overlap with the solid elements in the model, the displacement, velocity, and acceleration of the overlapping spheres are calculated using the unit interpolation method, and the load is transferred to the discrete model range through the force transmission between the overlapping sphere particles; Step S5: Calibrate the microscopic mechanical parameters of the overlap area, transfer the states between the discrete model and the continuous numerical model according to Newton's second law, gradually calculate the stress wave propagation process, record the displacement, velocity, and acceleration peaks of the nodes and discrete spheres, and obtain the response spectrum data; Step S5 includes the following steps: Step S51: When calibrating the microscopic parameters of the particle system, the following conditions should be met: the macroscopic properties of the particle parameters in the model overlap region are consistent or nearly consistent with the macroscopic properties of the continuous model, so as to ensure the stability of the stress wave in the model overlap region; Step S52: Set the calculation time step of the particle system to dt, with a value of dt of 1.0e-7. When the velocity of each node of the model at any time t is known, the velocity of the particle in the overlapping area is calculated by interpolation, and the obtained velocity is assigned to the ball i at the corresponding position, thereby realizing the transmission of the fluctuation generated by the stress wave to the particle system in the discontinuous area; Step S53: After the force exerted on the particle system at time t+dt is transmitted to the ball i, the force is transmitted as a concentrated force to the continuous numerical model for calculation. Steps S52 to S53 are repeated until the stress wave is continuously and smoothly transmitted on the boundary, thereby achieving stress wave transmission; Step S54: Compare the displacement, velocity, and acceleration peaks of the node and the discrete sphere at each time step, and record the maximum value to obtain response spectrum data.

2. The method for connecting continuous and discontinuous numerical models according to claim 1, characterized in that: Step S1 includes the following steps: Step S11: Treat the discontinuous domain as a continuous model and import the established continuous model into the FLAC3D 6.0 or higher version platform. The model consists of node coordinates and units, where the node coordinates are (X, Y, Z). Each unit contains 8 nodes, and the units are distinguished by unit group names. Step S12: In principle, the model is projected as a regular rectangle on the horizontal plane, so the minimum value of the model range in the x direction is set to x min , initial value 100000, maximum value is x max , the initial value is -100000; the model range is: the minimum value in the y direction is y min , initial value 100000, maximum value is ymax, initial value -100000; the model has a minimum value z in the z direction min , initial value 100000, traverse all nodes of the model and update x min 、x max 、y min 、y max 、z min .

3. The method for connecting continuous and discontinuous numerical models according to claim 2, characterized in that: The update steps are: if the coordinate of node x is less than x min , then replace x with x min ; If the node coordinate x value is greater than x max , then replace x with x max ; If the node coordinate y value is less than y min , then replace y with y min ; If the node coordinate y is greater than ymax, replace ymax with y; if the node coordinate z is less than zmin, replace zmin with z.

4. The method for connecting continuous and discontinuous numerical models according to claim 1, wherein: Step S12 further includes: setting a boundary tolerance error = 0.05, setting the x-coordinate of the left boundary of the model to [xmin-error, xmin+error], the x-coordinate of the right boundary to [xmax-error, xmax+error], the y-coordinate of the front boundary to [ymin-error, ymin+error], the y-coordinate to [ymax-error, ymax+error], and the z-coordinate of the bottom boundary to [zmin-error, zmin+error]; and setting independent dampers in the normal direction and horizontal direction of the model to absorb incident waves generated by internal vibration of the model; Normal viscous force t provided by the damper n and tangential viscosity t s As shown in formula (1): Where: v n ,v s are the normal and tangential velocity components on the model boundary, ρ is the medium density, C p ,C s are the p-wave and s-wave velocities of the model, respectively.

5. The method for connecting continuous and discontinuous numerical models according to claim 1, wherein: The group boundary enlargement operation in step S2 includes the following steps: Step S21: group the continuous models that are at risk of damage, output their outer boundary geometry as a discontinuous model boundary, which is composed of a series of triangles, and estimate the value of the average particle radius rave0 of the continuous model; Step S22: Perform normal judgment on the boundary geometry of the discontinuous model, judgment conditions and corresponding operations: 1) If the normal vector of the model boundary is positive, move the lower half of the boundary down by at least 2 times the average particle radius; 2) If the upper boundary of the model is a slope surface, the upper boundary cannot be changed; 3) If the upper boundary of the model is not the outer surface of the model but the middle layer, the upper boundary is translated upward by at least 2 times the average particle radius.

6. The method for connecting continuous and discontinuous numerical models according to claim 1, wherein: Step S4 includes the following steps: Step S41: Traverse all entity unit zones and particle balls, and divide the model into the following areas according to the geometric positions of the entity units and particles: continuous model area, discrete model area, and model overlap area. The particles in the rock bridge area of ​​the model are boundary control particles, and this area is the overlapping area where the entity units and particles overlap; Step S42: Since a solid grid cell may overlap with one or more particles, the displacement transfer in the overlapping area is calculated as follows: Where: α j and β i are the action coefficients of discrete particles and continuous elements, ranging from 0.0 to 1.0; in the discrete region α j Equal to 1, β i Equal to 0; in the continuous region, α j Equal to 0, β i Equal to 1; for the ball-zone coupling region, α j and β i In the overlapping area, it changes linearly from 0.00 to 1.0, and the sum of the two is 1; m j is the mass of discrete particles, d j is the displacement of discrete particles, F j tot is the external force acting on the discrete particles, λ j is the Lagrange multiplier of the discrete particle, n j is the number of discrete particles; m i is the mass of the continuous unit, u i is the displacement of the continuous element, F i tot is the external force acting on the continuous element, λ k is the Lagrange multiplier on the continuous element nodes, n i is the number of continuous units; K is the motion matrix, k jk is the unit motion matrix defined by the classical interpolation function, u k is the displacement function of the eight nodes of the continuous element i surrounding the discrete particle j.

Images