A debris flow confluence simulation method and system based on multi-grid method

Abstract

The present invention provides a method and system for simulating debris flow confluence areas based on a multi-grid method, relating to the technical field of multiphase flow simulation. Based on the existing MAC method for simulating debris flow confluence areas, the present invention introduces multi-grid technology to accelerate the solution process. The simulation results can well simulate the morphology and water level distribution of the intersection of the debris flow and the mainstream, and thus be used in the field of debris flow disaster prevention and control. The multi-grid technology makes the algorithm more computationally efficient and converges faster, making it more suitable for large-scale grid computing solutions. Compared with traditional MAC solution technology, the simulation of debris flow confluence areas is more efficient and less time-consuming, which is beneficial to large-scale grid computing solutions.

Description

Technical Field

[0001] The present invention relates to the technical field of multiphase flow simulation, and in particular to a method and system for simulating a debris flow intersection area based on a multi-grid method. Background Art

[0002] Debris flows are a common natural disaster in mountainous areas, causing significant casualties and property damage worldwide each year. Debris flows can travel great distances, spreading over vast areas, and possess high speed and propagation capabilities even on very gentle slopes. Debris flows have become a pressing issue in disaster prevention and mitigation efforts. Existing technologies, using multiphase flow simulation, can simulate debris flow confluence zones, thereby predicting flow trends and other conditions.

[0003] The challenge in multiphase flow simulation lies in tracking and identifying material interfaces. Existing multiphase flow simulation techniques can be categorized into two main types: the Euler-Lagrangian method and the Euler-Euler method. Among these, the Marked Grid Method (MAC) is a preferred method that belongs to the Euler-Lagrangian method.

[0004] Existing research, "MAC Method for Coupled Analysis of Debris Flow and Water Flow Confluence, Chen Chunguang, Yao Lingkan, et al., 2001(02)," shows that the MAC method can effectively simulate the morphology and water level distribution of the intersection of debris flow and mainstream, and is suitable for the field of debris flow disaster prevention and control. However, this method requires an iterative method to solve the pressure Poisson equation during the solution process. This process is usually computationally intensive and time-consuming, making it unsuitable for large-scale grid calculation solutions, limiting its application in the simulation of large-scale debris flow intersection areas.

[0005] Therefore, it is necessary to provide a debris flow confluence area simulation method and system based on the multi-grid method to solve the above technical problem. Summary of the Invention

[0006] To solve one of the above technical problems, the present invention proposes a multigrid-based debris flow confluence simulation method. A semi-staggered grid is established, which includes a number of discrete points. The dynamic equations of the confluence are established in the semi-staggered grid, and the debris flow / main river are simulated using a number of marker points. During the simulation, initialization parameters are set, the positions of the marker points at each time are solved, and the parameter values of each parameter are iterated to obtain the simulation results of the debris flow confluence. The specific steps are as follows:

[0007] Step S1: Parameter initialization and assignment: Set each initialization parameter and assign the parameter value;

[0008] Step S2: Discretization and construction of the grid: construct grid cells and marker points by initializing parameters and assigning corresponding parameter values, construct a semi-staggered grid, and set boundary conditions;

[0009] Step S3: Time advancement process and iteration: according to the initialization parameters, time advancement is performed, and an iterative method is used in the advancement process. When the iteration converges within the error, the position of the mark point and the parameter value at the next moment are output; the position of the mark point and the parameter value at the next moment are used in the iterative time advancement process to obtain the position of the mark point and the parameter value at each moment; wherein, the time advancement process includes the calculation of the intermediate speed, the iteration of the density and viscosity, and the solution of the pressure Poisson equation;

[0010] Step S4 post-processes the display image: obtains the marker point position and parameter value, and post-processes the marker point position and parameter value according to the display settings set by the user to obtain a post-processed display image, thereby completing the simulation of the debris flow intersection area.

[0011] As a further solution, the initialization parameters include intersection area simulation parameters and post-processing display parameters; the intersection area simulation parameters are used to set the initialization parameters required for debris flow intersection area simulation, including grid encryption coefficient, main river length, main river width, debris flow width, main river lateral initial velocity, main river longitudinal initial velocity, debris flow lateral initial velocity, debris flow longitudinal initial velocity, main river initial density, debris flow initial density, main river initial viscosity, debris flow initial viscosity and riverbed bottom friction coefficient; the post-processing display parameters are used to set the user's display settings, including debris flow marker outer contour line parameters, real-time cloud map parameters, post-processing display parameters and post-processing display variable parameters.

[0012] As a further solution, the boundary of the semi-staggered grid is composed of the main river inflow boundary, the main river outflow boundary, the river bank solid wall boundary and the debris flow inflow boundary;

[0013] The semi-staggered grid is composed of several grid cell units, and each marked point is numbered by a two-dimensional coordinate axis; the left discrete point of the grid cell unit stores the velocity u with a number (i-1 / 2, j), the grid cell unit discrete point stores the right velocity u with a number (i+1 / 2, j), the upper velocity v stored by the grid cell unit discrete point is numbered (i, j-1 / 2), the lower velocity v stored by the grid cell unit discrete point is numbered (i, j+1 / 2), and the grid center discrete point stores the velocity (i, j); the main river transverse initial velocity and the main river longitudinal initial velocity are stored at the cell boundary center discrete point of the main river inflow boundary; the debris flow transverse initial velocity and the debris flow longitudinal initial velocity are stored at the cell boundary center discrete point of the debris flow inflow boundary.

[0014] As a further solution, the main river inflow boundary and the main river outflow boundary adopt fixed boundary conditions, the river bank solid wall boundary adopts sliding boundary conditions, the middle inflow area of the debris flow inflow boundary adopts fixed boundary conditions, and the edge areas on both sides adopt sliding boundary conditions.

[0015] As a further solution, the time advancement process is iteratively performed through the following steps:

[0016] Step A1 calculates the parameter corresponding to the current moment n and And read the fluid density and fluid viscosity at the current moment;

[0017] Step A2: Based on the parameters obtained in step A1, the pressure field is further solved iteratively using the multi-grid technique;

[0018] Step A3: Calculate the intermediate step speed according to the preset formula based on the pressure field, where the intermediate step speed includes the intermediate step lateral speed and the intermediate step longitudinal velocity

[0019] Step A4: Calculate the movement position of the marker point within the time dt according to the intermediate step speed, and update the distribution of the marker point in the semi-staggered grid;

[0020] Step A5: updating the fluid density and fluid viscosity of each marking point according to the distribution of the marking points in the semi-staggered grid;

[0021] Step A6: Repeat steps A2 to A5 until the density change of the marked points in each grid cell is minimal, and set the current fluid density and fluid viscosity values as the initial values of the next iteration at time n+1; wherein the minimal change is determined by whether the density change is less than a preset minimum value;

[0022] Step A7: Repeat steps A1 to A6, and further iterate to obtain the position of the marking point and the parameter value at each moment.

[0023] As a further solution, the parameter corresponding to the current moment n and Calculated by the following formula 1:

[0024]

[0025] Where the superscript represents the time at which the parameter is located, the subscript represents the index of the semi-staggered grid position where the parameter is located, ρ represents the fluid density at the corresponding time and position, u represents the transverse velocity of the fluid at the corresponding time and position, v represents the longitudinal velocity of the fluid at the corresponding time and position, and Vis represents the discrete format of the viscosity term in the x / y direction.

[0026] As a further solution, the pressure field is iteratively solved by the following formula 2:

[0027]

[0028] Among them, the superscript represents the time at which the parameter is located, and the subscript represents the number of the semi-staggered grid position where the parameter is located. Represents the convergence judgment factor, the parameter at each moment and All are solved by the results of the previous moment; ρ represents the fluid density at the corresponding moment and position; P represents the fluid pressure at the corresponding moment and position;

[0029] The pressure field expression is solved by the multigrid method of V cycle, and the steps are as follows:

[0030] Step B1: Pre-smoothing: Formula 2 is rewritten to obtain Formula 3, which is:

[0031]

[0032] Among them, at time n The value at each grid point is known, and the It is calculated by formula 2, which requires Substitute into formula 3 and iterate to get the pressure at the next moment

[0033] The left side of Formula 3 is replaced by the coefficient matrix A h and the solution matrix u of pressure P h The product of is expressed as f on the right side. h Matrix representation, we get Formula 4:

[0034]

[0035] A h u h =f h

[0036] The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid. 2h is a grid twice as thick as the h layer. 4h is a grid twice as thick as the 2h layer. And so on, up to Nh. N is the maximum value that n can take. The default minimum number of grids in the coarsest grid layer Nh is not less than 8.

[0037] Perform 3 to 5 Gauss-Seidel iterations according to Formula 4 to obtain the solution matrix u h Estimated value of v h ;

[0038] Step B2 defines: through the f of the previous grid h Calculate the residual value f of the next coarse grid 2h , and insert it into the coarse grid 2h to get the estimated value v on the coarse grid 2h The specific process is:

[0039] First calculate f using formula 5 2h , the formula 5:

[0040]

[0041] Then perform 3 to 5 Gauss-Seidel iterations to solve Equation 6 and obtain the solution matrix u 2h Estimated value of v 2h ; Formula 6:

[0042] A 2h u 2h =f 2h

[0043] The process of defining the fine grid (n-1)h to the coarse grid nh is completed by performing multiple iterations and is described by Formula 7:

[0044] A (n-1)h v (n-1)h (fine grid) → A nh v nh (coarse grid);

[0045] Step B3: Grid coarsening: The coarse grid is solved by performing a fully weighted operation using Formula 8, and step B2 is repeated to coarsen the grid to the coarsest grid Nh at the bottom layer, and then Gauss-Seidel iteration is performed until convergence, Formula 8;

[0046]

[0047] The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid. 2h is a grid twice as thick as the h layer. 4h is a grid twice as thick as the 2h layer. And so on, up to Nh. N is the maximum value that n can take. The default minimum number of grids in the coarsest grid layer Nh is not less than 8.

[0048] Step B4: According to the interpolation operator formula, the approximate solution of the residual value ε is extended from the coarsest grid. nh Interpolate to the fine grid of the previous layer to get ε (n-1)h , and use the revised formula to estimate the value v of the previous layer (n-1)h Make corrections and get a new estimated value v(new) (n-1)h; Then continue to interpolate to the finer grid (n-2)h of the next layer, repeat until you return to the finest grid h, and get the final corrected solution v(new) h =v h +ε h , at this time v(new) h That is the solution u h ; The interpolation operator formula:

[0049]

[0050] Among them, the residual value Calculated using the residual value calculation formula; the superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid, 2h is a grid twice as thick as the h layer, 4h is a grid twice as thick as the 2h layer, and so on, all the way to Nh, where N is the maximum value that n can take. The default minimum number of grids in the coarsest grid layer Nh is not less than 8;

[0051] The correction formula:

[0052] v(new) (n-1)h =v (n-1)h +ε (n-1)h

[0053] The residual value calculation formula is:

[0054]

[0055] Among them, ε nh is the residual value on the n-th grid, A nh v nh It can be obtained by formula 4, is the value calculated by formula 4 for the previous grid The initial value on the nth grid layer after interpolation.

[0056] As a further solution, the intermediate step lateral speed Calculated by the following formula:

[0057]

[0058] Among them, the superscript represents the time at which the parameter is located, the subscript represents the number of the semi-staggered grid position where the parameter is located, and the parameter is the quantity obtained at this moment, ρ represents the fluid density at the corresponding moment and position, and is obtained through pressure field iteration; the longitudinal velocity of the intermediate step is Calculated by the following formula:

[0059]

[0060] Among them, the superscript represents the time at which the parameter is located, the subscript represents the number of the semi-staggered grid position where the parameter is located, and the parameter is the quantity that has been calculated at that moment, ρ represents the fluid density at the corresponding moment and position, and is obtained through pressure field iteration.

[0061] As a further solution, the fluid density ρ is iterated by the following formula:

[0062]

[0063] The marking points include debris flow marking points and main river marking points, ρ1 main river initial density, ρ2 debris flow initial density, m i,j Indicates the number of main river markers in the grid cell corresponding to the coordinates, n i,j Indicates the number of debris flow markers in the grid cell corresponding to the coordinates;

[0064] The fluid viscosity μ is iterated by the following steps:

[0065]

[0066] The marking points include debris flow marking points and main river marking points, μ1 initial viscosity of the main river, μ2 initial viscosity of the debris flow, m i,j Indicates the number of main river markers in the grid cell corresponding to the coordinates, n i,j Indicates the number of debris flow markers in the grid cell corresponding to the coordinates.

[0067] A debris flow confluence area simulation system based on a multi-grid method comprises an input parameter terminal, an intersection area simulation terminal and a simulation display terminal; the input parameter terminal is used to set various initialization parameters and perform parameter assignment; the intersection area simulation terminal performs initialization settings according to the intersection area simulation parameters, and uses a debris flow confluence area simulation method based on a multi-grid method as described in any of the above items to perform debris flow confluence area simulation; the simulation display terminal performs post-processing display settings according to the post-processing display parameters, and outputs images and display variables corresponding to the post-processing display settings.

[0068] Compared with related technologies, the method and system for simulating debris flow confluence areas based on the multi-grid method provided by the present invention have the following beneficial effects:

[0069] The present invention simulates the debris flow confluence area through the MAC method, which can well simulate the morphology and water level distribution of the intersection of the debris flow and the mainstream, and can be well applied to the field of debris flow disaster prevention and control. In addition, when solving the pressure Poisson equation, the MAC method using multi-grid technology is used to accelerate the iteration of this solution process. The multi-grid method has higher computational efficiency and faster convergence, making the algorithm more suitable for large-scale grid computing and solving. Compared with the traditional MAC solution method, it can more effectively simulate the debris flow confluence area, with less computational complexity and shorter time consumption. Compared with existing methods, it is more suitable for large-scale grid computing and solving. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 A flowchart of an embodiment of a debris flow confluence area simulation method based on a multi-grid method provided by the present invention;

[0071] Figure 2 A schematic diagram of setting boundary conditions provided by the present invention;

[0072] Figure 3 A diagram of a semi-staggered grid structure provided by the present invention;

[0073] Figure 4 A flow map of the marking points in the debris flow area provided by the present invention;

[0074] Figure 5 A distribution diagram of marking points provided by the present invention;

[0075] Figure 6 A simulation diagram of the intersection area of debris flows with different bulk densities provided by the present invention;

[0076] Figure 7 A schematic diagram of a preferred embodiment of the input parameter end provided by the present invention;

[0077] Figure 8 A simulated image of a debris flow confluence area of the debris flow confluence area simulation system provided by the present invention;

[0078] Figure 9 This is a time-consuming comparison chart between the method provided by the present invention and the traditional method. DETAILED DESCRIPTION

[0079] Various exemplary embodiments, features, and aspects of the present disclosure will be described in detail below with reference to the accompanying drawings. The same reference numerals in the accompanying drawings represent elements with the same or similar functions. Although various aspects of the embodiments are shown in the accompanying drawings, the drawings are not necessarily drawn to scale unless otherwise indicated.

[0080] The word “exemplary” is used exclusively herein to mean “serving as an example, example, or illustration.” Any embodiment described herein as “exemplary” is not necessarily to be construed as preferred or advantageous over other embodiments.

[0081] In addition, numerous specific details are provided in the following detailed description to better illustrate the present disclosure. Those skilled in the art will appreciate that the present disclosure can be practiced without certain specific details. In some instances, methods, means, components, and circuits well known to those skilled in the art are not described in detail in order to highlight the main points of the present disclosure.

[0082] like Figures 1 to 9 As shown, the present invention provides a debris flow confluence area simulation method based on the multigrid method, establishing a semi-staggered grid, the semi-staggered grid including a plurality of marking points, establishing a confluence area dynamic equation in the semi-staggered grid, and simulating the debris flow / main river through the plurality of marking points; during the simulation: setting initialization parameters, solving the positions of the marking points at each time, and iterating the parameter values of each parameter, thereby obtaining the simulation result of the debris flow confluence area; the specific steps are as follows:

[0083] Step S1: Parameter initialization and assignment: Set each initialization parameter and assign the parameter value;

[0084] Step S2: Discretization and construction of the grid: construct grid cells and marker points by initializing parameters and assigning corresponding parameter values, construct a semi-staggered grid, and set boundary conditions;

[0085] Step S3: Time advancement process and iteration: according to the initialization parameters, time advancement is performed, and an iterative method is used in the advancement process. When the iteration converges within the error, the position of the mark point and the parameter value at the next moment are output; the position of the mark point and the parameter value at the next moment are used in the iterative time advancement process to obtain the position of the mark point and the parameter value at each moment; wherein, the time advancement process includes the calculation of the intermediate speed, the iteration of the density and viscosity, and the solution of the pressure Poisson equation;

[0086] Step S4 post-processes the display image: obtains the marker point position and parameter value, and post-processes the marker point position and parameter value according to the display settings set by the user to obtain a post-processed display image, thereby completing the simulation of the debris flow intersection area.

[0087] It should be noted that since the pressure Poisson equation obtained by the intersection coupling equation has the properties of an elliptic equation, the multi-grid technology is used to accelerate the calculation of the intersection dynamic equation, which is suitable for solving this process. The intersection equation is as follows:

[0088]

[0089] Formula (1) is the two-dimensional flow equation in the intersection area, where is the single-width tangential force vector of the fluid, is the bottom friction force vector, is the static pressure of the mixed flow.

[0090]

[0091] Formula (2) is the vector representation of the dynamic equation of the intersection zone. In order to describe its calculation process more specifically, the discrete method is first derived based on the shallow water equation of the intersection zone. The shallow water equation is expressed as:

[0092]

[0093] Formula (3) assumes that H≈h, and H does not change with time and space, and has the relationship:

[0094]

[0095] Divide both sides of equation (3) by H(h) and substitute into equation (4), we have:

[0096]

[0097] in,

[0098]

[0099] Substituting formula (6) into formula (5), we have:

[0100]

[0101] Among them, C f is the bottom friction coefficient, η m is the generalized eddy viscosity coefficient of the intersection zone, and H(h) is the water depth. By deducing and solving the above equation, we can simulate the debris flow intersection zone. After obtaining the dynamic equation of the intersection zone, we consider the inviscid term of the equation:

[0102]

[0103] Discretize Equation (A-1), we have:

[0104]

[0105] in,

[0106]

[0107] According to the incompressible condition, we have:

[0108]

[0109] In the above formula, the velocities u and v are stored at the midpoints of the grid boundaries; ρ and η m , P is stored at the center of the grid.

[0110] Substituting formula (A-2) into (A-4), we have:

[0111]

[0112] The above formula is the variable coefficient pressure Poisson equation. This equation can be solved by iteratively solving the density during the calculation time, and then solving the pressure and then solving the velocity.

[0113] As a further solution, the initialization parameters include intersection area simulation parameters and post-processing display parameters; the intersection area simulation parameters are used to set the initialization parameters required for debris flow intersection area simulation, including grid encryption coefficient, main river length, main river width, debris flow width, main river lateral initial velocity, main river longitudinal initial velocity, debris flow lateral initial velocity, debris flow longitudinal initial velocity, main river initial density, debris flow initial density, main river initial viscosity, debris flow initial viscosity and riverbed bottom friction coefficient; the post-processing display parameters are used to set the user's display settings, including debris flow marker outer contour line parameters, real-time cloud map parameters, post-processing display parameters and post-processing display variable parameters.

[0114] As a further solution, Figure 2 As shown, the boundary of the semi-staggered grid consists of the main river inflow boundary, the main river outflow boundary, the river bank solid wall boundary and the debris flow inflow boundary; the main river inflow boundary and the main river outflow boundary adopt fixed boundary conditions, the river bank solid wall boundary adopts a sliding boundary condition, the middle inflow area of the debris flow inflow boundary adopts a fixed boundary condition, and the edge areas on both sides adopt a sliding boundary condition.

[0115] It should be noted that: Figure 2 As shown, fixed boundary conditions (Dirichlet boundary conditions) are used for the inflow and outflow boundaries of the main river and debris flow, and sliding boundary conditions are used for the river bank boundaries.

[0116] As a further solution, Figure 3As shown, the semi-staggered grid is composed of several grid cell units, and each marked point is labeled by a two-dimensional coordinate axis; the left discrete point of the grid cell unit stores the velocity u with a label of (i-1 / 2, j), the grid cell unit discrete point stores the right velocity u with a label of (i+1 / 2, j), the upper velocity v stored by the grid cell unit discrete point is labeled as (i, j-1 / 2), the lower velocity v stored by the grid cell unit discrete point is labeled as (i, j+1 / 2), and the grid center discrete point stores the label (i, j); the main river transverse initial velocity and the main river longitudinal initial velocity are stored at the cell boundary center discrete point of the main river inflow boundary; the debris flow transverse initial velocity and the debris flow longitudinal initial velocity are stored at the cell boundary center discrete point of the debris flow inflow boundary.

[0117] It should be noted that: Figure 4 As shown, the marking points are numbered, and the last row of marking points is somewhere outside the boundary. When the last row is greater than the boundary value, the marking points in the debris flow area are refreshed and the flow continues.

[0118] As a further solution, the time advancement process is iteratively performed through the following steps:

[0119] Step A1 calculates the parameter corresponding to the current moment n and And read the fluid density and fluid viscosity at the current moment;

[0120] Step A2: Based on the parameters obtained in step A1, the pressure field is further solved iteratively using the multi-grid technique;

[0121] Step A3: Calculate the intermediate step speed according to the preset formula based on the pressure field, where the intermediate step speed includes the intermediate step lateral speed and the intermediate step longitudinal velocity

[0122] Step A4: Calculate the movement position of the marker point within the time dt according to the intermediate step speed, and update the distribution of the marker point in the semi-staggered grid;

[0123] Step A5: updating the fluid density and fluid viscosity of each marking point according to the distribution of the marking points in the semi-staggered grid;

[0124] Step A6: Repeat steps A2 to A5 until the density change of the marked points in each grid cell is minimal, and set the current fluid density and fluid viscosity values as the initial values of the next iteration at time n+1; wherein the minimal change is determined by whether the density change is less than a preset minimum value;

[0125] Step A7: Repeat steps A1 to A6, and further iterate to obtain the position of the marking point and the parameter value at each moment.

[0126] As a further solution, the parameter corresponding to the current moment n and Calculated by the following formula:

[0127]

[0128] Where the superscript represents the time at which the parameter is located, the subscript represents the index of the semi-staggered grid position where the parameter is located, ρ represents the fluid density at the corresponding time and position, u represents the transverse velocity of the fluid at the corresponding time and position, v represents the longitudinal velocity of the fluid at the corresponding time and position, and Vis represents the discrete format of the viscosity term in the x / y direction.

[0129] As a further solution, the pressure field is iteratively solved by the following expression 2:

[0130]

[0131] Among them, the superscript represents the time at which the parameter is located, and the subscript represents the number of the semi-staggered grid position where the parameter is located. Represents the convergence judgment factor, the parameter at each moment and All are solved by the results of the previous moment; ρ represents the fluid density at the corresponding moment and position; P represents the fluid pressure at the corresponding moment and position;

[0132] The pressure field expression is solved by the multigrid method of "V" cycle, the steps are as follows:

[0133] Step B1: Pre-smoothing: Rewrite Formula 2 to obtain Formula 3, which is:

[0134]

[0135] Among them, at time n The value at each grid point is known, and the It can be calculated by formula 2 by Substitute into formula 3 and iterate to get the pressure at the next moment

[0136] The left side of Formula 3 is replaced by the coefficient matrix A h and the solution matrix u of pressure P h The product of is expressed as f on the right side. h Matrix representation:

[0137]

[0138] The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid (for example, a 64×64 grid), 2h is a grid twice as thick as the h layer (for example, a 32×32 grid), 4h is a grid twice as thick as the 2h layer (for example, a 16×16 grid), and so on... all the way to Nh, where N is the maximum value that n can take. By default, the minimum number of grids in the coarsest grid layer Nh is no less than 8.

[0139] Perform 3 to 5 Gauss-Seidel iterations according to Formula 4 to obtain the solution matrix u h Estimated value of v h .

[0140] Step B2 defines: through the f of the previous grid h Calculate the residual value f of the next coarse grid 2h , and insert it into the coarse grid 2h to get the estimated value v on the coarse grid 2h The specific process is:

[0141] First calculate f using formula 5 2h ,have:

[0142]

[0143] Then, by solving formula 4 in a similar way, perform 3 to 5 Gauss-Seidel iterations to solve formula 6) to obtain the solution matrix u 2h Estimated value of v 2h .

[0144] A 2h u 2h =f 2h (6)

[0145] By performing the above operation multiple times, the process of limiting the fine grid (n-1)h to the coarse grid nh is completed, which can be described by the following formula 7.

[0146] A (n-1)h v (n-1)h (fine grid) → A nh v nh (coarse grid)(7)

[0147] Step B3: Grid coarsening: The fully weighted equation 8 is used to solve the coarse grid, and step B2 is repeated to coarsen the grid to the coarsest grid Nh at the bottom layer (the default minimum number of grids in the coarsest grid layer is not less than 8), and then Gauss-Seidel iteration is performed until convergence;

[0148]

[0149] The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid (for example, a 64×64 grid), 2h is a grid twice as thick as the h layer (for example, a 32×32 grid), 4h is a grid twice as thick as the 2h layer (for example, a 16×16 grid), and so on... all the way to Nh, where N is the maximum value that n can take. By default, the minimum number of grids in the coarsest grid layer Nh is no less than 8.

[0150] Step B4: According to the interpolation operator formula 9, the approximate solution of the residual value ε is taken from the coarsest grid. nh Interpolate to the fine grid of the previous layer to get ε (n-1)h , and use the modified formula 10 to estimate the value v of the previous layer (n-1)h Make corrections and get a new estimated value v(new) (n-1)h Then continue to interpolate to the finer grid (n-2)h of the next layer, and repeat this process until you return to the finest grid h and get the final corrected solution v(new) h =v h +ε h , at this time v(new) h That is the solution u h .

[0151] Interpolation operator formula:

[0152]

[0153] Among them, the residual value This can be calculated using Equation 11. The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it indicates the initial finest grid (e.g., a 64×64 grid). 2h is a grid twice as thick as the h layer (e.g., a 32×32 grid). 4h is a grid twice as thick as the 2h layer (e.g., a 16×16 grid). And so on, up to Nh, where N is the maximum value that n can take. By default, the minimum number of grid cells in the coarsest grid layer Nh is no less than 8.

[0154] Correction formula:

[0155] v(new) (n-1)h =v (n-1)h +ε (n-1)h (10)

[0156] The residual value calculation formula is:

[0157]

[0158] Among them, ε nh is the residual value on the n-th grid, A nh v nh It can be obtained by formula 4, is the value calculated by formula 4) for the previous grid The initial value on the nth grid layer after interpolation.

[0159] As a further solution, the intermediate step lateral speed Calculated by the following formula:

[0160]

[0161] Among them, the superscript represents the time at which the parameter is located, the subscript represents the number of the semi-staggered grid position where the parameter is located, and the parameter is the quantity obtained at this moment, ρ represents the fluid density at the corresponding moment and position, and is obtained through pressure field iteration; the longitudinal velocity of the intermediate step is Calculated by the following formula:

[0162]

[0163] Among them, the superscript represents the time at which the parameter is located, the subscript represents the number of the semi-staggered grid position where the parameter is located, and the parameter is the quantity that has been calculated at that moment, ρ represents the fluid density at the corresponding moment and position, and is obtained through pressure field iteration.

[0164] As a further solution, the fluid density ρ is iterated by the following formula:

[0165]

[0166] The marking points include debris flow marking points and main river marking points, ρ1 main river initial density, ρ2 debris flow initial density, m i,j Indicates the number of main river markers in the grid cell corresponding to the coordinates, n i,j Indicates the number of debris flow markers in the grid cell corresponding to the coordinates;

[0167] The fluid viscosity μ is iterated by the following steps:

[0168]

[0169] The marking points include debris flow marking points and main river marking points, μ1 initial viscosity of the main river, μ2 initial viscosity of the debris flow, m i,j Indicates the number of main river markers in the grid cell corresponding to the coordinates, n i,j Indicates the number of debris flow markers in the grid cell corresponding to the coordinates.

[0170] It should be noted that: Figure 5As shown, the iteration of density and viscosity is related to the marking points. Assume that there are two kinds of marking points in a grid and their numbers are m and m, respectively. i,j and n i,j , the density and viscosity are ρ1, ρ2, μ1, μ2 respectively, then we can calculate the density and viscosity of the intersection area by their number ratio.

[0171] Figure 6 This is a cloud diagram of a debris flow confluence area simulation from a specific embodiment. Experimental verification shows that the proposed simulation method can effectively simulate debris flow confluence areas. Compared to existing methods, this method is less computationally intensive and time-consuming, making it more suitable for large-scale grid computing solutions.

[0172] A debris flow confluence area simulation system based on a multi-grid method comprises an input parameter terminal, an intersection area simulation terminal and a simulation display terminal; the input parameter terminal is used to set various initialization parameters and perform parameter assignment; the intersection area simulation terminal performs initialization settings according to the intersection area simulation parameters, and uses a debris flow confluence area simulation method based on a multi-grid method as described in any of the above items to perform debris flow confluence area simulation; the simulation display terminal performs post-processing display settings according to the post-processing display parameters, and outputs images and display variables corresponding to the post-processing display settings.

[0173] It should be noted that: this implementation uses the above method to construct a debris flow intersection area simulation system based on the multi-grid method. The system only requires the user to input some initial parameters to simulate the debris flow intersection area, without the need for human intervention in the simulation process.

[0174] The specific setting interface can be referred to as Figure 7 The input parameter terminal shown is Figure 8 The system can quickly and easily simulate the debris flow confluence area based on the input parameters to facilitate subsequent analysis and development. Figure 9 As shown (Note: the test is based on my personal thin and light notebook Lenovo Xiaoxin 14 Pro), we can clearly see that with the increase in the number of grids, the efficiency of this embodiment is significantly improved.

[0175] The above are merely embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention description and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A debris flow confluence simulation method based on a multi-grid method, characterized in that: A semi-staggered grid is established, which includes several discrete points. The dynamic equation of the intersection area is established in the semi-staggered grid, and the debris flow / main river is simulated through several marked points. During the simulation, initialization parameters are set, the positions of the marked points at each time are solved, and the parameter values of each parameter are iterated to obtain the simulation results of the debris flow intersection area. The specific steps are as follows: Step S1: Parameter initialization and assignment: Set each initialization parameter and assign the parameter value; Step S2: Discretization and construction of the grid: construct grid cells and marker points by initializing parameters and assigning corresponding parameter values, construct a semi-staggered grid, and set boundary conditions; Step S3: Time advancement process and iteration: according to the initialization parameters, time advancement is performed, and an iterative method is used in the advancement process. When the iteration converges within the error, the position of the mark point and the parameter value at the next moment are output; The position and parameter value of the marker point at the next moment are used in an iterative time advancement process to obtain the position and parameter value of the marker point at each moment; wherein the time advancement process includes calculation of intermediate speed, iteration of density and viscosity, and solution of the pressure Poisson equation; The time advancement process is iteratively performed through the following steps: Step A1 calculates the parameter corresponding to the current moment n and And read the fluid density and fluid viscosity at the current moment; Step A2: Based on the parameters obtained in step A1, the multi-grid technique is further used to iteratively solve the pressure field; Step A3: Calculate the intermediate step speed according to the preset formula based on the pressure field, where the intermediate step speed includes the intermediate step lateral speed and the intermediate step longitudinal velocity Step A4: Calculate the movement position of the marker point within the time dt according to the intermediate step speed, and update the distribution of the marker point in the semi-staggered grid; Step A5: updating the fluid density and fluid viscosity of each marking point according to the distribution of the marking points in the semi-staggered grid; Step A6: Repeat steps A2 to A5 until the density change of the marked points in each grid cell is minimal, and set the current fluid density and fluid viscosity values as the initial values of the next iteration at time n+1; wherein the minimal change is determined by whether the density change is less than a preset minimum value; Step A7: Repeat steps A1 to A6, and further iterate to obtain the position of the marking point and the parameter value at each moment; The pressure field is iteratively solved by the following formula 2: Among them, the superscript represents the time at which the parameter is located, and the subscript represents the number of the semi-staggered grid position where the parameter is located. Represents the convergence judgment factor, the parameter at each moment and All are solved by the results of the previous moment; ρ represents the fluid density at the corresponding moment and position; P represents the fluid pressure at the corresponding moment and position; The pressure field expression is solved by the multigrid method of V cycle, and the steps are as follows: Step B1: Pre-smoothing: Formula 2 is rewritten to obtain Formula 3, which is: Among them, at time n The value at each grid point is known, and the It is calculated by formula 2, which requires Substitute into formula 3 and iterate to get the pressure at the next moment The left side of Formula 3 is replaced by the coefficient matrix A h and the solution matrix u of pressure P h The product of is expressed as f on the right side. h Matrix representation, we get Formula 4: A h u h =f h The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid. 2h is a grid twice as thick as the h layer. 4h is a grid twice as thick as the 2h layer. And so on, up to Nh. N is the maximum value that n can take. The default minimum number of grids in the coarsest grid layer Nh is not less than 8. Perform 3 to 5 Gauss-Seidel iterations according to Formula 4 to obtain the solution matrix u h Estimated value of v h ; Step B2 defines: through the f of the previous grid h Calculate the residual value f of the next coarse grid 2h , and insert it into the coarse grid 2h to get the estimated value v on the coarse grid 2h The specific process is: First calculate f using formula 5 2h , the formula 5: Then perform 3 to 5 Gauss-Seidel iterations to solve Equation 6 and obtain the solution matrix u 2h Estimated value of v 2h ; Formula 6: A 2h u 2h =f 2h The process of defining the fine grid (n-1)h to the coarse grid nh is completed by performing multiple iterations and is described by Formula 7: A (n-1)h v (n-1)h (fine grid) → A nh v nh (coarse grid); Step B3: Grid coarsening: The coarse grid is solved by performing a fully weighted operation using Formula 8, and step B2 is repeated to coarsen the grid to the coarsest grid Nh at the bottom layer, and then Gauss-Seidel iteration is performed until convergence, Formula 8; The superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid. 2h is a grid twice as thick as the h layer. 4h is a grid twice as thick as the 2h layer. And so on, up to Nh. N is the maximum value that n can take. The default minimum number of grids in the coarsest grid layer Nh is not less than 8. Step B4: According to the interpolation operator formula, the approximate solution of the residual value ε is extended from the coarsest grid. nh Interpolate to the fine grid of the previous layer to get ε (n-1)h , and use the revised formula to estimate the value v of the previous layer (n-1)h Make corrections and get a new estimated value v(new) (n-1)h ; Then continue to interpolate to the finer grid (n-2)h of the next layer, repeat until you return to the finest grid h, and get the final corrected solution v(new) h =v h +ε h , at this time v(new) h That is the solution u h ; The interpolation operator formula: Among them, the residual value Calculated using the residual value calculation formula; the superscript nh indicates the thickness of the current grid layer. When the superscript is h, it is the initial finest grid, 2h is a grid twice as thick as the h layer, 4h is a grid twice as thick as the 2h layer, and so on, all the way to Nh, where N is the maximum value that n can take. The default minimum number of grids in the coarsest grid layer Nh is not less than 8; The correction formula: in(new) (n-1)h =in (n-1)h +ε (n-1)h The residual value calculation formula is: Among them, ε nh is the residual value on the n-th grid, A nh v nh It can be obtained by formula 4, is the value calculated by formula 4 for the previous grid The initial value on the nth grid after interpolation; Step S4 post-processes the display image: obtains the marker point position and parameter value, and post-processes the marker point position and parameter value according to the display settings set by the user to obtain a post-processed display image, thereby completing the simulation of the debris flow intersection area.

