GPU (Graphics Processing Unit) accelerated water environment simulation method considering inorganic state substances

Abstract

The invention discloses a GPU (Graphics Processing Unit) accelerated water environment simulation method considering inorganic state substances, which comprises the following steps of: firstly, processing terrain scatter data by utilizing Arcgis to form a high-resolution DEM (Digital Elevation Model) terrain to represent a real terrain of a river, lake and sea area, and taking the real terrain as an input condition of terrain data; and then a hydrodynamic model is coupled with the hydrodynamic-substance transport coupling equation, and finally a transport process distribution diagram of each environmental factor is output to reflect a transport rule of each environmental factor in the water body, and the eutrophication condition of the water body is evaluated. According to the method, refined simulation of the transport rule of each environment variable and the mutual transformation relation between the variables in the complex river, lake and sea area water can be realized, the eutrophication condition of the river, lake and sea area water can be accurately and effectively evaluated, and an important simulation tool can be provided for effective prediction and prevention of water environment pollution.

Description

Technical Field

[0001] The invention belongs to the technical field of water environment simulation in rivers, lakes and seas, and in particular relates to a GPU-accelerated water environment simulation method taking inorganic substances into consideration. Background Art

[0002] Due to the rapid development of industrialization and urbanization, a large amount of industrial and agricultural wastewater and domestic sewage has been generated. These sewages contain rich nitrogen, phosphorus and other substances. If they are not treated adequately, they will be directly or indirectly discharged into rivers, lakes, seas, reservoirs and canal water bodies. In the long run, water pollution is difficult to recover, aquatic organisms die, and even pose a direct threat to human health. In order to effectively carry out water pollution prevention and control, water environment models are favored by many scholars. However, the early simple BOD-DO water quality model and the modified multi-component water quality model cannot meet the simulation needs of river water environment. It is necessary to construct a water environment model characterized by multidimensional, multi-media and dynamic simulation to cope with complex river pollution. However, there are complex transformation relationships between the environmental factors in the model, which leads to complex calculations and long simulation time. Based on this, a GPU-accelerated water environment simulation method considering inorganic substances is constructed. Based on the transformation equilibrium equation of the state variables of each substance, the transport law of various environmental variables in the water body of rivers, lakes and seas and their mutual transformation relationship are accurately and efficiently simulated, so as to achieve accurate simulation and effective evaluation of the eutrophication of complex water bodies in rivers, lakes and seas. Summary of the invention

[0003] The purpose of this invention is to provide a GPU-accelerated water environment simulation method that takes inorganic matter into consideration, which can realize the refined simulation of the transport laws of various environmental variables in complex river, lake and sea water bodies and the mutual transformation relationship between variables, so as to accurately and effectively evaluate the eutrophication status of river, lake and sea water bodies, and provide an important simulation tool for the effective prediction and prevention of water environment pollution.

[0004] The technical solution adopted by the present invention is a GPU accelerated water environment simulation method considering inorganic substances, which is specifically implemented according to the following steps:

[0005] Step 1: Use ArcGIS to process the scattered terrain data to form a high-resolution DEM terrain representation of the real terrain of rivers, lakes and sea areas as the input condition of terrain data;

[0006] Step 2: discretize the above DEM terrain data through rectangular grids, and assign values ​​to relevant parameters according to known boundary conditions;

[0007] Step 3: Construct the hydrodynamic-material transport coupling equation;

[0008] Step 4, coupling the hydrodynamic model with the above hydrodynamic-material transport coupling equation;

[0009] Step 5: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h, using the split operator method to transform the equilibrium equation of each material state variable and express them one by one;

[0010] Step 6: Use the improved hydrostatic reconstruction method to meet the balance between the flux term and the bottom slope source term in hydrostatic simulation on complex grid terrain;

[0011] Step 7: Use the MUSCL second-order accuracy and the second-order Runge-Kutta method to achieve the time-space second-order accuracy of the model;

[0012] Step 8: Output the distribution map of the transport process of each environmental factor to reflect the transport law of each environmental factor in the water body and evaluate the eutrophication status of the water body.

[0013] The present invention is also characterized in that:

[0014] The hydrodynamic-mass transport coupling equation in step 3 is expressed as follows:

[0015]

[0016] In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h can be expressed as a biochemical reaction term, with the unit of kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

[0017] The vector form of the two-dimensional hydrodynamic-mass transport coupling equations in step 4 is as follows:

[0018]

[0019] in,

[0020]

[0021] Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, temperature diffusion term and temperature source term; z b is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average concentration of the substance along the water depth, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s;

[0022] S in source and sink terms sc h is the biochemical reaction term, and its unit is kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

[0023] Step 4 is implemented according to the following steps:

[0024] Step 4.1: In the control unit grid i, the central finite volume method is used to integrate the two-dimensional hydrodynamic-material transport coupling control equations. Based on the Gaussian divergence theorem, the area integral of the flux term is converted into a line integral, which is:

[0025]

[0026] Where Ω is the volume of control volume i; Г is the boundary of control volume i; n is the unit vector in the direction of the external normal corresponding to the control volume boundary Г; F(Q)·n is the external normal convection term flux vector on the control volume boundary Г; S is the source term vector, including the bottom slope source term, friction source term, diffusion term and biochemical reaction term;

[0027] The flux F(Q)·n line integral is:

[0028]

[0029] Where Ω is the volume of the i-th control unit grid, in m 3 ; Г is the boundary of the i-th control unit grid; n is the unit vector in the direction of the external normal corresponding to the boundary Г; S b represents the bottom slope source flux; S f represents the friction source flux; k is the number of the control unit grid edge; lk is the length of the kth side of the i-th control unit grid;

[0030] The flux term is solved using the HLLC Riemann solver:

[0031] The wave speed is calculated as follows:

[0032]

[0033] Among them, the intermediate variable h * and According to the calculation formula,

[0034]

[0035]

[0036] The calculation result of F(Q)·n is:

[0037]

[0038] Among them, the control volume unit grid interface intermediate flux F * for,

[0039]

[0040] Control volume unit grid interface flux intermediate wave sides, and for,

[0041]

