Pollutant transport prediction method based on smoothed particle hydrodynamics

Abstract

The present invention relates to a pollutant transport prediction method based on smoothed particle hydrodynamics, which is implemented through the following steps: Step 1, system initialization: initializing parameter information; Step 2, particle generation: including generating fluid particles, inflow particles, and outflow particles; Step 3, listing the discrete format of the governing equations and the boundary condition application method and iterative calculation: using the SPH method to discretize the governing equations describing pollutant transport, and iterative calculation starting from time 0 based on the initial parameter conditions; Step 4, outputting the results: including the position, velocity, and concentration of all types of particles at the final moment. The present invention can accurately predict the changes in concentration when the pollutants are initially discontinuously distributed, without generating non-physical oscillation phenomena.

Description

Technical Field

[0001] The present invention relates to the technical field of SPH fluid numerical simulation, in particular to a new method for using SPH to predict concentration changes during the transport process of pollutants under the condition of initial discontinuous distribution. Background Art

[0002] Smooth Particle Hydrodynamics (SPH) is a numerical method and an interpolation method. It was first introduced in the field of astrophysics[1][2] and is now widely used in fluid numerical simulation. SPH fluid numerical simulation is based on the SPH numerical method to solve mathematical models. The continuous mathematical model is discretized in time and space, so that the computer can automatically solve the changing trends of various physical properties during the fluid motion process based on the initial parameter conditions.

[0003] In the SPH numerical simulation process, the area of ​​fluid flow studied is called the problem domain. In the case of a large problem domain, the key fluid flow problems in a small area of ​​the problem domain can be studied. This small area is called the computational domain. For example, for the problem of pollutant transport in a river, we only focus on the changes in pollutants in the polluted area of ​​the river. The entire river is the problem domain, and the polluted area is the computational domain. The size of the computational domain directly determines the amount of calculation, so the computational domain cannot be expanded indefinitely and needs to be controlled within a certain range. The fluid flows into the computational domain from the inlet boundary and flows out of the computational domain from the outflow boundary. The inlet boundary and the outflow boundary are also the intersection of the computational domain and the problem domain.

[0004] SPH considers the fluid as a group of discrete particles, which include inflow particles, fluid particles and outflow particles in the pollutant transport process. Each particle usually has five physical properties: position, velocity, density, mass and concentration. As we all know, the new position of the fluid particle can be calculated by the original position and velocity of the fluid particle. The inflow particle and the outflow particle can be regarded as special fluid particles and also need to have the physical properties of fluid particles. Determining the new concentration of the fluid particle at the next time step by the original concentration and diffusion coefficient of the fluid particle is the key to ensuring the correctness of the simulation results. The initial discontinuous distribution of pollutants is a relatively special type of pollutant transport problem. For this problem, the previous method [3] is to calculate the diffusion flux of the current fluid particle by the concentration of other fluid particles searched in the support domain, and then calculate the new concentration of the current fluid particle at the next moment by the diffusion flux of other fluid particles searched in the support domain, the velocity, concentration and diffusion coefficient of the current fluid particle. However, this method will produce non-physical oscillation problems during the calculation process. The support domain refers to the local area around the current fluid particle.

[0005] References:

[0006] [1]Lucy L BA numerical approach to the testing of the fissionhypothesis.The Astrophysical Journal,1977,82(12):1013-1024.

[0007] [2]Gingold RA, Monaghan J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375-389.

[0008] [3]Fatehi R,Fayazbakhsh MA,Manzari M T.On discretization of second-order derivatives in Smoothed Particle Hydrodynamics[J].Proceedings of World Academy of Science Engineering&Technolog,2008. Summary of the Invention

[0009] The purpose of the present invention is to provide a pollutant transport prediction method based on smoothed particle hydrodynamics. This method can predict the concentration changes of pollutants during the flow of water, obtaining a realistic pollutant propagation process. This method can be applied to practical engineering problems, such as predicting the propagation process of oil leaks, and has the characteristics of low error and wide applicability. The technical solution of the present invention is as follows:

[0010] A pollutant transport prediction method based on smoothed particle hydrodynamics includes the following steps:

[0011] Step 1: Initialize the system parameters. The parameter information includes the vertex coordinates of the computational domain, the vertex coordinates of the inflow boundary, the vertex coordinates of the outflow boundary, the fluid particle spacing, the smoothing length, the total simulation time, the time step, and the diffusion coefficient.

[0012] Step 2: Generate particles;

[0013] 1) 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 x-direction velocity, fluid particle y-direction velocity, fluid particle type identification number 1, and fluid particle concentration;

[0014] 2) 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 x-direction velocity, inflow particle y-direction velocity, inflow particle type identification number 2, and inflow particle concentration;

[0015] 3) Generate outflow particles, evenly distribute the outflow particles in the fluid domain, and add initial outflow particle information; the initial outflow particle information includes outflow particle density, outflow particle mass, outflow particle x-direction velocity, outflow particle y-direction velocity, outflow particle type identification number 3, and outflow particle concentration;

[0016] Step 3: List the discrete format of the control equations and the boundary conditions and perform iterative calculations:

[0017] 1) Discrete format of the control equation: The control equation describing the pollutant transport problem is discretized using the SPH numerical method into the following format:

[0018]

[0019]

[0020]

[0021]

[0022]

[0023] Among them, W(x i -x j ,h) is called the smooth kernel function (hereinafter referred to as the "kernel function"), is the gradient of the kernel function, h is called the smoothing length, which is equal to 1.33dx, dx is the initial regular distribution of particle spacing, the kernel function usually selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |x i -x j | / h, then the kernel function W(x i -x j ,h)The expression is:

[0024]

[0025] i represents the current fluid particle, and its adjacent 2h range is its support domain. j is one of the fluid particles in the support domain of the current fluid particle. i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain, N is the number of all fluid particles in the support domain of the current particle i, n represents the current moment, and n+1 represents the next moment; 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, The intermediate variables in the calculation process have no actual meaning, C i and C j are the concentrations of the current fluid particle i and the fluid particle j in its support domain, D represents the diffusion coefficient, and t is the time step;

[0026] In formulas (1) to (4), is the characteristic function, and its expression is as follows:

[0027]

[0028]

[0029]

[0030]

[0031] x i 、y i Represents the position components of the current fluid particle i in the x direction and y direction, x j 、y j represent the position components of the fluid particle j in the support domain in the x-direction and y-direction respectively;

[0032] In formulas (1) to (4),

[0033]

[0034]

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

[0036] 2) Boundary condition application method;

[0037] 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.

[0038] 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. Other properties remain unchanged.

[0039] 3) Iterative calculation; the current time starts from 0, and the total time is given in step 1; the specific calculation process is as follows:

[0040] ① Determine whether the current moment has reached the total time; if yes, end the calculation; if not, repeat the calculation process ② to ④;

[0041] ② All particles in the computational domain are assigned coordinates during initialization. Here, the new position of each particle is equal to the original position plus the particle's own velocity 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 to become a new fluid particle, and the inlet particle is supplemented at the corresponding position of 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 become a new outflow particle.

[0042] ③ Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses the other fluid particles searched in the support domain by the current fluid particle to calculate the current fluid particle according to equations (1) to (2). and

[0043] ④ Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other fluid particles searched in the support domain by the current fluid particle to calculate the new concentration of the current fluid particle after one time step according to equations (3) to (5);

[0044] ⑤Current moment = current moment + 1 time step; repeat calculation process ①;

[0045] Step 4: Output the results. When the total time is up, end the loop and output the position, velocity, and concentration of all fluid particles, the position, velocity, and concentration of all inflow particles, and the position, velocity, and concentration of all outflow particles at this moment as the final results.

[0046] Furthermore, the generated particles described in step 2: the concentration of the inflow particles and the outflow particles is 0.

[0047] Furthermore, the method for supplementing the incoming particles at the corresponding position of the outermost layer of the inlet boundary described in the calculation process ② in step 3 is as follows: when the incoming particle flows into the fluid domain, a new fluid particle is generated in the fluid domain, and the new fluid particle inherits the position, velocity, concentration and other properties of the original incoming particle; 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.

[0048] The method of the present invention has the following characteristics: it can accurately predict the change of concentration of pollutants under the condition of initial discontinuous distribution, and will not produce non-physical oscillation phenomenon. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] In order to illustrate the reliability and accuracy of the technical solutions in the embodiments of the present invention in more detail and intuitively, the following briefly introduces the drawings required for use in the embodiments. The drawings described below are only some embodiments of the present invention. Through some simple modifications, the present invention can be applied to other embodiments without paying any creative work.

[0050] Figure 1 Physical Model of Pollutant Transport Problem

[0051] Figure 2 Flowchart of the new method

[0052] Figure 3 Initial distribution map of pollutants

[0053] Figure 4 Final distribution map of pollutants

[0054] Figure 5 On the y=2m section, the results of the pollutant concentration change prediction using the new SPH-based method of the present invention are compared with the exact solution. DETAILED DESCRIPTION

[0055] The specific solution process of this method is illustrated by the example of discontinuous initial distribution of pollutants. This problem can be described as follows: the initial pollutants are only distributed in a part of the fluid and the initial concentration values ​​in the polluted area are the same. The physical model of this problem is as follows Figure 1 As shown, 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 initial distribution area of ​​the pollutants is a square area A″B″C″D″. The length of the inlet boundary A′D′ is 4m, the length of the outlet boundary B′C′ is 4m, the width of the fluid domain is 4m, and the width of the pollutant area is 1m, located in the center of the fluid domain. In this problem, all types of particles have a flow velocity of 0.1m / s in the x-direction and no flow velocity in the y-direction. The diffusion coefficient (how fast the pollutant diffuses) is 0.02m2 The present invention will be described in further detail below with reference to the accompanying drawings. Figure 2 The program flow chart of the proposed method is given. It includes:

[0056] Step 1: Initialize the system parameters. The vertex coordinates of the computational domain are A(-0.04, 4.0), B(4.04, 4.0), C(4.04, 0), and D(-0.04, 0). The vertex coordinates of the pollutant distribution area are A″(1.5, 2.5), B″(2.5, 2.5), C″(2.5, 1.5), and D″(1.5, 1.5). The vertex coordinates of the inflow boundary are A′D′(0, 4.0), (0, 0), and the vertex coordinates of the outflow boundary are B′C′(4.0, 4.0), (4.0, 0). The fluid particle spacing Δx = Δy = 0.02m, the smoothing length h = 1.33Δx = 1.33Δy, the total simulation time is 5s, using an equal time step, a time step Δt = 0.01s, and a diffusion coefficient of 0.02m. 2 / s.

[0057] Step 2: Generate particles, including

[0058] 1) Generate fluid particles, distribute them evenly in the fluid domain, and add initial fluid particle information. Evenly distribute 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 , x-direction velocity u = 0.1 m / s, y-direction velocity v = 0, and fluid particle type identification number is 1. Among them, the fluid particle concentration in the pollutant distribution area A″B″C″D″ is C = 10 mg / L, and the concentration of other fluid particles is 0.

[0059] 2) Generate inflow particles: Evenly distribute inflow particles 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. The inflow particle mass m = 0.1 kg and density ρ = 1000 kg / m 3 , x-direction velocity u = 0.1 m / s, y-direction velocity v = 0, inflow particle type identification number is 2, and inflow particle concentration is 0.

[0060] 3) Generate outflow particles: Distribute outflow particles evenly at the outlet boundary and add initial outflow particle information. Arrange two layers of inlet particles outward from the outflow boundary B′C′ with the particle spacing as the unit. The outflow particle mass m = 0.1 kg and density ρ = 1000 kg / m 3 , x-direction velocity u = 0.1 m / s, y-direction velocity v = 0, outflow particle type identification number is 3, and outflow particle concentration is 0.

[0061] Step 3: List the discrete format of the control equations and the boundary conditions and perform iterative calculations:

[0062] 2) Discrete format of the control equation: The control equation describing the pollutant transport problem is discretized using the SPH numerical method into the following format:

[0063]

[0064]

[0065]

[0066]

[0067]

