CN112257313A - Pollutant transport high-resolution numerical simulation method based on GPU acceleration - Google Patents
Pollutant transport high-resolution numerical simulation method based on GPU acceleration Download PDFInfo
- Publication number
- CN112257313A CN112257313A CN202011131148.0A CN202011131148A CN112257313A CN 112257313 A CN112257313 A CN 112257313A CN 202011131148 A CN202011131148 A CN 202011131148A CN 112257313 A CN112257313 A CN 112257313A
- Authority
- CN
- China
- Prior art keywords
- gpu
- pollutant
- vector
- flow
- unit grid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000003344 environmental pollutant Substances 0.000 title claims abstract description 77
- 231100000719 pollutant Toxicity 0.000 title claims abstract description 77
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000004088 simulation Methods 0.000 title claims abstract description 41
- 230000001133 acceleration Effects 0.000 title claims abstract description 26
- 230000004907 flux Effects 0.000 claims abstract description 34
- 238000004364 calculation method Methods 0.000 claims abstract description 20
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 61
- 239000013598 vector Substances 0.000 claims description 47
- 238000009792 diffusion process Methods 0.000 claims description 21
- 230000008569 process Effects 0.000 claims description 15
- 239000000126 substance Substances 0.000 claims description 7
- 239000000463 material Substances 0.000 claims description 4
- 230000010355 oscillation Effects 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- ZLSWBLPERHFHIS-UHFFFAOYSA-N Fenoprop Chemical compound OC(=O)C(C)OC1=CC(Cl)=C(Cl)C=C1Cl ZLSWBLPERHFHIS-UHFFFAOYSA-N 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 3
- 210000003205 muscle Anatomy 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 239000003403 water pollutant Substances 0.000 abstract description 4
- 238000013508 migration Methods 0.000 abstract description 2
- 230000005012 migration Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 12
- 239000000356 contaminant Substances 0.000 description 10
- 230000003068 static effect Effects 0.000 description 6
- 238000003911 water pollution Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012821 model calculation Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000006424 Flood reaction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 230000034994 death Effects 0.000 description 1
- 231100000517 death Toxicity 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000007716 flux method Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 230000005180 public health Effects 0.000 description 1
- 239000010865 sewage Substances 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Geometry (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Computer Graphics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a GPU acceleration-based pollutant transport high-resolution numerical simulation method, which comprises the steps of reading a DEM terrain file consisting of high-resolution unit grids, and obtaining an assigned variable; reading initial conditions and boundary conditions, and replacing variables from a CPU to a GPU by adopting a finite volume method space discrete control equation; dividing dry and wet calculation areas, and calculating flux on a unit grid interface on a GPU by adopting an HLLC approximate Riemann solver; calculating a source item of a control equation on a GPU; calculating the time step length of iterative update according to the Curian number; acquiring hydraulic element and concentration information of all unit grids; and copying the hydraulic element and the concentration information from the GPU video memory to a host memory and outputting. The problem of among the prior art study water pollutant migration motion law's numerical model operating efficiency and precision relatively poor is solved.
Description
Technical Field
The invention belongs to the technical field of numerical simulation acceleration, and particularly relates to a pollutant transport high-resolution numerical simulation method based on GPU acceleration.
Background
In recent years, a number of serious mountain torrent outbreaks have been reported around the world. For example, in 3 months in 2019, severe mountain torrent outbreaks occurred in different areas of babuca province in indonesia, and 63 lives were lost in large-scale flood disasters. In 7 months 2012, a catastrophic rainstorm caused 79 deaths in Beijing, China. Due to the violence and unpredictability of mountain torrents, it poses a great threat to human life and property. Flood water also has a profound impact on cities. In cities, large floods may damage sewage treatment plants or other pollutant-releasing facilities, and may also pose a significant risk to public health and local economy. Therefore, the understanding of the transportation process of the pollutants has important guiding significance for water conservancy and environmental engineering.
In order to deal with the unpredictability and the acuteness of water pollution events, researchers use a two-dimensional hydrodynamics equation and a transport equation to establish a large number of numerical models to analyze and research the water pollutant migration motion law. However, different numerical methods cause differences in the calculation accuracy of the model, and the model runs in a long time and even cannot run on a computer; once large-scale water body pollution occurs, emergency decision and corresponding measures are needed, corresponding judgment is made on the transportation of pollutants in the water body, and the model is operated on a computer without time at all. Therefore, the rapid and efficient prediction of the sudden water pollution event is still the key point of research of scholars at home and abroad.
Disclosure of Invention
The invention aims to provide a GPU acceleration-based pollutant transport high-resolution numerical simulation method, which solves the problem that a numerical model for researching the water pollutant transport motion law in the prior art is poor in operation efficiency and accuracy.
The invention adopts the technical scheme that a pollutant transport high-resolution numerical simulation method based on GPU acceleration is implemented according to the following steps:
step 2, reading the initial condition and the boundary condition, and replacing variables from the CPU to the GPU by adopting a finite volume method space discrete control equation;
step 3, dividing dry and wet calculation areas, and calculating the flux on a unit grid interface on a GPU by adopting an HLLC approximate Riemann solver;
step 4, calculating a source item of the control equation on the GPU;
step 5, calculating the time step length of iterative update according to the Curian number;
and 7, copying the hydraulic power element and the concentration information from the GPU video memory to a host memory and outputting.
The invention is also characterized in that:
the finite volume method space discrete control equation comprises a two-dimensional shallow water equation and a pollutant transport equation, and the vector form of the coupled conservation format is expressed as follows:
wherein the content of the first and second substances,
in formulas (1) and (2): t is time; x and y are horizontal and vertical coordinates, respectively; q is a variable vector; f and G are flux vectors in the x and y directions, respectively; q. q.sxSingle wide flow in the x direction, qx＝uh；qySingle wide flow in the y direction, qyVh; s is a source item vector; g is the gravitational acceleration, u and v are the flow velocities in the x and y directions, respectively; h is the water depth; z is a radical ofbIs the elevation of the bottom surface of the riverbed; cfIs the coefficient of friction of the bed surface, Cf＝gn2/h1/3N is a Manning coefficient; c is the average concentration of the vertical line of the pollutants; dxAnd DyRepresents the diffusion coefficients in the x and y directions; q. q.sinThe flow intensity of the point source discharge; cinIs the average concentration of the perpendicular line of material of the point source.
The initial conditions comprise flow data, Manning data and initial distribution positions of pollutants; the boundary conditions include a closed boundary and an open boundary.
The calculation process of step 3 is as follows:
step 3.1, integrating the formula (1) on any control body omega by adopting a finite volume method:
in the formula (3), omega is a control body; t is time; x and y are horizontal and vertical coordinates, respectively; q is a variable vector; f and G are flux vectors in the x and y directions, respectively; s is a source item vector;
step 3.2, by Gauss-Green formula will
The volume fraction of (a) is converted into a surface integral along the control volume boundary:
in the formula (4), Γ is a boundary of the control body, and n is a unit vector of the outer normal direction corresponding to the boundary Γ; f (q) is a flux vector; s (q) is a source item vector;
step 3.3, in a cartesian coordinate system, setting the line integral of f (q) n on the unit grid as an algebraic expression, specifically:
in the formula (5), k is the side length number of the unit grid i; lkThe length of the k-th edge of the unit grid; fk(qn) Interface fluxes for 4 interfaces east, west, south and north of the cell grid; n isk＝(nx，ny) An outer normal unit vector that is a side; q. q.snAnd the variables of the time n comprise water depth, flow rate and pollutant concentration.
When the HLLC approximate Riemann solver is adopted to calculate the interface flux of the unit grid, firstly, spatial second-order precision reconstruction is carried out on the interface variable of the unit grid to form a Riemann problem; and then, carrying out numerical reconstruction on variables on the left side and the right side of a unit grid interface by adopting an MUSCL state interpolation method, and limiting the gradient by combining a Min-Mod limiter so as to effectively inhibit numerical oscillation.
The calculation process of step 4 is as follows:
step 4.1, processing by using finite difference method
Item, update q of cell grid i from n to n +1 period to:
wherein in the formula (6), q ═ h, qx,qy,hC]TIs a variable vector comprising water depth, flow velocity, flow and pollutant concentration; n is the time; k is the side length number of the unit grid i; lkThe length of the k-th edge of the unit grid; delta t is a time step, needs to meet the CFL condition limit, and is 0.5; fk(qn) Is a flux vector; s is source term approximation and comprises a bottom slope source term and a friction source term;
a variable vector at the nth moment for the unit grid i;
a variable vector at the nth moment for the unit grid i;
step 4.2, obtaining unit grid i of second-order time precision in new time step by adopting two-step Runge-Kutta method
Are assigned as:
wherein the content of the first and second substances,
in the formula (7), in the formula,
is an intermediate variable value;
is a correction factor.
In step 4.1, the bottom slope source item is processed as follows; decomposing the integrals on the unit grids into integrals of the subunits, and converting the bottom slope source terms of the unit grids into the flux of the grid surfaces by assuming that the bed bottom height and the water level in the grids are linearly changed;
the calculation formula of the friction source term is as follows:
in formulas (8) and (9): q. q.sxSingle wide flow in the x direction, qx＝uh；qyIs a single wide flow in the y direction, qyVh; h is the water depth; u is the flow velocity in the x direction; v is the flow velocity in the y direction; sfxIs the friction resistance source item in the x direction; sfyThe friction resistance source item in the y direction; cfIs the coefficient of roughness of the bed surface, Cf＝gn2/h1/3(ii) a n is the Manning coefficient.
In step 5, the condition limit to be satisfied by the time step is specifically:
in the formula (10), Δ t is a time step; CFL is the Kuran number; delta xiAnd Δ yiThe length of the left side, the right side, the upper side and the lower side of the unit grid i in a Cartesian coordinate system respectively; u. ofiAnd viFlow velocities of the cell grid i in the x and y directions, respectively; h isiIs the water depth at the center of the unit grid i.
In step 6, the hydraulic power factor and the concentration information are water depth, flow velocity, flow and pollutant concentration respectively.
The invention has the beneficial effects that:
the invention relates to a pollutant transport high-resolution numerical simulation method based on GPU acceleration, which adopts high-resolution DEM data to describe the terrain of a simulation area, applies a finite volume method numerical solution of Godunov format to solve a two-dimensional shallow water equation and a convection diffusion equation in the surface runoff evolution process, adopts a second-order explicit Runge-Kutta method to update the time step length, calculates the water flux and the material flux by adopting an HLLC approximate Riemann solver capable of effectively capturing shock wave capability at a unit interface, processes a bottom slope source item by adopting a new bottom slope flux method, calculates frictional resistance by using a semi-implicit method with better stability to realize high-precision and steady confluence evolution process simulation, and simultaneously introduces a GPU parallel acceleration technology to improve the model simulation efficiency; the pollutant transport high-resolution numerical simulation method based on GPU acceleration has the advantages of stable calculation, high precision and simulation efficiency, capability of realizing accurate simulation of a water pollutant transport process, obtaining a more detailed water depth distribution map and a pollutant concentration distribution map, and capability of providing reliable theoretical basis and powerful data support for rapid evaluation and early warning of sudden water pollution accidents; according to the pollutant transport high-resolution numerical simulation method based on GPU acceleration, a GPU acceleration technology is applied to model calculation, the model calculation efficiency is improved, and early warning and evaluation can be rapidly and accurately carried out on water pollution events.
Drawings
FIG. 1 is a flow chart of GPU computation in a method for simulating pollutant transport high resolution numerical value based on GPU acceleration according to the invention;
FIG. 2 is a diagram showing the result of pure convection transport simulation of pollutants by the GPU-based accelerated pollutant transport high-resolution numerical simulation method of the present invention;
FIG. 2(a) is a distribution diagram of the calculated domain extent, initial conditions, boundary conditions and initial position and shape of contaminants
FIG. 2(b) is the concentration and profile distribution of the contaminant at t 5s
FIG. 2(c) is the concentration and profile distribution of the contaminant at t 10s
FIG. 3 is a graph showing the result of a pollutant concentration peak diffusion simulation in a still water flow field by using a pollutant transport high-resolution numerical simulation method based on GPU acceleration;
FIG. 3(a) shows a diffusion coefficient of 0m2Simulation result diagram of pollutant concentration peak diffusion in static water flow field at/s
FIG. 3(b) is a graph showing a diffusion coefficient of 0.1m2Simulation result diagram of pollutant concentration peak diffusion in static water flow field at/s
FIG. 3(c) is a graph showing a diffusion coefficient of 0.3m2Simulation result diagram of pollutant concentration peak diffusion in static water flow field at/s
FIG. 3(d) is a graph showing a diffusion coefficient of 0.5m2Simulation result diagram of pollutant concentration peak diffusion in static water flow field at/s
FIG. 4 is a graph showing the simulation result of the pollutant transport high-resolution numerical simulation method based on GPU acceleration on the Malpasset dam break flood evolution process;
FIG. 4(a) is a water depth distribution diagram of Malpaset dam break flood evolution at t-1000 s
FIG. 4(b) is a water depth distribution diagram of Malpaset dam break flood evolution at t 1800s
FIG. 5 is a graph showing the results of point source pollution source transport simulation driven by Malpaset dam break by the GPU acceleration-based high-resolution numerical simulation method for pollutant transport.
FIG. 5(a) is a pollutant distribution diagram of point source pollution source transportation driven by Malpaset dam break flood at t 900s
FIG. 5(b) is a pollutant distribution diagram of point source pollution source transportation driven by Malpaset dam break flood at t 1000s
FIG. 5(c) is a pollutant distribution diagram of point source pollution source transportation driven by Malpaset dam break flood at t 1800s
FIG. 5(d) is a pollutant distribution diagram of point source pollution source transportation driven by Malpaset dam break flood at t-3000 s
Detailed Description
The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.
A pollutant transport high-resolution numerical simulation method based on GPU acceleration is characterized by comprising the following steps:
step 2, reading the initial condition and the boundary condition, and replacing variables from the CPU to the GPU by adopting a finite volume method space discrete control equation;
the finite volume method space discrete control equation comprises a two-dimensional shallow water equation and a pollutant transport equation, and the vector form of the coupled conservation format is expressed as follows:
wherein the content of the first and second substances,
in formulas (1) and (2): t is time; x and y are horizontal and vertical coordinates, respectively; q is a variable vector; f and G are flux vectors in the x and y directions, respectively; q. q.sxSingle wide flow in the x direction, qx＝uh；qySingle wide flow in the y direction, qyVh; s is a source item vector; g is the gravitational acceleration, u and v are the flow velocities in the x and y directions, respectively; h is the water depth; z is a radical ofbIs the elevation of the bottom surface of the riverbed; cfIs the coefficient of friction of the bed surface, Cf＝gn2/h1/3N is a Manning coefficient; c is the average concentration of the vertical line of the pollutants; dxAnd DyRepresents the diffusion coefficients in the x and y directions; q. q.sinThe flow intensity of the point source discharge; cinThe average concentration of the perpendicular line of the material which is a point source;
the initial conditions comprise flow data, Manning data and initial distribution positions of pollutants; the boundary conditions include a closed boundary and an open boundary;
step 3, dividing dry and wet calculation areas, and calculating the flux on a unit grid interface on a GPU by adopting an HLLC approximate Riemann solver; as shown in fig. 1;
the calculation process of step 3 is as follows:
step 3.1, integrating the formula (1) on any control body omega by adopting a finite volume method:
in the formula (3), omega is a control body; t is time; x and y are horizontal and vertical coordinates, respectively; q is a variable vector; f and G are flux vectors in the x and y directions, respectively; s is a source item vector;
step 3.2, by Gauss-Green formula will
The volume fraction of (a) is converted into a surface integral along the control volume boundary:
in the formula (4), Γ is a boundary of the control body, and n is a unit vector of the outer normal direction corresponding to the boundary Γ; f (q) is a flux vector; s (q) is a source item vector;
step 3.3, in a cartesian coordinate system, setting the line integral of f (q) n on the unit grid as an algebraic expression, specifically:
in the formula (5), k is the side length number of the unit grid i; lkThe length of the k-th edge of the unit grid; fk(qn) Interface fluxes for 4 interfaces east, west, south and north of the cell grid; n isk＝(nx，ny) An outer normal unit vector that is a side; q. q.snVariables at time n, including depth, flow rate and contaminant concentration;
step 4, calculating a source item of the control equation on the GPU;
the calculation process of step 4 is as follows:
step 4.1, processing by using finite difference method
Item, update q of cell grid i from n to n +1 period to:
wherein in the formula (6), q ═ h, qx,qy,hC]TIs a variable vector comprising water depth, flow velocity, flow and pollutant concentration; n is the time; k is the side length number of the unit grid i; lkThe length of the k-th edge of the unit grid; delta t is a time step, needs to meet the CFL condition limit, and is 0.5; fk(qn) Is a flux vector; s is source term approximation and comprises a bottom slope source term and a friction source term;
a variable vector at the nth moment for the unit grid i;
a variable vector at the nth moment for the unit grid i;
step 4.2, obtaining unit grid i of second-order time precision in new time step by adopting two-step Runge-Kutta method
Are assigned as:
wherein the content of the first and second substances,
in the formula (7), in the formula,
is an intermediate variable value;
is a correction factor;
in step 4.1, the bottom slope source item is processed as follows; decomposing the integrals on the unit grids into integrals of the subunits, and converting the bottom slope source terms of the unit grids into the flux of the grid surfaces by assuming that the bed bottom height and the water level in the grids are linearly changed;
the calculation formula of the friction source term is as follows:
in formulas (8) and (9): q. q.sxSingle wide flow in the x direction, qx＝uh；qyIs a single wide flow in the y direction, qyVh; h is the water depth; u is the flow velocity in the x direction; v is the flow velocity in the y direction; sfxIs the friction resistance source item in the x direction; sfyThe friction resistance source item in the y direction; cfIs the coefficient of roughness of the bed surface, Cf＝gn2/h1/3(ii) a n is a Manning coefficient;
setting a minimum water depth of 1 × 10 when updating the pollutant concentration-6m, when the water depth of the unit grid is less than the minimum water depth, the concentration of the pollutants is 0;
step 5, calculating the time step length of iterative update according to the Curian number;
the condition limit that the time step needs to meet is specifically as follows:
in the formula (10), Δ t is a time step; CFL is the Kuran number; delta xiAnd Δ yiRespectively cartesian coordinate systemThe left and right sides and the upper and lower sides of the middle unit grid i are long; u. ofiAnd viFlow velocities of the cell grid i in the x and y directions, respectively; h isiIs the water depth at the center of the unit grid i;
the hydraulic element and the concentration information are water depth, flow velocity, flow and pollutant concentration respectively;
and 7, copying the hydraulic power element and the concentration information from the GPU video memory to a host memory and outputting.
When the HLLC approximate Riemann solver is adopted to calculate the interface flux of the unit grid, firstly, spatial second-order precision reconstruction is carried out on the interface variable of the unit grid to form a Riemann problem; then, carrying out numerical reconstruction on variables on the left side and the right side of a unit grid interface by adopting an MUSCL state interpolation method, and limiting the gradient by combining a Min-Mod limiter so as to effectively inhibit numerical oscillation; the method specifically comprises the following steps:
introducing a non-negative water depth and concentration reconstruction method into an HLLC approximate Riemann solver, and calculating the water quantity, momentum and pollutant transport flux on a unit grid interface;
the calculation formula is as follows:
in the formula (11), F*LAnd F*RAre respectively F*L＝(F*1,F*2,vLF*1,CLF*1)TAnd F*L＝(F*1,F*2,vLF*1,CLF*1)TWherein v and C are the tangential velocity and the concentration of the substance at the Riemann problem, respectively, and remain unchanged at the left and right waves; f*LAnd F*RRespectively the flux in the middle of the left and right grids of the interface; fLAnd FRRespectively the flux of the left and right sides of the grid unit interface; wave velocity SL、SMAnd SRRespectively the left, middle and right wave velocities; subscripts "L" and "R" are the left and right states of the cell grid interface, respectively;
FL＝F(qL) And FR＝F(qR) Calculating according to Riemann states at two sides of the boundary;
in the formula (12), qRAnd q isLVariable vectors which are respectively left and right of a unit grid interface comprise water depth, flow velocity, flow and pollutant concentration;
in the formula (13), hLAnd hRThe water depth of the left side and the right side of the unit grid interface respectively; g is the gravitational acceleration generated by the gravity of the earth; u. ofRAnd uLFlow rates on the left and right sides of the cell grid interface respectively; h is*And u*Respectively the water depth and flow speed in the middle of the grid unit and the intermediate variable h*And u*Calculated from the following formula:
using the gaussian divergence theorem, the pollutant diffusion flux is expressed as:
in the formula (15), FdifDiffusing flux for the contaminant; h isLAnd hRThe water depth of the left side and the right side of the unit grid interface respectively;
and
the original slopes of the variables, respectively; n isxAnd nyUnit vectors representing x and y directions in a coordinate system, respectively;
realizing the coupling solution of water flow and pollutant transport flux by adopting an HLLC approximate Riemann solver, and calculating flux Fk(q)·nk＝(F1,F2,F3,F4+Fdif)TIterative calculations were performed for equation (6).
And (3) experimental verification:
(1) pure convection transport simulation of pollutants
The calculated domain size of the contaminants was 100m × 20m, the grid size was set to 0.5m × 0.5m, the initial water depth was 0.2m, the Manning coefficient was assumed to be 0.018, the contaminants were distributed in the region of [45m, 50m ] × [7m, 13m ] at the initial moment, the concentration value was 1mg/L, the right boundary was the open boundary of the free outflow, and the remaining boundary conditions were the solid boundary, as shown in FIG. 2 (a); fig. 2(b) and (c) show the calculated results (profile distribution) of the contaminant concentration and morphology at different times.
As can be seen from the graphs (a), (b) and (c) of FIG. 2, the pollutants do not have numerical damping or false numerical oscillation phenomenon in the process of moving with the water flow, and the numerical dissipation is small. The method has higher precision on the treatment of the pollutants by the convection item.
(2) Simulation of pollutant concentration peak diffusion in static water flow field
Simulating a pollutant diffusion process by adopting a Gaussian concentration distribution peak in a calculation range of 75m multiplied by 30m, wherein the periphery of the calculation range is a closed boundary, and the water level is constant at 1 m; respectively calculating the diffusion coefficient to be 0m by simulation2/s、0.1m2/s、0.3m2/s、0.5m2The diffusion process of the contaminant concentration peak at/s is shown in FIGS. 3(a), (b), (c) and (d).
the comparison between the numerical solution and the analytic solution when t is 10s is shown in fig. 3, and it can be seen from fig. 3 that, in the whole process of the finite volume method space discrete control equation calculation, the pollutant flow rate is always kept at 0, no spurious flow occurs, and the pollutant peak concentration is kept at x037.5m, showing that the finite volume method space discrete control equation has better still water harmony. In addition, as the diffusion coefficient increases, the peak concentration of the contaminant goes to the next stepThe concentration peak value tends to be flat, and the numerical solution and the analytic solution are well matched, which shows that the invention has higher precision for the pollutant concentration peak diffusion treatment in the static water flow field.
(3) The pollutant transport high-resolution numerical simulation method based on GPU acceleration is adopted to simulate the Malpasset dam break flood evolution process
The upstream water depth of the dam is 100m, the downstream initial water depth is 0m, the downstream outflow port is an open boundary, and the rest boundaries are solid closed boundaries. Assuming that a pollutant source discharge port at (x-9660 m, y-3074 m) releases pollutants, the discharge concentration and intensity are 1mg/L and 1m/s respectively, and the diffusion coefficients in the x and y directions are both 0.5m2And s. Taking the coefficient of the Manning roughness of the river channel to be 0.033s/m1 /3. The terrain data adopts DEM data with the resolution of 10m, and the total number of the DEM data is 155 ten thousand square grid units.
When t is 0s, the flood rapidly progresses to the lower part or the flood area of the downstream terrain, and fig. 4(a) and (b) show the water depth distribution of the flood at the time t is 1000s and 1800s, and the water flow velocity is fast and the flood is distributed in the riverbed. Fig. 5(a), (b), (c), (d) show the pollutant distribution at the time t 900s, 1000s, 1800s and 3000s, respectively, and when the water flow reaches the point source position at the time t 900s, the point source starts to release pollutants, and the water flow rapidly and pollutants are mixed to evolve downstream. As can be seen from fig. 5, the pollutants always migrate downstream along with the water flow in the river and are matched with the dry-wet interface, which indicates that the point source pollutant discharge can cause large-scale water pollution and rapidly affect the downstream water, and conforms to the physical law. The practicability of the pollutant transport high-resolution numerical simulation method based on GPU acceleration is verified.
Claims (9)
1. A pollutant transport high-resolution numerical simulation method based on GPU acceleration is characterized by comprising the following steps:
step 1, reading a DEM terrain file consisting of high-resolution unit grids to obtain an assigned variable;
step 2, reading the initial condition and the boundary condition, and replacing the variable from the CPU to the GPU by adopting a finite volume method space discrete control equation;
step 3, dividing dry and wet calculation areas, and calculating the flux on a unit grid interface on a GPU by adopting an HLLC approximate Riemann solver;
step 4, calculating a source item of the control equation on the GPU;
step 5, calculating the time step length of iterative update according to the Curian number;
step 6, repeating the step 3-5 to obtain the hydraulic elements and concentration information of all unit grids;
and 7, copying the hydraulic element and the concentration information from the GPU video memory to a host memory and outputting.
2. The GPU-acceleration-based pollutant transport high-resolution numerical simulation method according to claim 1, wherein the finite volume method space discrete control equation comprises a two-dimensional shallow water equation and a pollutant transport equation, and the vector form of a coupled conservation format is represented as follows:
wherein the content of the first and second substances,
in formulas (1) and (2): t is time; x and y are horizontal and vertical coordinates, respectively; q is a variable vector; f and G are flux vectors in the x and y directions, respectively; q. q.sxSingle wide flow in the x direction, qx＝uh；qySingle wide flow in the y direction, qyVh; s is a source item vector; g is the gravitational acceleration, u and v are the flow velocities in the x and y directions, respectively; h is the water depth; z is a radical ofbIs the elevation of the bottom surface of the riverbed; cfIs the coefficient of friction of the bed surface, Cf＝gn2/h1/3N is a Manning coefficient; c is the average concentration of the vertical line of the pollutants; dxAnd DyIndicating diffusion systems in the x and y directionsCounting; q. q.sinThe flow intensity of the point source discharge; cinIs the average concentration of the perpendicular line of material of the point source.
3. The GPU-acceleration-based pollutant transport high-resolution numerical simulation method according to claim 2, wherein the initial conditions comprise flow data, Manning data and initial pollutant distribution positions; the boundary conditions include a closed boundary and an open boundary.
4. The method for high-resolution numerical simulation of pollutant transport based on GPU acceleration as claimed in claim 2, characterized in that the calculation process of step 3 is as follows:
step 3.1, integrating the formula (1) on any control body omega by adopting a finite volume method:
in the formula (3), omega is a control body; t is time; x and y are horizontal and vertical coordinates, respectively; q is a variable vector; f and G are flux vectors in the x and y directions, respectively; s is a source item vector;
step 3.2, by Gauss-Green formula will
The volume fraction of (a) is converted into a surface integral along the control volume boundary:
in the formula (4), Γ is a boundary of the control body, and n is a unit vector of the outer normal direction corresponding to the boundary Γ; f (q) is a flux vector; s (q) is a source item vector;
step 3.3, in a cartesian coordinate system, setting the line integral of f (q) n on the unit grid as an algebraic expression, specifically:
in the formula (5), k is the side length number of the unit grid i; lkThe length of the k-th edge of the unit grid; fk(qn) Interface fluxes for 4 interfaces east, west, south and north of the cell grid; n isk＝(nx，ny) An outer normal unit vector that is a side; q. q.snAnd the variables of the time n comprise water depth, flow rate and pollutant concentration.
5. The GPU acceleration-based pollutant transport high-resolution numerical simulation method is characterized in that when an HLLC approximation Riemann solver is used for calculating unit grid interface flux, spatial second-order precision reconstruction is carried out on unit grid interface variables to form a Riemann problem; and then, carrying out numerical reconstruction on variables on the left side and the right side of a unit grid interface by adopting an MUSCL state interpolation method, and limiting the gradient by combining a Min-Mod limiter so as to effectively inhibit numerical oscillation.
6. The GPU-acceleration-based pollutant transport high-resolution numerical simulation method according to claim 4, characterized in that the calculation process of step 4 is as follows:
step 4.1, processing by using finite difference method
Item, update q of cell grid i from n to n +1 period to:
wherein in the formula (6), q ═ h, qx,qy,hC]TIs a variable vector comprising water depth, flow velocity, flow and pollutant concentration; n is the time; k is the side length number of the unit grid i; lkLength of edge of k-th edge of unit grid(ii) a Delta t is a time step, needs to meet the CFL condition limit, and is 0.5; fk(qn) Is a flux vector; s is source term approximation and comprises a bottom slope source term and a friction source term;
a variable vector at the nth moment for the unit grid i;
a variable vector at the nth moment for the unit grid i;
step 4.2, obtaining unit grid i of second-order time precision in new time step by adopting two-step Runge-Kutta method
Are assigned as:
wherein the content of the first and second substances,
in the formula (7), in the formula,
is an intermediate variable value;
is a correction factor.
7. The GPU-acceleration-based pollutant transport high-resolution numerical simulation method according to claim 6, characterized in that in step 4.1, the bottom slope source term is processed as follows; decomposing the integrals on the unit grids into integrals of the subunits, and converting the bottom slope source terms of the unit grids into the flux of the grid surfaces by assuming that the bed bottom height and the water level in the grids are linearly changed;
the calculation formula of the friction source term is as follows:
in formulas (8) and (9): q. q.sxSingle wide flow in the x direction, qx＝uh；qyIs a single wide flow in the y direction, qyVh; h is the water depth; u is the flow velocity in the x direction; v is the flow velocity in the y direction; sfxIs the friction resistance source item in the x direction; sfyThe friction resistance source item in the y direction; cfIs the coefficient of roughness of the bed surface, Cf＝gn2/h1/3(ii) a n is the Manning coefficient.
8. The method according to claim 6, wherein in step 5, the time step satisfies the following condition:
in the formula (10), Δ t is a time step; CFL is the Kuran number; delta xiAnd Δ yiThe length of the left side, the right side, the upper side and the lower side of the unit grid i in a Cartesian coordinate system respectively; u. ofiAnd viFlow velocities of the cell grid i in the x and y directions, respectively; h isiIs the water depth at the center of the unit grid i.
9. The GPU-acceleration-based pollutant transport high-resolution numerical simulation method according to claim 1, wherein in the step 6, the hydraulic element and the concentration information are water depth, flow velocity, flow rate and pollutant concentration respectively.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011131148.0A CN112257313B (en) | 2020-10-21 | 2020-10-21 | GPU acceleration-based high-resolution numerical simulation method for pollutant transportation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011131148.0A CN112257313B (en) | 2020-10-21 | 2020-10-21 | GPU acceleration-based high-resolution numerical simulation method for pollutant transportation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112257313A true CN112257313A (en) | 2021-01-22 |
| CN112257313B CN112257313B (en) | 2024-05-14 |
Family
ID=74263877
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011131148.0A Active CN112257313B (en) | 2020-10-21 | 2020-10-21 | GPU acceleration-based high-resolution numerical simulation method for pollutant transportation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112257313B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112836872A (en) * | 2021-01-29 | 2021-05-25 | 西安理工大学 | Multi-GPU-based pollutant convection diffusion equation high-performance numerical solution method |
| CN117010232A (en) * | 2023-06-27 | 2023-11-07 | 西安理工大学 | Urban non-point source pollution whole process high-resolution simulation method based on GPU (graphic processing Unit) acceleration technology |
| CN117933148A (en) * | 2024-03-25 | 2024-04-26 | 中国空气动力研究与发展中心计算空气动力研究所 | Method and system for determining iteration times of volume fraction equation based on free interface |
Citations (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009044282A2 (en) * | 2007-10-04 | 2009-04-09 | Mental Images Gmbh | Quasi-monte carlo light transport simulation by efficient ray tracing |
| KR20110103242A (en) * | 2010-03-12 | 2011-09-20 | 서울대학교산학협력단 | Method for analyzing a longitudinal dispersion mechanism of pollutant in a river using time-splitting method |
| TW201201035A (en) * | 2010-06-24 | 2012-01-01 | jin-song Lai | Two-dimensional hydraulic sediment transport simulation calculation method and computer program product |
| KR20120002206A (en) * | 2010-06-30 | 2012-01-05 | 경북대학교 산학협력단 | Numerical analysis method of 2-dimentional flow and inundation |
| KR20130022800A (en) * | 2011-08-26 | 2013-03-07 | 경북대학교 산학협력단 | Numerical analysis method using 2-dimensional finite volume model applicable to mixed meshes |
| CN106815448A (en) * | 2017-02-07 | 2017-06-09 | 长江水资源保护科学研究所 | A kind of river attenuation type pollutant analogy method |
| CN107451372A (en) * | 2017-08-09 | 2017-12-08 | 中国水利水电科学研究院 | A kind of flood of a mountain area numerical simulation method that kinematic wave is combined with dynamic wave |
| CN108021780A (en) * | 2018-01-24 | 2018-05-11 | 国电南瑞科技股份有限公司 | A kind of mountain torrents dynamic emulation method based on random unstrctured grid model |
| KR101885350B1 (en) * | 2017-07-25 | 2018-08-03 | 주식회사 비전아이티 | System for integrated managment of air pollution |
| KR101912627B1 (en) * | 2017-05-30 | 2018-10-30 | 에스지에이블록체인 주식회사 | Method for Integration Visualizing GIS based Runoff-Hydraulic Model Analysis result |
| CN109190261A (en) * | 2018-09-07 | 2019-01-11 | 中国水利水电科学研究院 | A kind of flood risk analysis method that one-dimensional river network generalization is coupled with one two dimension of high-performance |
| CN109857543A (en) * | 2018-12-21 | 2019-06-07 | 中国地质大学（北京） | A kind of streamline simulation accelerated method calculated based on the more GPU of multinode |
| US20190196058A1 (en) * | 2017-12-21 | 2019-06-27 | Amit Kumar | Method and System for Modeling in a Subsurface Region |
| CN110275733A (en) * | 2019-06-27 | 2019-09-24 | 上海交通大学 | The GPU parallel acceleration method of phonon Boltzmann equation is solved based on finite volume method |
| CN110765694A (en) * | 2019-11-21 | 2020-02-07 | 华南理工大学 | Urban surface water flow numerical simulation method based on simplified shallow water equation set |
| CN111222240A (en) * | 2020-01-06 | 2020-06-02 | 中国人民解放军国防科技大学 | Thermochemical unbalanced flow field data calculation method and device accelerated by GPU |
| CN111767684A (en) * | 2020-06-30 | 2020-10-13 | 西安理工大学 | Optimized friction resistance source term implicit format two-dimensional shallow water equation modeling method |
| CN111768502A (en) * | 2020-07-08 | 2020-10-13 | 西安理工大学 | Non-structural grid two-dimensional flood simulation system based on GPU acceleration technology |
-
2020
- 2020-10-21 CN CN202011131148.0A patent/CN112257313B/en active Active
Patent Citations (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2009044282A2 (en) * | 2007-10-04 | 2009-04-09 | Mental Images Gmbh | Quasi-monte carlo light transport simulation by efficient ray tracing |
| KR20110103242A (en) * | 2010-03-12 | 2011-09-20 | 서울대학교산학협력단 | Method for analyzing a longitudinal dispersion mechanism of pollutant in a river using time-splitting method |
| TW201201035A (en) * | 2010-06-24 | 2012-01-01 | jin-song Lai | Two-dimensional hydraulic sediment transport simulation calculation method and computer program product |
| KR20120002206A (en) * | 2010-06-30 | 2012-01-05 | 경북대학교 산학협력단 | Numerical analysis method of 2-dimentional flow and inundation |
| KR20130022800A (en) * | 2011-08-26 | 2013-03-07 | 경북대학교 산학협력단 | Numerical analysis method using 2-dimensional finite volume model applicable to mixed meshes |
| CN106815448A (en) * | 2017-02-07 | 2017-06-09 | 长江水资源保护科学研究所 | A kind of river attenuation type pollutant analogy method |
| KR101912627B1 (en) * | 2017-05-30 | 2018-10-30 | 에스지에이블록체인 주식회사 | Method for Integration Visualizing GIS based Runoff-Hydraulic Model Analysis result |
| KR101885350B1 (en) * | 2017-07-25 | 2018-08-03 | 주식회사 비전아이티 | System for integrated managment of air pollution |
| CN107451372A (en) * | 2017-08-09 | 2017-12-08 | 中国水利水电科学研究院 | A kind of flood of a mountain area numerical simulation method that kinematic wave is combined with dynamic wave |
| US20190196058A1 (en) * | 2017-12-21 | 2019-06-27 | Amit Kumar | Method and System for Modeling in a Subsurface Region |
| CN108021780A (en) * | 2018-01-24 | 2018-05-11 | 国电南瑞科技股份有限公司 | A kind of mountain torrents dynamic emulation method based on random unstrctured grid model |
| CN109190261A (en) * | 2018-09-07 | 2019-01-11 | 中国水利水电科学研究院 | A kind of flood risk analysis method that one-dimensional river network generalization is coupled with one two dimension of high-performance |
| CN109857543A (en) * | 2018-12-21 | 2019-06-07 | 中国地质大学（北京） | A kind of streamline simulation accelerated method calculated based on the more GPU of multinode |
| CN110275733A (en) * | 2019-06-27 | 2019-09-24 | 上海交通大学 | The GPU parallel acceleration method of phonon Boltzmann equation is solved based on finite volume method |
| CN110765694A (en) * | 2019-11-21 | 2020-02-07 | 华南理工大学 | Urban surface water flow numerical simulation method based on simplified shallow water equation set |
| CN111222240A (en) * | 2020-01-06 | 2020-06-02 | 中国人民解放军国防科技大学 | Thermochemical unbalanced flow field data calculation method and device accelerated by GPU |
| CN111767684A (en) * | 2020-06-30 | 2020-10-13 | 西安理工大学 | Optimized friction resistance source term implicit format two-dimensional shallow water equation modeling method |
| CN111768502A (en) * | 2020-07-08 | 2020-10-13 | 西安理工大学 | Non-structural grid two-dimensional flood simulation system based on GPU acceleration technology |
Non-Patent Citations (2)
| Title |
|---|
| 侯精明;李钰茜;同玉;刘菲菲;马利平;梁行行;: "植草沟径流调控效果对关键设计参数的响应规律模拟", 水科学进展, no. 01, 31 January 2020 (2020-01-31) * |
| 郭敏鹏;侯精明;付德宇;康永德;石宝山;李昌镐;: "面源污染输移过程高性能数值模拟方法", 环境工程, no. 08, 15 August 2020 (2020-08-15) * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112836872A (en) * | 2021-01-29 | 2021-05-25 | 西安理工大学 | Multi-GPU-based pollutant convection diffusion equation high-performance numerical solution method |
| CN112836872B (en) * | 2021-01-29 | 2023-08-18 | 西安理工大学 | Multi-GPU-based high-performance numerical solution method for pollutant convection diffusion equation |
| CN117010232A (en) * | 2023-06-27 | 2023-11-07 | 西安理工大学 | Urban non-point source pollution whole process high-resolution simulation method based on GPU (graphic processing Unit) acceleration technology |
| CN117010232B (en) * | 2023-06-27 | 2024-06-14 | 西安理工大学 | Urban non-point source pollution whole process high-resolution simulation method based on GPU (graphic processing Unit) acceleration technology |
| CN117933148A (en) * | 2024-03-25 | 2024-04-26 | 中国空气动力研究与发展中心计算空气动力研究所 | Method and system for determining iteration times of volume fraction equation based on free interface |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112257313B (en) | 2024-05-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112257313A (en) | Pollutant transport high-resolution numerical simulation method based on GPU acceleration | |
| Smyth | A review of Computational Fluid Dynamics (CFD) airflow modelling over aeolian landforms | |
| Shen et al. | Modeling study of the influences of tide and stratification on age of water in the tidal James River | |
| CN110008509B (en) | Method for analyzing internal solitary wave action force characteristics under consideration of background flow field | |
| Xiao et al. | Numerical modeling of wave–current forces acting on horizontal cylinder of marine structures by VOF method | |
| Ravensbergen et al. | A variational multiscale framework for atmospheric turbulent flows over complex environmental terrains | |
| Wei et al. | Impacts of rainfall intensity and urbanization on water environment of urban lakes | |
| CN108021780A (en) | A kind of mountain torrents dynamic emulation method based on random unstrctured grid model | |
| CN115691709B (en) | Simulation method for researching transportation rule of micro plastic particles in water flow intersection area | |
| Balas et al. | An implicit three‐dimensional numerical model to simulate transport processes in coastal water bodies | |
| CN110955998A (en) | GIS-based large-range debris flow numerical simulation and numerical processing method | |
| CN111428401A (en) | Method for simulating damming process of barrier lake | |
| CN111881596A (en) | Oil spill pollution source reverse-time tracking simulation method based on Lagrange interpolation | |
| Khosronejad et al. | Large-eddy simulation of flash flood propagation and sediment transport in a dry-bed desert stream | |
| CN110990926B (en) | Urban surface building hydrodynamic simulation method based on area correction rate | |
| Audusse et al. | Multilayer Saint-Venant equations over movable beds | |
| US11868690B1 (en) | Method, device, electronic equipment and medium for analyzing disaster prevention and mitigation effectiveness of ecological seawall | |
| Sheng et al. | Modeling coastal currents and sediment transport | |
| Ding et al. | Spatial-temporal variation of permanganate index under accident conditions during different water periods in Zigui drinking water source area, Three Gorges reservoir area, China | |
| CN115130935B (en) | Ecological seawall disaster prevention and reduction efficiency analysis method and device, electronic equipment and medium | |
| Zhang et al. | Numerical study of hydrodynamic and solute transport with discontinuous flows in coastal water | |
| Dong et al. | Simulation of unidirectional propagating wave trains in deep water using a fully non-hydrostatic model | |
| Zeng et al. | Modeling flows over gravel beds by a drag force method and a modified S–A turbulence closure | |
| Tang et al. | Two-dimensional water environment numerical simulation research based on EFDC in Mudan River, Northeast China | |
| Sun et al. | Analysis of numerical factors affecting large eddy simulation of pollutant diffusion around buildings |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |