Underground water solute transport prediction method based on FVM-SPH coupling algorithm

Abstract

The invention discloses a groundwater solute transport prediction method based on an FVM-SPH coupling algorithm. The method comprises the following steps: initializing system parameters; generating a grid and initializing water head information of grid nodes; listing a discrete format of a control equation for describing the underground water seepage problem and an inlet and outlet water head boundary based on a finite volume method (FVM), and carrying out iterative calculation; outputting a steady-state water head result and taking steady-state water head information as a background grid for solute transportation; arranging particles on the background grid; carrying out interpolation calculation on the speed and the diffusion coefficient of the position of the particle; listing a discrete format and a boundary condition applying mode of a control equation for describing a solute transport problem based on a smoothed particle hydrodynamics (SPH) numerical method, and performing iterative calculation; and outputting a result. The method can accurately simulate and predict the change of the concentration in the solute transport process of the coupled groundwater flow field, can effectively adapt to the nonlinear characteristic of hydraulic parameters in groundwater, and can provide instructive suggestions for salt lake brine pumping work.

Description

Technical Field

[0001] The present invention belongs to the technical field of computational fluid dynamics simulation research, and specifically relates to a method for solving and predicting groundwater solute transport problems by using a coupling algorithm of a finite volume method (FVM) and a smoothed particle hydrodynamics method (SPH). Background Art

[0002] In addition to the visible transport processes above the surface, solute transport in salt lake mining areas also involves a significant portion of solute transport in groundwater. However, long-term brine mining activities have severely damaged the hydrogeological structure of the underground mining area. For example, changes in the porosity and permeability of the underground rocks have caused the groundwater flow structure to continuously change, severely impacting the solute transport process. Therefore, when studying the solute transport process in salt lake mining areas, the impact of the groundwater seepage field on the solute transport process cannot be ignored. It is necessary to clarify the changing patterns of the groundwater seepage field and understand the groundwater flow structure in order to accurately simulate and predict the solute transport process in salt lake mining areas.

[0003] Smoothed Particle Hydrodynamics (SPH) is a Lagrangian-based numerical method and interpolation technique. First introduced in astrophysics, it is now widely used in numerical simulations of fluid mechanics. The key concept behind SPH is to discretize the fluid into a series of particle points. These particles are classified as inflow particles, fluid particles, and outflow particles. Each particle typically has five physical properties: position, velocity, density, mass, and concentration. Determining the new concentration of a fluid particle at the next time step based on its original concentration and diffusion coefficient is crucial to ensuring the accuracy of simulation results.

[0004] The finite volume method (FVM) is a numerical method based on Euler's ideas. It was originally a grid method proposed to improve the finite difference method (FDM). At the same time, it borrowed some advantages of the finite element method (FEM). It is currently mainly used to solve fluid mechanics and thermodynamics problems. As a traditional Euler method, FVM discretizes the computational domain into grid cells. Each grid cell is called a control volume. The integral region of the control equation is these independent control volumes, which represents the flux conservation within the control volume. This also ensures that each physical quantity satisfies the integral form of conservation relations within each control volume. Then, FVM discretizes the integral form of conservation control equations through a specific interpolation algorithm, and finally converts it into a discrete form of conservation equations.

[0005] Previous methods [1] have successfully verified the accuracy of SPH in solving the convection-diffusion equation. However, they do not consider the effect of the flow field on the concentration field transport. The seepage field in salt lake groundwater exhibits heterogeneity, so the effect of the groundwater seepage field on solute transport must be considered. Currently popular software for solving groundwater flow fields, such as MODFLOW, all use traditional grid methods with high accuracy and precision.

[0006] References:

[0007] [1] Hou Qingzhi, Chen Shaokang, Miao Chunfu, et al. Pollutant transport prediction method based on smoothed particle hydrodynamics: CN202210560376.2[P]. CN202210560376.2[2024-04-12]. Summary of the Invention

[0008] The purpose of the present invention is to overcome the deficiencies in the prior art and provide a groundwater solute transport prediction method based on the FVM-SPH coupling algorithm, so as to predict the changes in solute concentration during the flow process of the groundwater seepage field, obtain a groundwater solute transport process that conforms to reality, play a guiding role in brine extraction work in salt lakes, and have the characteristics of small error and high computational efficiency.

[0009] The purpose of the present invention is achieved through the following technical solutions:

[0010] This will be completed after the claim is determined

[0011] Compared with the prior art, the technical solution of the present invention has the following beneficial effects:

[0012] 1. The present invention combines the finite volume method (FVM) to solve the groundwater seepage field with the smoothed particle hydrodynamics (SPH) to solve the solute transport process. At the same time, it takes advantage of the high precision of FVM in heterogeneous seepage problems and the advantages of SPH in tracking fluid interfaces and concentration changes, significantly improving the overall simulation accuracy and avoiding the error accumulation of a single method under complex geological conditions. FVM strictly handles the conservation of flux within the grid unit, ensuring the accuracy of the head and velocity; SPH is based on Lagrangian particle motion, does not require gridding of the concentration field, and naturally adapts to strong nonlinear seepage areas to achieve high-fidelity coupled simulation.

[0013] 2. This method constructs a uniform quadrilateral grid within the fluid domain and iteratively solves the seepage equation until steady state. This method captures heterogeneous features such as formation porosity and permeability with high resolution, providing a stable and accurate background hydraulic head field, laying a solid foundation for subsequent solute transport. The gridded FVM ensures flux conservation within each unit cell, and the fine grid improves spatial resolution, accurately reflecting the distribution of local hydraulic gradients.

[0014] 3. This paper uses the cubic Hermite interpolation method to calculate the velocity and diffusion coefficient at the particle position. This significantly reduces the interpolation error of the velocity and diffusion coefficient, ensuring high accuracy of the velocity field and diffusion coefficient input in the SPH calculation, thereby improving the accuracy of solute transport simulation.

[0015] 4. The present invention arranges three types of particles on the background grid: fluid particles (type 1), inflow particles (type 2), and outflow particles (type 3), and dynamically allocates or converts them at the boundary. This method can automatically maintain boundary conditions to ensure flux conservation throughout the simulation; continuous replenishment at the inlet and automatic removal at the outlet achieve a truly seamless boundary flow process. Specifically, through particle type identification and position determination, inflow particles are converted into fluid particles once they enter the fluid domain, maintaining unchanged concentration and mass; outflow particles are converted into outflow particles and removed from outside the domain to ensure precise balance of the system.

[0016] 5. Using a cubic spline smoothing kernel function and its gradient, the solute transport equation is discretized and the concentration is updated iteratively. This method of tracking solutes significantly reduces numerical diffusion and is adaptable to complex seepage paths and heterogeneous media.

[0017] 6. A small time step (0.001 d) was used for the FVM-based steady-state head calculation, while a larger time step (0.01 d) was used for the SPH transport calculation. This balances the flow field's stable convergence with the computational efficiency of the transport process, shortening the overall simulation time while ensuring the numerical stability and convergence of each subprocess.

[0018] 7. The FVM-SPH coupling method naturally supports nonlinear inputs for parameters such as spatial permeability, porosity, and diffusion coefficient. Even in the face of severely damaged hydrogeological structures during salt lake brine mining, it can still accurately reflect the nonlinear characteristics of the medium, ensuring the reliability of simulation results.

[0019] 8. After the simulation is complete, multi-dimensional information such as the position, velocity, diffusion coefficient, and concentration of all particles at the final moment is output. This provides rich and intuitive data support for analyzing solute migration patterns in the salt lake area and optimizing brine extraction plans, facilitating subsequent decision-making.

[0020] In summary, the present invention achieves high-precision and high-efficiency simulation of heterogeneous underground water solute transport through the coupling of refined grid FVM and high-fidelity SPH, supplemented by high-order interpolation, dynamic particle management and step-by-step time advancement, with small error and fast convergence, which can provide powerful guidance for brine mining and environmental protection in salt lake mining areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 It is a physical model for groundwater solute transport problems.

[0022] Figure 2Schematic diagram of the process of the present invention.

[0023] Figure 3 Schematic diagram of particle distribution in the background grid at a certain moment.

[0024] Figure 4 Schematic diagram of particle flux calculation within the calculation domain in the present invention.

[0025] Figure 5 Schematic diagram of particle concentration calculation inside the calculation domain in the present invention.

[0026] Figure 6 Comparison between the results of the FVM prediction of the steady-state hydraulic head change of the groundwater seepage field in the present invention and the reference solution under high-resolution grid on the y=50m section.

[0027] Figure 7 This is the result of the change in x-direction velocity on the y=50m section.

[0028] Figure 8 This is the change in the diffusion coefficient in the x direction on the y=50m section.

[0029] Figure 9 Comparison of the solute transport concentration change prediction results of the present invention using FVM-SPH on the y=50m section with the reference solution under high-resolution grid. DETAILED DESCRIPTION

[0030] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0031] This embodiment illustrates the specific solution process of the method of the present invention through a specific example. The problem can be described as follows: in a phreatic aquifer, the left and right hydraulic boundaries are constant Dirichlet boundaries, and the upper and lower hydraulic boundaries are free outflow boundaries. A certain mass of solute is introduced into the aquifer from a well at position P. The physical model of this problem is as follows: Figure 1 As shown in the figure, the computational domain is a rectangular area ABCD. The area enclosed by A′B′C′D′ is the area where the fluid is located, referred to as the "fluid domain". The length of the inlet boundary A′D′ is 100m, the length of the outlet boundary B′C′ is 100m, and the width of the fluid domain is 500m. In this problem, the permeability coefficient K = 100m / d, the longitudinal diffusion α L =1.0m, lateral dispersion α T =0.1m, the initial solute mass is M=4000g, and the coordinates of point P are (50,50). The initial head distribution, concentration distribution, head boundary, and concentration boundary conditions are as follows:

[0032]

[0033] The present invention will be described in further detail below with reference to the accompanying drawings. Figure 2 A schematic flow chart of the method of the present invention is given.

[0034] include:

[0035] Step 1: Initialize the system. The vertex coordinates of the computational domain are A(-10,100), B(510,100), C(510,0), and D(-10,0). The vertex coordinates of the inflow boundary are A′(0,100), D′(0,0), and the vertex coordinates of the outflow boundary are B′(500,100), C′(500,0). The grid size and fluid particle spacing are both Δx=Δy=5.0m, the smoothing length h=1.33Δx=1.33Δy, the total simulation time is 1 year, and the equal time step is used. The time step is set to Δt=0.001d when FVM solves the seepage equation, and the time step is set to Δt=0.01d when SPH solves the solute transport equation. The permeability coefficient K=100m / d, and the longitudinal diffusion α L =1.0m, lateral dispersion α T =0.1m.

[0036] Step 2: Generate the grid and initialize the water head information of the grid nodes. The initial water head distribution is obtained according to formula (20). The fluid domain A′B′C′D′ is evenly divided into quadrilateral grids with a grid size of Δx = Δy = 5.0m.

[0037] Step 3: Based on the finite volume method (FVM), the discrete format of the governing equations describing the groundwater seepage problem and the inlet and outlet head boundaries are listed, and iterative calculations are performed:

[0038] (301) The governing equation describing the groundwater seepage problem is discretized using FVM in the following format:

[0039]

[0040]

[0041] Among them, μ d is the gravity water supply of the phreatic aquifer, K is the permeability coefficient, H i'+1 / 2,j' 、H i'-1 / 2,j' 、H i',j'+1 / 2 and H i',j'-1 / 2 They are intermediate variables in the calculation process, the superscript n represents the current moment, n+1 represents the next moment, Δt is the time step, H i',j' 、H i'-1,j' 、H i'+1,j' 、H i',j'-1 、H i',j'+1are the water heads at each grid node, is the hydraulic head at the grid node (i', j') at time n+1, Δx and Δy are the grid sizes in the x and y directions, respectively;

[0042] (302) Inlet and outlet head boundary settings: The inlet head boundary is set to a constant head boundary h = 10m, and the outlet head boundary is set to a constant head boundary h = 6m;

[0043] (303) Iterative calculation:

[0044] The current time starts from 0 until the head distribution reaches the steady-state head; the specific calculation process is as follows:

[0045] (3031) Determine whether the current head distribution reaches the steady-state head by comparing the current head distribution with the previous head distribution; if yes, end the calculation; if not, recalculate the process (3032) to (3033);

[0046] (3032) For the internal grid in the computational domain, calculate the H of each node according to equations (1) to (4). i'+1 / 2,j' 、H i'-1 / 2,j' 、H i',j'+1 / 2 and H i',j'-1 / 2 ;

[0047] (3033) Calculate the new hydraulic head at the internal grid nodes in the computational domain after 0.001s according to equation (5);

[0048] (3034) Current time = current time + 1 time step; repeat the calculation process (3031);

[0049] Step 4: Output the steady-state hydraulic head result: When the simulation time reaches the steady-state hydraulic head, the cycle ends and the steady-state hydraulic head information of the grid node at this moment is output as the background grid for solute transport;

[0050] Step 5: Arrange particles on the output background grid, including

[0051] 501) Generate fluid particles, evenly distribute the fluid particles in the fluid domain, and add initial fluid particle information. Evenly distribute the fluid particles in the fluid domain A'B'C'D' with the fluid particle spacing as the unit. The fluid particle mass m = 0.1 kg, density ρ = 1000 kg / m 3 , the fluid particle type identification number is 1. Among them, the fluid particle concentration is obtained according to formula (21), and the concentrations of other fluid particles are 0.

[0052] 502) Generate inflow particles, evenly distribute them at the inflow boundary, and add initial inflow particle information. Arrange two layers of inflow particles outward from the inflow boundary A′D′ with the particle spacing as the unit. Inflow mass m = 0.1 kg, density ρ = 1000 kg / m 3 , the inflow particle type identification number is 2, and the inflow particle concentration is 0.

[0053] 503) Generate outflow particles, evenly distribute outflow particles at the outflow boundary, and add initial outflow particle information. Arrange two layers of inflow particles outward at the outflow boundary B′C′ with the particle spacing as the unit. The outflow particles have a mass m = 0.1 kg and a density ρ = 1000 kg / m 3 , the outflow particle type identification number is 3, and the outflow particle concentration is 0.

