A method for simulating pore water pressure dissipation of debris flow using porous filter grids based on SPH solid-liquid mixed flow model

Abstract

The present invention discloses a method for simulating the pore water pressure dissipation of debris flow using a porous filtration grid based on an SPH solid-liquid mixed flow model, comprising the following steps: 1) reading the basic terrain coordinates of the disaster area according to the debris flow geological survey report, and constructing a regional node spatial information array; 2) calculating the index in the terrain grid in combination with the grid spacing according to the node spatial information array read in step 1; 3) traversing all grid nodes in the terrain area; 4) initializing the required calculation variables and memory space. This method effectively identifies the control area of ​​the permeable grille. Appropriate boundary conditions and physical equations are applied to the nodes in the control area of ​​the permeable grille, thereby successfully simulating the effect of the permeable grille on the pore water pressure dissipation of the debris flow. It not only provides a scientific simulation method for the prevention and control of debris flow disasters, but also significantly improves the application effect of engineering measures in practice.

Description

Technical Field

[0001] The invention relates to the technical field of geological disaster simulation, in particular to a method for simulating the pore water pressure dissipation effect of debris flow in a porous filtering grid based on an SPH solid-liquid mixed flow model. Background Art

[0002] As a cutting-edge meshless numerical simulation technology, the smooth particle hydrodynamic method has shown its wide application value in the analysis and simulation of geological disasters such as debris flows. Given that these disasters are often caused by sudden heavy rainfall, and water is the main driving force for the long-distance movement of debris flows, how to effectively drain water to reduce the fluidity and destructive power of debris flows is an urgent problem to be solved. At present, most debris flow disaster prevention and control measures are not ideal in terms of blocking effect. As a new type of prevention method, porous filter grids use their unique structural design to achieve the discharge of water in debris flows, thereby reducing the flow velocity of debris flows. However, existing numerical simulation methods are difficult to accurately depict the changes in pore water pressure during debris flow movement. Summary of the invention

[0003] In view of this, the purpose of the present invention is to provide a method for simulating the dissipation of pore water pressure of debris flow in a porous filter grid based on the SPH solid-liquid mixed flow model, which significantly enhances the effect of debris flow prevention and control projects and provides key technical support for reducing disaster risks and losses.

[0004] To achieve the above object, the present invention adopts the following technical solution: a method for simulating the pore water pressure dissipation of debris flow in a porous filter grid based on the SPH solid-liquid mixed flow model, characterized in that it includes the following steps:

[0005] Step 1: Read the basic terrain coordinates of the disaster area according to the debris flow geological survey report and construct the regional node spatial information array; at the same time, create SPH particles based on the scale of debris flow disasters triggered by the source area;

[0006] Step 2: Based on the node spatial information array read in step 1, the index in the terrain grid is calculated in combination with the grid spacing; using the index, the four nearest terrain points are determined, and the elevation of the target point is calculated by applying the bilinear interpolation formula to construct the terrain grid of the study area;

[0007] Step 3: Traverse all grid nodes in the terrain area, and determine whether each grid node belongs to the porous filter grid control area according to the base type and elevation conditions of each grid node; once a point is determined to belong to the control area, immediately stop the loop processing of the point, and record the index of the point in the list of porous filter grid control areas;

[0008] Step 4: Initialize the required calculation variables and memory space; according to the results determined in step 3, for the nodes belonging to the porous filter grid control area, select and apply the corresponding boundary conditions, combine the diffusion and consolidation equations, calculate and simulate the changes in the excess pore water pressure of the debris flow when it flows in this area, and deduce the dissipation effect of the porous filter grid on the total pore water pressure during the movement of the debris flow.

[0009] 2. A method for simulating the pore water pressure dissipation of debris flow in a porous filter grid based on the SPH solid-liquid mixed flow model according to claim 1, characterized in that the regional node spatial information array topol(n, :) in step 1 retrieves different terrain spatial information stored therein according to different indexes in the first dimension, in particular:

[0010] 1) topol(1, ipoin): When the first dimension index of topol is 1, the terrain elevation value Z(ipoin) of each grid point is stored;

