Pore solute migration numerical simulation method based on discrete pore density flow method

Abstract

A particle pore-solute migration adsorption equilibrium efficient solving method based on a discrete element pore density flow method comprises the following steps: 1) constructing a porous medium skeleton based on a discrete element method and a pore network model, and establishing a topological mapping relation between a pore space and adsorbent particles; based on a solid particle accumulation skeleton and a pore-pore throat channel thought of a pore network model, performing matrix storage on pore concentration and particle adsorption capacity; 2) through decoupling of a solute migration process, (1) constructing a solute convection term based on a pore network model so as to describe solute convection guided by fluid movement; (2) introducing a solute diffusion term based on a Fick second law; and 3) based on an isothermal adsorption equation, adopting a quasi-kinetic model for driving, evolving the adsorption equilibrium state of each central pore and the associated particle system by the system, representing the solid-liquid distribution rule of the solute in the complex porous medium, and adopting explicit time step iterative calculation to realize particle pore-solute adsorption equilibrium.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of discrete element porous medium seepage and solute transport simulation, and particularly relates to a numerical simulation method for particle pore-solute adsorption equilibrium based on a discrete element pore density flow method. BACKGROUND

[0002] With the rapid development of China's economy, industrial and agricultural activities intensify, and farmland soil and industrial land pollution problems are increasingly prominent. By studying the interfacial equilibrium behavior of pollutants in porous media, it is helpful to reveal its environmental chemical behavior, predict its morphological transformation and migration trend in the soil-plant-water system, provide a scientific basis for risk assessment and pollution prevention, and has important significance for protecting agricultural and ecological safety.

[0003] In the process of studying solute migration, both convection and diffusion transport mechanisms are considered, and the pore solute concentration is updated through convection-diffusion coupling. Adsorption behavior is a major process of pollutant accumulation in soil, which determines the mobility and bioavailability of pollutants in soil. The adsorption effect in a narrow sense is a physical and chemical process in which a solute transfers from a liquid phase to a solid phase, and in a broad sense, it is a process including precipitation, dissolution, complexation and other chemical actions. Based on the above analysis, accurate simulation of solute transport and adsorption behavior is the core of scientifically predicting the migration and transformation rules and long-term environmental risks of pollutants in rock-soil bodies. The coupling model provides a key theoretical tool for in-depth understanding of the spatiotemporal distribution characteristics and ecological effects of pollutants in porous media. SUMMARY

[0004] The present application aims to provide a high-efficiency solution method for particle pore-solute migration and adsorption equilibrium based on a discrete element pore density flow method, to realize fine and efficient coupling simulation of the migration, diffusion and adsorption process of solutes in porous media in rock-soil bodies. The method constructs the topological structure of the pore-particle system, decouples the convection and diffusion mechanisms in solute transport, and simultaneously solves the adsorption equilibrium state and adsorption rate by combining the isothermal adsorption and quasi-kinetic models, so as to accurately depict the distribution behavior of pollutants between the solid and liquid phases, and provide a reliable numerical simulation tool for pollutant migration prediction and risk assessment in rock-soil environments.

[0005] The present application provides a high-efficiency solution method for particle pore-solute migration and adsorption equilibrium based on a discrete element pore density flow method. The method comprises the following steps: (1) Constructing a porous medium skeleton based on the discrete element method and the pore network model, establishing the topological mapping relationship between the pore space and the adsorbent particles, and initializing the adsorption characteristic parameters of each pore and its adjacent particles. Based on the pore-pore throat channel idea of the solid particle accumulation skeleton and the pore network model, the pore concentration and the particle adsorption capacity are stored in matrix form.

[0006] (2) By decoupling the solute migration process, two core control terms were established, which mainly include the following two points: ① Constructing a solute convection term based on a pore network model to describe solute convection dominated by fluid motion; ② Introducing a solute diffusion term based on Fick's second law to characterize the molecular diffusion process driven by concentration gradient.

[0007] (3) Based on the isothermal adsorption equation, a quasi-kinetic model is used to drive the system to evolve the adsorption equilibrium state of each pore and its associated particle system, thereby characterizing the solid-liquid distribution law of solute in complex porous media.