[0054] Step 6: Interpolate and calculate the velocity and diffusion coefficient at the particle's location. Based on the output steady-state hydraulic head information for the grid nodes, Darcy's law and the diffusion coefficient calculation formula are used to calculate the x- and y-direction velocities and diffusion coefficient of the grid nodes. Interpolation is then used to calculate the x- and y-direction velocities and diffusion coefficient at the particle's location. In this embodiment, the high-precision Hermite Cubic Interpolation in Space Line (HCSL) method is used.

[0055] (601) Darcy’s law can be expressed as:

[0056]

[0057] Where I represents the hydraulic gradient, is the total flow rate, v x 、v y are the flow velocities in the x and y directions, respectively;

[0058] (602) Diffusion coefficient calculation:

[0059]

[0060] Among them, D x 、D y are the diffusion coefficients in the x and y directions, α L and α T are the longitudinal and transverse diffusivities, respectively.

[0061] Step 7: Based on the smoothed particle hydrodynamics (SPH) numerical method, the discrete format of the governing equations describing the solute transport problem and the boundary condition application method are listed and iteratively calculated:

[0062] (701) Discrete format of the solute transport governing equation:

[0063]

[0064]

[0065] Among them, W(x i -x j ,h) is called the smooth kernel function, is the gradient of the smooth kernel function, h = 1.33dx is called the smooth length, where dx is the initial regular distribution of particle spacing, the smooth kernel function selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |xi-xj| / h, then the smooth kernel function W(x i -x j ,h) is:

[0066]

[0067] Among them, i represents the current fluid particle, its adjacent 2h range is its support domain, j is one of the fluid particles in the support domain of the current fluid particle, x i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain, respectively; N is the number of all fluid particles in the support domain of the current particle i; m j is the mass of fluid particle j, ρ i and ρ j are the densities of the current fluid particle i and the fluid particle j in its support domain, are the intermediate variables in the calculation process, and are the concentrations of the current fluid particle i and the fluid particle j in its support domain, D i is the diffusion coefficient at the current position of fluid particle i;

[0068]

[0069] I y(-) (j) is the characteristic function, and its expression is as follows:

[0070]

[0071] x i and y i They are the position components of the current fluid particle i in the x and y directions, x j and y j represent the position components of the fluid particle j in the support domain in the x-direction and y-direction respectively;

[0072] In formulas (8) to (11),

[0073]

[0074] in, and Represent the gradient of the smooth kernel function of the current fluid particle i in the x-direction and y-direction respectively;

[0075] (702) Boundary condition application method:

[0076] A free outflow boundary condition is applied to the outlet boundary. When a fluid particle flows into the outlet domain, its fluid particle type identification number 1 is changed to outflow particle type identification number 3, becoming a new outflow particle. Other properties remain unchanged.

[0077] The velocity condition is applied to the inlet boundary: when the inflow particle flows into the fluid domain, its inflow particle type identification number 2 is changed to fluid particle type identification number 1, becoming a new fluid particle, and other properties remain unchanged;

[0078] (703) Iterative calculation:

[0079] The current time starts from 0, and the total simulation time has been given in step 1; the specific calculation process is as follows:

[0080] (7031) Determine whether the current time has reached 1 year; if yes, end the calculation; if not, recalculate the process (7032) to (7035);

[0081] (7032) All particles are initially assigned coordinates. The new position of each particle is equal to the original position plus the velocity at the particle's location multiplied by the time step of 0.01d. For inflow particles that enter the fluid domain after the position is updated, the type identification number in their particle properties is changed to 1, and they become new fluid particles. The inflow particles are supplemented at the corresponding positions on the outermost layer of the inlet boundary. For fluid particles that flow out of the fluid domain after the position is updated, the type identification number in their particle properties is changed to 3, and they become new outflow particles.

[0082] (7033) For all particles distributed in the background grid, the HCSL interpolation algorithm is used to calculate the flow velocity and diffusion coefficient at the location of each particle according to the coordinate position of each particle;

[0083] (7034) Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other particles searched in the support domain by the current fluid particle to calculate the current fluid particle according to equations (8) to (9). and

[0084] (7035) Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other particles searched by the current fluid particle in the support domain to calculate the new concentration of the current fluid particle after 0.01d according to equations (10) to (12);

[0085] (7036) Current time = current time + 1 time step; repeat calculation process (7031).

[0086] Step 8: Output the results. After the total time reaches 1 year, the loop ends and the calculation stops. The position, velocity, diffusion coefficient, and concentration information of all particles at the current moment are output as the final calculation results.

[0087] Furthermore, the system initialization in step 1 is as follows:

[0088] 101) The grid size is 5 m in both x and y directions;

