A non-saturated transient fluid-structure interaction calculation method based on FDM-DDA

Abstract

The present invention relates to the field of data processing, and specifically to an unsaturated transient fluid-solid coupling calculation method based on FDM-DDA, comprising: using MATLAB to generate cracks to establish a rock and soil model; identifying rock cracks in the rock and soil model to form a fluid network; calculating the rock crack aperture in the fluid network; performing DDA unsaturated steady-state seepage field calculation, using FDM to perform unsaturated transient seepage field calculation, and obtaining steady-state calculation results; writing the steady-state calculation results into DDA as a fluid calculation module and assigning unsaturated transient seepage field calculation parameters, solving the water head and volumetric water content at each node as transient calculation results; and performing fluid-solid coupling analysis based on the transient calculation results. The present invention adds the function of rock crack identification, uses FDM to write the Richards equation into DDA, and realizes DDA simulation of the unsaturated transient seepage process of rock and soil, which can be used in various working conditions such as rainfall and reservoir water level changes.

Description

Technical Field

[0001] The present invention relates to the field of data processing, and in particular to a non-saturated transient fluid-solid coupling calculation method based on FDM-DDA. Background Art

[0002] When analyzing the safety and stability of slope engineering projects, the slope's condition under heavy rain conditions is often considered. Rainfall weakens the strength of rock and soil, and is therefore often a significant factor in instability and poses a safety hazard.

[0003] Currently, the finite element method (FEM) is widely used in geotechnical modeling for slope engineering. However, the FEM is based on the continuity assumption, which does not conform to the nature of geotechnical materials as discrete media. Therefore, DDA, a numerical method with large deformation and displacement characteristics, is widely used in the numerical simulation of geotechnical materials. Simulations based on DDA have played a positive role in engineering design and construction stability assessment.

[0004] However, it has the following disadvantages:

[0005] 1. DDA numerical simulation of rock and soil mostly focuses on the solid aspect, that is, simple rock and soil, ignoring the impact of seepage on rock and soil;

[0006] 2. When considering water-rock-soil interactions, or fluid-solid coupling simulations, DDA often only considers steady-state seepage. Simulating transient processes such as rainfall and reservoir level fluctuations using steady-state seepage cannot fully capture the impact of these processes on the rock and soil.

[0007] In view of the above problems, a transient fluid-solid coupling calculation method suitable for DDA is urgently needed to realize DDA simulation of non-stable seepage rock and soil masses such as rainfall and reservoir water level fluctuations. Summary of the Invention

[0008] The present invention aims to provide an FDM-DDA-based method for calculating unsaturated transient fluid-solid coupling. This method utilizes DDA for preprocessing and stress field calculation. FDM is used to calculate the unsaturated transient seepage field, and the results are then imported into DDA for fluid-solid coupling analysis. This allows DDA to simulate unstable seepage rock and soil conditions, such as rainfall and reservoir level fluctuations.

[0009] To achieve the above object, the present invention provides an unsaturated transient fluid-solid coupling calculation method based on FDM-DDA, comprising generating cracks using MATLAB to establish a rock and soil model;

[0010] Identifying rock fractures in the rock and soil model to form a fluid network;

[0011] Calculation of rock fracture aperture in fluid networks;

[0012] Perform DDA unsaturated steady-state seepage field calculations and use FDM to calculate unsaturated transient seepage field to obtain steady-state calculation results;

[0013] The steady-state calculation results are written into DDA as a fluid calculation module and assigned with unsaturated transient seepage field calculation parameters. The hydraulic head and volumetric water content at each node are solved as transient calculation results.

[0014] Fluid-structure coupling analysis is performed based on transient calculation results.

[0015] The specific method of using matlab to generate cracks is: using matlab to write a Monte Carlo algorithm to generate the coordinates of the two endpoints of a random crack, and then importing them into the DDA pre-processing dc program to generate cracks.

[0016] The specific steps of identifying rock fractures in the rock and soil model and forming a fluid network are:

[0017] Obtain a line segment formed by two adjacent points, label the endpoints sequentially, and store the labels in the first column of the first two-dimensional data;

[0018] Traversing the line segments, searching for two line segments whose endpoint coordinate differences are within a preset range and do not belong to the same block as the two sides of the short crack;

[0019] Storing the line segment numbers on both sides of the short crack in the second column of the first two-dimensional data;

[0020] identifying long cracks in the line segments based on parallel lines;