[0008] Based on the calculated equilibrium deviation, a kinetic model is used to drive the system evolution. Building upon this, explicit time-step iterative calculations are employed to achieve the particle pore-solute adsorption equilibrium.

[0009] Step 10 discretizes the soil and rock mass into a set of solid particles to generate a random stacking model; based on the spatial distance between particles, a geometric threshold is set to identify and segment the pore space between the particle skeleton, and systematically divide it into an interconnected pore unit network, thereby constructing a pore-particle coupling model.

[0010] Step 11: Based on the fluid mass flux in the porous network, construct a solute convection migration model dominated by fluid motion; and based on Fick's second law, construct a molecular diffusion migration model driven by concentration gradient.

[0011] Step 12, based on isothermal adsorption theory and combined with a quasi-kinetic adsorption model, systematically describes the adsorption-desorption dynamic process between each pore unit and its associated particles.

[0012] Step 13 Iterative calculations continue until solute adsorption equilibrium is reached and pore fluid convection and diffusion are stable.

[0013] Step 14: Output the simulation results.

[0014] Steps 20-21 explain the relevant methods for calculating solute convection and diffusion in step 11.

[0015] Steps 30-32 explain the relevant program modules for calculating solute adsorption in particle pores in step 12: This invention proposes an adsorption equilibrium solution architecture based on the discrete element pore density flow method, and constructs a solute convection module, diffusion module, and adsorption equilibrium module under the discrete element framework, realizing efficient and stable coupled simulation of solute transport-adsorption reaction at the pore scale, making up for the shortcomings of existing methods in the field of solute transport in porous media.

[0016] The beneficial effects of this invention are: while realizing the coupling of mechanisms such as convection, diffusion and adsorption reaction, it transforms the calculation of complex two-dimensional porous media systems into highly vectorized matrix operations, so that the model has both the accuracy of physical processes and the efficiency of large-scale calculations, which significantly improves the practical value of pore-scale solute migration simulation and provides a numerical simulation tool that can be applied to engineering practice for pollutant migration prediction and risk assessment in geotechnical environmental engineering. Attached Figure Description

[0017] Figure 1 Flowchart for solving particle pore-solute adsorption equilibrium; Figure 2 This is a flow chart of solute convection and diffusion. Figure 3 Flowchart for calculating adsorption rate; Figure 4 Schematic diagram of the central angle corresponding to the contact surface of particle pores; Figure 5 A schematic diagram of the model box modeling for verification and debugging of this invention; Figure 6 application Figure 5 This diagram illustrates the adsorption capacity of soil particles in the initial stage. At this point, the rate of increase in Cd heavy metal ion adsorption capacity in each soil layer is basically consistent. Figure 6 b represents the mid-stage of the simulation. When heavy metal ions migrate downwards and come into contact with the top of the fracture, it is evident that the adsorption amount of Cd heavy metal ions on the soil particles around the fracture channel shows a significant accelerated trend compared to particles at the same level. Figure 6 c represents the later stage of the simulation. After the seepage stabilizes, there are significant differences in the adsorption of Cd heavy metal ions on the soil particles. The heavy metals preferentially pass through the fracture channels and migrate and adsorb with the surrounding pores and particles. The simulation results are basically consistent with the theoretical situation.

[0018] Figure 7 This invention selects particles from the center of the model surface and plots a graph showing the change in the amount of adsorption of soil and rock particles over time. Detailed Implementation

[0019] Figure 1 Flowchart for solving particle pore-solute adsorption equilibrium Step 10: Discretize the soil and rock mass into a set of contacting solid particles using the discrete element method (DEM), generating a random packing model to reproduce its macroscopic structure. Based on the contact relationships between particles, interconnected fluid domains are delineated by connecting the center points of adjacent particles. A pore fluid topology differential network is established with pore fluid domain units as nodes and pore throat channels connected, resulting in a pore fluid network-solid particle skeleton model. Geometric parameters such as the effective fluid area and actual seepage channel width in the grid cells are calculated based on particle radius and particle coordinates.