[0089] 102) The distance between fluid particles, the distance between inflow particles and the distance between outflow particles are all 5m;

[0090] 103) The inlet and outlet boundaries are both located outside the fluid domain and are fixed areas. The coordinates of the inlet and outlet vertices do not change.

[0091] Furthermore, when arranging particles on the background grid in step five: the number of layers of inflow particles is two, and the number of layers of outflow particles is two;

[0092] Furthermore, the calculation process (7032) in step seven supplements the inflow particle at the corresponding position on the outermost layer of the inlet boundary: when the abscissa of the inflow particle position is greater than or equal to 0, a new fluid particle is generated in the fluid domain. This new fluid particle inherits the position, velocity, mass, density, concentration, and other properties of the original inflow particle. Then, the position of the inflow particle with an abscissa greater than or equal to 0 is changed to the abscissa minus 10m.

[0093] Figure 5 For the section y = 50m, the results of the FVM prediction of the steady-state head change of the groundwater seepage field in this invention are consistent with the reference solution under high-resolution grid.

[0094] Figure 6 and Figure 7 The results of the changes in velocity and diffusion coefficient in the x-direction on the y=50m section show that both the velocity and diffusion coefficient change with the coordinate position, showing nonlinearity.

[0095] Figure 8 For the section y=50m, the results of the FVM-SPH prediction of solute transport concentration changes in the present invention are consistent with the reference solution under high-resolution grids.

[0096] Based on the same inventive concept, the embodiments of the present application also provide a groundwater solute transport prediction device based on the FVM-SPH coupling algorithm, which can be used to implement the method described in the above embodiments, as described in the following embodiments. Since the principle of solving the problem by the groundwater solute transport prediction device is similar to that of the groundwater solute transport prediction method, the implementation of the groundwater solute transport prediction device can refer to the implementation of the groundwater solute transport prediction method, and the repeated parts will not be repeated. As used below, the term "unit" or "module" can be a combination of software and / or hardware that implements a predetermined function. Although the system described in the following embodiments is preferably implemented in software, the implementation of hardware, or a combination of software and hardware, is also possible and conceived.

[0097] An embodiment of the present invention provides a groundwater solute transport prediction device based on the FVM-SPH coupling algorithm, comprising:

[0098] Initialization unit, used to initialize system parameter information, including the coordinates of the vertex of the computational domain, the vertex coordinates of the inflow boundary, the vertex coordinates of the outflow boundary, the grid size, the particle spacing, the smoothing length, the permeability coefficient, the lateral dispersion, the longitudinal dispersion, the total simulation time and the time step;

[0099] The grid generation unit is used to initialize the water head information of all nodes in the grid and evenly divide the fluid domain into several quadrilateral grids with the grid size as the size;

[0100] FVM calculation unit, used to list the discrete format and inlet and outlet head boundaries of the control equation describing the groundwater seepage problem based on the finite volume method FVM, and perform iterative calculations;

[0101] The background grid generation unit is used to end the cycle when the simulation time reaches the steady-state hydraulic head, and output the steady-state hydraulic head information of the grid node at the current moment as the background grid for solute transport;

[0102] The particle placement unit is used to place particles on the background grid. The particle placement unit includes:

[0103] Generate fluid particles module, used to evenly distribute fluid particles in the fluid domain and add initial fluid particle information; the initial fluid particle information includes fluid particle density, fluid particle mass, fluid particle type identification number 1, and fluid particle concentration;

[0104] Generate inflow particle module, which is used to evenly distribute inflow particles in the fluid domain and add initial inflow particle information; the initial inflow particle information includes inflow particle density, inflow particle mass, inflow particle type identification number 2, and inflow particle concentration;

[0105] Generate outflow particle module, which is used to evenly distribute outflow particles in the fluid domain and add initial outflow particle information; the initialization particle information includes outflow particle density, outflow particle mass, outflow particle type identification number 3, and outflow particle concentration.

[0106] The interpolation calculation unit is used to interpolate and calculate the velocity and diffusion coefficient at the location of the particle. Based on the output steady-state water head information of the grid node, Darcy's law and the diffusion coefficient calculation formula are used to calculate the x-direction velocity, y-direction velocity and diffusion coefficient of the grid node. Then, the x-direction velocity, y-direction velocity and diffusion coefficient at the location of the particle are calculated using the interpolation method. The cubic Hermite interpolation algorithm HCSL on the spatial line is used.

[0107] The SPH calculation unit is used to list the discrete format of the control equations describing the solute transport problem and the boundary condition application method based on the smoothed particle hydrodynamics (SPH) numerical method and perform iterative calculations;

[0108] The output unit is used to output the position, velocity, concentration of all fluid particles, the position, velocity, concentration of all inflow particles, and the position, velocity and concentration of outflow particles at the end of the simulation cycle as the final result.