[0021] Split the long cracks into short cracks, and store the line segment numbers of the short cracks in the first two columns of the second two-dimensional array;

[0022] The same points in the first two-dimensional array and the second two-dimensional array are merged to form a fluid network.

[0023] The specific steps of calculating the rock mass fracture aperture in the fluid network are:

[0024] Get the endpoint numbers of the line segments on both sides of the crack. The numbers on one side are the first and second points, and the numbers on the other side are the third and fourth points.

[0025] Calculate a first distance to a third point based on the first point and the second point;

[0026] Calculate a second distance to a fourth point based on the first point and the second point;

[0027] The equivalent hydraulic opening is calculated based on the first distance and the second distance.

[0028] The specific steps of calculating the DDA unsaturated steady-state seepage field are:

[0029] The cubic law is used to calculate the seepage of a single fracture;

[0030] The Gauss–Seidel iteration was used to obtain the hydraulic head at each node, and then the volumetric water content at each node was calculated;

[0031] The relevant program is written into the DDA program as a steady-state seepage module.

[0032] The specific steps of calculating the unsaturated transient seepage field using FDM are:

[0033] The rock mass fractures are used as the fluid network required for FDM seepage field calculation;

[0034] FDM was used to transform the Richards equation into a volumetric water content matrix that conforms to the DDA format.

[0035] The specific steps of performing fluid-solid coupling analysis based on transient calculation results are:

[0036] The transient calculation results are converted into water load matrix according to the minimum potential energy principle;

[0037] The water load matrix is added to the DDA load matrix to perform DDA stress field calculation. When the block displacement is less than the specified value, the stress field calculation reaches equilibrium.

[0038] The present invention provides an unsaturated transient fluid-solid coupling calculation method based on FDM-DDA, comprising: generating fractures using MATLAB to establish a rock and soil model; identifying rock fractures in the rock and soil model to form a fluid network; calculating the rock fracture aperture in the fluid network; performing DDA unsaturated steady-state seepage field calculation, performing unsaturated transient seepage field calculation using FDM to obtain steady-state calculation results; writing the steady-state calculation results into DDA as a fluid calculation module and assigning unsaturated transient seepage field calculation parameters, solving for the hydraulic head and volumetric water content at each node as transient calculation results; and performing fluid-solid coupling analysis based on the transient calculation results. The present invention addresses the shortcoming of DDA, which is its inability to simulate unsaturated transient seepage in rock and soil. The function of identifying rock cracks has been added, which can identify more complex situations such as concave polygons and long cracks, calculate the crack opening, and distinguish the outer boundary of the model; the present invention uses FDM to write the Richards equation into DDA, realizing DDA simulation of the unsaturated transient seepage process of rock and soil, and can be used in various working conditions such as rainfall and reservoir water level changes; the present invention can consider the interaction between rock and soil and water, that is, fluid-solid coupling calculation, which is more in line with the engineering practice of the influence of fluids such as rainfall on geotechnical engineering in engineering, and further provides benefits for engineering design and safety assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0040] Figure 1 It is a structural diagram of the DDA random fractured rock and soil calculation model of the present invention.

[0041] Figure 2 This is a schematic diagram of an output file of a two-dimensional array obtained by identifying rock mass fractures according to the present invention;

[0042] Figure 3 It is an effective saturation cloud diagram of the fluid-solid coupling calculation result of the present invention; Figure 4 yes Figure 1 Schematic cross-section diagram.

[0043] Figure 4 This is a flow chart of a non-saturated transient fluid-solid coupling calculation method based on FDM-DDA of the present invention.

[0044] Figure 5 It is a flow chart of the present invention for identifying rock mass fractures in the rock mass model and forming a fluid network.

[0045] Figure 6 It is a flow chart of the present invention for calculating the aperture of rock mass fractures in a fluid network. DETAILED DESCRIPTION

[0046] The following describes embodiments of the present invention in detail, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and are not to be construed as limiting the present invention.

[0047] See also Figures 1 to 6 The present invention provides a non-saturated transient fluid-solid coupling calculation method based on FDM-DDA, comprising:

[0048] S101 uses Matlab to generate cracks to build a rock and soil model;

[0049] The specific method of using matlab to generate cracks is: using matlab to write a Monte Carlo algorithm to generate the coordinates of the two endpoints of a random crack, and then importing them into the DDA pre-processing dc program to generate cracks.