[0011] 2) topol(2, ipoin) and topol(3, ipoin): When the if_Gradient flag is true, the code calculates the terrain gradient in the x and y directions and stores them in Zpx and Zpy respectively; that is, the second and third dimensions of topol store the x and y derivatives of each grid point.

[0012] 3. A method for simulating the dissipation of pore water pressure of debris flow in a porous filter grid based on the SPH solid-liquid mixed flow model according to claim 1, characterized in that step 2 comprises the following steps:

[0013] Step 21: According to the two-dimensional coordinates xx and yy of each node combined with the grid spacing in the x direction deltxg , grid spacing in the y direction deltyg Calculate mesh element index ipoig :

[0014] (1)

[0015] In the formula, ipoigx, ipoigy are the grid element indices in the x and y directions respectively, npoigx is the number of grid nodes in the x direction, x min g and y min g are the minimum values ​​of the two-dimensional coordinates of the grid in the x and y directions respectively;

[0016] Step 22: Indexing Based on Mesh Elements ipoig, determine the four grid nodes n1, n2, n3, n4 closest to the research node:

[0017] (2)

[0018] Step 23: Call topol(1,ipoin) and grid element index ipoigx and ipoigy , determine the node reference coordinates x ref and y ref , and analyze the grid interpolation weights xi and η :

[0019] (3)

[0020] Step 24: Using Shape Functions sh_func And the node spatial information data topol(1, :), the elevation is calculated by bilinear interpolation Z :

[0021] (4)

[0022] in, topol(1,n1),topol(1,n2),topol(1,n3),topol(1,n4) They are the node space information data corresponding to the four grid nodes n1, n2, n3, and n4; sh_func(1), sh_func(2), sh_func(3), sh_ func(4) They correspond to the shape functions of the four mesh nodes based on bilinear interpolation.

[0023] In a preferred embodiment, step 3 comprises the following steps:

[0024] Step 31: Traverse all nodes of the terrain grid and record their regional geometry type typ_geom; make the following assumptions based on the corresponding geometry type:

[0025] Assumption a: When typ_geom=1, check whether the point is within a certain range of the Z axis;

[0026] Assumption b: When typ_geom=2, check whether the point is within the rectangular area;

[0027] Assumption c: When typ_geom=3, check whether the point is within the inclined rectangular area;

[0028] Step 32: Select the judgment logic according to the geometry type typ_geom to filter the ownership of each point, where:

[0029] Referring to hypothesis a, read the maximum value Z in the Z-axis direction of the terrain coordinates of the area max With the minimum value Z min, when the Z coordinate of the point is determined to be [Z min , Z max ], then the point belongs to this area;

[0030] Referring to hypothesis b, read the x and y coordinates of the lower left and upper right endpoints of the region, recorded as (x sw , y sw )、(x ne ,y ne ) and connect the two points to construct the diagonal of the rectangle to generate the rectangular area; when the x and y axis coordinates of the point fall within the range [x sw ,x ne ] and [y sw , y ne ], then the point belongs to this area;

[0031] Referring to hypothesis c, read the x and y coordinates of the four endpoints of the lower left, upper left, upper right, and lower right of the region, and record them as (x sw ,y sw )、(x nw , y nw )、(x ne , y ne ), (x se , y se ), connect the four endpoints in sequence to form an inclined rectangle. If a point falls in this inclined rectangle, the point belongs to this area;

[0032] When a point does not belong to the current determination area, it will proceed to the determination of the next area until the point is determined to belong to.

[0033] In a preferred embodiment, step 4 comprises the following steps:

[0034] Step 41: Allocate the memory space of the array by using the allocate statement, and initialize the excess pore water pressure value Δp of the current time step w rel , used to store information related to excess pore water pressure;

[0035] Step 42: Based on the region belonging to each point, its boundary condition LBc1 is determined. If the node boundary condition LBc1=1, the node belongs to the porous filtration grid control area and the bottom permeable boundary condition needs to be considered; otherwise, it does not belong to the porous filtration grid control area and the permeable boundary condition does not need to be considered.