[0068] Among them, W(x i -x j ,h) is called the smooth kernel function (hereinafter referred to as the "kernel function"), is the gradient of the kernel function, h is called the smoothing length, which is equal to 1.33dx, dx is the initial regular distribution of particle spacing, the kernel function usually selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |x i -x j | / h, then the kernel function W(x i -x j ,h)The expression is:

[0069]

[0070] i represents the current fluid particle, and its adjacent 2h range is its support domain. j is one of the fluid particles in the support domain of the current fluid particle. i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain, N is the number of all fluid particles in the support domain of the current particle i, n represents the current moment, and n+1 represents the next moment. 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, The intermediate variables in the calculation process have no actual meaning, C i and C j are the concentrations of the current fluid particle i and the fluid particle j in its support domain, D represents the diffusion coefficient, and t is the time step.

[0071] In formulas (1) to (4), is the characteristic function, and its expression is as follows:

[0072]

[0073]

[0074]

[0075]

[0076] x i 、y i Represents the position components of the current fluid particle i in the x direction and y direction, x j 、y j represent the position components of the fluid particle j in the support domain in the x and y directions, respectively.

[0077] In formulas (1) to (4),

[0078]

[0079]

[0080] in, and Represent the gradient of the kernel function of the current fluid particle i in the x direction and y direction respectively.

[0081] 2) How boundary conditions are applied.

[0082] 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.

[0083] The velocity condition is applied to the inlet boundary. When an inflow particle flows into the fluid domain, its inflow particle type identification number is changed from inflow particle type identification number 2 to fluid particle type identification number 1, becoming a new fluid particle. Other properties remain unchanged.

[0084] 3) Iterative calculation. The current time starts at 0 and the simulation is performed for 5 seconds. The specific calculation process is as follows:

[0085] ① Determine whether the current time has reached 5 seconds. If yes, end the calculation; if not, repeat the calculation process ② to ④;

[0086] ② All particles in the computational domain are assigned coordinates during initialization. Here, each particle's new position is equal to its original position plus the particle's velocity multiplied by the time step of 0.01s. For inflow particles that enter the fluid domain after the position update, the type identification number in their particle properties is changed to 1, becoming a new fluid particle. These inflow particles are then added to the corresponding positions at the outermost layer of the inlet boundary. For fluid particles that exit the fluid domain after the position update, the type identification number in their particle properties is changed to 3, becoming a new outflow particle.

[0087] ③ Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses the other fluid particles searched in the support domain by the current fluid particle to calculate the current fluid particle according to equations (1) to (2). and

[0088] ④ Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other fluid particles searched in the support domain by the current fluid particle to calculate the new concentration of the current fluid particle after 0.01s according to equations (3) to (5);

[0089] ⑤Current time = current time + 1 time step. Repeat the calculation process ①.

[0090] Step 4: Output the results. When the total time is up, end the loop and output the position, velocity, and concentration of all fluid particles, the position, velocity, and concentration of all inflow particles, and the position, velocity, and concentration of all outflow particles at that moment to a file as the final result.

[0091] Step 4: Output the results. After the total time reaches 5 seconds, the loop ends and the calculation stops. The position, velocity, and concentration information of all particles at the current moment are output as the final calculation results.

[0092] Furthermore, the system initialization described in step 1:

[0093] 1) The fluid particle spacing, inflow particle spacing, and outflow particle spacing are all 0.02m;

[0094] 2) 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.

[0095] Furthermore, the particles generated in step 2 have two layers of inflow particles and two layers of outflow particles;

[0096] Furthermore, the inflow particle is supplemented at the corresponding position on the outermost layer of the inlet boundary as described in step 3 (2). When the abscissa of the inflow particle position is greater than or equal to 0, a new fluid particle is generated within the fluid domain. This new fluid particle inherits the position, velocity, mass, density, concentration, and other properties of the original inflow particle. The position of the inflow particle with an abscissa greater than or equal to 0 is then changed to the abscissa minus 0.04 m.

[0097] Figure 5 On the section y=2m, the results of the present invention's prediction of pollutant concentration changes using the new SPH-based method are consistent with the exact solution.

[0098] Although the present invention has been described above with reference to the drawings, the present invention is not limited to the above-mentioned specific embodiments. The above-mentioned specific embodiments are merely illustrative and not restrictive. Under the guidance of the present invention, ordinary technicians in this field can make many variations without departing from the purpose of the present invention, and these are all protected by the present invention.