[0020] Step 11 Solute convection-diffusion equilibrium Step 21) Initialize pore solute concentration If this is the starting time step for the calculation, then the initial value of the pore solute concentration is defined according to the initial conditions; otherwise, the pore solute concentration after the calculation in the previous time step is inherited as the initial value.

[0021] Step 22) Calculation of solute convection and diffusion in pores Based on the fluid mass exchange flow rate between pore units, the change in solute concentration in the pore units due to convection is calculated, representing the change in solute concentration per unit time caused by solute convection. for: ; In the formula, It is the mass of seepage from the surrounding pores to the central pore. The concentration difference between the surrounding pores and the central pores. The mass of the fluid in the central pore.

[0022] Based on the solute concentration gradient and conductivity coefficient between pore units, the change in solute concentration in the central pore unit due to diffusion is calculated, according to Fick's second law: ; In the formula, J is the diffusion flux. D Where is the diffusion coefficient. For the concentration gradient.

[0023] When applied to solute diffusion in porous networks, the change in solute concentration per unit time due to solute diffusion is... for: ; In the formula, The density of the fluid within the pores. The mass of the fluid within the pores. D Where is the diffusion coefficient. This represents the solute concentration gradient between adjacent fluid domains. This represents the throat area between adjacent pores.

[0024] Step 23) Update the pore solute concentration Solute concentration in the fluid domain after convection-diffusion at one time step for: ; In the formula, The solute concentration in the fluid domain is in the state.

[0025] Step 12 Solute Adsorption Equilibrium Step 31) Initialize particle adsorption capacity If this is the starting time step for the calculation, then the initial value of the particle adsorption amount is defined according to the initial conditions; otherwise, the particle adsorption amount after the calculation in the previous time step is inherited as the initial value.

[0026] Step 32) Adsorption equilibrium model calculation Isothermal adsorption models are used to describe the adsorption and desorption capacity and characteristics of pollutants in soil. Currently, commonly used models include the nonlinear Langmuir adsorption model, the Freundlich adsorption model, and the linear Henry adsorption model.

[0027] The Langmuir model is suitable for describing ideal adsorption processes with uniform adsorption sites in a monolayer: ; In the formula, This represents the equilibrium adsorption capacity of the particles. Qmax For maximum adsorption capacity, K L It is the Langmuir constant. Ce This is the pore equilibrium concentration.

[0028] The Freundlich model is suitable for describing non-ideal adsorption processes involving multilayered molecules and non-uniform adsorption surfaces. ; In the formula Qe This represents the equilibrium adsorption capacity of the particles. K F Where is the adsorption capacity constant. It is a heterogeneity index. Ce This is the pore equilibrium concentration.

[0029] The Henry model is suitable for describing linear adsorption processes with constant adsorption affinity at low concentrations. ; In the formula, Q e This represents the equilibrium adsorption capacity of the particles. K H For allocation coefficients, Ce This is the pore equilibrium concentration.

[0030] Step 33) Calculate the equilibrium adsorption capacity of the particles. The adsorption process occurs at the contact surface between the particle and the adjacent pores, and the law of conservation of mass also applies during this process. Its governing equation can be expressed as: ; In the formula, Ct This represents the current pore equilibrium concentration. Ce For pore equilibrium concentration, A For pore volume, Qe This represents the equilibrium adsorption capacity of the particles. Qt This represents the current particle adsorption capacity. M This refers to particle mass.

[0031] In summary, the equilibrium adsorption capacity of particles can be determined by the current pore concentration and the particle adsorption capacity: ; In the formula, Qe This represents the equilibrium adsorption capacity of the particles. This represents the current pore concentration. Qt This represents the current particle adsorption capacity.

[0032] Step 34) Calculate the adsorption rate Adsorption rates are often described by quasi-kinetic adsorption processes.

[0033] When the adsorption process is dominated by physical adsorption, the pseudo-first-order kinetic adsorption equation is as follows: ; In the formula, k 1 The adsorption coefficient is the pseudo-first-order kinetic coefficient. Q e This represents the equilibrium adsorption capacity of the particles. Qt For pore equilibrium concentration, It is the ratio of the central angle to the circumferential angle corresponding to the solid-liquid interface.

[0034] When the adsorption process is dominated by chemisorption, the pseudo-second-order kinetic adsorption equation is as follows: ; In the formula, k 2 The adsorption coefficient is the pseudo-second-order kinetic coefficient. Q e This represents the equilibrium adsorption capacity of the particles. Qt For pore equilibrium concentration, It is the ratio of the central angle to the circumferential angle corresponding to the solid-liquid interface.

[0035] Step 35) Adsorption capacity of new particles and concentration of solute in pores For each soil particle unit, the change in particle adsorption amount per unit time due to solute adsorption in the neighboring M pores is as follows: ; In the formula, This represents the change in particle adsorption capacity. For particles to the first The adsorption rate of each pore.

[0036] For each pore grid cell, the change in pore solute concentration per unit time due to the adsorption of N neighboring particles is: ; In the formula, The change in pore concentration Let be the adsorption rate of the j-th particle on this pore. Let the mass of the j-th particle be... This represents the area of ​​the pore.

[0037] Step 13: Repeat steps 11-12 until solute adsorption equilibrium is reached and pore fluid convection and diffusion are stable. Iteratively calculate the solute convection-diffusion equilibrium and solute adsorption equilibrium processes until the following convergence conditions are met simultaneously.

[0038] The pore concentration change converges, and the solute concentration change in all pore units within a continuous time step is less than a set threshold: ; In the formula, The pore concentration at the next time step. The pore concentration at the current time step. This represents the concentration convergence tolerance.

[0039] The changes in particle adsorption capacity converge, and the changes in adsorption capacity of all particle units within a continuous time step should be less than a set threshold: ; In the formula, This represents the particle adsorption amount at the next time step. This represents the particle adsorption amount at the current time step. This represents the convergence tolerance of adsorption capacity.

[0040] Step 14: Output the simulation results. (See attached image.) Figure 5 , Figure 6 a, Figure 6 b、 Figure 6 c. Figure 7 .

[0041] A discrete element model of a model box was designed to verify the solute adsorption equilibrium model of the pore-fracture dual medium. Figure 5A schematic diagram of the model box is provided. The central part of the model represents the soil and rock model to be tested. Initially, it is set as clean soil and rock particles. The upper and lower voids are defined as pores, which are set as fluid boundaries during simulation. A vertical crack is defined in the central region of the model. The upper boundary void is set to a constant high pressure, and the lower boundary void is given a constant low pressure. Cd heavy metal ions are selected as the characteristic pollutant. A Cd heavy metal solution with a concentration of 1.6 mg / L is constantly injected into the upper boundary void. The Langmuir isothermal adsorption model is selected, with Qmax = 50 mg / g and K... L =0.62 L / mg, pseudo-second-order kinetic adsorption process, K2=0.016 h⁻¹ -1 Under the influence of a vertically downward pressure field, pollutants begin to migrate and be adsorbed.

[0042] The above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. An efficient solution method for particle pore-solute migration and adsorption equilibrium based on discrete element method (DEM), characterized in that, Includes the following steps: (1) Based on the discrete element method and pore network model, a porous medium skeleton is constructed, and the topological mapping relationship between the pore space and the adsorbent particles is established. The adsorption characteristic parameters of each pore and its neighboring particles are initialized. Based on the pore-throat channel idea of ​​the solid particle stacking skeleton and pore network model, the pore concentration and particle adsorption amount are stored in a matrix. (2) By decoupling the solute migration process, two core control terms were established, including the following two points: ① Constructing a solute convection term based on a pore network model to describe solute convection dominated by fluid motion; ② Introduce a solute diffusion term based on Fick's second law to characterize the molecular diffusion process driven by the concentration gradient; (3) Based on the isothermal adsorption equation, a quasi-kinetic model is used to drive the system to evolve the adsorption equilibrium state of each central pore and its associated particle system, thereby characterizing the solid-liquid partitioning law of solute in complex porous media. Based on the calculated equilibrium deviation, a kinetic model is used to drive the system evolution; on this basis, explicit time-step iterative calculation is used to realize the particle pore-solute adsorption equilibrium.

