CN111931263A - 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
CN111931263A
CN111931263A CN202010424842.5A CN202010424842A CN111931263A CN 111931263 A CN111931263 A CN 111931263A CN 202010424842 A CN202010424842 A CN 202010424842A CN 111931263 A CN111931263 A CN 111931263A
Authority
CN
China
Prior art keywords
grid
flow velocity
interface
equation
water level
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010424842.5A
Other languages
Chinese (zh)
Other versions
CN111931263B (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)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Civil Engineering (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Architecture (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 BDA0002498266510000031
Figure BDA0002498266510000032
Figure BDA0002498266510000033
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 BDA0002498266510000034
is the friction term in the x direction,
Figure BDA0002498266510000035
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 BDA0002498266510000036
Figure BDA0002498266510000037
in the formula
Figure BDA0002498266510000038
Figure BDA0002498266510000039
Figure BDA00024982665100000310
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 BDA0002498266510000041
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 BDA0002498266510000042
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 BDA0002498266510000043
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
Figure BDA0002498266510000044
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in the current iteration are respectively,
Figure BDA0002498266510000045
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
Figure BDA0002498266510000046
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 BDA0002498266510000047
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 BDA0002498266510000048
is the main diagonal coefficient of the momentum equation,
Figure BDA0002498266510000049
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 BDA0002498266510000051
Figure BDA0002498266510000052
in the formula:
Figure BDA0002498266510000053
Figure BDA0002498266510000054
Figure BDA0002498266510000055
respectively correcting the flow velocity in the x direction and the flow velocity in the y direction of the grid P;
Figure BDA0002498266510000056
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 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 BDA0002498266510000057
Figure BDA0002498266510000058
in the formula:
Figure BDA0002498266510000059
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
Figure BDA00024982665100000510
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 BDA0002498266510000061
in the formula:
Figure BDA0002498266510000062
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 BDA0002498266510000063
Figure BDA0002498266510000064
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 BDA0002498266510000065
in the formula
Figure BDA0002498266510000066
Figure BDA0002498266510000067
In the formula:
Figure BDA0002498266510000068
is the main diagonal coefficient of the water level correction equation,
Figure BDA0002498266510000069
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 BDA00024982665100000610
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 BDA0002498266510000071
in the formula:
Figure BDA0002498266510000072
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 BDA0002498266510000091
Figure BDA0002498266510000092
Figure BDA0002498266510000093
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 BDA0002498266510000094
is the friction term in the x direction,
Figure BDA0002498266510000095
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 BDA0002498266510000096
Figure BDA0002498266510000097
in the formula
Figure BDA0002498266510000098
Figure BDA0002498266510000099
Figure BDA0002498266510000101
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 BDA0002498266510000102
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 BDA0002498266510000103
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 BDA0002498266510000104
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
Figure BDA0002498266510000105
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in the current iteration are respectively,
Figure BDA0002498266510000106
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
Figure BDA0002498266510000107
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 BDA0002498266510000108
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 BDA0002498266510000109
is the main diagonal coefficient of the momentum equation,
Figure BDA00024982665100001010
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 BDA0002498266510000111
Figure BDA0002498266510000112
in the formula:
Figure BDA0002498266510000113
Figure BDA0002498266510000114
Figure BDA0002498266510000115
respectively corrected x-direction flow for mesh PVelocity and y-direction flow rate;
Figure BDA0002498266510000116
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 BDA0002498266510000117
Figure BDA0002498266510000118
in the formula:
Figure BDA0002498266510000121
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
Figure BDA0002498266510000122
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 BDA0002498266510000123
in the formula:
Figure BDA0002498266510000124
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 BDA0002498266510000125
Figure BDA0002498266510000126
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 BDA0002498266510000127
in the formula
Figure BDA0002498266510000131
Figure BDA0002498266510000132
In the formula:
Figure BDA0002498266510000133
is the main diagonal coefficient of the water level correction equation,
Figure BDA0002498266510000134
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 BDA0002498266510000135
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 BDA0002498266510000136
in the formula:
Figure BDA0002498266510000137
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 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 BDA0002498266510000141
Figure BDA0002498266510000151
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)
Figure BDA0002498266510000161
Figure BDA0002498266510000171
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 (6)

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;
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 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. According toThe steady flow simulation method based on the optimal solution of the planar two-dimensional shallow water equation of claim 1, wherein in step 3, the allowable value in the calculation convergence condition is 10-3～10-6。
3. The method for simulating constant flow based on the optimized solution of the planar two-dimensional shallow water equation as claimed in claim 1, wherein 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 FDA0002498266500000021
Figure FDA0002498266500000022
Figure FDA0002498266500000023
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 FDA0002498266500000024
is the friction term in the x direction,
Figure FDA0002498266500000025
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 FDA0002498266500000031
Figure FDA0002498266500000032
in the formula
Figure FDA0002498266500000033
Figure FDA0002498266500000034
Figure FDA0002498266500000035
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 FDA0002498266500000036
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 FDA0002498266500000037
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 FDA0002498266500000038
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid P in the current iteration are respectively,
Figure FDA0002498266500000039
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid E in the current iteration are respectively,
Figure FDA00024982665000000310
the x-direction temporary flow velocity and the y-direction temporary flow velocity of the grid W in the current iteration are respectively,
Figure FDA00024982665000000311
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 FDA0002498266500000041
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 FDA0002498266500000042
is the main diagonal coefficient of the momentum equation,
Figure FDA0002498266500000043
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.
4. The constant flow simulation method based on the optimized solution of the planar two-dimensional shallow water equation according to claim 1, wherein in step 5, the interfacial flow rate is obtained by:
step 51, the optimized flow rate correction formula is:
Figure FDA0002498266500000044
Figure FDA0002498266500000045
in the formula:
Figure FDA0002498266500000046
Figure FDA0002498266500000047
Figure FDA0002498266500000048
respectively correcting the flow velocity in the x direction and the flow velocity in the y direction of the grid P;
Figure FDA0002498266500000049
are respectively asGrid P pseudo-flow velocity in x-direction and y-direction, 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 FDA0002498266500000051
Figure FDA0002498266500000059
in the formula:
Figure FDA0002498266500000052
interface flow rates at an interface e, an interface w, an interface n and an interface s are respectively;
Figure FDA0002498266500000053
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.
5. 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 6 is specifically as follows:
step 61, the discrete continuous equation of the finite volume method is as follows:
Figure FDA0002498266500000054
in the formula:
Figure FDA0002498266500000055
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 FDA0002498266500000056
Figure FDA0002498266500000057
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 FDA0002498266500000058
Figure FDA0002498266500000061
in the formula
Figure FDA0002498266500000062
Figure FDA0002498266500000063
In the formula:
Figure FDA0002498266500000064
is the main diagonal coefficient of the water level correction equation,
Figure FDA0002498266500000065
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 FDA0002498266500000066
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.
6. 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 FDA0002498266500000067
in the formula:
Figure FDA0002498266500000068
corrected water level, omega, for the grid PpIs the water levelA 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 true CN111931263A (en) 2020-11-13
CN111931263B 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)