[0036] Step 43: Use the corresponding excess pore water pressure equation to calculate the relative value Δp of the excess pore water pressure of debris flow in the porous filter grid control area at any time step. w rel , the calculation steps are as follows:

[0037] (5)

[0038] in, VBc1 is the ratio of liquid density to effective average density, hh is the relative height of debris flow mixture, alphaS is the volume fraction of the saturated layer of debris flow mixture, densw is the liquid density, dens_bar is the average density of the mixture, dens_bar_eff is the effective average density;

[0039] According to the total pore water pressure p w The excess pore water pressure Δp w and hydrostatic pressure p h The total pore water pressure is calculated based on the concept of the sum of the values ​​to determine the effectiveness of the porous filter grid.

[0040] Compared with the prior art, the present invention has the following beneficial effects: by adopting the SPH solid-liquid mixed flow model, the present invention finely constructs the terrain grid and accurately identifies the action area of ​​the porous filter grid, and successfully simulates the dissipation effect of the porous filter grid on the pore water pressure of debris flow. This achievement provides a scientific theoretical basis and numerical simulation tool for the design and application of porous filter grids, significantly enhances the effect of debris flow prevention and control projects, and provides key technical support for reducing disaster risks and losses. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 A method flow chart of a preferred embodiment of the present invention;

[0042] Figure 2 A schematic diagram of the region determination logic of a preferred embodiment of the present invention;

[0043] Figure 3 Schematic diagram of a method for calculating total pore water pressure according to a preferred embodiment of the present invention. DETAILED DESCRIPTION

[0044] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0045] It should be noted that the following detailed descriptions are illustrative and are intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present application belongs.

[0046] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application; as used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or their combinations.

[0047] At present, most debris flow disaster prevention and control measures are not ideal in terms of blocking effect. As a new type of prevention method, porous filter grids use their unique structural design to achieve the discharge of water in debris flows, thereby reducing the flow velocity of debris flows. However, existing numerical simulation methods are difficult to accurately depict the changes in pore water pressure during debris flow movement. The present invention successfully simulates the dissipation effect of porous filter grids on debris flow pore water pressure by adopting the SPH solid-liquid mixed flow model, finely constructing the terrain grid, and accurately identifying the action area of ​​the porous filter grid. This achievement provides a scientific theoretical basis and numerical simulation tools for the design and application of porous filter grids, significantly enhances the effectiveness of debris flow prevention and control projects, and provides key technical support for reducing disaster risks and losses. The specific implementation flow chart is shown in the attached figure. Figure 1 As shown. It mainly includes the following 4 steps:

[0048] Step 1: Read the basic terrain coordinates of the disaster area according to the debris flow geological survey report and construct the regional node spatial information array. At the same time, create SPH particles based on the scale of debris flow disasters triggered by the source area;

[0049] Step 2: Based on the node spatial information array read in step 1, the index in the terrain grid is calculated in combination with the grid spacing. Using these indexes, the four nearest terrain points are determined, and the elevation of the target point is calculated using the bilinear interpolation formula to construct the terrain grid of the study area;

[0050] Step 3: Traverse all grid nodes in the terrain area and determine whether they belong to the porous filter grid control area based on the base type and elevation conditions of each point. Once a point is determined to belong to the control area, stop the loop processing of the point immediately and record the index of the point in the list of porous filter grid control areas;

[0051] Step 4: Initialize the required calculation variables and memory space. According to the results determined in step 3, for the nodes belonging to the porous filter grid control area, select and apply the corresponding boundary conditions, combine the diffusion and consolidation equations, calculate and simulate the changes in the excess pore water pressure of the debris flow when it flows in this area, and deduce the dissipation effect of the porous filter grid on the total pore water pressure during the movement of the debris flow.

[0052] According to the three different types of regions, apply its judgment logic (see Figure 2 ), according to the method flow ( Figure 1 ) for a detailed introduction:

[0053] Step 1: Read the basic terrain coordinates of the disaster area according to the debris flow geological survey report. The coordinates of the lower left, upper left, upper right, and lower right endpoints of the area are (0,0), (0,100), (100,100), (100,0), and the range in the Z direction is between [0,100]. This is used to construct the regional node spatial information array topol(n,:). At the same time, 2000 debris flow SPH solid phase particles and liquid phase particles are generated based on the scale of debris flow disasters triggered by the source area.

[0054] Step 2: Based on the node spatial information array topol(n, :) read in step 1, the index ipoig in the terrain grid is calculated in combination with the grid spacing deltxg and deltyg in the x and y directions. Using these indexes, the four nearest terrain points n1, n2, n3, and n4 are determined, and the elevation Zp of the target point is calculated using the bilinear interpolation formula to construct the terrain grid of the study area. The following steps are included:

[0055] S2.1: Calculate the grid element index ipoig based on the two-dimensional coordinates xx and yy of each node combined with the grid spacing deltxg and deltyg. Take the node coordinate A (50,50) and the grid spacing of 2 and 2 as an example:

[0056] (1)

[0057] Where ipoigx and ipoigy are the grid element indices in the x and y directions respectively, npoigx is the number of grid nodes in the x direction, and xming and yming are the minimum values ​​of the two-dimensional coordinates of the grid in the x and y directions respectively.

[0058] S2.2: Based on the grid element index ipoig, determine the four grid nodes closest to the study node A:

[0059] (2)

[0060] S2.3: Call topol(1,ipoin) and grid indexes ipoigx and ipoigy to determine the node reference coordinate x ref and ref , and analyze the grid interpolation weights xi and η:

[0061] (3)

[0062] S2.4: Using the shape function sh_func and the node spatial information data topol(1, :), the elevation Zp is calculated by bilinear interpolation, where the elevations of the four nodes topol(1,n1), topol(1,n2), topol(1,n3), and topol(1,n4) are 48, 49, 51, and 52 respectively:

[0063] (4)

[0064] Therefore, the elevation of node A is calculated by bilinear interpolation to be 48. In the same way, the elevation of each node can be calculated.

[0065] Step 3: Traverse all grid nodes in the terrain area and determine whether they belong to the current control area based on the base type and elevation conditions of each point. Once a point is determined to belong to the control area, the loop processing of the point is immediately stopped and the index of the point is recorded in the list of porous filter grid control areas. The following steps are included:

[0066] S3.1: Traverse all nodes of the terrain grid and record their regional geometry type typ_geom. The following will make assumptions around three different geometry types:

[0067] Assumption a: When typ_geom=1, a porous filtering grid is planned to be arranged within the elevation range [45,50].

[0068] Assumption b: When typ_geom=2, a porous filter grid is arranged in a rectangular area surrounded by four two-dimensional coordinates: (45,45), (45,55), (55,55), and (55,45).

[0069] Assumption c: When typ_geom=3, a porous filter grid is arranged within the inclined rectangular range surrounded by four two-dimensional coordinates (40,45), (45,55), (65,55), and (60,45).

[0070] S3.2: Select the judgment logic according to the geometry type typ_geom to filter the belonging of each point:

[0071] Referring to hypothesis a, the elevation of node A is Zp=48, which is within the range controlled by the porous filtration grid, so it can be determined that it belongs to this area.

[0072] Referring to hypothesis b, the two-dimensional coordinates of node A are (50,50). The position of this point is within the rectangular range controlled by the porous filter grid, so it can be determined that it belongs to this area.

[0073] Referring to assumption c, the two-dimensional coordinates of node A are (50,50). The position of this point is within the inclined rectangular range controlled by the porous filter grid, so it can be determined that it belongs to this area.

[0074] Step 4: Initialize the required calculation variables and memory space. According to the results determined in step 3, for the nodes belonging to the porous filter grid control area, select and apply the corresponding boundary conditions, combine the diffusion and consolidation equations, calculate and simulate the changes in the excess pore water pressure of the debris flow when it flows in this area, and deduce the dissipation effect of the porous filter grid on the total pore water pressure during the movement of the debris flow. It includes the following steps:

[0075] S4.1: Initialize the calculation variables and memory space. Use the allocate statement to allocate the memory space of the array and initialize the excess pore water pressure value Δp of the current time step. w rel , used to store information related to excess pore water pressure.

[0076] S4.2: Based on the region to which each point belongs, determine its boundary condition LBc1:

[0077] Since the node boundary condition LBc1=1, the node belongs to the porous filtration grid control area, and the bottom permeable boundary condition needs to be considered.

[0078] S4.3: Use the corresponding excess pore water pressure equation to calculate the relative value of excess pore water pressure Δp of debris flow in the porous filter grid control area at any time step w rel , the calculation steps are as follows:

[0079] (5)

[0080] Wherein, VBc1 is the ratio of liquid density to effective average density, hh is the relative height of debris flow mixture, which is 1; alphaS is the volume fraction of saturated layer of debris flow mixture, which is 1; densw is the liquid density, which is 1000 kg / m 3 ; dens_bar is the average density of the mixture, which is 2000kg / m 3 ;dens_bar_eff is the effective average density.

[0081] According to the total pore water pressure p w The excess pore water pressure Δp w and hydrostatic pressure p h The concept of the sum of Figure 3 ), calculate the relative value of excess pore water pressure Δp w rel = -1, relative to the hydrostatic pressure p hrel The sum is the relative value of the total pore water pressure p w rel is 0, which means that the pore water pressure is completely dissipated on the porous filter grid.

[0082] It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above and that the invention can be implemented in other specific forms without departing from the spirit or essential features of the invention. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description, and it is intended that all variations within the meaning and scope of the elements of the claims be included in the invention. Any reference numeral in a claim should not be considered as limiting the claim to which it relates.

[0083] As described above, the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features thereof may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for simulating the pore water pressure dissipation of debris flow in porous filter grids based on the SPH solid-liquid mixed flow model, characterized in that: The following steps are involved: Step 1: Read the basic terrain coordinates of the disaster area according to the debris flow geological survey report and construct the regional node spatial information array; at the same time, create SPH particles based on the scale of debris flow disasters triggered by the source area; Step 2: Based on the node spatial information array read in step 1, the index in the terrain grid is calculated in combination with the grid spacing; using the index, the four nearest terrain points are determined, and the elevation of the target point is calculated by applying the bilinear interpolation formula to construct the terrain grid of the study area; Step 3: Traverse all grid nodes in the terrain area and determine whether each grid node belongs to the porous filtration grid control area according to its base type and elevation conditions; Once a point is determined to belong to the control area, the loop processing of the point is immediately stopped, and the index of the point is recorded in the list of the porous filter grid control area; Step 4: Initialize the required calculation variables and memory space; according to the results determined in step 3, for the nodes belonging to the porous filter grid control area, select and apply the corresponding boundary conditions, combine the diffusion and consolidation equations, calculate and simulate the change of excess pore water pressure of the debris flow when it flows in this area, and deduce the dissipation effect of the porous filter grid on the total pore water pressure during the movement of the debris flow; Step 4 includes the following steps: Step 41: Allocate the memory space of the array by using the allocate statement, and initialize the relative value of the excess pore water pressure Δp of the current time step w rel , used to store information related to excess pore water pressure; Step 42: Based on the region belonging to each point, its boundary condition LBc1 is determined. If the node boundary condition LBc1=1, the node belongs to the porous filtration grid control area and the bottom permeable boundary condition needs to be considered; otherwise, it does not belong to the porous filtration grid control area and the permeable boundary condition does not need to be considered. Step 43: Use the corresponding excess pore water pressure equation to calculate the relative value Δp of the excess pore water pressure of the debris flow in the porous filter grid control area at any time step. w rel , the calculation steps are as follows: （5） in, VBc1 is the ratio of liquid density to effective average density, hh is the relative height of debris flow mixture, alphaS is the volume fraction of the saturated layer of debris flow mixture, densw is the liquid density, dens_bar is the average density of the mixture, dens_bar_ eff is the effective average density; According to the total pore water pressure p w The excess pore water pressure Δp w and hydrostatic pressure p h The total pore water pressure is calculated based on the concept of the sum of the values ​​to determine the effectiveness of the porous filter grid.