2. The efficient solution method for particle pore-solute migration and adsorption equilibrium based on discrete element density flow method according to claim 1, characterized in that, The specific steps are as follows: Step 10 discretizes the soil and rock mass into a set of solid particles to generate a random stacking model; based on the spatial distance between particles, a geometric threshold is set to identify and segment the pore space between the particle skeleton, and systematically divide it into an interconnected pore unit network, thereby constructing a pore-particle coupling model. Step 11: Based on the fluid mass flux in the porous network, construct a solute convection migration model dominated by fluid motion; and based on Fick's second law, construct a molecular diffusion migration model driven by concentration gradient. Step 12, based on the isothermal adsorption theory and combined with the quasi-kinetic adsorption model, systematically describes the adsorption-desorption dynamic process between each pore unit and its associated particles. Step 13 Iterative calculations continue until solute adsorption equilibrium is reached and pore fluid convection and diffusion are stable; Step 14: Output the simulation results.

3. The efficient solution method for particle pore-solute migration and adsorption equilibrium based on discrete element pore density flow method according to claim 1, characterized in that, In step 10, the rock and soil mass is discretized into a set of solid particles in contact with each other using the discrete element method, and a random stacking model is generated to reproduce its macroscopic structure. Based on the contact relationship between particles, the interconnected fluid domains are decomposed by connecting the center points of adjacent particles. A pore fluid topology differential network is established with pore fluid domain units as nodes and pore throat channels connected to each other, to obtain the pore fluid network-solid particle skeleton model. Based on the pore-throat channel model of solid particle packing skeleton and pore network, the computational fluid domain is meshed; starting directly from the "point" represented by the particle center, a three-dimensional tetrahedral fluid computational mesh is established based on the Delaunay triangulation algorithm; at the same time, the mesh re-division caused by particle motion during the calculation is considered. (2) Calculate the geometric parameters of the pore grid: Based on the particle radius and particle coordinates, calculate the geometric parameters such as the effective fluid volume and the actual seepage channel area in the grid cell; (3) Establish a three-dimensional discrete element fluid-structure interaction numerical simulation model; construct program modules to handle fluid seepage, fluid-solid interaction, and solid-fluid interaction; and on the basis of the above, use an explicit time-step iterative algorithm to realize the three-dimensional discrete element fluid-structure interaction calculation model.

4. The efficient solution method for particle pore-solute migration and adsorption based on discrete element pore density flow method according to claim 1, characterized in that, Step 11: Method for solute convection-diffusion equilibrium Step 21) Initialize pore solute concentration If this is the starting time step for the calculation, then the initial value of the pore solute concentration is defined according to the initial conditions; otherwise, the pore solute concentration after the calculation in the previous time step is inherited as the initial value. Step 22) Calculation of solute convection and diffusion in pores Based on the fluid mass exchange flow rate between pore units, the change in solute concentration in the pore units due to convection is calculated, representing the change in solute concentration per unit time caused by solute convection. for: ； In the formula, It is the mass of seepage from the surrounding pores to the central pore. The concentration difference between the surrounding pores and the central pores. For the fluid mass in the central pore; Based on the solute concentration gradient and conductivity coefficient between pore units, the change in solute concentration in the central pore unit due to diffusion is calculated, according to Fick's second law: ； In the formula, J is the diffusion flux. D Where is the diffusion coefficient. For concentration gradient; When applied to solute diffusion in porous networks, the change in solute concentration per unit time due to solute diffusion is... for: ； In the formula, The density of the fluid within the pores. The mass of the fluid within the pores. D Where is the diffusion coefficient. This represents the solute concentration gradient between adjacent fluid domains. The throat area between adjacent pores; Step 23) Update the pore solute concentration Solute concentration in the fluid domain after convection-diffusion at one time step for: ； In the formula, The solute concentration in the fluid domain is in the state.