Citations (3)

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

Patent Citations (3)

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

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
FEI WANG 等: "Fourier Analysis of the Effect of Momentum Interpolation on the SIMPLE Series Formulated on a Collocated Grid", 《NUMERICAL HEAT TRANSFER, PART B: FUNDAMENTALS》 *
汪飞: "SIMPLE系列算法的Fourier分析及其改进的初步研究", 《中国博士学位论文全文数据库基础科学辑》 *
董炳江 等: "基于非结构混合网格求解二维浅水方程的一种数值方法", 《水科学进展》 *

Also Published As

Publication number Publication date
CN111931263B (en) 2021-09-14

Similar Documents

Publication Publication Date Title
CN105138787B (en) The supersonic flow field design method of feature based line tracking
CN102722611B (en) Method for carrying out parallelization numerical simulation on hydrodynamic force conditions of river provided with cascade hydropower station
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
CN113158451A (en) Large-area river three-dimensional simulation method based on one-dimensional flood routing model
CN114117877B (en) Topological optimization method based on isogeometric particle description
CN111931263B (en) Constant flow simulation method based on optimized solution of planar two-dimensional shallow water equation
US6405142B1 (en) Fluid analyzer and program recording medium
CN104572575B (en) A kind of especially big deformation dynamics mess generation method
CN105224726B (en) The method that structured grid Dynamic mesh is used for unstrctured grid flow field calculation device
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
JP2011065360A (en) Flow numerical analysis method for setting boundary condition of unramified and non-orthogonal structural grid and using iterative computation
CN115618635A (en) Method and equipment for calculating working parameters of water cooling system of energy storage container
CN114818531A (en) Motion interface tracking numerical dissipation calculation method
CN114417675A (en) Finite element calculation method for continuous casting, solidification and heat transfer of special-shaped blank
CN103699748B (en) Skid chemical plant frame hoisting decorates method for determining position
CN109224796B (en) Flue flow equalizer device and design method
CN109408927B (en) Two-dimensional static magnetic field parallel finite element acceleration method based on black box transmission line model
CN108682046B (en) Suspended nozzle calculation method with flexible base material
CN112560247B (en) Automatic generation method of similar-appearance structural grid based on reference grid
CN108595726B (en) Riverbed adjusting method based on underwater repose angle of sediment
Avvari et al. Flow apportionment algorithm for optimization of power plant ducting
CN110427700A (en) Interpolation fitting method is adsorbed for the grid of tail flow field irregular three-D scatterplot coordinate
CN111144008A (en) Streamline tracking method based on high-order velocity field fitting

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