[0109] Preferably, the embodiments of the present application further provide a specific implementation of an electronic device capable of implementing all steps of the groundwater solute transport prediction method based on the FVM-SPH coupling algorithm in the above embodiment, and the electronic device specifically includes the following contents:

[0110] Processor, memory, communications interface, and bus;

[0111] Among them, the processor, memory, and communication interface communicate with each other through the bus; the communication interface is used to realize information transmission between related devices such as server-side devices, metering devices, and user-side devices.

[0112] The processor is used to call the computer program in the memory, and when the processor executes the computer program, it implements all the steps in the groundwater solute transport prediction method based on the FVM-SPH coupling algorithm in the above embodiment.

[0113] An embodiment of the present application also provides a computer-readable storage medium capable of implementing all steps of the groundwater solute transport prediction method based on the FVM-SPH coupling algorithm in the above-mentioned embodiment. The computer-readable storage medium stores a computer program, which, when executed by a processor, implements all steps of the groundwater solute transport prediction method based on the FVM-SPH coupling algorithm in the above-mentioned embodiment.

[0114] The various embodiments in this specification are described in a progressive manner. Similar portions between the various embodiments can be referenced to each other. Each embodiment focuses on the differences between the other embodiments. In particular, the hardware + program embodiments are generally similar to the method embodiments, so their description is relatively simple. For relevant portions, refer to the description of the method embodiments.

[0115] The foregoing description of this specification describes specific embodiments. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recited in the claims can be performed in an order different from that described in the embodiments and still achieve the desired results. Furthermore, the processes depicted in the accompanying drawings do not necessarily require the specific order shown or the sequential order to achieve the desired results. In certain embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0116] Although the present application provides method operation steps such as embodiments or flowcharts, more or fewer operation steps may be included based on conventional or non-creative work. The order of steps listed in the embodiments is only one way of executing the steps among many steps and does not represent the only execution order. When an actual device or client product is executed, it can be executed in the order shown in the embodiments or the drawings or in parallel (for example, in a parallel processor or multi-threaded processing environment).

[0117] It will be understood by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Thus, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0118] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0119] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 The steps for the function specified in one or more boxes.

[0120] The present invention is not limited to the embodiments described above. The above description of the specific embodiments is intended to describe and illustrate the technical solutions of the present invention. The above specific embodiments are merely illustrative and not restrictive. Without departing from the scope of the present invention and the scope of protection of the claims, those skilled in the art may make various specific modifications based on the teachings of the present invention, all of which fall within the scope of protection of the present invention.

Claims

1. A groundwater solute transport prediction method based on FVM-SPH coupling algorithm, characterized in that: include: Step 1: Initialize the system parameters. The parameter information includes the coordinates of the vertex of the computational domain, the vertex coordinates of the inflow boundary, the vertex coordinates of the outflow boundary, the grid size, the particle spacing, the smoothing length, the permeability coefficient, the lateral dispersion, the longitudinal dispersion, the total simulation time, and the time step. Step 2: Generate a grid and initialize the water head information of the grid nodes; divide the fluid domain into several quadrilateral grids evenly with the grid size as the size; Step 3: Based on the finite volume method (FVM), the discrete format of the governing equations describing the groundwater seepage problem and the inlet and outlet head boundaries are listed, and iterative calculations are performed. Step 4: Output the steady-state hydraulic head result: When the simulation time reaches the steady-state hydraulic head, the cycle ends and the steady-state hydraulic head information of the grid node at the current moment is output as the background grid for solute transport; Step 5: Arrange particles on the output background grid, including: 501) Generate fluid particles, evenly distribute the fluid particles in the fluid domain, and add initial fluid particle information; the initial fluid particle information includes fluid particle density, fluid particle mass, fluid particle type identification number 1, and fluid particle concentration; 502) Generate inflow particles, evenly distribute the inflow particles in the fluid domain, and add initial inflow particle information; the initial inflow particle information includes inflow particle density, inflow particle mass, inflow particle type identification number 2, and inflow particle concentration; 503) Generate outflow particles, evenly distribute the outflow particles in the fluid domain, and add initial outflow particle information; the initialization particle information includes outflow particle density, outflow particle mass, outflow particle type identification number 3, and outflow particle concentration; Step 6: Interpolate and calculate the velocity and diffusion coefficient at the particle's location. Based on the output steady-state water head information of the grid node, Darcy's law and the diffusion coefficient calculation formula are used to calculate the x-direction velocity, y-direction velocity, and diffusion coefficient of the grid node. Then, the x-direction velocity, y-direction velocity, and diffusion coefficient at the particle's location are calculated using the interpolation method. The cubic Hermite interpolation algorithm HCSL on the spatial line is used. Step 7: Based on the smoothed particle hydrodynamics (SPH) numerical method, the discrete format of the governing equations describing the solute transport problem and the boundary condition application method are listed and iterative calculations are performed; Step 8. Output results: When the total simulation time is up, the loop ends and the position, velocity, and concentration of all fluid particles at the corresponding moment, the position, velocity, and concentration of all inflow particles, and the position, velocity, and concentration of outflow particles are output as the final results.

