Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation
Technical Field
The invention belongs to the technical field of numerical simulation of hydraulic engineering, and particularly relates to a constant flow simulation method based on an optimized solution plane two-dimensional shallow water equation.
Background
In engineering design of reservoir power stations, river channel treatment, channel improvement, water environment treatment and the like, a plane two-dimensional mathematical model is generally required to be applied to analyze influences of engineering construction on river channel water level, flow velocity, riverbed erosion and deposition, water quality and the like, and the two-dimensional numerical simulation is based on plane two-dimensional constant flow simulation. At present, a SIMPLE algorithm (a semi-implicit method for solving a pressure coupling equation) is most widely applied in a plane two-dimensional constant flow simulation method, but the SIMPLE algorithm is subjected to a plurality of assumptions and simplifications in the calculation process, so that the coupling relation between the flow rate and the water level is weakened, the robustness and the calculation efficiency of the SIMPLE algorithm are influenced, and the simulation efficiency is relatively low.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a constant flow simulation method based on an optimized solving plane two-dimensional shallow water equation, which can effectively solve the problems.
The technical scheme adopted by the invention is as follows:
the invention provides a constant flow simulation method based on an optimized solution plane two-dimensional shallow water equation, which comprises the following steps of:
step 1, obtaining the range of a research area;
step 2, mesh subdivision is carried out on the research area by adopting uniform rectangular meshes, and M multiplied by N meshes are obtained through common division; wherein, the number of the grids in the x direction is M, the number of the grids in the y direction is N, and the grid number is marked as (i, j); wherein i is more than or equal to 1 and less than or equal to M, j is more than or equal to 1 and less than or equal to N, the size of each grid is delta x delta y, delta x is the length of the grid in the x direction, delta y is the length of the grid in the y direction, and the grid center of each grid stores a variable to be solved; the variables to be solved comprise flow rate and water level;
step 3, setting boundary conditions, initial conditions and calculation convergence conditions:
the research area is provided with an inlet flow at an inlet boundary and an outlet boundary water level at an outlet boundary;
the initial conditions include: setting the initial water level of each grid as the outlet boundary water level, setting the initial water depth of each grid as the initial water level minus the river bottom elevation, and uniformly setting the initial flow rate of each grid as the section average flow rate of the section where the grid is located;
the calculation convergence condition is set as: the ratio of the continuous equation allowance two norm to the inlet flow is smaller than an allowable value;
step 4, dispersing an x-direction momentum equation and a y-direction momentum equation in a planar two-dimensional shallow water equation by adopting a finite volume method, and performing sub-relaxation treatment to obtain a dispersed momentum equation; solving the discrete momentum equation to obtain the temporary flow velocity of each grid;
step 5, constructing an optimized interface flow velocity interpolation formula, and calculating the interface flow velocity by adopting the optimized interface flow velocity interpolation formula;
step 6, calculating a two-norm of the continuous equation allowance by adopting a continuous equation in a finite volume method discrete plane two-dimensional shallow water equation, obtaining a water level correction equation by utilizing an optimized correction relation between the interface flow velocity and the water level, and solving the water level correction equation to obtain a water level correction value of each grid;
step 7, updating the temporary flow rate obtained in the step 4 and the interface flow rate obtained in the step 5 by adopting the water level correction value obtained in the step 6; updating the corresponding grid water level by adopting the water level correction value obtained in the step 6;
step 8, judging whether the ratio of the continuous equation allowance two norm to the inlet flow meets the calculation convergence condition set in the step 3, if so, iteratively converging, and outputting the flow velocity and water level distribution under the constant condition;
if not, the updated grid water level and the updated grid flow rate are used as initial values, the step 4 to the step 7 are repeated, and the next iteration is carried out until the calculation convergence condition is met.
Preferably, in step 3, the allowable value in the calculation convergence condition is 10-3~10-6。
Preferably, in step 4, the temporary flow rate of each grid is obtained by:
step 41, a constant plane two-dimensional shallow water equation, an x-direction momentum equation and a y-direction momentum equation in a cartesian coordinate system are respectively as follows:
in the formula:
u is the flow velocity in the x direction, v is the flow velocity in the y direction, h is the water depth, z is the water level, g is the gravity acceleration, r is the turbulent diffusion coefficient,
is the friction term in the x direction,
is a friction resistance term in the y direction, and n is a roughness coefficient;
step 42, discretizing an x-direction momentum equation (2) and a y-direction momentum equation (3) by using a finite volume method, discretizing a convection diffusion term by adopting a power format, discretizing a water level partial derivative term by adopting a central difference format, processing a friction term by adopting a semi-implicit format, and performing sub-relaxation processing to obtain a discrete momentum equation:
in the formula
Wherein:
p, E, W, N, S are grids with different relative positions, grid P is a grid with a central position, the adjacent grid in the north direction of grid P is grid N, the adjacent grid in the south direction of grid P is grid S, the adjacent grid in the west direction of grid P is grid W, and the adjacent grid in the east direction of grid P is grid E;
e. w, N and S are interfaces of grid P and grid E, grid P and grid W, grid P and grid N, grid P and grid S respectively, and the number of grid P is (i, j), so that the numbers of grid E, grid W, grid N and grid S are (i +1, j), (i-1, j), (i, j +1) and (i, j-1), the numbers of interface E, interface W, interface N and interface S are (i +1/2, j), (i-1/2, j), (i, j +1/2) and (i, j-1/2), respectively;
the superscript "0" represents the calculated value of the previous iteration,
the water depth, the flow velocity in the x direction and the flow velocity in the y direction at the grid P in the previous iteration,
respectively the water levels of the grid E, the grid W, the grid N and the grid S in the previous iteration;
the superscript "+" represents a temporary value,
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in the current iteration are respectively,
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
the temporary flow velocity in the x direction and the temporary flow velocity in the y direction of the grid N in the iteration are respectively,
respectively obtaining the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid S in the iteration of the current round;
Δ x and Δ y are the lengths of the grids in the x and y directions respectively; omega
uIs a velocity relaxation factor; d
e、D
w、D
n、D
sDiffusion flux at interface e, interface w, interface n, and interface s, respectively, F
e、F
w、F
n、F
sConvection fluxes at an interface e, an interface w, an interface n and an interface s are respectively;
is the main diagonal coefficient of the momentum equation,
respectively are the influence coefficients of grid E, grid W, grid N and grid S in the momentum equation on grid P;
and 43, respectively solving a discrete momentum equation (4) and a discrete momentum equation (5) by adopting an alternating direction iteration method to obtain the temporary flow velocity of each grid.
Preferably, step 5, the interfacial flow rate is obtained by:
step 51, the optimized flow rate correction formula is:
in the formula:
respectively correcting the flow velocity in the x direction and the flow velocity in the y direction of the grid P;
pseudo-x-direction flow velocity and pseudo-y-direction for grid P, respectivelyFlow rate, z 'of formula'
E、z'
W、z'
N、z'
SRespectively grid E, grid W, grid N and grid S water level correction values; b is
P、D
PRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the x direction; c
P、E
PRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the y direction;
step 52, according to the optimized flow rate correction formula, constructing an optimized interface flow rate interpolation formula as follows:
in the formula:
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
respectively simulating flow velocities at a grid E, a grid W, a grid N and a grid S; b is
E、B
W、C
N、C
SRespectively is a grid E, a grid W, a grid N and a grid S flow velocity coefficient;
and obtaining the interface flow velocity by adopting the optimized interface flow velocity interpolation formula.
Preferably, step 6 specifically comprises:
step 61, the discrete continuous equation of the finite volume method is as follows:
in the formula:
the interface flow velocity after the correction at the interface e, the interface w, the interface n and the interface s respectively has an optimized correction relationship with the water level as follows:
in the formula:
z'Pas grid P water level correction value, DE、DW、EN、ESRespectively taking flow velocity correction coefficients of a grid E, a grid W, a grid N and a grid S;
step 62, integrating the formulas (10) to (12), wherein the water level correction equation is as follows:
in the formula
In the formula:
is the main diagonal coefficient of the water level correction equation,
influence coefficients of a grid E, a grid W, a grid N and a grid S in the water level correction equation on a grid P are respectively;
step 62, the number of the grid P is recorded as (i, j), and the calculation formula of the margin two-norm a of the continuous equation is as follows:
and step 63, solving a water level correction equation (13) by adopting an alternating direction iteration method to obtain a water level correction value of each grid.
Preferably, step 7 specifically comprises:
step 71, updating the flow rate of each grid by adopting optimized flow rate correction formulas (6) to (7);
step 72, updating the interface flow rate by adopting the optimized interface flow rate correction formulas (11) to (12);
and step 73, introducing a water level relaxation factor, and correcting the water level as follows:
in the formula:
corrected water level, omega, for the grid P
pIs a water level relaxation factor.
The constant flow simulation method based on the optimized solution of the planar two-dimensional shallow water equation has the following advantages:
the constant flow simulation method based on the optimized solution of the planar two-dimensional shallow water equation has better robustness and simulation efficiency, and the simulation calculation amount of the method is less than 50% of that of the common SIMPLE algorithm when the speed relaxation factor is small.
Drawings
FIG. 1 is a schematic flow chart of a constant flow simulation method based on an optimized solution of a planar two-dimensional shallow water equation provided by the invention;
FIG. 2 is a schematic diagram of a grid arrangement provided by the present invention;
fig. 3 is a schematic diagram of a grid arrangement provided in embodiment 1 of the present invention;
FIG. 4 is a schematic view of a calculated flow field provided in example 1 of the present invention;
fig. 5 is a schematic diagram of a grid arrangement provided in embodiment 2 of the present invention;
fig. 6 is a schematic view of a calculated flow field provided in embodiment 2 of the present invention.
Detailed Description
In order to make the technical problems, technical solutions and advantageous effects solved by the present invention more clearly apparent, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Compared with the traditional SIMPLE simulation method, the method adopts the optimized flow rate interpolation and correction formula, namely the influence of adjacent point speed correction is considered when the flow rate is corrected, the momentum conservation is promoted on the premise of mass conservation, the coupling relation between the flow rate and the water level is improved, and the method has better robustness and calculation efficiency when the two-dimensional shallow water equation is solved, so that the simulation efficiency of the constant flow is improved.
Referring to the attached drawings, the constant flow simulation method based on the optimization solution of the plane two-dimensional shallow water equation comprises the following steps:
step 1, obtaining the range of a research area;
step 2, mesh subdivision is carried out on the research area by adopting uniform rectangular meshes, and M multiplied by N meshes are obtained through common division; wherein, the number of the grids in the x direction is M, the number of the grids in the y direction is N, and the grid number is marked as (i, j); wherein i is more than or equal to 1 and less than or equal to M, j is more than or equal to 1 and less than or equal to N, the size of each grid is delta x delta y, delta x is the length of the grid in the x direction, delta y is the length of the grid in the y direction, and the grid center of each grid stores a variable to be solved; the variables to be solved comprise flow rate and water level;
step 3, setting boundary conditions, initial conditions and calculation convergence conditions:
the research area is provided with an inlet flow at an inlet boundary and an outlet boundary water level at an outlet boundary;
the initial conditions include: setting the initial water level of each grid as the outlet boundary water level, setting the initial water depth of each grid as the initial water level minus the river bottom elevation, and uniformly setting the initial flow rate of each grid as the section average flow rate of the section where the grid is located;
the calculation convergence condition is set as: the ratio of the two-norm residual of the equation of continuity to the inlet flow is less than the allowable value, typically 10-3~10-6。
Step 4, dispersing an x-direction momentum equation and a y-direction momentum equation in a planar two-dimensional shallow water equation by adopting a finite volume method, and performing sub-relaxation treatment to obtain a dispersed momentum equation; solving the discrete momentum equation to obtain the temporary flow velocity of each grid;
in this step, the temporary flow rate of each grid is obtained specifically by:
step 41, a constant plane two-dimensional shallow water equation, an x-direction momentum equation and a y-direction momentum equation in a cartesian coordinate system are respectively as follows:
in the formula:
u is the flow velocity in the x direction, v is the flow velocity in the y direction, h is the water depth, z is the water level, g is the gravity acceleration, r is the turbulent diffusion coefficient,
is the friction term in the x direction,
is a friction resistance term in the y direction, and n is a roughness coefficient;
step 42, discretizing an x-direction momentum equation (2) and a y-direction momentum equation (3) by using a finite volume method, discretizing a convection diffusion term by adopting a power format, discretizing a water level partial derivative term by adopting a central difference format, processing a friction term by adopting a semi-implicit format, and performing sub-relaxation processing to obtain a discrete momentum equation:
in the formula
Wherein:
p, E, W, N, S are grids with different relative positions, grid P is a grid with a central position, the adjacent grid in the north direction of grid P is grid N, the adjacent grid in the south direction of grid P is grid S, the adjacent grid in the west direction of grid P is grid W, and the adjacent grid in the east direction of grid P is grid E;
e. w, N and S are interfaces of grid P and grid E, grid P and grid W, grid P and grid N, grid P and grid S respectively, and the number of grid P is (i, j), so that the numbers of grid E, grid W, grid N and grid S are (i +1, j), (i-1, j), (i, j +1) and (i, j-1), the numbers of interface E, interface W, interface N and interface S are (i +1/2, j), (i-1/2, j), (i, j +1/2) and (i, j-1/2), respectively;
the superscript "0" represents the calculated value of the previous iteration,
the water depth, the flow velocity in the x direction and the flow velocity in the y direction at the grid P in the previous iteration,
respectively the water levels of the grid E, the grid W, the grid N and the grid S in the previous iteration;
the superscript "+" represents a temporary value,
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in the current iteration are respectively,
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
the temporary flow velocity in the x direction and the temporary flow velocity in the y direction of the grid N in the iteration are respectively,
respectively obtaining the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid S in the iteration of the current round;
Δ x and Δ y are the lengths of the grids in the x and y directions respectively; omega
uIs a velocity relaxation factor; d
e、D
w、D
n、D
sDiffusion flux at interface e, interface w, interface n, and interface s, respectively, F
e、F
w、F
n、F
sConvection fluxes at an interface e, an interface w, an interface n and an interface s are respectively;
is the main diagonal coefficient of the momentum equation,
respectively are the influence coefficients of grid E, grid W, grid N and grid S in the momentum equation on grid P;
and 43, respectively solving the discrete momentum equation (4) and the discrete momentum equation (5) by adopting an alternating direction iterative method (ADI method) to obtain the temporary flow velocity of each grid.
Step 5, constructing an optimized interface flow velocity interpolation formula, and calculating the interface flow velocity by adopting the optimized interface flow velocity interpolation formula;
in this step, the interfacial flow rate is obtained specifically by:
step 51, under the condition of considering the adjacent point velocity correction influence and promoting momentum conservation, the optimized flow velocity correction formula is as follows:
in the formula:
respectively corrected x-direction flow for mesh PVelocity and y-direction flow rate;
respectively, a pseudo-flow velocity in the x-direction and a pseudo-flow velocity in the y-direction of grid P, where z'
E、z'
W、z'
N、z'
SRespectively grid E, grid W, grid N and grid S water level correction values; b is
P、D
PRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the x direction; c
P、E
PRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the y direction;
step 52, according to the optimized flow rate correction formula, constructing an optimized interface flow rate interpolation formula as follows:
in the formula:
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
respectively simulating flow velocities at a grid E, a grid W, a grid N and a grid S; b is
E、B
W、C
N、C
SRespectively is a grid E, a grid W, a grid N and a grid S flow velocity coefficient;
and obtaining the interface flow velocity by adopting the optimized interface flow velocity interpolation formula.
Step 6, calculating a two-norm residual quantity of the continuous equation by adopting a continuous equation in a finite volume method discrete plane two-dimensional shallow water equation, obtaining a water level correction equation by utilizing an optimized correction relation between the interface flow rate and the water level, and solving the water level correction equation to obtain a water level correction value of each grid;
the step 6 specifically comprises the following steps:
step 61, the discrete continuous equation of the finite volume method is as follows:
in the formula:
the interface flow velocity after the correction at the interface e, the interface w, the interface n and the interface s respectively has an optimized correction relationship with the water level as follows:
in the formula:
z'Pas grid P water level correction value, DE、DW、EN、ESRespectively taking flow velocity correction coefficients of a grid E, a grid W, a grid N and a grid S;
step 62, integrating the formulas (10) to (12), wherein the water level correction equation is as follows:
in the formula
In the formula:
is the main diagonal coefficient of the water level correction equation,
influence coefficients of a grid E, a grid W, a grid N and a grid S in the water level correction equation on a grid P are respectively;
step 62, the number of the grid P is recorded as (i, j), and the calculation formula of the margin two-norm a of the continuous equation is as follows:
and step 63, solving a water level correction equation (13) by adopting an alternating direction iteration method to obtain a water level correction value of each grid.
Step 7, updating the temporary flow rate obtained in the step 4 and the interface flow rate obtained in the step 5 by adopting the water level correction value obtained in the step 6; updating the water level of the corresponding grid by adopting the water level correction value obtained in the step 6;
the step 7 specifically comprises the following steps:
step 71, updating the flow rate of each grid by adopting optimized flow rate correction formulas (6) to (7);
step 72, updating the interface flow rate by adopting the optimized interface flow rate correction formulas (11) to (12);
and step 73, introducing a water level relaxation factor, and correcting the water level as follows:
in the formula:
corrected water level, omega, for the grid P
pIs a water level relaxation factor.
Step 8, judging whether the ratio of the continuous equation allowance two norm to the inlet flow meets the calculation convergence condition set in the step 3, if so, iteratively converging, and outputting the flow velocity and water level distribution under the constant condition;
if not, the updated grid water level and flow rate are used as initial values, the steps 4 to 7 are repeated, and the next iteration is carried out until the calculation convergence condition is met.
Therefore, through the optimized flow velocity interpolation and correction formula, the coupling relation between the flow velocity and the water level is enhanced, and the simulation efficiency and the robustness are improved.
Two examples are presented below:
the first embodiment is as follows:
fig. 3 and 4 are schematic diagrams of a single-sided relaxation flow algorithm. The inlet section is 8m long and 20m wide, the width of the widening section is 92m long and 30m wide, the elevation of the inlet river bottom is 0.05m, the elevation of the outlet river bottom is 0m, and the vertical gradient of the river reach is 0.0005. The grid is arranged as follows: 100 pieces in the x direction, 30 pieces in the y direction, and 1m in Δ x ═ Δ y. Turbulent diffusion coefficient f 0.02m2And s, roughness n is 0.03. The boundary conditions of the inlet and the outlet are as follows: inlet flow rate of 10m3And/s, the outlet water level is 1 m. Calculating initial conditions: the initial water level of each grid is the outlet water level, the flow velocity in the x direction is the average flow velocity of the cross section, and the flow velocity in the y direction is 0 m/s. Calculating convergence condition as that the ratio of two norms of the discrete continuous equation residual error to inlet flow is less than 2.0 multiplied by 10-6。
Table 1 compares the computational efficiency and robustness of the method of the present invention and the SIMPLE algorithm in simulating a single-sided relaxation flow example. Relaxation factor omega at the same speeduNext, the water level relaxation factor omega converged by the method of the inventionpThe range is wider than that of the SIMPLE algorithm, and the robustness of the method is better than that of the SIMPLE algorithm; under the condition of the same speed and water level relaxation factors, the number of iteration steps required by convergence of the method is generally smaller than that of a SIMPLE algorithm, the calculation efficiency is higher, and particularly when the speed relaxation factors are smaller, the number of the convergence steps calculated by the method is less than half of that of the SIMPLE algorithm.
TABLE 1 comparison of computational efficiency and robustness (one-sided relaxation example, grid number 100X 30)
Note "-" represents a calculation of no convergence
Example two:
fig. 5 and 6 are schematic diagrams of two-sided symmetrical relaxation flow examples. The inlet section is 10m long and 28m wide, the width of the relaxation section is 70m long and 40m wide, and the river bottom elevation is 0 m. The grid is arranged as follows: the number of x-direction is 40, the number of y-direction is 40, Δ x is 2m, and Δ y is 1 m. Turbulent diffusion coefficient f 0.02m2And s, roughness n is 0.025. The boundary conditions of the inlet and the outlet are as follows: the inlet flow rate is 14m3And/s, the outlet water level is 1 m. Calculating initial conditions: the initial water level of each grid is the outlet water level, the flow velocity in the x direction is the average flow velocity of the cross section, and the flow velocity in the y direction is 0 m/s. Calculating convergence condition as that the ratio of two norms of the discrete continuous equation residual error to inlet flow is less than 1.0 multiplied by 10-6。
Table 2 compares the computational efficiency and robustness of the method of the present invention and the SIMPLE algorithm in simulating a bilateral symmetry relaxation flow example. Relaxation factor omega at the same speeduNext, the water level relaxation factor omega converged by the method of the inventionpThe range is wider than that of a SIMPLE algorithm, and the robustness is relatively better; under the condition of the same speed and water level relaxation factor combination, the calculated convergence step number of the method is generally smaller than that of a SIMPLE algorithm, and especially when the speed relaxation factor is smaller, the calculation efficiency of the method is obviously superior to that of the SIMPLE algorithm.
TABLE 2 comparison of computational efficiency and robustness (bilateral relaxation example, grid number 40X 40)
Note "-" represents a calculation of no convergence
Therefore, the invention provides a constant flow simulation method based on an optimized solution plane two-dimensional shallow water equation, which comprises the steps of firstly carrying out mesh subdivision on a research area, and setting inlet and outlet boundary conditions, initial conditions and calculation convergence conditions; then, a finite volume method is adopted to disperse a momentum equation and solve the equation; calculating the interface flow velocity by adopting an optimized interface flow velocity interpolation formula; then, a finite volume method discrete continuous equation is adopted, a water level correction equation is deduced and solved; correcting the flow rate by adopting an optimized flow rate correction formula, and synchronously correcting the water level; and repeating iteration until the convergence condition is met, thereby obtaining the flow velocity, water level and water depth distribution of each grid under the constant condition. Compared with the traditional SIMPLE algorithm (a semi-implicit method for solving a pressure coupling equation), the simulation method provided by the invention is better in robustness and simulation efficiency, and when the speed relaxation factor is small, the simulation calculation amount of the method provided by the invention is less than 50% of that of the common SIMPLE algorithm.
The foregoing is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, various modifications and improvements can be made without departing from the principle of the present invention, and such modifications and improvements should also be considered within the scope of the present invention.