2. The method for simulating debris flow confluence area based on multi-grid method according to claim 1, characterized in that: The initialization parameters include intersection area simulation parameters and post-processing display parameters; the intersection area simulation parameters are used to set the initialization parameters required for debris flow intersection area simulation, including grid encryption coefficient, main river length, main river width, debris flow width, main river lateral initial velocity, main river longitudinal initial velocity, debris flow lateral initial velocity, debris flow longitudinal initial velocity, main river initial density, debris flow initial density, main river initial viscosity, debris flow initial viscosity and riverbed bottom friction coefficient; the post-processing display parameters are used to set the user's display settings, including debris flow marker outer contour line parameters, real-time cloud map parameters, post-processing display parameters and post-processing display variable parameters.

3. The method for simulating debris flow confluence area based on multi-grid method according to claim 2, characterized in that: The boundary of the semi-staggered grid consists of the main river inflow boundary, the main river outflow boundary, the river bank solid wall boundary and the debris flow inflow boundary; The semi-staggered grid is composed of several grid cell units, and each marked point is numbered by a two-dimensional coordinate axis; the left discrete point of the grid cell unit stores the velocity u with a number (i-1 / 2, j), the grid cell unit discrete point stores the right velocity u with a number (i+1 / 2, j), the upper velocity v stored by the grid cell unit discrete point is numbered (i, j-1 / 2), the lower velocity v stored by the grid cell unit discrete point is numbered (i, j+1 / 2), and the grid center discrete point stores the velocity (i, j); the main river transverse initial velocity and the main river longitudinal initial velocity are stored at the cell boundary center discrete point of the main river inflow boundary; the debris flow transverse initial velocity and the debris flow longitudinal initial velocity are stored at the cell boundary center discrete point of the debris flow inflow boundary.

4. The method for simulating debris flow confluence area based on multi-grid method according to claim 3, characterized in that: The main river inflow boundary and the main river outflow boundary adopt fixed boundary conditions, the river bank solid wall boundary adopts sliding boundary conditions, the middle inflow area of the debris flow inflow boundary adopts fixed boundary conditions, and the edge areas on both sides adopt sliding boundary conditions.

5. The method for simulating debris flow confluence area based on multi-grid method according to claim 1, characterized in that: The parameter corresponding to the current moment n and Calculated by the following formula 1: Where the superscript represents the time at which the parameter is located, the subscript represents the index of the semi-staggered grid position where the parameter is located, ρ represents the fluid density at the corresponding time and position, u represents the transverse velocity of the fluid at the corresponding time and position, v represents the longitudinal velocity of the fluid at the corresponding time and position, and Vis represents the discrete format of the viscosity term in the x / y direction.

6. The method for simulating debris flow confluence area based on multi-grid method according to claim 1, characterized in that: The intermediate step lateral speed Calculated by the following formula 12: Among them, the superscript represents the time at which the parameter is located, the subscript represents the number of the semi-staggered grid position where the parameter is located, and the parameter is the quantity obtained at this moment, ρ represents the fluid density at the corresponding moment and position, and is obtained through pressure field iteration; the longitudinal velocity of the intermediate step is Calculated by the following formula 13: Among them, the superscript represents the time at which the parameter is located, the subscript represents the number of the semi-staggered grid position where the parameter is located, and the parameter is the quantity that has been calculated at that moment, ρ represents the fluid density at the corresponding moment and position, and is obtained through pressure field iteration.

7. The method for simulating debris flow confluence area based on multi-grid method according to claim 6, characterized in that: The fluid density ρ is iterated by the following formula 14: The marking points include debris flow marking points and main river marking points, ρ1 main river initial density, ρ2 debris flow initial density, m i,j Indicates the number of main river markers in the grid cell corresponding to the coordinates, n i,j Indicates the number of debris flow markers in the grid cell corresponding to the coordinates; The fluid viscosity μ is iterated by the following formula 15: The marking points include debris flow marking points and main river marking points, μ1 initial viscosity of the main river, μ2 initial viscosity of the debris flow, m i,j Indicates the number of main river markers in the grid cell corresponding to the coordinates, n i,j Indicates the number of debris flow markers in the grid cell corresponding to the coordinates.

8. A debris flow confluence simulation system based on multi-grid method, characterized by: It includes an input parameter terminal, an intersection area simulation terminal and a simulation display terminal; the input parameter terminal is used to set various initialization parameters and perform parameter assignment; the intersection area simulation terminal performs initialization settings according to the intersection area simulation parameters, and uses a debris flow intersection area simulation method based on a multi-grid method as described in any one of claims 1 to 7 to perform debris flow intersection area simulation; the simulation display terminal performs post-processing display settings according to the post-processing display parameters, and outputs images and display variables corresponding to the post-processing display settings.

Images