2. The method for predicting groundwater solute transport based on the FVM-SPH coupling algorithm according to claim 1, characterized in that: Step 3 is as follows: (301) The governing equation describing the groundwater seepage problem is discretized using FVM in the following format: Among them, μ d is the gravity water supply of the phreatic aquifer, K is the permeability coefficient, H i'+1 / 2,j' 、H i'-1 / 2,j' 、H i',j'+1 / 2 and H i',j'-1 / 2 They are intermediate variables in the calculation process, the superscript n represents the current moment, n+1 represents the next moment, Δt is the time step, H i',j' 、H i'-1,j' 、H i'+1,j' 、H i',j'-1 、H i',j'+1 are the water heads at each grid node, is the hydraulic head at the grid node (i', j') at time n+1, Δx and Δy are the grid sizes in the x and y directions, respectively; (302) Inlet and outlet head boundary setting: The head of the inlet and outlet head boundaries are all subject to a constant head boundary; (303) Iterative calculation: The current time starts from 0 until the head distribution reaches the steady-state head; the specific calculation process is as follows: (3031) By comparing the head distribution at the current moment with the head distribution at the previous moment, it is determined whether the head distribution at the current moment reaches the steady-state head; if yes, the calculation is terminated; if not, the calculation process (3032) to (3033) is repeated; (3032) For the internal grid in the computational domain, calculate the H of each node according to equations (1) to (4). i'+1 / 2,j' 、H i'-1 / 2,j' 、H i',j'+1 / 2 and H i',j'-1 / 2 ; (3033) Calculate the new hydraulic head of the internal grid nodes in the computational domain after one time step according to equation (5); (3034) Current time = current time + 1 time step; repeat calculation process (3031).

3. The method for predicting groundwater solute transport based on the FVM-SPH coupling algorithm according to claim 1, characterized in that: Step 6 is as follows: 601) Darcy's law is expressed as: Where K is the permeability coefficient, I represents the hydraulic gradient, is the total flow rate, v x 、v y are the flow velocities in the x and y directions, respectively; 602) Diffusion coefficient calculation: Among them, D x 、D y are the diffusion coefficients in the x and y directions, α L and α T are the longitudinal and transverse diffusivities, respectively.

4. The method for predicting groundwater solute transport based on the FVM-SPH coupling algorithm according to claim 1, characterized in that: Step seven is as follows: (701) The governing equations describing the solute transport problem are discretized using the SPH numerical method in the following discretization format: Among them, W(x i -x j ,h) is called the smooth kernel function, is the gradient of the smooth kernel function, h = 1.33dx is called the smooth length, where dx is the initial regular distribution of particle spacing, the smooth kernel function selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |xi-xj| / h, then the smooth kernel function W(x i -x j ,h) is: Among them, i represents the current fluid particle, its adjacent 2h range is its support domain, j is one of the fluid particles in the support domain of the current fluid particle, x i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain, respectively; N is the number of all fluid particles in the support domain of the current particle i; m j is the mass of fluid particle j, ρ i and ρ j are the densities of the current fluid particle i and the fluid particle j in its support domain, are intermediate variables in the calculation process, C i and C j are the concentrations of the current fluid particle i and the fluid particle j in its support domain, D i is the diffusion coefficient at the current position of fluid particle i; In formulas (8) to (11), I x(+) (j), I x(-) (j), I y(+) (j), I y(-) (j) is the characteristic function, and its expression is as follows: x i and y i They are the position components of the current fluid particle i in the x and y directions, x j and y j represent the position components of the fluid particle j in the support domain in the x-direction and y-direction respectively; In formulas (8) to (11), in, and Represent the gradient of the smooth kernel function of the current fluid particle i in the x-direction and y-direction respectively; (702) Boundary condition application method: A free outflow boundary condition is applied to the outlet boundary. When a fluid particle flows into the outlet domain, its fluid particle type identification number 1 is changed to outflow particle type identification number 3, becoming a new outflow particle. Other properties remain unchanged. The velocity condition is applied to the inlet boundary: when the inflow particle flows into the fluid domain, its inflow particle type identification number 2 is changed to fluid particle type identification number 1, becoming a new fluid particle, and other properties remain unchanged; (703) Iterative calculation: The current time starts at 0, and the total simulation time has been given in step 1. The specific calculation process is as follows: (7031) Determine whether the current moment has reached the simulation time; if yes, end the calculation; if not, recalculate the process (7032) to (7035); (7032) All particles are initially assigned coordinates. The new position of each particle is equal to the original position plus the velocity at the particle's location multiplied by the time step. For inflow particles that enter the fluid domain after the position is updated, the type identification number in their particle properties is changed, and they become new fluid particles. The inflow particles are supplemented at the corresponding positions on the outermost layer of the inlet boundary. For fluid particles that flow out of the fluid domain after the position is updated, the type identification number in their particle properties is changed, and they become new outflow particles. (7033) For all particles distributed in the background grid, the HCSL interpolation algorithm is used to calculate the flow velocity and diffusion coefficient at the location of each particle according to the coordinate position of each particle; (7034) Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other particles searched in the support domain by the current fluid particle to calculate the current fluid particle according to equations (8) to (9). and (7035) Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other particles searched by the current fluid particle in the support domain to calculate the new concentration of the current fluid particle after one time step according to equations (10) to (12); (7036) Current time = current time + 1 time step; repeat calculation process (7031).