[0042] In the formula, F L and F R are the fluxes on the left and right sides of the unit grid interface respectively; F * , and are the left wave flux, middle wave flux and right wave flux on both sides of the intermediate wave of the control volume unit grid interface flux respectively; u || is the tangential velocity along the interface, u || =-un y +vn x ,u ||L and u ||R are the left and right tangential velocities of the control volume unit grid interface, in m 3 / s;n x and n y are the indicated directions along x and y, respectively; u and v are the single-width flow rates along x and y, respectively, in m / s; variable Q ⊥ =[h,q x n x +qy n y ];q x and q y are the single width flow along x and y, in m 2 / s;S L , S M , S R They are left wave speed, middle wave speed, and right wave speed, in m / s; u ⊥ is the velocity perpendicular to the grid unit interface, in m / s; u ⊥ =un x +vn y ;h * and is the intermediate variable; g is the acceleration due to gravity, which is 9.81 and the unit is m / s 2 ; C s is the average concentration of the substance along the water depth, kg / m 3 ;

[0043] Step 4.2: Treat the diffusion term as a source term and use the second-order central difference method to discretize the diffusion term. The material transport diffusion term is expressed as:

[0044]

[0045] According to the second-order central difference method, the first-order derivative and second-order derivative of the material transport diffusion term are solved, and for any two-dimensional unit grid center (i, j) we get:

[0046]

[0047] Thus, the discrete form of the material transport diffusion term is obtained as follows:

[0048]

[0049] Step 4.3: The bottom slope source term is solved using the bottom slope flux method, and the vector form of the bottom slope flux at a certain side f is obtained as:

[0050]

[0051] In the formula, F sf (Q) is the bottom slope flux at a certain side f; n f is the unit external normal vector of the surface f; n fx and n fy Yes f and (the component of the unit external normal vector of face f), h L and z bL are the water depth and bottom elevation at the center of the grid cell, respectively. With z MLThe water depth and bottom elevation at the edge of the grid. g is the acceleration due to gravity in m / s 2 ;

[0052] Step 4.4: The friction source term is solved using the explicit and implicit friction method. The calculation formula is:

[0053]

[0054] Among them, the calculation formula for single-width flow update to time step n+1 is:

[0055]

[0056] In the formula, S fx and S fy S f Friction source term in x and y directions; C f is the bed surface roughness coefficient, C f =gn 2 / h 1 / 3 , n is the Manning coefficient, s / m 1 / 3 ; h is the water depth, in meters; q x ,q y are the average flow velocities in the x and y directions, in m 2 / s;q x n+1 ,q y n+1 They are the single-width flow rates in the x and y directions at time step n+1, in m 2 / s;q x n ,q y n They are the single-width flow rates in the x and y directions at time step n, in m 2 / s; Δt is the time step, in seconds; Δq x is the difference in single-width flow in the x direction at time step n and time step n+1, that is, Δq x =q x n+1 -q x n , unit is m 2 / s; In order to avoid the denominator being zero, θ is taken as 1.0e -12 .

[0057] Step 5 is implemented according to the following steps:

[0058] Step 5.1, dissolved oxygen balance includes dissolved oxygen DO, biochemical oxygen demand BOD and chemical oxygen demand COD; dissolved oxygen DO concentration includes aquatic plant respiration, nitrification, biodegradation and sediment or silt oxygen consumption and photosynthesis oxygen production;

[0059] Step 5.2: The nitrogen cycle process includes ammonia nitrogen NH 3 、Nitrate nitrogen NO 3 and nitrite nitrogen NO 2 Three forms;

[0060] Step 5.3: Phosphate PO in the Phosphorus Cycle 4 The concentration includes the consumption by aquatic plants and bacterial microorganisms, the release from organic matter degradation reactions and anoxic biological respiration;

[0061] Step 5.4: Taking the growth dynamics of phytoplankton as the main line, the ratio of the three main elements C:N:P reflects the conversion relationship between phytoplankton, i.e. chlorophyll α, and nutrients in the environment.

[0062] Step 6 is implemented according to the following steps:

[0063] Step 6.1: Based on the hydrostatic reconstruction method, the bottom slope elevation at the center point M of the dry-wet unit boundary edge f is taken as:

[0064]

[0065] In the formula, are the bottom slope elevations on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m;

[0066] Reconstruct the water depth on both sides and make the water level on the right side of the dry-wet interface the same as the elevation of the reconstructed bottom slope, that is:

[0067]

[0068] In the formula, They are the water levels on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; are the water depths on the left and right sides after reconstruction, in meters; z bM is the bottom slope elevation after reconstruction of the center point M of the dry-wet unit boundary edge f, in meters;

[0069] If expressed in terms of water level, the reconstructed water level is:

[0070]

[0071] In the formula, are the water levels after reconstruction, in m;

[0072] Step 6.2: Based on the above, reconstruct the traffic on both sides of the interface:

[0073]

[0074] In the formula, They are the single-width flows on the left and right sides in the x and y directions after reconstruction, in m 2 / s; are the flow velocities on the left and right sides of the center point M of the dry and wet unit boundary edge f in the x and y directions, respectively, in m / s;

[0075] Step 6.3: The bottom slope elevation value at the center point M of the dry-wet boundary cell edge f should be the lowest value of the bottom elevation and the left water level after reconstruction, that is:

[0076]

[0077] Step 6.4: Further correct u M 、v M , to ensure that the boundary edge of the dry and wet unit changes monotonically, the correction formula is as follows:

[0078]

[0079] In the formula, q xM ,q yM They are the single-width flow rates in the x and y directions after reconstruction of the center point M of the dry and wet unit boundary edge f, in m 3 / s;u M 、v M are the corresponding flow velocities in the x and y directions, in m / s; h M is the corresponding water depth after reconstruction, in meters; ε is the allowable tolerance, which is given by ε=T 2 It is obtained that T is the area of ​​the largest unit;

[0080] On complex grids, the water depth change value and water depth at the dry-wet boundary interface are used as the discrimination condition to solve the dry-wet boundary problem in the second-order or higher-order accuracy format that satisfies the complex grid terrain. The discrimination formula is as follows:

[0081]

[0082] Where: is the water depth on the left side of the center point M of the dry-wet unit boundary edge f, in meters; is the slope elevation to the left of the center point M of the dry-wet unit boundary edge f, in meters; h i is the water depth of the i-th unit grid, in meters; is the bottom slope elevation on the left side of the i-unit grid, in meters; γ is the ratio of the water depth reconstruction value at the edge of the grid unit to the water depth reconstruction value at the centroid of the unit, γ < 1, the higher the γ value, the more grid units meet this standard, ε wd is the water depth tolerance, usually ε wd =1×10 -6 .

[0083] Step 7 is implemented according to the following steps:

[0084] Step 7.1, use the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtain the restricted variable gradient by restricting the variable gradient, and interpolate the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction.

[0085] Step 7.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The i-th control unit grid in the new time step is The values ​​are:

[0086]

[0087] Intermediate variable value for

[0088]

[0089] Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency;

[0090] Step 7.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows:

[0091]

[0092] In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

[0093] The beneficial effect of the present invention is that the GPU-accelerated water environment simulation method considering inorganic substances is proposed to address the shortcomings of the existing water environment model, such as complex calculations and long simulation time due to the complex conversion relationship between various environmental factors. The method can accurately and efficiently simulate the transport laws of various environmental variables in the water bodies of rivers, lakes and seas and the conversion relationship between them, thereby achieving accurate simulation and effective evaluation of the eutrophication status of complex water bodies in rivers, lakes and seas. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 It is a schematic diagram of the research area in the GPU accelerated water environment simulation method considering inorganic substances of the present invention and its high-resolution DEM topographic map;

[0095] Figure 2 It is a schematic diagram of a GPU-accelerated water environment model considering inorganic substances in a GPU-accelerated water environment simulation method considering inorganic substances of the present invention;

[0096] FIG3( a ) is a comparison diagram of simulated and measured water depths at the lake outlet in the GPU-accelerated water environment simulation method considering inorganic substances of the present invention;

[0097] FIG3( b ) is a comparison diagram of simulated values ​​and measured values ​​of various environmental factor variables at monitoring points in the GPU-accelerated water environment simulation method considering inorganic substances of the present invention;

[0098] Figure 4 It is a quantitative evaluation result diagram of the performance of the GPU-accelerated water environment model considering inorganic substances in the GPU-accelerated water environment simulation method considering inorganic substances of the present invention;

[0099] FIG5(a) is a GPU accelerated water environment simulation method of the present invention considering inorganic substances with a throughput of 0.2 m 3 / s distribution map of hydrodynamic transfer.

[0100] FIG5( b ) is a GPU accelerated water environment simulation method of the present invention considering inorganic substances with a throughput of 0.2 m 3 / sTransport distribution diagram of various environmental factors. DETAILED DESCRIPTION

[0101] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0102] The present invention considers a GPU-accelerated water environment simulation method for inorganic substances, which is specifically implemented according to the following steps:

[0103] Step 1: Use ArcGIS to process the scattered terrain data to form a high-resolution DEM terrain representation of the real terrain of rivers, lakes and sea areas as the input condition of terrain data;

[0104] Step 2: discretize the above DEM terrain data through rectangular grids, and assign values ​​to relevant parameters according to known boundary conditions;

[0105] Step 3: Based on the biochemical reaction principle between substances and the convection diffusion theory, construct the hydrodynamic-material transport coupling equation;

[0106] The hydrodynamic-mass transport coupling equation in step 3 is expressed as follows:

[0107]

[0108] In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h can be expressed as a biochemical reaction term, with the unit of kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

[0109] Step 4: Couple the hydrodynamic model with the above hydrodynamic-material transport coupling equation to form a two-dimensional hydrodynamic-material transport coupling equation group, and represent it in vector form; use the central finite volume method of Godunov format for spatial discretization to solve the large deformation and large gradient water flow; at the same time, solve the flux term, diffusion term, bottom slope source term and friction source term of the coupling control equation group;

[0110] The vector form of the two-dimensional hydrodynamic-mass transport coupling equations in step 4 is as follows:

[0111]

[0112] in,

[0113]

[0114] Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, temperature diffusion term and temperature source term; z b is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average concentration of the substance along the water depth, in kg / m 3 ;D x , D yare the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s;

[0115] S in source and sink terms sc h is the biochemical reaction term, and its unit is kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

[0116] Step 4 is implemented according to the following steps:

[0117] Step 4.1: In the control unit grid i, the central finite volume method is used to integrate the two-dimensional hydrodynamic-material transport coupling control equations. Based on the Gaussian divergence theorem, the area integral of the flux term is converted into a line integral, which is:

[0118]

[0119] Where Ω is the volume of control volume i; Г is the boundary of control volume i; n is the unit vector in the direction of the external normal corresponding to the control volume boundary Г; F(Q)·n is the external normal convection term flux vector on the control volume boundary Г; S is the source term vector, including the bottom slope source term, friction source term, diffusion term and biochemical reaction term;

[0120] The flux F(Q)·n line integral is:

[0121]

[0122] Where Ω is the volume of the i-th control unit grid, in m 3 ; Г is the boundary of the i-th control unit grid; n is the unit vector in the direction of the external normal corresponding to the boundary Г; S b represents the bottom slope source flux; S f represents the friction source flux; k is the number of the control unit grid edge; l k is the length of the kth side of the i-th control unit grid;

[0123] The flux term is solved using the HLLC Riemann solver:

[0124] The wave speed is calculated as follows:

[0125]

[0126] Among them, the intermediate variable h * and According to the calculation formula,

[0127]

[0128] The calculation result of F(Q)·n is:

[0129]

[0130] Among them, the control volume unit grid interface intermediate flux F * for,

[0131]

[0132] Control volume unit grid interface flux intermediate wave sides, and for,

[0133]

[0134] In the formula, F L and F R are the fluxes on the left and right sides of the unit grid interface respectively; F * , and are the left wave flux, middle wave flux and right wave flux on both sides of the intermediate wave of the control volume unit grid interface flux respectively; u || is the tangential velocity along the interface, u || =-un y +vn x ,u ||L and u ||R are the left and right tangential velocities of the control volume unit grid interface, in m 3 / s;n x and n y are the indicated directions along x and y, respectively; u and v are the single-width flow rates along x and y, respectively, in m / s; variable Q ⊥ =[h,q x n x +q y n y ];q x and q y are the single width flow along x and y, in m 2 / s;S L , S M , S R They are left wave speed, middle wave speed, and right wave speed, in m / s; u ⊥ is the velocity perpendicular to the grid unit interface, in m / s; u ⊥ =un x +vn y ;h * and u * ⊥ is the intermediate variable; g is the acceleration due to gravity, which is 9.81 and the unit is m / s 2 ; Cs is the average concentration of the substance along the water depth, kg / m 3 ;

[0135] Step 4.2: Treat the diffusion term as a source term and use the second-order central difference method to discretize the diffusion term. The material transport diffusion term is expressed as:

[0136]

[0137] According to the second-order central difference method, the first-order derivative and second-order derivative of the material transport diffusion term are solved, and for any two-dimensional unit grid center (i, j) we get:

[0138]

[0139] Thus, the discrete form of the material transport diffusion term is obtained as follows:

[0140]

[0141] Step 4.3: The bottom slope source term is solved using the bottom slope flux method, and the vector form of the bottom slope flux at a certain side f is obtained as:

[0142]

[0143] In the formula, F sf (Q) is the bottom slope flux at a certain side f; n f is the unit external normal vector of the surface f; n fx and n fy Yes f and (the component of the unit external normal vector of face f), h L and z bL are the water depth and bottom elevation at the center of the grid cell, respectively. With z ML The water depth and bottom elevation at the edge of the grid. g is the acceleration due to gravity in m / s 2 ;

[0144] Step 4.4: The friction source term is solved using the explicit and implicit friction method. The calculation formula is:

[0145]

[0146] Among them, the calculation formula for single-width flow update to time step n+1 is:

[0147]

[0148] In the formula, S fx and S fy S f Friction source term in x and y directions; Cf is the bed surface roughness coefficient, C f =gn 2 / h 1 / 3 , n is the Manning coefficient, s / m 1 / 3 ; h is the water depth, in meters; q x ,q y are the average flow velocities in the x and y directions, in m 2 / s;q x n+1 ,q y n+1 They are the single-width flow rates in the x and y directions at time step n+1, in m 2 / s;q x n ,q y n They are the single-width flow rates in the x and y directions at time step n, in m 2 / s; Δt is the time step, in seconds; Δq x is the difference in single-width flow in the x direction at time step n and time step n+1, that is, Δq x =q x n+1 -q x n , unit is m 2 / s; In order to avoid the denominator being zero, θ is taken as 1.0e -12 .

[0149] Step 5: Substitute the biochemical reaction term S in the source term of the above hydrodynamic-material transport coupling equations sc h, and the dissection operator method is used to transform the equilibrium equation of each substance state variable to express it one by one. Among them, each substance state variable only considers the inorganic form of nutrients such as nitrogen and phosphorus in the water body. That is, the state variables of each substance are inorganic nitrogen (ammonia nitrogen NH 3 、Nitrate nitrogen NO 3 and nitrite nitrogen NO 2 ), inorganic phosphorus (phosphate PO 4 ), biochemical oxygen demand BOD, chemical oxygen demand COD, dissolved oxygen DO and aquatic plant chlorophyll α, etc.

[0150] Step 5 is implemented according to the following steps:

[0151] Step 5.1, dissolved oxygen balance includes dissolved oxygen DO, biochemical oxygen demand BOD and chemical oxygen demand COD; dissolved oxygen DO concentration includes aquatic plant respiration, nitrification, biodegradation and sediment or sediment oxygen consumption and photosynthetic oxygen production, biochemical oxygen demand BOD and chemical oxygen demand COD concentration are mainly related to their own degradation, and degradation is carried out with first-order reaction kinetics; at the same time, the Michaelis-Menten formula is used to describe the reaction process under low oxygen conditions;

[0152] Step 5.2: Nitrogen is mainly considered to exist in the form of inorganic nitrogen. Its cycle process includes ammonia nitrogen NH 3 、Nitrate nitrogen NO 3 and nitrite nitrogen NO 2 Three forms;

[0153] Step 5.3: Phosphorus in the phosphorus cycle mainly exists in the form of inorganic phosphorus, and only phosphate PO is considered. 4 Its existing form. Phosphate PO 4 The concentration includes the consumption by aquatic plants and bacterial microorganisms, the release from organic matter degradation reactions and anoxic biological respiration;

[0154] Step 5.4: Chlorophyll α concentration is an important indicator for evaluating the eutrophication of rivers and lakes. It takes the growth dynamics of phytoplankton as the main line and the ratio of the three main elements C:N:P to reflect the conversion relationship between phytoplankton, i.e. chlorophyll α, and nutrients in the environment.

[0155] Step 6: To further ensure the full stability condition and mass conservation of the above-mentioned coupling model, the improved hydrostatic reconstruction method is used to reconstruct the bottom slope elevation on the unit grid boundary surface by defining it as the maximum value at the boundary between two adjacent units; at the same time, the negative water depth value reconstructed at the interface is defined as zero to maintain non-negative water depth, and the corresponding water level and bottom slope elevation are modified to meet the balance between the flux term and the bottom slope source term in the hydrostatic simulation on complex grid terrain;

[0156] Step 6 is implemented according to the following steps:

[0157] Step 6.1: Based on the hydrostatic reconstruction method, the bottom slope elevation at the center point M of the dry-wet unit boundary edge f is taken as:

[0158]

[0159] In the formula, are the bottom slope elevations on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m;

[0160] Reconstruct the water depth on both sides, and at the same time, make the water level on the right side of the dry-wet interface the same as the elevation of the reconstructed bottom slope, that is:

[0161]

[0162] In the formula, They are the water levels on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; are the water depths on the left and right sides after reconstruction, in meters; z bM is the bottom slope elevation after reconstruction of the center point M of the dry-wet unit boundary edge f, in meters;

[0163] If expressed in terms of water level, the reconstructed water level is:

[0164]

[0165] In the formula, are the water levels after reconstruction, in m;

[0166] Step 6.2: Based on the above, reconstruct the traffic on both sides of the interface:

[0167]

[0168] In the formula, They are the single-width flows on the left and right sides in the x and y directions after reconstruction, in m 2 / s; are the flow velocities on the left and right sides of the center point M of the dry and wet unit boundary edge f in the x and y directions, respectively, in m / s;

[0169] Step 6.3: To ensure that the reconstructed bottom slope elevation is not higher than the water level on the left, the bottom slope elevation value at the center point M of the dry-wet boundary cell edge f should be the lowest value of the bottom elevation after reconstruction and the water level on the left, that is:

[0170]

[0171] Step 6.4: Further correct u M 、v M , to ensure that the boundary edge of the dry and wet unit changes monotonically, the correction formula is as follows:

[0172]

[0173] In the formula, q xM ,q yM They are the single-width flow rates in the x and y directions after reconstruction of the center point M of the dry and wet unit boundary edge f, in m 3 / s;u M 、v M are the corresponding flow velocities in the x and y directions, in m / s; h M is the corresponding water depth after reconstruction, in meters; ε is the allowable tolerance, which is given by ε=T 2 It is obtained that T is the area of ​​the largest unit;

[0174] On complex grids, the water depth change value and water depth at the dry-wet boundary interface are used as the discrimination condition to solve the dry-wet boundary problem in the second-order or higher-order accuracy format that satisfies the complex grid terrain. The discrimination formula is as follows:

[0175]

[0176] Where: is the water depth on the left side of the center point M of the dry-wet unit boundary edge f, in meters; is the slope elevation to the left of the center point M of the dry-wet unit boundary edge f, in meters; h i is the water depth of the i-th unit grid, in meters; is the bottom slope elevation on the left side of the i-unit grid, in meters; γ is the ratio of the water depth reconstruction value at the edge of the grid unit to the water depth reconstruction value at the centroid of the unit, γ < 1, the higher the γ value, the more grid units meet this standard, and it is taken as 0.25 in this paper; ε wd is the water depth tolerance, usually ε wd =1×10 -6 .

[0177] Step 7: Based on the above, the second-order precision of MUSCL and the second-order Runge-Kutta method are further used to achieve the second-order spatiotemporal precision of the model. The temporal stability is restricted by the CFL condition. GPU parallel computing is introduced in the whole model simulation process to improve the simulation efficiency of the model.

[0178] Step 7 is implemented according to the following steps:

[0179] Step 7.1, use the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtain the restricted variable gradient by restricting the variable gradient, and interpolate the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction.

[0180] Step 7.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The i-th control unit grid in the new time step is The values ​​are:

[0181]

[0182] Intermediate variable value for

[0183]

[0184] Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency;

[0185] Step 7.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows:

[0186]

[0187] In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

[0188] Step 8: Based on the simulation results of various environmental factors obtained from the above water environment simulation, use MATLAB program to output the distribution diagram of the transport process of various environmental factors to reflect the transport law of various environmental factors in the water body and evaluate the eutrophication status of the water body.

[0189] Example 1

[0190] The present invention considers a GPU-accelerated water environment simulation method for inorganic substances, which is specifically implemented according to the following steps:

[0191] Step 1: Use ArcGIS to process the scattered terrain data to form a high-resolution DEM terrain representation of the real terrain of rivers, lakes and sea areas as the input condition of terrain data;

[0192] Step 2: discretize the above DEM terrain data through rectangular grids, and assign values ​​to relevant parameters according to known boundary conditions;

[0193] Step 3: Based on the biochemical reaction principle between substances and the convection diffusion theory, construct the hydrodynamic-material transport coupling equation;

[0194] Step 4: Couple the hydrodynamic model with the above hydrodynamic-material transport coupling equation to form a two-dimensional hydrodynamic-material transport coupling equation group, and represent it in vector form; use the central finite volume method of Godunov format for spatial discretization to solve the large deformation and large gradient water flow; at the same time, solve the flux term, diffusion term, bottom slope source term and friction source term of the coupling control equation group;

[0195] Step 5: Substitute the biochemical reaction term S in the source term of the above hydrodynamic-material transport coupling equations sch, and the dissection operator method is used to transform the equilibrium equation of each substance state variable to express it one by one. Among them, each substance state variable only considers the inorganic form of nutrients such as nitrogen and phosphorus in the water body. That is, the state variables of each substance are inorganic nitrogen (ammonia nitrogen NH 3 、Nitrate nitrogen NO 3 and nitrite nitrogen NO 2 ), inorganic phosphorus (phosphate PO 4 ), biochemical oxygen demand BOD, chemical oxygen demand COD, dissolved oxygen DO and aquatic plant chlorophyll α, etc.

[0196] Step 6: To further ensure the full stability condition and mass conservation of the above-mentioned coupling model, the improved hydrostatic reconstruction method is used to reconstruct the bottom slope elevation on the unit grid boundary surface by defining it as the maximum value at the boundary between two adjacent units; at the same time, the negative water depth value reconstructed at the interface is defined as zero to maintain non-negative water depth, and the corresponding water level and bottom slope elevation are modified to meet the balance between the flux term and the bottom slope source term in the hydrostatic simulation on complex grid terrain;

[0197] Step 7: Based on the above, the second-order precision of MUSCL and the second-order Runge-Kutta method are further used to achieve the second-order spatiotemporal precision of the model. The temporal stability is restricted by the CFL condition. GPU parallel computing is introduced in the whole model simulation process to improve the simulation efficiency of the model.

[0198] Step 8: Based on the simulation results of various environmental factors obtained from the above water environment simulation, use MATLAB program to output the distribution diagram of the transport process of various environmental factors to reflect the transport law of various environmental factors in the water body and evaluate the eutrophication status of the water body.

[0199] Example 2

[0200] The present invention considers a GPU-accelerated water environment simulation method for inorganic substances, which is specifically implemented according to the following steps:

[0201] Step 1: Use ArcGIS to process the scattered terrain data to form a high-resolution DEM terrain representation of the real terrain of rivers, lakes and sea areas as the input condition of terrain data;

[0202] Step 2: discretize the above DEM terrain data through rectangular grids, and assign values ​​to relevant parameters according to known boundary conditions;

[0203] Step 3: Based on the biochemical reaction principle between substances and the convection diffusion theory, construct the hydrodynamic-material transport coupling equation;

[0204] The hydrodynamic-mass transport coupling equation in step 3 is expressed as follows:

[0205]

[0206] In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h can be expressed as a biochemical reaction term, with the unit of kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

[0207] Step 4: Couple the hydrodynamic model with the above hydrodynamic-material transport coupling equation to form a two-dimensional hydrodynamic-material transport coupling equation group, and represent it in vector form; use the central finite volume method of Godunov format for spatial discretization to solve the large deformation and large gradient water flow; at the same time, solve the flux term, diffusion term, bottom slope source term and friction source term of the coupling control equation group;

[0208] The vector form of the two-dimensional hydrodynamic-mass transport coupling equations in step 4 is as follows:

[0209]

[0210] in,

[0211]

[0212] Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, temperature diffusion term and temperature source term; z b is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average concentration of the substance along the water depth, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s;

[0213] S in source and sink terms sc h is the biochemical reaction term, and its unit is kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

[0214] Step 5: Substitute the biochemical reaction term S in the source term of the above hydrodynamic-material transport coupling equations sc h, and the dissection operator method is used to transform the equilibrium equation of each substance state variable to express it one by one. Among them, each substance state variable only considers the inorganic form of nutrients such as nitrogen and phosphorus in the water body. That is, the state variables of each substance are inorganic nitrogen (ammonia nitrogen NH 3 、Nitrate nitrogen NO 3 and nitrite nitrogen NO 2 ), inorganic phosphorus (phosphate PO 4 ), biochemical oxygen demand BOD, chemical oxygen demand COD, dissolved oxygen DO and aquatic plant chlorophyll α, etc.

[0215] Step 6: To further ensure the full stability condition and mass conservation of the above-mentioned coupling model, the improved hydrostatic reconstruction method is used to reconstruct the bottom slope elevation on the unit grid boundary surface by defining it as the maximum value at the boundary between two adjacent units; at the same time, the negative water depth value reconstructed at the interface is defined as zero to maintain non-negative water depth, and the corresponding water level and bottom slope elevation are modified to meet the balance between the flux term and the bottom slope source term in the hydrostatic simulation on complex grid terrain;

[0216] Step 7: Based on the above, the second-order precision of MUSCL and the second-order Runge-Kutta method are further used to achieve the second-order spatiotemporal precision of the model. The temporal stability is restricted by the CFL condition. GPU parallel computing is introduced in the whole model simulation process to improve the simulation efficiency of the model.

[0217] Step 7 is implemented according to the following steps:

[0218] Step 7.1, use the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtain the restricted variable gradient by restricting the variable gradient, and interpolate the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction.

[0219] Step 7.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The i-th control unit grid in the new time step is The values ​​are:

[0220]

[0221] Intermediate variable value for

[0222]

[0223]

[0224] Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency;

[0225] Step 7.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows:

[0226]

[0227] In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the gravitational acceleration, with a value of 9.81 meters per second 2 ; h is the water depth, in meters.

[0228] Step 8: Based on the simulation results of various environmental factors obtained from the above water environment simulation, use MATLAB program to output the distribution diagram of the transport process of various environmental factors to reflect the transport law of various environmental factors in the water body and evaluate the eutrophication status of the water body.

[0229] Example 3

[0230] This method selects Yanming Lake, located on the west bank of the Chan River in Xi'an, as an example for model application. The total area of ​​the lake is about 740,000 m 2 The water depth is about 1.5 to 2.5 meters, the maximum width of the lake is about 300 meters, the normal water level and the lowest lake level are 434.5 meters and 433.0 meters respectively, and the lake water mainly comes from the Chan River. Figure 1 The terrain of Yanming Lake No. 5 was discretized using rectangular grids with a grid size of 3.0 m, totaling 102,659 grids.

[0231] Based on the above numerical model of coupling hydrodynamics and material transport, the transformation balance equation of each material state variable is used as the basic theoretical basis, and only the inorganic form of nutrients such as nitrogen and phosphorus in the water body is considered, that is, inorganic nitrogen (ammonia nitrogen NH 3 、Nitrate nitrogen NO3 and nitrite nitrogen NO 2 ) and inorganic phosphorus (phosphate PO 4 ) form. At the same time, the material state variables such as biochemical oxygen demand BOD, chemical oxygen demand COD, dissolved oxygen DO and aquatic plant chlorophyll α are considered. The biochemical reaction term S in formula (3-1) is converted into iC h is expressed using each state variable equation. GPU parallel computing is introduced to improve simulation efficiency, and an efficient and high-resolution water environment model considering inorganic substances is constructed to simulate the nitrogen cycle, phosphorus cycle, dissolved oxygen balance, and the transport process of water environmental factors such as chlorophyll α. The relationship between the state variables in the water environment model is shown in Figure 2 .

[0232] In the simulation, the measured water depth data at the outlet of Yanming Lake No. 5 and the measured COD, TN, TP and chlorophyll α concentration data at the monitoring points were selected to verify the model and calibrate the parameters. The comparison chart of the verification results of the measured water depth data at the outlet of Yanming Lake No. 5 is shown in Figure 3(a). It can be seen from the figure that the change trend of the measured water depth is consistent with that of the simulated value, and they are basically distributed near the simulated value curve. The comparison chart of the verification results of the measured COD, TN, TP and chlorophyll α concentration data at the monitoring point of Yanming Lake No. 5 is shown in Figure 3(b). At the same time, in order to illustrate the simulation effect of the model, based on the measured and simulated values ​​of each state variable at the above-mentioned monitoring points, three statistical variables, MAE, RMSE and NSE, are used to quantitatively evaluate the simulation results, see. Figure 4 . According to the evaluation results, the error variables MAE, RMSE and NSE of the measured water depth at the outlet of No. 5 Yanming Lake are compared with the simulated value, which are 0.0221, 0.0273 and 0.9619 respectively. At the same time, the error results of the simulated and measured values ​​of water environmental factors such as COD, TN, TP and chlorophyll α at the monitoring point are further calculated, and the error statistical variables MAE, RMSE and NSE are 0.0020~0.8222, 0.0024~0.9203 and 0.7298~0.9671 respectively. Combined with the above comparison chart, it can be seen that the model has good simulation accuracy and high credibility, and can truly reflect the changes in the hydrodynamics and water environment in the water body of No. 5 Yanming Lake. On the basis of the calibration of these relevant parameters, a throughput flow of 0.2m 3 The hydrodynamic distribution of Yanming Lake under the condition of / s is shown in Figure 5(a); and the distribution of various environmental factors is shown in Figure 5(b), which reflects the distribution law of the lake with various environmental factor variables.

Claims

1. GPU-accelerated water environment simulation method considering inorganic substances, It is characterized in that Follow the steps below to implement it: Step 1: Use ArcGIS to process the scattered terrain data to form a high-resolution DEM terrain representation of the real terrain of rivers, lakes and sea areas as the input condition of terrain data; Step 2: discretize the above DEM terrain data through rectangular grids, and assign values ​​to relevant parameters according to known boundary conditions; Step 3: Construct the hydrodynamic-material transport coupling equation; Step 4, coupling the hydrodynamic model with the above hydrodynamic-material transport coupling equation; Step 5: Substitute the biochemical reaction term S in the source term of the hydrodynamic-material transport coupling equation sc h, using the split operator method to transform the equilibrium equation of each material state variable and express them one by one; Step 6: Use the improved hydrostatic reconstruction method to meet the balance between the flux term and the bottom slope source term in hydrostatic simulation on complex grid terrain; Step 7: Use the MUSCL second-order accuracy and the second-order Runge-Kutta method to achieve the time-space second-order accuracy of the model; Step 8: Output the distribution map of the transport process of each environmental factor to reflect the transport law of each environmental factor in the water body and evaluate the eutrophication status of the water body.

2. The GPU-accelerated water environment simulation method considering inorganic substances according to claim 1, It is characterized in that The hydrodynamic-material transport coupling equation in step 3 is expressed as follows: In the formula, C s is the average substance concentration along the water depth direction, in kg / m 3 ; D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink items sc h can be expressed as a biochemical reaction term, with the unit of kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

3. The GPU-accelerated water environment simulation method considering inorganic substances according to claim 2, It is characterized in that The vector form of the two-dimensional hydrodynamic-material transport coupling equations in step 4 is expressed as follows: in, Where t is time in seconds; Q is a variable vector, where h is water depth in meters, and q is x and q y They are the single width flow in the x and y directions, in m 2 / s; F, G are flux vectors in the x and y directions respectively; u, v are flow velocities in the x and y directions respectively, in units of m / s; S is the source term vector, including bottom slope source term, friction source term, temperature diffusion term and temperature source term; z b is the riverbed bottom elevation, in m; C f is the friction coefficient of the bed surface, calculated as C f =gn 2 / h 1 / 3 , where n is the Manning coefficient; g is the acceleration due to gravity, in m / s 2 ; C s is the average concentration of the substance along the water depth, in kg / m 3 ;D x , D y are the turbulent diffusion coefficients in the x and y directions, respectively, in m 2 / s; S in source and sink terms sc h is the biochemical reaction term, and its unit is kg / m 3 ; C 0 is the average substance concentration along the water depth direction, in kg / m 3 .

4. The GPU-accelerated water environment simulation method considering inorganic substances according to claim 3, It is characterized in that The step 4 is specifically implemented according to the following steps: Step 4.1: In the control unit grid i, the central finite volume method is used to integrate the two-dimensional hydrodynamic-material transport coupling control equations. Based on the Gaussian divergence theorem, the area integral of the flux term is converted into a line integral, which is: Where Ω is the volume of control volume i; Г is the boundary of control volume i; n is the unit vector in the direction of the external normal corresponding to the control volume boundary Г; F(Q)·n is the external normal convection term flux vector on the control volume boundary Г; S is the source term vector, including the bottom slope source term, friction source term, diffusion term and biochemical reaction term; The flux F(Q)·n line integral is: Where Ω is the volume of the i-th control unit grid, in m 3 ; Г is the boundary of the i-th control unit grid; n is the unit vector in the direction of the external normal corresponding to the boundary Г; S b represents the bottom slope source flux; S f represents the friction source flux; k is the number of the control unit grid edge; l k is the length of the kth side of the i-th control unit grid; The flux term is solved using the HLLC Riemann solver: The wave speed is calculated as follows: Among them, the intermediate variable h * and According to the calculation formula, The calculation result of F(Q)·n is: Among them, the control volume unit grid interface intermediate flux F * for, Control volume unit grid interface flux intermediate wave sides, and for, In the formula, F L and F R are the fluxes on the left and right sides of the unit grid interface respectively; F * , and are the left wave flux, middle wave flux and right wave flux on both sides of the intermediate wave of the control volume unit grid interface flux respectively; u || is the tangential velocity along the interface, u || =-un y +vn x ,u ||L and u ||R are the left and right tangential velocities of the control volume unit grid interface, in m 3 / s;n x and n y are the indicated directions along x and y, respectively; u and v are the single-width flow rates along x and y, respectively, in m / s; variable Q ⊥ =[h,q x n x +q y n y ];q x and q y are the single width flow along x and y, in m 2 / s;S L , S M , S R They are left wave speed, middle wave speed, and right wave speed, in m / s; u ⊥ is the velocity perpendicular to the grid unit interface, in m / s; u ⊥ =un x +vn y ;h * and is the intermediate variable; g is the acceleration due to gravity, which is 9.81 and the unit is m / s 2 ; C s is the average concentration of the substance along the water depth, kg / m 3 ; Step 4.2: Treat the diffusion term as a source term and use the second-order central difference method to discretize the diffusion term. The material transport diffusion term is expressed as: According to the second-order central difference method, the first-order derivative and second-order derivative of the material transport diffusion term are solved, and for any two-dimensional unit grid center (i, j) we get: Thus, the discrete form of the material transport diffusion term is obtained as follows: Step 4.3: The bottom slope source term is solved using the bottom slope flux method, and the vector form of the bottom slope flux at a certain side f is obtained as: In the formula, F sf (Q) is the bottom slope flux at a certain side f; n f is the unit external normal vector of the surface f; n fx and n fy Yes f and (the component of the unit external normal vector of face f), h L and z bL are the water depth and bottom elevation at the center of the grid cell, respectively. With z ML The water depth and bottom elevation at the edge of the grid. g is the acceleration due to gravity in m / s 2 ; Step 4.4: The friction source term is solved using the explicit and implicit friction method. The calculation formula is: Among them, the calculation formula for single-width flow update to time step n+1 is: In the formula, S fx and S fy S f Friction source term in x and y directions; C f is the bed surface roughness coefficient, C f =gn 2 / h 1 / 3 , n is the Manning coefficient, s / m 1 / 3 ; h is the water depth, in meters; q x ,q y are the average flow velocities in the x and y directions, in m 2 / s;q x n+1 ,q y n+1 They are the single-width flow rates in the x and y directions at time step n+1, in m 2 / s;q x n ,q y n They are the single-width flow rates in the x and y directions at time step n, in m 2 / s; Δt is the time step, in seconds; Δq x is the difference in single-width flow in the x direction at time step n and time step n+1, that is, Δq x =q x n+1 -q x n , unit is m 2 / s; In order to avoid the denominator being zero, θ is taken as 1.0e -12 .

