CN111931263B - Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation - Google Patents

Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation Download PDF

Info

Publication number
CN111931263B
CN111931263B CN202010424842.5A CN202010424842A CN111931263B CN 111931263 B CN111931263 B CN 111931263B CN 202010424842 A CN202010424842 A CN 202010424842A CN 111931263 B CN111931263 B CN 111931263B
Authority
CN
China
Prior art keywords
grid
flow velocity
interface
water level
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010424842.5A
Other languages
Chinese (zh)
Other versions
CN111931263A (en
Inventor
汪飞
徐高洪
邴建平
徐长江
郭海晋
贾建伟
邓鹏鑫
李林娟
王栋
何康洁
郭熙望
邹宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Bureau of Hydrology Changjiang Water Resources Commission
Original Assignee
Bureau of Hydrology Changjiang Water Resources Commission
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Bureau of Hydrology Changjiang Water Resources Commission filed Critical Bureau of Hydrology Changjiang Water Resources Commission
Priority to CN202010424842.5A priority Critical patent/CN111931263B/en
Publication of CN111931263A publication Critical patent/CN111931263A/en
Application granted granted Critical
Publication of CN111931263B publication Critical patent/CN111931263B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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: mesh generation is carried out on the research area; setting boundary conditions, initial conditions and calculation convergence conditions; solving a discrete momentum equation to obtain the temporary flow velocity of each grid; calculating the interface flow velocity by adopting the optimized interface flow velocity interpolation formula; solving a water level correction equation to obtain a water level correction value of each grid; updating the flow rate, the interface flow rate and the water level of the grid; and if the iterative calculation meets the convergence condition, outputting the flow velocity and water level distribution under the constant condition of the research area. The invention has the following advantages: the method is based on the optimization solution of the plane two-dimensional shallow water equation, the robustness and the simulation efficiency are better, 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.

Description

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:
Figure GDA0003197956740000031
Figure GDA0003197956740000032
Figure GDA0003197956740000033
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,
Figure GDA0003197956740000034
is the friction term in the x direction,
Figure GDA0003197956740000035
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:
Figure GDA0003197956740000036
Figure GDA0003197956740000037
in the formula
Figure GDA0003197956740000038
Figure GDA0003197956740000039
Figure GDA00031979567400000310
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,
Figure GDA0003197956740000041
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,
Figure GDA0003197956740000042
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,
Figure GDA0003197956740000043
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
Figure GDA0003197956740000044
respectively a local iteration middle netThe x-direction temporary flow rate and the y-direction temporary flow rate of the cell E,
Figure GDA0003197956740000045
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
Figure GDA0003197956740000046
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,
Figure GDA0003197956740000047
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; omegauIs a velocity relaxation factor; de、Dw、Dn、DsDiffusion flux at interface e, interface w, interface n, and interface s, respectively, Fe、Fw、Fn、FsConvection fluxes at an interface e, an interface w, an interface n and an interface s are respectively;
Figure GDA0003197956740000048
is the main diagonal coefficient of the momentum equation,
Figure GDA0003197956740000049
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:
Figure GDA0003197956740000051
Figure GDA0003197956740000052
in the formula:
Figure GDA0003197956740000053
Figure GDA0003197956740000054
Figure GDA0003197956740000055
respectively correcting the flow velocity in the x direction and the flow velocity in the y direction of the grid P;
Figure GDA0003197956740000056
mesh P pseudo-x-direction flow velocity and pseudo-y-direction flow velocity, respectively, where z'E、z'W、z'N、z'SRespectively grid E, grid W, grid N and grid S water level correction values; b isP、DPRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the x direction; cP、EPRespectively 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:
Figure GDA0003197956740000057
Figure GDA0003197956740000058
in the formula:
Figure GDA0003197956740000059
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
Figure GDA00031979567400000510
respectively simulating flow velocities at a grid E, a grid W, a grid N and a grid S; b isE、BW、CN、CSRespectively 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:
Figure GDA0003197956740000061
in the formula:
Figure GDA0003197956740000062
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:
Figure GDA0003197956740000063
Figure GDA0003197956740000064
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:
Figure GDA0003197956740000065
in the formula
Figure GDA0003197956740000066
Figure GDA0003197956740000067
In the formula:
Figure GDA0003197956740000068
is the main diagonal coefficient of the water level correction equation,
Figure GDA0003197956740000069
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:
Figure GDA00031979567400000610
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:
Figure GDA0003197956740000071
in the formula:
Figure GDA0003197956740000072
corrected water level, omega, for the grid PpIs 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:
Figure GDA0003197956740000091
Figure GDA0003197956740000092
Figure GDA0003197956740000093
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,
Figure GDA0003197956740000094
is the friction term in the x direction,
Figure GDA0003197956740000095
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:
Figure GDA0003197956740000096
Figure GDA0003197956740000097
in the formula
Figure GDA0003197956740000098
Figure GDA0003197956740000099
Figure GDA0003197956740000101
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,
Figure GDA0003197956740000102
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,
Figure GDA0003197956740000103
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,
Figure GDA0003197956740000104
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
Figure GDA0003197956740000105
are respectively asThe x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in this iteration,
Figure GDA0003197956740000106
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
Figure GDA0003197956740000107
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,
Figure GDA0003197956740000108
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; omegauIs a velocity relaxation factor; de、Dw、Dn、DsDiffusion flux at interface e, interface w, interface n, and interface s, respectively, Fe、Fw、Fn、FsConvection fluxes at an interface e, an interface w, an interface n and an interface s are respectively;
Figure GDA0003197956740000109
is the main diagonal coefficient of the momentum equation,
Figure GDA00031979567400001010
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:
Figure GDA0003197956740000111
Figure GDA0003197956740000112
in the formula:
Figure GDA0003197956740000113
Figure GDA0003197956740000114
Figure GDA0003197956740000115
respectively correcting the flow velocity in the x direction and the flow velocity in the y direction of the grid P;
Figure GDA0003197956740000116
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 isP、DPRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the x direction; cP、EPRespectively 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:
Figure GDA0003197956740000117
Figure GDA0003197956740000118
in the formula:
Figure GDA0003197956740000121
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
Figure GDA0003197956740000122
respectively simulating flow velocities at a grid E, a grid W, a grid N and a grid S; b isE、BW、CN、CSRespectively 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:
Figure GDA0003197956740000123
in the formula:
Figure GDA0003197956740000124
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:
Figure GDA0003197956740000125
Figure GDA0003197956740000126
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:
Figure GDA0003197956740000127
in the formula
Figure GDA0003197956740000131
Figure GDA0003197956740000132
In the formula:
Figure GDA0003197956740000133
is the main diagonal coefficient of the water level correction equation,
Figure GDA0003197956740000134
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:
Figure GDA0003197956740000135
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:
Figure GDA0003197956740000136
in the formula:
Figure GDA0003197956740000137
corrected water level, omega, for the grid PpIs 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 flowDynamic 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)
Figure GDA0003197956740000141
Figure GDA0003197956740000151
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 conditionsThe ratio of two norm of the discrete continuous equation residual to the 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)
Figure GDA0003197956740000161
Figure GDA0003197956740000171
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.