[0050] Based on statistical principles, the probability distribution of the rock mass fracture's midpoint, inclination, dip, and length is determined. Using MATLAB's random number seed, the random command is used to generate random numbers that match the specified probability distribution. The coordinates of the two end points of the fracture are calculated using the straight-line point-slope formula of the randomly generated fracture. The coordinate data is exported and plotted using the .scr AutoCAD script file. The resulting plot is saved as a .dxf file and imported into DDA, split into line and point layers, to create the rock mass model.

[0051] S102 identifies rock mass fractures in the rock mass model to form a fluid network;

[0052] The specific steps are:

[0053] S201 obtains a line segment formed by two adjacent points, labels the endpoints sequentially, and stores the labels in the first column of the first two-dimensional data;

[0054] A DDA block is formed by connecting the coordinates of its vertices. Starting from a point in the block, the line is connected counterclockwise back to the starting point. Due to the characteristics of DDA data, the line segment formed by point i and the next point i+1 is represented by the sequence number of point i and stored in the first column of the first two-dimensional array lem.

[0055] S202 traverses the line segments to find two line segments whose endpoint coordinate differences are within a preset range and do not belong to the same block as two sides of the short crack;

[0056] Traverse the line segments and find two line segments whose endpoint coordinates are close enough and do not belong to the same block as the two sides of the rock mass fracture. Store the found line segment serial numbers in the second column of the two-dimensional array.

[0057] S203 stores the line segment numbers on both sides of the short crack in the second column of the first two-dimensional data;

[0058] Store the found line segment numbers in the second column of the two-dimensional array.

[0059] S204 confirms long cracks in the line segment based on parallel lines;

[0060] In the case of long fractures, the DDA block may not have a block vertex on one side of the fracture, so the above method cannot identify this type of fracture. Apart from this type of fracture, all other fractures are short fractures. Because each sequence number i uniquely corresponds to the line segment consisting of points i and i+1, long fractures only need to be screened among the nodes not stored in the first two-dimensional array.

[0061] S205 splits the long crack into short cracks, and stores the line segment numbers of the short cracks into the first two columns of the second two-dimensional array;

[0062] For this type of node, the other side of the fissure is determined by searching for parallel lines. Once a parallel line is found, the long fissure is split into several shorter fissures based on the containment relationship of the line segments. For the short fissures formed by splitting, the cosine of the angle between the vectors is used to determine whether any points fall on the short fissure line segment. If so, the fissure is further split; if not, the endpoint numbers of the line segment are stored in the first two columns of the second two-dimensional array lemm.

[0063] S206: Merge the same points in the first two-dimensional array and the second two-dimensional array to form a fluid network.

[0064] Finally, we merge the identical points in the two 2D arrays above. This merging occurs within the lem and lemm arrays, and between the lem and lemm arrays. When the coordinates of two points are close enough, they are considered the same point, representing two different line segments passing through the same node. By merging these line segments, we find the connection between each node and its surrounding nodes, and store the merged results in the first two columns of the 2D array leq.

[0065] S103 calculates the aperture of rock mass fractures in the fluid network;

[0066] The specific steps are:

[0067] S301 obtains the endpoint numbers of the line segments on both sides of the crack, the first and second points on one side, and the third and fourth points on the other side;

[0068] S302 calculates a first distance to a third point based on the first point and the second point;

[0069] According to the line segment numbers on both sides of the rock mass fracture stored in the above steps, the first distance e1 from one side of the fracture to one of the endpoints on the other side is calculated using formula (1).

[0070]

[0071] Among them, subscripts 1 and 2 represent the coordinates of the two endpoints on the same side of the crack (x1, y1) and (x2, y2), and subscript 3 represents the coordinates of the endpoint on the other side of the crack (x3, y3).

[0072] S303 calculates a second distance to a fourth point based on the first point and the second point;

[0073] Similarly, the second distance e2 between endpoints 1 and 2 and the endpoint 4 on the other side can be calculated.

[0074] S304 calculates the equivalent hydraulic opening based on the first distance and the second distance.

[0075] Calculate the equivalent hydraulic opening according to the empirical formula.

[0076]

[0077] In formula (2), r is the ratio of the distances e1 to e2. Since the initial cracks in the DDA model overlap on both sides, the initial crack opening needs to be manually input according to the project situation. When the result calculated by formula (2) is greater than the initial opening, the crack opening is updated.

[0078] S104 performs DDA unsaturated steady-state seepage field calculation, uses FDM to perform unsaturated transient seepage field calculation, and obtains steady-state calculation results;