5. The GPU-accelerated water environment simulation method considering inorganic substances according to claim 4, It is characterized in that The step 5 is specifically implemented according to the following steps: Step 5.1, dissolved oxygen balance includes dissolved oxygen DO, biochemical oxygen demand BOD and chemical oxygen demand COD; dissolved oxygen DO concentration includes aquatic plant respiration, nitrification, biodegradation and sediment or silt oxygen consumption and photosynthesis oxygen production; Step 5.2: The nitrogen cycle process includes ammonia nitrogen NH 3 、Nitrate nitrogen NO 3 and nitrite nitrogen NO 2 Three forms; Step 5.3: Phosphate PO in the Phosphorus Cycle 4 The concentration includes the consumption by aquatic plants and bacterial microorganisms, the release from organic matter degradation reactions and anoxic biological respiration; Step 5.4: Taking the growth dynamics of phytoplankton as the main line, the ratio of the three main elements C:N:P reflects the conversion relationship between phytoplankton, i.e. chlorophyll α, and nutrients in the environment.

6. The GPU-accelerated water environment simulation method considering inorganic substances according to claim 5, It is characterized in that The step 6 is specifically implemented according to the following steps: Step 6.1: Based on the hydrostatic reconstruction method, the bottom slope elevation at the center point M of the dry-wet unit boundary edge f is taken as: In the formula, are the bottom slope elevations on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; Reconstruct the water depth on both sides, and at the same time, make the water level on the right side of the dry-wet interface the same as the elevation of the reconstructed bottom slope, that is: In the formula, They are the water levels on the left and right sides of the center point M of the dry and wet unit boundary edge f, in m; are the water depths on the left and right sides after reconstruction, in meters; z bM is the bottom slope elevation after reconstruction of the center point M of the dry-wet unit boundary edge f, in meters; If expressed in terms of water level, the reconstructed water level is: In the formula, are the water levels after reconstruction, in m; Step 6.2: Based on the above, reconstruct the traffic on both sides of the interface: In the formula, They are the single-width flows on the left and right sides in the x and y directions after reconstruction, in m 2 / s; are the flow velocities on the left and right sides of the center point M of the dry and wet unit boundary edge f in the x and y directions, respectively, in m / s; Step 6.3: The bottom slope elevation value at the center point M of the dry-wet boundary cell edge f should be the lowest value of the bottom elevation and the left water level after reconstruction, that is: Step 6.4: Further correct u M 、v M , to ensure that the boundary edge of the dry and wet unit changes monotonically, the correction formula is as follows: In the formula, q xM ,q yM They are the single-width flow rates in the x and y directions after reconstruction of the center point M of the dry and wet unit boundary edge f, in m 3 / s;u M 、v M are the corresponding flow velocities in the x and y directions, in m / s; h M is the corresponding water depth after reconstruction, in meters; ε is the allowable tolerance, which is given by ε=T 2 It is obtained that T is the area of ​​the largest unit; On complex grids, the water depth change value and water depth at the dry-wet boundary interface are used as the discrimination condition to solve the dry-wet boundary problem in the second-order or higher-order accuracy format that satisfies the complex grid terrain. The discrimination formula is as follows: Where: is the water depth on the left side of the center point M of the dry-wet unit boundary edge f, in meters; is the slope elevation to the left of the center point M of the dry-wet unit boundary edge f, in meters; h i is the water depth of the i-th unit grid, in meters; is the bottom slope elevation on the left side of the i-unit grid, in meters; γ is the ratio of the water depth reconstruction value at the edge of the grid unit to the water depth reconstruction value at the centroid of the unit, γ < 1, the higher the γ value, the more grid units meet this standard, ε wd is the water depth tolerance, usually ε wd =1×10 -6 .

7. The GPU-accelerated water environment simulation method considering inorganic substances according to claim 6, It is characterized in that The step 7 is specifically implemented according to the following steps: Step 7.1, using the MUSCL reconstruction method to calculate the spatial variation gradient of each variable in the unit grid; obtaining the restricted variable gradient by restricting the variable gradient, and interpolating the variable value at the midpoint of the unit boundary based on the restricted gradient to achieve MUSCL second-order reconstruction; Step 7.2: Obtain the second-order time accuracy by the two-step Runge-Kutta method. The i-th control unit grid in the new time step is The values ​​are: Intermediate variable value for Where i is the number of control unit grids; Ω is the volume of the i-th control unit grid, in m 3 ;q n+1 is the single-width flow at time step n+1, in m 2 / s,q n is the single-width flow at time step n, in m 2 / s; Δt is the time step, in seconds; S(q n ) is the source term at time step n; during the entire solution process, GPU parallel computing technology is introduced to improve the computing efficiency; Step 7.3: Under the spatial step length Δx, the calculation formula for the value of the time step length Δt is as follows: In the formula, Δx min It represents the minimum distance from the center of the grid unit to the corresponding interface, in meters; u is the flow velocity, in meters per second; g is the acceleration due to gravity; and h is the water depth, in meters.

Images