5. The method for predicting groundwater solute transport based on the FVM-SPH coupling algorithm according to claim 4, characterized in that: The method for supplementing the incoming particles at the corresponding position of the outermost layer of the inlet boundary described in the calculation process (7032) is as follows: when the incoming particle flows into the fluid domain, a new fluid particle is generated in the fluid domain. The new fluid particle inherits the position and concentration of the original incoming particle, and the velocity and diffusion coefficient are calculated using the HCSL interpolation algorithm in step (7033); then, the coordinate value of the incoming particle flowing into the fluid domain in the normal direction of the inlet boundary is changed to the coordinate value in the normal direction of the inlet boundary minus the length of the inlet boundary, and a new incoming particle is generated.

6. The method for predicting groundwater solute transport based on the FVM-SPH coupling algorithm according to claim 1, characterized in that: In the process of arranging particles described in step 5, the concentrations of inflow particles and outflow particles need to be set to 0.

7. A groundwater solute transport prediction device based on FVM-SPH coupling algorithm, characterized in that: include: Initialization unit, used to initialize system parameter information, including the coordinates of the vertex of the computational domain, the vertex coordinates of the inflow boundary, the vertex coordinates of the outflow boundary, the grid size, the particle spacing, the smoothing length, the permeability coefficient, the lateral dispersion, the longitudinal dispersion, the total simulation time and the time step; The grid generation unit is used to initialize the water head information of all nodes in the grid and evenly divide the fluid domain into several quadrilateral grids with the grid size as the size; FVM calculation unit, used to list the discrete format and inlet and outlet head boundaries of the control equation describing the groundwater seepage problem based on the finite volume method FVM, and perform iterative calculations; The background grid generation unit is used to end the cycle when the simulation time reaches the steady-state hydraulic head, and output the steady-state hydraulic head information of the grid node at the current moment as the background grid for solute transport; Arrange particle unit, used to arrange particles on the background grid: The interpolation calculation unit is used to interpolate and calculate the velocity and diffusion coefficient at the location of the particle. Based on the output steady-state water head information of the grid node, Darcy's law and the diffusion coefficient calculation formula are used to calculate the x-direction velocity, y-direction velocity and diffusion coefficient of the grid node. Then, the x-direction velocity, y-direction velocity and diffusion coefficient at the location of the particle are calculated using the interpolation method. The cubic Hermite interpolation algorithm HCSL on the spatial line is used. The SPH calculation unit is used to list the discrete format of the control equations describing the solute transport problem and the boundary condition application method based on the smoothed particle hydrodynamics (SPH) numerical method and perform iterative calculations; The output unit is used to output the position, velocity, concentration of all fluid particles, the position, velocity, concentration of all inflow particles, and the position, velocity and concentration of outflow particles at the end of the simulation cycle as the final result.

8. The groundwater solute transport prediction device based on the FVM-SPH coupling algorithm according to claim 7, characterized in that: The particle arrangement unit includes: Generate fluid particles module, used to evenly distribute fluid particles in the fluid domain and add initial fluid particle information; the initial fluid particle information includes fluid particle density, fluid particle mass, fluid particle type identification number 1, and fluid particle concentration; Generate inflow particle module, which is used to evenly distribute inflow particles in the fluid domain and add initial inflow particle information; the initial inflow particle information includes inflow particle density, inflow particle mass, inflow particle type identification number 2, and inflow particle concentration; Generate outflow particle module, which is used to evenly distribute outflow particles in the fluid domain and add initial outflow particle information; the initialization particle information includes outflow particle density, outflow particle mass, outflow particle type identification number 3, and outflow particle concentration.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the steps of the groundwater solute transport prediction method based on the FVM-SPH coupling algorithm according to any one of claims 1 to 6 are implemented.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the groundwater solute transport prediction method based on the FVM-SPH coupling algorithm described in any one of claims 1 to 6 are implemented.

Images