[0079] The specific steps of calculating the DDA unsaturated steady-state seepage field are:

[0080] 401 The cubic law is used to calculate the seepage of a single fracture;

[0081] The cubic law is used to calculate the seepage of a single fracture. The seepage direction is determined according to the position relationship of the nodes and the flow of the single fracture is superimposed to obtain formula (3):

[0082]

[0083] In formula (3), T ij represents the seepage coefficient of fractures i and j connected to point i, T ii All T ij The algebraic sum of H i , Qi Represent the head and flow rate at point i respectively.

[0084] S402 uses Gauss–Seidel iteration to obtain the hydraulic head at each node and then calculates the volumetric water content at each node;

[0085] The Gauss-Seidel iteration is used to obtain the hydraulic head at each node, and then the volumetric water content of each node is calculated.

[0086] S e =[1+(αP c ) n ] -m (4)

[0087] In formula (4), P c is the pressure head at the node, α, n, and m are empirical coefficients in the van Genuchten model, which are related to the crack width and air intake value, and S e is the effective saturation of the node.

[0088] S403 writes the relevant program into the DDA program as a steady-state seepage module.

[0089] The specific steps of using FDM to calculate the unsaturated transient seepage field are:

[0090] S501 uses rock mass fractures as the fluid network required for FDM seepage field calculation;

[0091] S502 uses FDM to convert the Richards equation into a volumetric water content matrix that conforms to the DDA format.

[0092] The Richards equation in two dimensions is as follows:

[0093]

[0094] In formula (5), K(θ) is the unsaturated seepage coefficient, which is a function related to the volumetric water content θ of the rock mass. Rewrite it in the form of volumetric water content.

[0095]

[0096] In formula (6) Formula (6) is solved using FDM. First, the partial derivatives are written into a discretized format using central differences.

[0097]

[0098] In Equation (7), the superscript θ represents the time difference, and the subscript θ represents the spatial difference. The same method is used to calculate the y direction. Continue to discretize the partial derivatives in Equation (7).

[0099]

[0100] Since there is a certain hysteresis in the unsaturated seepage, the finite difference method D(θ j ) to approximate D(θ j+1 ). Arranging equations (7), (8), and (9) yields:

[0101]

[0102] There is a similar calculation process for the y direction. In order to simplify the expression, let

[0103] Combining equations (10) and (6), we can get

[0104]

[0105] In Equation (11), the superscript j represents the state at the current time step, which is a known quantity. j+1 represents the state at the next time step, which is solved using the state at the current time step. Δt is the time step of the transient seepage. Equation (11) is organized into a matrix form.

[0106]

[0107] In formula (12), A ij is the coefficient of the crack formed by points i and j, that is, α or β in formula (11).

[0108] A ii It's A ij The inverse of the sum of the algebraic sum of and the reciprocal of the time steps.

[0109] S105 writes the steady-state calculation results into the DDA as a fluid calculation module and assigns unsaturated transient seepage field calculation parameters, and solves the hydraulic head and volumetric water content at each node as the transient calculation results;

[0110] The coefficient matrix A in Equation (12) is a strictly diagonally dominant sparse matrix and is solved using SOR iteration similar to DDA. Solving Equation (12) requires finding the volumetric water content corresponding to time step j.

[0111]

[0112] In formula (13), S e is the effective saturation of the node, θ r is the residual water content, θ sis the saturated water content. The steady-state volumetric water content can be calculated based on the effective saturation, which serves as the initial and boundary conditions for the transient calculation. Boundary conditions are applied based on the actual project situation, and the parameters required for the unsaturated transient seepage field calculation are assigned. This is then incorporated into the DDA program as an unsaturated transient seepage module.

[0113] S106 performs fluid-solid coupling analysis based on transient calculation results.

[0114] The specific steps are:

[0115] S601 converts the transient calculation results into a water load matrix according to the minimum potential energy principle;

[0116] S602 adds the water load matrix to the DDA load matrix to perform DDA stress field calculation. When the block displacement is less than the specified value, the stress field calculation reaches equilibrium.

[0117] The overall equilibrium equation of DDA is as follows:

[0118]