Claims (3)

1. A constant flow simulation method based on an optimized solving plane two-dimensional shallow water equation is characterized by comprising 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 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;
in step 4, the temporary flow rate of each grid is obtained by the following method:
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:
Figure FDA0003197956730000011
Figure FDA0003197956730000021
Figure FDA0003197956730000022
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,
Figure FDA0003197956730000023
is the friction term in the x direction,
Figure FDA0003197956730000024
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:
Figure FDA0003197956730000025
Figure FDA0003197956730000026
in the formula
Figure FDA0003197956730000027
Figure FDA0003197956730000028
Figure FDA0003197956730000029
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,
Figure FDA0003197956730000031
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,
Figure FDA0003197956730000032
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,
Figure FDA0003197956730000033
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
Figure FDA0003197956730000034
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in the current iteration are respectively,
Figure FDA0003197956730000035
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
Figure FDA0003197956730000036
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,
Figure FDA0003197956730000037
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; omegauIs a velocity relaxation factor; de、Dw、Dn、DsDiffusion flux at interface e, interface w, interface n, and interface s, respectively, Fe、Fw、Fn、FsConvection fluxes at an interface e, an interface w, an interface n and an interface s are respectively;
Figure FDA0003197956730000038
is the main diagonal coefficient of the momentum equation,
Figure FDA0003197956730000039
respectively are the influence coefficients of grid E, grid W, grid N and grid S in the momentum equation on grid P;
step 43, 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;
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;
and 5, specifically obtaining the interfacial flow rate in the following mode:
step 51, the optimized flow rate correction formula is:
Figure FDA00031979567300000310
Figure FDA0003197956730000041
in the formula:
Figure FDA0003197956730000042
Figure FDA0003197956730000043
Figure FDA0003197956730000044
respectively correcting the flow velocity in the x direction and the flow velocity in the y direction of the grid P;
Figure FDA0003197956730000045
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 isP、DPRespectively is a flow velocity coefficient and a flow velocity correction coefficient of the grid P in the x direction; cP、EPRespectively 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:
Figure FDA0003197956730000046
Figure FDA0003197956730000047
in the formula:
Figure FDA0003197956730000048
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
Figure FDA0003197956730000049
respectively simulating flow velocities at a grid E, a grid W, a grid N and a grid S; b isE、BW、CN、CSRespectively is a grid E, a grid W, a grid N and a grid S flow velocity coefficient;
obtaining 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;
the step 6 specifically comprises the following steps:
step 61, the discrete continuous equation of the finite volume method is as follows:
Figure FDA0003197956730000051
in the formula:
Figure FDA0003197956730000052
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:
Figure FDA0003197956730000053
Figure FDA0003197956730000054
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:
Figure FDA0003197956730000055
in the formula
Figure FDA0003197956730000056
Figure FDA0003197956730000057
In the formula:
Figure FDA0003197956730000058
is the main diagonal coefficient of the water level correction equation,
Figure FDA0003197956730000059
grid E and grid W in water level correction equation respectivelyInfluence coefficients of the grid N and the grid S on the grid P;
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:
Figure FDA0003197956730000061
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 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 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.
2. The method for constant flow simulation based on optimized solution of planar two-dimensional shallow water equation as claimed in claim 1, wherein in step 3, the allowable value in the calculation convergence condition is 10-3～10-6。
3. The constant flow simulation method based on the optimal solution of the planar two-dimensional shallow water equation according to claim 1, wherein the step 7 is specifically as follows:
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:
Figure FDA0003197956730000062
in the formula:
Figure FDA0003197956730000063
corrected water level, omega, for the grid PpIs a water level relaxation factor.
CN202010424842.5A 2020-05-19 2020-05-19 Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation Active CN111931263B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010424842.5A CN111931263B (en) 2020-05-19 2020-05-19 Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010424842.5A CN111931263B (en) 2020-05-19 2020-05-19 Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation

Publications (2)

Publication Number Publication Date
CN111931263A CN111931263A (en) 2020-11-13
CN111931263B true CN111931263B (en) 2021-09-14

Family

ID=73316362

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010424842.5A Active CN111931263B (en) 2020-05-19 2020-05-19 Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation

Country Status (1)

Country Link
CN (1) CN111931263B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113094878A (en) * 2021-03-23 2021-07-09 华中科技大学 Parallel computing method and system for casting flow field numerical simulation

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106503396B (en) * 2016-11-14 2019-04-12 中国电建集团昆明勘测设计研究院有限公司 Multidimensional hydraulic system transient simulation method based on finite difference method and finite volume method coupling
CN109543275B (en) * 2018-11-15 2019-07-09 中国水利水电科学研究院 A kind of city rainwash Two-dimensional numerical simulation method
CN110765694B (en) * 2019-11-21 2021-11-05 华南理工大学 Urban surface water flow numerical simulation method based on simplified shallow water equation set

Also Published As

Publication number Publication date
CN111931263A (en) 2020-11-13

Similar Documents

Publication Publication Date Title
CN109086534B (en) Wind farm wake correction method and system based on CFD hydrodynamic model
CN102722611B (en) Method for carrying out parallelization numerical simulation on hydrodynamic force conditions of river provided with cascade hydropower station
CN111931263B (en) Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation
CN113158451A (en) Large-area river three-dimensional simulation method based on one-dimensional flood routing model
CN102890732B (en) Hydrodynamic condition parallel numerical simulation method for channel with various structures
CN115329689A (en) High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion
CN114117877B (en) Topological optimization method based on isogeometric particle description
CN109408927B (en) Two-dimensional static magnetic field parallel finite element acceleration method based on black box transmission line model
CN108694299B (en) ICEM-CFD-based two-dimensional finite element neutronics steady-state calculation method
CN114048558B (en) Modeling method for blade profile of air compressor with non-uniform contour error
CN114741798B (en) Motor rotor structure topology optimization method considering electromagnetic and mechanical properties
CN115618635A (en) Method and equipment for calculating working parameters of water cooling system of energy storage container
JP2011065360A (en) Flow numerical analysis method for setting boundary condition of unramified and non-orthogonal structural grid and using iterative computation
CN114818531A (en) Motion interface tracking numerical dissipation calculation method
CN111539153B (en) Water-sand joint optimization scheduling method based on preconfigured sediment information base
CN114417675A (en) Finite element calculation method for continuous casting, solidification and heat transfer of special-shaped blank
CN112560247B (en) Automatic generation method of similar-appearance structural grid based on reference grid
CN109224796B (en) Flue flow equalizer device and design method
CN108682046B (en) Suspended nozzle calculation method with flexible base material
CN108595726B (en) Riverbed adjusting method based on underwater repose angle of sediment
CN118212377B (en) CAD curved surface high-order triangular mesh generation method based on extremely small curved surface theory
Avvari et al. Flow apportionment algorithm for optimization of power plant ducting
CN115788908B (en) Bidirectional axial flow pump blade space coordinate design and construction method thereof
Papadakis et al. A local grid refinement method for three‐dimensional turbulent recirculating flows
CN113128039B (en) Method, device and storage medium for calculating vertical joint type fishway step type water surface line

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Wang Fei

Inventor after: He Kangjie

Inventor after: Guo Xiwang

Inventor after: Zou Ning

Inventor after: Bing Jianping

Inventor after: Xu Gaohong

Inventor after: Jia Jianwei

Inventor after: Xu Changjiang

Inventor after: Guo Haijin

Inventor after: Deng Pengxin

Inventor after: Li Linjuan

Inventor after: Wang Dong

Inventor before: Wang Fei

Inventor before: He Kangjie

Inventor before: Guo Xiwang

Inventor before: Zou Ning

Inventor before: Xu Gaohong

Inventor before: Bing Jianping

Inventor before: Xu Changjiang

Inventor before: Guo Haijin

Inventor before: Jia Jianwei

Inventor before: Deng Pengxin

Inventor before: Li Linjuan

Inventor before: Wang Dong