Claims

1. A pollutant transport prediction method based on smoothed particle hydrodynamics, comprising the following steps: Step 1: Initialize the system parameters. The parameter information includes the vertex coordinates of the computational domain, the vertex coordinates of the inflow boundary, the vertex coordinates of the outflow boundary, the fluid particle spacing, the smoothing length, the total simulation time, the time step, and the diffusion coefficient. Step 2: Generate particles; 1) 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 x-direction velocity, fluid particle y-direction velocity, fluid particle type identification number 1, and fluid particle concentration; 2) 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 x-direction velocity, inflow particle y-direction velocity, inflow particle type identification number 2, and inflow particle concentration; 3) Generate outflow particles, evenly distribute the outflow particles in the fluid domain, and add initial outflow particle information; the initial outflow particle information includes outflow particle density, outflow particle mass, outflow particle x-direction velocity, outflow particle y-direction velocity, outflow particle type identification number 3, and outflow particle concentration; Step 3: List the discrete format of the control equations and the boundary conditions and perform iterative calculations: 1) Discrete format of the control equation: The control equation describing the pollutant transport problem is discretized using the SPH numerical method into the following format: Among them, W(x i -x j ,h) is called the smooth kernel function, is the gradient of the kernel function, h is called the smoothing length, which is equal to 1.33dx, dx is the initial regular distribution of particle spacing, the kernel function selects the cubic spline function, and R represents x i and x j The relative distance between them, R = |x i -x j | / h, then the kernel function W(x i -x j ,h)The expression is: i represents the current fluid particle, and its adjacent 2h range is its support domain. j is one of the fluid particles in the support domain of the current fluid particle. i and x j are the positions of the current fluid particle i and the fluid particle j in its support domain, N is the number of all fluid particles in the support domain of the current particle i, n represents the current moment, and n+1 represents the next moment; 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, The intermediate variables in the calculation process have no actual meaning, C i and C j are the concentrations of the current fluid particle i and the fluid particle j in its support domain, D represents the diffusion coefficient, and t is the time step; In formulas (1) to (4), is the characteristic function, and its expression is as follows: x i 、y i Represents the position components of the current fluid particle i in the x direction and y direction, x j 、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 (1) to (4), in, and Represent the gradient of the kernel function of the current fluid particle i in the x direction and y direction respectively; 2) 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. Other properties remain unchanged. 3) Iterative calculation; the current time starts from 0, and the total time is given in step 1; the specific calculation process is as follows: ① Determine whether the current moment has reached the total time; if yes, end the calculation; if not, repeat the calculation process ② to ④; ② All particles in the computational domain are assigned coordinates during initialization. Here, the new position of each particle is equal to the original position plus the particle's own velocity 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 to become a new fluid particle, and the inlet particles are supplemented at the corresponding position of 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 become a new outflow particle. ③ Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses the other fluid particles searched in the support domain by the current fluid particle to calculate the current fluid particle according to equations (1) to (2). and ④ Each fluid particle takes turns to participate in the calculation as the current fluid particle, and uses other fluid particles searched in the support domain by the current fluid particle to calculate the new concentration of the current fluid particle after one time step according to equations (3) to (5); ⑤Current moment = current moment + 1 time step; repeat calculation process ①; Step 4: Output the results. When the total time is up, end the loop and output the position, velocity, and concentration of all fluid particles, the position, velocity, and concentration of all inflow particles, and the position, velocity, and concentration of all outflow particles at this moment as the final results.

2. The pollutant transport prediction method according to claim 1, wherein: The generated particles described in step 2: the concentration of inflow particles and outflow particles is 0.

3. The pollutant transport prediction method according to claim 1, wherein: The method for supplementing the incoming particles at the corresponding position of the outermost layer of the inlet boundary described in the calculation process ② in step 3 is as follows: when the incoming particle flows into the fluid domain, a new fluid particle will be generated in the fluid domain. The new fluid particle inherits the position, velocity, and concentration properties of the original inflow particle; then, the coordinate value of the inflow 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 inflow particle is generated.

Images