2. The method for simulating the pore water pressure dissipation of debris flow using a porous filter grid based on the SPH solid-liquid mixed flow model according to claim 1, characterized in that: In step 1, the regional node spatial information array topol(n, :) retrieves different terrain spatial information stored therein according to different indexes in the first dimension, where: 1) topol(1, ipoin): When the first dimension index of topol is 1, the terrain elevation value Z(ipoin) of each grid point is stored; 2) topol(2, ipoin) and topol(3, ipoin): When the if_Gradient flag is true, the code calculates the terrain gradient in the x and y directions and stores them in Zpx and Zpy respectively; that is, the second and third dimensions of topol store the x and y derivatives of each grid point.

3. The method for simulating the pore water pressure dissipation of debris flow using a porous filter grid based on the SPH solid-liquid mixed flow model according to claim 1, characterized in that: Step 2 contains the following steps: Step 21: According to the two-dimensional coordinates xx and yy of each node combined with the grid spacing in the x direction deltxg , grid spacing in the y direction deltyg Calculate mesh element index ipoig : （1） In the formula, ipoigx, ipoigy are the grid element indices in the x and y directions respectively, npoigx is the number of grid nodes in the x direction, x min g and y min g are the minimum values ​​of the two-dimensional coordinates of the grid in the x and y directions respectively; Step 22: Indexing Based on Mesh Elements ipoig , determine the four grid nodes n1, n2, n3, n4 closest to the research node: （2） Step 23: Call topol(1,ipoin) and grid element index ipoigx and ipoigy , determine the node reference coordinates x ref and y ref , and analyze the grid interpolation weights xi and η : （3） Step 24: Using Shape Functions sh_func And the node spatial information data topol(1, :), the elevation is calculated by bilinear interpolation Z : （4） in, topol(1,n1),topol(1,n2),topol(1,n3),topol(1,n4) They are the node space information data corresponding to the four grid nodes n1, n2, n3, and n4; sh_func(1), sh_func(2), sh_func(3), sh_func(4) They correspond to the shape functions of the four mesh nodes based on bilinear interpolation.

4. The method for simulating the pore water pressure dissipation of debris flow using a porous filter grid based on the SPH solid-liquid mixed flow model according to claim 1, characterized in that: Step 3 contains the following steps: Step 31: Traverse all nodes of the terrain grid and record their regional geometry type typ_geom; make the following assumptions based on the corresponding geometry type: Assumption a: When typ_geom=1, check whether the point is within a certain range of the Z axis; Assumption b: When typ_geom=2, check whether the point is within the rectangular area; Assumption c: When typ_geom=3, check whether the point is within the inclined rectangular area; Step 32: Select the judgment logic according to the geometry type typ_geom to filter the ownership of each point, where: Referring to hypothesis a, read the maximum value Z in the Z-axis direction of the terrain coordinates of the area max With the minimum value Z min , when the Z coordinate of the point is determined to be [Z min , Z max ], the point belongs to this area; Referring to hypothesis b, read the x and y coordinates of the lower left and upper right endpoints of the region, recorded as (x sw , y sw )、(x ne ,y ne ) and connect the two points to construct the diagonal of the rectangle to generate the rectangular area; when the x and y axis coordinates of the point fall within the range [x sw ,x ne ] and [y sw , y ne ], then the point belongs to this area; Referring to hypothesis c, read the x and y coordinates of the four endpoints of the lower left, upper left, upper right, and lower right of the region, and record them as (x sw ,y sw )、(x nw , y nw )、(x ne , y ne ), (x se , y se ), connect the four endpoints in sequence to form an inclined rectangle. If a point falls in this inclined rectangle, the point belongs to this area; When a point does not belong to the current determination area, it will proceed to the determination of the next area until the point is determined to belong to.

Images