5. The efficient solution method for particle pore-solute migration and adsorption based on discrete element pore density flow method according to claim 1, characterized in that, Step 12 involves the following steps in the procedure for calculating the adsorption of solutes in particulate pores. Step 31) Initialization of particle adsorption capacity If this is the starting time step for the calculation, then the initial value of the particle adsorption amount is defined according to the initial conditions; otherwise, the particle adsorption amount after the calculation in the previous time step is inherited as the initial value. Step 32) Adsorption equilibrium model calculation Isothermal adsorption models are used to describe the adsorption and desorption capacity and characteristics of pollutants in soil. Currently, commonly used models include the nonlinear Langmuir adsorption model, the Freundlich adsorption model, and the linear Henry adsorption model. The Langmuir model is suitable for describing ideal adsorption processes with uniform adsorption sites in a monolayer: ； In the formula, This represents the equilibrium adsorption capacity of the particles. Qmax For maximum adsorption capacity, K L It is the Langmuir constant. Ce The pore equilibrium concentration; The Freundlich model is suitable for describing non-ideal adsorption processes involving multilayered molecules and non-uniform adsorption surfaces. ； In the formula Qe This represents the equilibrium adsorption capacity of the particles. K F Where is the adsorption capacity constant. It is a heterogeneity index. Ce The pore equilibrium concentration; The Henry model is suitable for describing linear adsorption processes with constant adsorption affinity at low concentrations. ； In the formula, Q e This represents the equilibrium adsorption capacity of the particles. K H For allocation coefficients, Ce The pore equilibrium concentration; Step 33) Calculate the equilibrium adsorption capacity of the particles. The adsorption process occurs at the contact surface between the particle and the adjacent pores, and the law of conservation of mass also applies during this process. Its governing equation can be expressed as: ； In the formula, Ct This represents the current pore equilibrium concentration. Ce For pore equilibrium concentration, A For pore volume, Qe This represents the equilibrium adsorption capacity of the particles. Qt This represents the current particle adsorption capacity. M For particle mass; In summary, the equilibrium adsorption capacity of particles can be determined by the current pore concentration and the particle adsorption capacity: ； In the formula, Qe This represents the equilibrium adsorption capacity of the particles. This represents the current pore concentration. Qt This represents the current particle adsorption capacity. Step 34) Calculate the adsorption rate Adsorption rates are often described by quasi-kinetic adsorption processes; When the adsorption process is dominated by physical adsorption, the pseudo-first-order kinetic adsorption equation is as follows: ； In the formula, k 1 The adsorption coefficient is the pseudo-first-order kinetic coefficient. Q e This represents the equilibrium adsorption capacity of the particles. Qt For pore equilibrium concentration, It is the ratio of the central angle to the circumferential angle corresponding to the solid-liquid interface. When the adsorption process is dominated by chemisorption, the pseudo-second-order kinetic adsorption equation is as follows: ； In the formula, k 2 The adsorption coefficient is the pseudo-second-order kinetic coefficient. Q e This represents the equilibrium adsorption capacity of the particles. Qt For pore equilibrium concentration, It is the ratio of the central angle to the circumferential angle corresponding to the solid-liquid interface. Step 35) Adsorption capacity of new particles and concentration of solute in pores For each soil particle unit, the change in particle adsorption amount per unit time due to solute adsorption in the neighboring M pores is as follows: ； In the formula, This represents the change in particle adsorption capacity. For particles to the first The adsorption rate of each pore; For each pore grid cell, the change in pore solute concentration per unit time due to the adsorption of N neighboring particles is: ； In the formula, The change in pore concentration Let be the adsorption rate of the j-th particle on this pore. Let the mass of the j-th particle be... This represents the area of ​​the pore.

6. The efficient solution method for particle pore-solute migration and adsorption based on discrete element pore density flow method according to claim 1, characterized in that, Step 13: Repeat steps 11-12 until solute adsorption equilibrium is reached and pore fluid convection and diffusion are stable. Iteratively calculate the solute convection-diffusion equilibrium and solute adsorption equilibrium processes until the following convergence conditions are simultaneously met; The pore concentration change converges, and the solute concentration change in all pore units within a continuous time step is less than a set threshold: ； In the formula, The pore concentration at the next time step. The pore concentration at the current time step. For concentration convergence tolerance; The changes in particle adsorption capacity converge, and the changes in adsorption capacity of all particle units within a continuous time step should be less than a set threshold: ； In the formula, This represents the particle adsorption amount at the next time step. This represents the particle adsorption amount at the current time step. This represents the convergence tolerance of adsorption capacity.

Images