[0119] In Equation (14), K is the block stiffness matrix, D is the deformation matrix, and F is the load matrix. The nodal head obtained in step S105 is converted into a water load matrix according to the principle of minimum potential energy and added to the DDA load matrix F to perform the DDA stress field calculation. When the block displacement at a certain time step is less than the specified value, the stress field calculation is considered to have reached equilibrium. The fracture aperture in step S103 is updated and the seepage module calculation is performed. Then, steps S104 and S105 are repeated to implement the DDA fluid-solid coupling analysis.

[0120] The invention is further described below using a two-dimensional problem as an example. A rock and soil model is established and imported into DDA. Two groups of cracks are set, with inclinations of 100-128° and 180-185°, dip angles of 73-84° and 69-78°, and crack lengths of 3-5m and 2-30m. Random cracks are generated by a Matlab program and imported into DDA to obtain Figure 1 ; The array obtained after the rock mass fracture identification in step 2 is as follows Figure 2 As shown, the initial fracture opening of the rock mass is set to 0.0003 m. When the equivalent hydraulic opening calculated in step S103 is greater than this value, it is updated.

[0121] The above disclosure is only a preferred embodiment of the present invention, and certainly cannot be used to limit the scope of the rights of the present invention. Ordinary technicians in this field can understand that all or part of the processes of the above embodiment and equivalent changes made in accordance with the claims of the present invention are still within the scope of the invention.

Claims

1. A non-saturated transient fluid-structure interaction calculation method based on FDM-DDA, It is characterized by: Including: using Matlab to generate cracks to build rock and soil models; Identifying rock fractures in the rock and soil model to form a fluid network; Calculation of rock mass fracture aperture in fluid networks; Perform DDA unsaturated steady-state seepage field calculations and use FDM to calculate unsaturated transient seepage field to obtain steady-state calculation results; The steady-state calculation results are written into DDA as a fluid calculation module and assigned with unsaturated transient seepage field calculation parameters. The hydraulic head and volumetric water content at each node are solved as transient calculation results. Fluid-solid coupling analysis is performed based on transient calculation results; the specific steps for calculating the DDA unsaturated steady-state seepage field are: using the cubic law to calculate the seepage of a single fracture; using Gauss-Seidel iteration to obtain the head of each node, and then calculating the volumetric water content of each node; writing the relevant program into the DDA program as a steady-state seepage module; the specific steps for using FDM to calculate the unsaturated transient seepage field are: using rock fractures as the fluid network required for FDM seepage field calculation; using FDM to convert the Richards equation into a volumetric water content matrix that conforms to the DDA format.

2. The unsaturated transient fluid-structure interaction calculation method based on FDM-DDA according to claim 1, characterized in that: The specific method of using matlab to generate cracks is: using matlab to write a Monte Carlo algorithm to generate the coordinates of the two endpoints of a random crack, and then importing them into the DDA pre-processing dc program to generate cracks.

3. The unsaturated transient fluid-solid coupling calculation method based on FDM-DDA according to claim 1, characterized in that: The specific steps of identifying rock mass fractures in the rock and soil model and forming a fluid network are: Obtain a line segment formed by two adjacent points, label the endpoints sequentially, and store the labels in the first column of the first two-dimensional data; Traversing the line segments, searching for two line segments whose endpoint coordinate differences are within a preset range and do not belong to the same block as the two sides of the short crack; Storing the line segment numbers on both sides of the short crack in the second column of the first two-dimensional data; identifying long cracks in the line segments based on parallel lines; Split the long cracks into short cracks, and store the line segment numbers of the short cracks in the first two columns of the second two-dimensional array; The same points in the first two-dimensional array and the second two-dimensional array are merged to form a fluid network.

4. The unsaturated transient fluid-structure interaction calculation method based on FDM-DDA according to claim 1, characterized in that: The specific steps of calculating the rock mass fracture aperture in the fluid network are: Get the endpoint numbers of the line segments on both sides of the crack. The numbers on one side are the first and second points, and the numbers on the other side are the third and fourth points. Calculate a first distance to a third point based on the first point and the second point; Calculate a second distance to a fourth point based on the first point and the second point; The equivalent hydraulic opening is calculated based on the first distance and the second distance.

5. The unsaturated transient fluid-solid coupling calculation method based on FDM-DDA according to claim 1, characterized in that: The specific steps of performing fluid-solid coupling analysis based on transient calculation results are: The transient calculation results are converted into water load matrix according to the minimum potential energy principle; The water load matrix is added to the DDA load matrix to perform DDA stress field calculation. When the block displacement is less than the specified value, the stress field calculation reaches equilibrium.

Images