A FVM simulation method for free surface transient flow considering unsteady friction model

Abstract

The present invention discloses a FVM simulation method for transient free surface flow considering an unsteady friction model, comprising the following steps: constructing a transient free surface flow control equation for a pipeline including a Bruno unsteady friction model; meshing the calculation area according to the calculation characteristics of the finite volume method; using the finite volume method to calculate the free surface flow including the Bruno unsteady friction model; and performing model evaluation and verification. The present invention utilizes the FVM of the improved Bruno unsteady friction model to simulate the transient problem of the free surface of the pipeline, which can solve the problems of mass non-conservation, large simulation dissipation when the Courant number is less than 1.0, and inaccurate constant friction model in the existing characteristic line MOC format, and can also achieve higher accuracy, so as to achieve real-time and accurate simulation in the actual system, which is of great significance to the accuracy, stability and efficiency of real-time control simulation of the actual pipeline system.

Description

Technical Field

[0001] The present invention belongs to the field of hydraulic numerical simulation and calculation technology, and relates to hydraulics of hydropower stations (pumping stations), sewage pipes, and urban water supply (drainage) pipe networks. Specifically, it relates to an FVM simulation method for free surface transient flow considering an unsteady friction model. Background Art

[0002] Free surface flow may often occur in various piping systems, including sewage pipes, power plant pipes, and urban water supply pipes. An efficient and accurate free surface flow model is a key factor in achieving real-time control of piping systems.

[0003] Many scholars have made great contributions to the study and dissemination of free surface flows. The numerical solutions range from kinematic waves to dynamic waves in the governing equations. Most of the existing numerical models for free surface transient flows are based on the method of characteristic lines (MOC) in an implicit format. MOC is easy to solve complex boundary conditions. However, when the Courant number is less than 1.0, interpolation is inevitable, which often leads to numerical attenuation of the crest water depth and delayed response time of the water surface line. Importantly, the MOC scheme cannot guarantee the conservation of the water body when calculating free surface transient flows. At the same time, existing models generally only consider steady-state friction in the numerical simulation of free surface transient pipe flows. However, in fast transient flows, the steady-state friction approximation cannot accurately simulate the subsequent water depth attenuation. Therefore, it is crucial to develop alternative solution schemes that are easy to handle boundary conditions and have robustness and high resolution. Summary of the invention

[0004] Purpose of the invention: In order to overcome the problems in the prior art of numerical simulation of free surface transient pipelines, such as mass non-conservation in the existing characteristic line method, large simulation dissipation when the Courant number is less than 1.0, and inaccurate constant friction model, a FVM (finite volume method) simulation method for free surface transient flow considering an unsteady friction model is provided. The research considers the finite volume solution of the Bruno unsteady friction model to solve the free surface transient flow, which can solve the problem of mass non-conservation in the existing characteristic line MOC format, and can also achieve higher accuracy, so as to achieve real-time and accurate simulation in the actual system, which is of great significance to the accuracy, stability and efficiency of real-time control simulation of the actual pipeline system.

[0005] Technical solution: To achieve the above object, the present invention provides a FVM simulation method for free surface transient flow considering an unsteady friction model, comprising the following steps:

[0006] S1: Construct the governing equations of transient flow on the free surface of the pipe including the Bruno unsteady friction model;

[0007] S2: According to the calculation characteristics of the finite volume method, the calculation area is divided into grids;

[0008] S3: Based on the meshing in step S2, the finite volume method is used to calculate the free surface flow including the Bruno unsteady friction model;

[0009] S4: Based on the calculation results of step S3, perform model evaluation and verification.

[0010] Furthermore, the transient flow control equation of the pipeline free surface in which the Bruno unsteady friction model is introduced in step S1 is:

[0011]

[0012] The variables W, flux G, and source term S can be written as:

[0013]

[0014] In the formula, Area is the cross-sectional area of ​​the water; Q is the flow rate; is the average pressure of the water column on the cross section; ρ is the density of the liquid; g is the acceleration due to gravity; Slope is the slope of the pipe; S f is the pipe friction;

[0015] The improved Bruno unsteady friction model is introduced into the control equation. The specific model is as follows:

[0016] S f =S q +S b #(3)

[0017] In the formula, S q is the steady-state friction term; S b is the unsteady friction term. According to the Manning formula, the steady-state friction term S q It can be defined as:

[0018]

[0019] Where n is the Manning roughness coefficient of the pipe, and R is the hydraulic radius;

[0020] The improved Runeau unsteady friction model proposes the following method to calculate the unsteady friction term S b The formula is:

[0021]

[0022] Where c is the wave velocity; if v<0, SGN(v)=-1; if v≥0, SGN(v)=1; f is the friction coefficient.

[0023] The Waddi shear attenuation coefficient B* is used to predict f analytically:

[0024]

[0025] Furthermore, in step S2, the calculation area is gridded. Gridding is a necessary step for finite volume method calculation. The calculation domain is gridded to divide the calculation domain into calculation units. The specific method is as follows: the entire calculation domain is gridded by finite volume, and a calculation grid of the finite volume system of the free surface flow of the pipeline is established. The free surface flow pipeline space domain is discretized into N units, each unit length is Δx, the time domain is discretized into a time interval Δt, the midpoint of the control body is numbered i, the left and right boundaries of the i-th control body are i-1 / 2 and i+1 / 2 respectively, and the local discrete grid of the calculation model is used.

[0026] Furthermore, the integral when solving by the finite volume method in step S3 is specifically:

[0027] Integrate equation (1) along the x direction of the spatial domain from the boundary i-1 / 2 to the boundary i+1 / 2, and substitute get:

[0028]

[0029] In the formula: superscript n+1 is the time t+Δt; n is the time t.

[0030] Furthermore, the calculation of the free surface flow of the Bruno unsteady friction model in step S3 specifically includes the following steps:

[0031] A1: Second-order Accuracy Reconstruction

[0032] The center-upwind format method is used for second-order linear reconstruction, so as to obtain the second-order precision FVM format in space. At the same time, the MINMOD slope limiter function is introduced to avoid the false oscillation phenomenon during linear interpolation reconstruction.

[0033] The second-order precision reconstruction specifically includes the following steps:

[0034] A1-1: Reconstruction

[0035]

[0036]

[0037] Where:

[0038] A1-2: Time Advancement

[0039]

[0040]

[0041] A2: Solution for the double shock wave hypothesis

[0042] The double shock wave pseudo-solver is used to solve the problem by assuming that there are shock waves on both sides of the calculation area. Since there is no rarefaction wave involved, the calculation is simple to understand, but the result is reliable. The double shock wave pseudo-solver is used for approximate solution. Combined with the conservation law, the flux G in the double shock wave pseudo-solver is given. i+1 / 2 ;

[0043]

[0044] In the formula, O L , O R are the wave speeds measured on the left and right sides respectively; G L and G R are the fluxes on the left and right sides of the calculation unit interface, respectively. L and W R It is the variable of left and right measurement of the interface;

[0045] O L =v L -E L #(11)

[0046] O R =v R +E R #(12)

[0047] In the formula, E K (K = L, R) can be expressed as:

[0048]

[0049] Area * is the estimated value of the area of ​​the star region; To calculate the centroid pressure of the star region water section; To calculate the centroid pressure of the water-passing cross-section on the left and right sides of the star region; c K (K = L, R) is used to calculate the velocity of gravity waves on the left and right sides of the star region; for a circular channel, the variables c, Area and It can be expressed as:

[0050]

[0051]

[0052]

[0053] Where d is the pipe diameter; φ is the central angle of the water surface; h is the water depth;

[0054] For Area * The solution is based on the following formula for Area * The value of is estimated:

[0055] Area * =0.5[1+(V L -V R ) / (γ L +γ R )]Area L +Area R )#(17)

[0056] Where γ = ∫(c / A)dA; for a circular cross-section pipe, γ can be expressed as follows:

[0057]

[0058] A3: Solving boundary conditions

[0059] All computational domains are calculated using Equation (7), including the upstream and downstream boundaries G 1 / 2 and G N+1 / 2 ; For each G i+1 / 2 , linear reconstruction using the center-upwind format method requires data on the left and right sides of the interface, and two calculation units are added on both sides of the calculation domain, that is, a total of four calculation units are added, and all parameter values ​​of these four calculation units are the same as the values ​​at the original boundary; compared with the existing method, this method can also achieve second-order accuracy at the boundary. At the same time, the boundary processing method proposed in the present invention processes the boundary and internal nodes uniformly when calculating the boundary, making the processing calculation simpler.

[0060] Based on the boundary processing method, the computational domain is expanded from 1 to N to -1 to N+2. All variables in the computational domain -1 to N+2 at the previous moment are known. By reconstructing the elements of the computational domain 1 to N, the variable values ​​of the upper and lower boundaries can be obtained.

[0061]

[0062]

[0063] For a circular cross-section channel, if the inflow is specified at the upstream boundary, and For The unknown variables involved; therefore, there is only one unknown variable in equation (19): Directly find the and The downstream boundary can be found in the same way;

[0064] In order to improve the accuracy of the boundary, two calculation units are introduced on both sides of the boundary. Their parameters are consistent with the boundary parameters, so that other errors are not introduced while improving the accuracy. The specific expression of the boundary unit is given:

[0065]

[0066]

[0067] A4: Considering the time integration of Bruno's unsteady friction model

[0068] G i-1 / 2 and G i+1 / 2 The second-order format of has been obtained in the first three steps. In order to advance the solution from time n to time n+1, time integration is required. In the presence of Bruno unsteady friction, the following time integration is performed to introduce the Bruno unsteady friction model into the solution

[0069]

[0070]

[0071]

[0072] A5: Stability Constraints

[0073] Stability constraints include Courant number constraints, taking into account the source term constraints of the Bruno unsteady friction model

[0074]

[0075]

[0076] The overall stability constraint is:

[0077] Δt max =min i＝1…n [Δt S,max ,Δt Cr,max ]#(28).

[0078] When solving the free surface transient flow with the characteristic line method MOC, when the Courant number is less than 1.0, spatial interpolation is required to achieve the effect of smooth damping and diffusion of waves and wavefronts. Diffusion enables the wave to reach the boundary earlier, and damping reduces the wave crest, artificially delaying the time when overload occurs. In combination with the above scheme, the present invention adopts the finite volume method of conservation form for solution. The high-order accuracy of the finite volume method is insensitive to the Courant number and can maintain mass conservation. Therefore, even when the Courant number is small, high accuracy can be maintained. At the same time, the present invention also takes into account the Bruno unsteady friction resistance, which is more accurate when solving the free surface flow than only considering constant friction.

[0079] Beneficial effects: Compared with the prior art, the present invention applies the unsteady friction model to the FVM simulation method of free surface transient flow, which can solve the problems of inaccurate calculation, long calculation time, complex solution, and difficulty in real-time prediction in the existing one-dimensional transient flow model. In particular, it solves the problems of mass non-conservation in the characteristic line method, large simulation dissipation when the Courant number is less than 1.0, and inaccurate constant friction model. It can achieve higher simulation accuracy, and the simulation results are closer to the actual situation, so as to achieve real-time and accurate simulation in the actual system, which is of great significance to the accuracy, stability and efficiency of real-time control simulation of the actual pipeline system. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 It is a flow chart of the implementation of the present invention;

[0081] Figure 2 is a finite volume method meshing diagram of a free surface flow conduit of the present invention;

[0082] Figure 3 is a cross-sectional view of a circular channel of a free surface flow conduit of the present invention;

[0083] Figure 4 It is a comparison diagram of the water depth result after the water depth of Example 1 tends to be stable, the simulation result of the characteristic line MOC method, and the water depth of the pipeline simulated by the FVM of the present invention at x=2.50m;

[0084] Figure 5 It is the accurate solution of water depth in Example 1, the simulation result of the characteristic line MOC method, and the comparison diagram of water depth along the pipeline at the time of 197.0s simulated by the FVM format of the present invention;

[0085] Figure 6 This is the FVM format of Example 2 of the present invention considering Bruno's non-constant friction, only considering constant friction, and the water depth comparison diagram of Ackers 1964 experimental values ​​at the pipeline x=77.940m. DETAILED DESCRIPTION

[0086] The present invention is further explained below in conjunction with the accompanying drawings and specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention. After reading the present invention, various equivalent forms of modifications to the present invention by those skilled in the art all fall within the scope defined by the claims attached to this application.

[0087] The present invention provides a FVM simulation method for free surface transient flow considering an unsteady friction model, referring to Figure 1 , which comprises the following steps:

[0088] S1: Construct the governing equations of transient flow on the free surface of the pipe including the Bruno unsteady friction model:

[0089] The governing equation for transient flow on the free surface of a pipe with the Bruno unsteady friction model is:

[0090]

[0091] The variables W, flux G, and source term S can be written as:

[0092]

[0093] In the formula, Area is the cross-sectional area of ​​the water; Q is the flow rate; is the average pressure of the water column on the cross section; ρ is the density of the liquid; g is the acceleration due to gravity; Slope is the slope of the pipe; S f is the pipe friction;

[0094] The improved Bruno unsteady friction model is introduced into the control equation. The specific model is as follows:

[0095] S f =S q +S b #(3)

[0096] In the formula, S q is the steady-state friction term; S b is the unsteady friction term. According to the Manning formula, the steady-state friction term S q It can be defined as:

[0097]

[0098] Where n is the Manning roughness coefficient of the pipe, and R is the hydraulic radius;

[0099] The improved Runeau unsteady friction model proposes the following method to calculate the unsteady friction term S b The formula is:

[0100]

[0101] Where c is the wave velocity; if v<0, SGN(v)=-1; if v≥0, SGN(v)=1; f is the friction coefficient.

[0102] The Waddi shear attenuation coefficient B* is used to predict f analytically:

[0103]

[0104] S2: According to the calculation characteristics of the finite volume method, the calculation area is divided into grids:

[0105] Meshing is a necessary step for finite volume method calculation. The computational domain is meshed and divided into computational units. The specific method is as follows: the entire computational domain is meshed with finite volume, and the computational grid of the finite volume system of free surface flow of the pipeline is established. The spatial domain of the free surface flow pipeline is discretized into N units, each unit length is Δx, the time domain is discretized into the time interval Δt, the midpoint of the control body is numbered i, the left and right boundaries of the i-th control body are i-1 / 2 and i+1 / 2 respectively, and the local discrete grid of the computational model is calculated as follows: Figure 2 shown.

[0106] S3: Free surface flow calculation including Bruno unsteady friction model using finite volume method:

[0107] Integrate equation (1) along the x direction of the spatial domain from the boundary i-1 / 2 to the boundary i+1 / 2, and substitute get:

[0108]

[0109] In the formula: superscript n+1 is the time t+Δt; n is the time t.

[0110] The calculation of free surface flow including the Bruno unsteady friction model consists of the following steps:

[0111] A1: Second-order Accuracy Reconstruction

[0112] The center-upwind format method is used for second-order linear reconstruction, so as to obtain the second-order precision FVM format in space. At the same time, the MINMOD slope limiter function is introduced to avoid the false oscillation phenomenon during linear interpolation reconstruction.

[0113] The second-order precision reconstruction specifically includes the following steps:

[0114] A1-1: Reconstruction

[0115]

[0116]

[0117] Where:

[0118] A1-2: Time Advancement

[0119]

[0120]

[0121] A2: Solution for the double shock wave hypothesis

[0122] The double shock wave pseudo-solver is used to solve the problem by assuming that there are shock waves on both sides of the calculation area. Since there is no rarefaction wave involved, the calculation is simple to understand, but the result is reliable. The double shock wave pseudo-solver is used for approximate solution. Combined with the conservation law, the flux G in the double shock wave pseudo-solver is given. i+1 / 2 ;

[0123]

[0124] In the formula, O L , O R are the wave speeds measured on the left and right sides respectively; G L and G R are the fluxes on the left and right sides of the calculation unit interface, respectively. L and W R It is the variable of left and right measurement of the interface;

[0125] O L =v L -E L #(11)

[0126] O R =v R +E R #(12)

[0127] In the formula, E K (K = L, R) can be expressed as:

[0128]

[0129] Area * is the estimated value of the area of ​​the star region; To calculate the centroid pressure of the star region water section; To calculate the centroid pressure of the water-passing cross-section on the left and right sides of the star region; c K (K = L, R) is the velocity of gravity waves on the left and right sides of the star region; for a circular channel, Figure 3 As shown, variables c, Area, and It can be expressed as:

[0130]

[0131]

[0132]

[0133] Where d is the pipe diameter; φ is the central angle of the water surface; h is the water depth;

[0134] For Area * The solution is based on the following formula for Area * The value of is estimated:

[0135] Area * =0.5[1+(V L -V R ) / (γ L +γ R )](Area L +Area R )#(17)

[0136] Where γ = ∫(c / A)dA; for a circular cross-section pipe, γ can be expressed as follows:

[0137]

[0138] A3: Solving boundary conditions

[0139] All computational domains are calculated using Equation (7), including the upstream and downstream boundaries G 1 / 2 and G N+1 / 2 ; For each G i+1 / 2 , linear reconstruction using the center-upwind format method requires data on the left and right sides of the interface, and two calculation units are added on both sides of the calculation domain, that is, a total of four calculation units are added, and all parameter values ​​of these four calculation units are the same as the values ​​at the original boundary; compared with the existing method, this method can also achieve second-order accuracy at the boundary. At the same time, the boundary processing method proposed in the present invention processes the boundary and internal nodes uniformly when calculating the boundary, making the processing calculation simpler.

[0140] Based on the boundary processing method, the computational domain is expanded from 1 to N to -1 to N+2. All variables in the computational domain -1 to N+2 at the previous moment are known. By reconstructing the elements of the computational domain 1 to N, the variable values ​​of the upper and lower boundaries can be obtained.

[0141]

[0142]

[0143] For a circular cross-section channel, if the inflow is specified at the upstream boundary, and For Related unknown variables; therefore, there is only one unknown variable in equation (19) Directly find the and The downstream boundary can be found in the same way;

[0144] In order to improve the accuracy of the boundary, two calculation units are introduced on both sides of the boundary. Their parameters are consistent with the boundary parameters, so that other errors are not introduced while improving the accuracy. The specific expression of the boundary unit is given:

[0145]

[0146]

[0147] A4: Considering the time integration of Bruno's unsteady friction model

[0148] G i-1 / 2 and G i+1 / 2 The second-order format of has been obtained in the first three steps. In order to advance the solution from time n to time n+1, time integration is required. In the presence of Bruno unsteady friction, the following time integration is performed to introduce the Bruno unsteady friction model into the solution

[0149]

[0150]

[0151]

[0152] A5: Stability Constraints

[0153] Stability constraints include Courant number constraints, taking into account the source term constraints of the Bruno unsteady friction model

[0154]

[0155]

[0156] The overall stability constraint is:

[0157] Δt max =min i＝1…n [Δt s,max ,Δt Cr,max ]#(28).

[0158] S4: Conduct model evaluation and verification.

[0159] Based on the above scheme, in order to verify the actual simulation effect of the method of the present invention, Example 1 and Example 2 were carried out respectively, as follows:

[0160] Embodiment 1:

[0161] In this embodiment, the water depth value after transient stabilization is selected without considering friction, and the simulation value of the method proposed in the present invention and the simulation value of the MOC method are compared to verify the conservation of the FVM model of the present invention.

[0162] The exact solution is obtained without considering the friction, and the method of the present invention is compared with the characteristic line MOC method to verify the effectiveness and conservation of the method of the present invention.

[0163] In Example 1, the length of the pipeline is 1000.0m, the inner diameter is 15.0m, the slope is 0.00, and the Manning roughness coefficient is 0.00. There is a closed gate in the middle of the pipeline, the water depth upstream of the gate is 10.0m, and the water depth downstream is 3.0m. The change of the water depth of the free surface flow in the pipeline is simulated when the middle gate is opened instantly.

[0164] Finally, the solution obtained by the method of the present invention is compared with the result of the characteristic line MOC method, and the following results are obtained respectively: Figure 4 and Figure 5 The pressure comparison curve shown is Figure 4 It can be seen that the model proposed by the present invention is compared with the actual stable water depth and the calculation results of MOC. After the water depth tends to be stable, the water depth simulated by the MOC method is lower than the actual stable water depth, while the model proposed by the present invention has the correct stable water depth. It can be seen that the model established by the method of the present invention satisfies the law of conservation of mass. On the contrary, the characteristic line MOC method cannot guarantee the conservation of mass. Figure 5 It can be seen that when the Courant number is less than 1.0, MOC dissipates more and the simulation error in the time dimension is large.

[0165] Embodiment 2:

[0166] The results of the free surface pipeline transient simulation experiment by Ackers in 1964 were selected to verify the effectiveness of the FVM simulation method of the transient flow of the free surface of the pipeline based on the Bruno unsteady friction model. The pipeline length is 304.80m, the inner diameter is 0.3048m, the slope is 0.001, and the Manning roughness coefficient is 0.0115. Considering Bruno unsteady friction, the kinematic viscosity of the fluid is 0.8937 10-6m2 / s. The initial base flow water depth of the system is 0.0768m, and the flow rate is 0.004984m3 / s. The upstream boundary flow flow is in a trapezoidal form, with a total flow time of 132s, a maximum flow of 0.018689m3 / s, and a maximum flow duration of 12.0s. In this case, transients will occur in the free surface pipeline as the upstream inflow changes.

[0167] The calculation results of the method of the present invention are compared with the transient simulation results of free surface pipeline by Ackers in 1964, and the following results are obtained: Figure 6 The comparison curve shown is Figure 6 It can be seen that compared with the model without considering the Bruno unsteady friction, the Bruno unsteady friction model can accurately describe the transient pipeline pressure amplitude and attenuation, especially in the time dimension, which is better than only considering the Bruno unsteady friction.

[0168] According to the results of Examples 1 and 2, the method of the present invention can well solve the problems of mass non-conservation and numerical dissipation in the MOC format during the transient process of the free surface of the pipeline, and the simulation results are also more accurate.

Claims

1. A FVM simulation method for transient free surface flow considering an unsteady friction model, characterized in that: The steps include: S1: Construct the governing equations of transient flow on the free surface of the pipe including the Bruno unsteady friction model; S2: According to the calculation characteristics of the finite volume method, the calculation area is divided into grids; S3: Based on the meshing in step S2, the finite volume method is used to calculate the free surface flow including the Bruno unsteady friction model; S4: Perform model evaluation and verification according to the calculation results of step S3; The transient flow control equation of the pipeline free surface with the Bruno unsteady friction model introduced in step S1 is: The variables W, flux G, and source term S are written as: In the formula, Area is the cross-sectional area of ​​the water; Q is the flow rate; is the average pressure of the water column on the cross section; ρ is the density of the liquid; g is the acceleration due to gravity; Slope is the slope of the pipe; S f is the pipe friction; The improved Bruno unsteady friction model is introduced into the control equation. The specific model is as follows: S f =S q +S b (3) In the formula, S q is the steady-state friction term; S b is the unsteady friction term; In step S1, according to the Manning formula, the steady-state friction term S q Defined as: Where n is the Manning roughness coefficient of the pipe, and R is the hydraulic radius; The improved Runeau unsteady friction model proposes the following method to calculate the unsteady friction term S b The formula is: Where c is the wave velocity; if v<0, SGN(v)=-1, if v≥0, SGN(v)=1; f is the friction coefficient; In step S1, the Waddi shear attenuation coefficient B* is used to analytically predict f:

2. The FVM simulation method of free surface transient flow considering an unsteady friction model according to claim 1, characterized in that: The specific method of meshing the calculation area in step S2 is: divide the entire calculation domain into finite volume meshes, establish a calculation grid for the finite volume system of the free surface flow of the pipeline, discretize the free surface flow pipeline space domain into N units, each unit length is Δx, the time domain is discretized into a time interval Δt, the midpoint of the control body is numbered i, the left and right boundaries of the i-th control body are i-1 / 2 and i+1 / 2 respectively, and the local discrete grid of the calculation model is calculated.

3. The FVM simulation method of free surface transient flow considering an unsteady friction model according to claim 1, characterized in that: The integral when solving by the finite volume method in step S3 is specifically: Integrate equation (1) along the x direction of the spatial domain from the boundary i-1 / 2 to the boundary i+1 / 2, and substitute get: In the formula: superscript n+1 is the time t+Δt; n is the time t.

4. The FVM simulation method of free surface transient flow considering an unsteady friction model according to claim 3, characterized in that: The calculation of the free surface flow of the Bruno unsteady friction model in step S3 specifically includes the following steps: A1: Second-order Accuracy Reconstruction The center-upwind format method is used for second-order linear reconstruction, so as to obtain the second-order precision FVM format in space. At the same time, the MINMOD slope limiter function is introduced to avoid the false oscillation phenomenon during linear interpolation reconstruction. The second-order precision reconstruction specifically includes the following steps: A1-1: Reconstruction Where: A1-2: Time Advancement A2: Solution for the double shock wave hypothesis The double shock pseudo solver is used for approximate solution. Combined with the conservation law, the flux G in the double shock pseudo solver is given. i+1 / 2 ; In the formula, O L , O R are the wave speeds measured on the left and right sides respectively; G L and G R are the fluxes on the left and right sides of the calculation unit interface respectively; W L and W R It is the variable of left and right measurement of the interface; About L =v L -E L (11) About R =v R +E R (12) In the formula, E K (K = L, R) is expressed as: Where Area* is the estimated value of the star region area; To calculate the centroid pressure of the star region water section; To calculate the centroid pressure of the water-passing cross-section on the left and right sides of the star region; c K (K = L, R) is used to calculate the velocity of gravity waves on the left and right sides of the star region; for a circular channel, the variables c, Area and It is expressed as: Where d is the pipe diameter; φ is the central angle of the water surface; h is the water depth; For the solution of Area*, the value of Area* is estimated based on the following formula: Area * =0.5[1+(V L V R ) / (γ L +γ E )](Area L +Area R ) (17) Where γ = ∫(c / A)dA; for a circular cross-section pipe, γ is expressed as follows: A3: Solving boundary conditions All computational domains are calculated using formula (7), including the upstream and downstream boundaries G 1 / 2 and G N+1 / 2 ; For each G i+1 / 2 ,Using the center-upwind format method for linear reconstruction requires data on the left and right sides of the interface, and two calculation units are added on both sides of the calculation domain, that is, a total of four calculation units are added, and all parameter values ​​of these four calculation units are the same as those at the original boundary; Based on the boundary processing method, the computational domain is expanded from 1 to N to -1 to N+2. All variables in the computational domain -1 to N+2 at the previous moment are known. The elements of the computational domain 1 to N are reconstructed to obtain the variable values ​​of the upper and lower boundaries. For a circular cross-section channel, if the inflow is specified at the upstream boundary, and For The unknown variables involved; therefore, there is only one unknown variable in equation (19): Directly find the and The downstream boundary is found in the same way; In order to improve the accuracy of the boundary, two calculation units are introduced on both sides of the boundary. Their parameters are consistent with the boundary parameters, so that other errors are not introduced while improving the accuracy. The specific expression of the boundary unit is given: A4: Considering the time integration of Bruno's unsteady friction model G i-1 / 2 and G i+1 / 2 The second-order format of has been obtained in the first three steps. In order to advance the solution from time n to time n+1, time integration is required. In the presence of Bruno unsteady friction, the following time integration is performed to introduce the Bruno unsteady friction model into the solution A5: Stability Constraints Stability constraints include Courant number constraints, taking into account the source term constraints of the Bruno unsteady friction model The overall stability constraint is: Δt max =min i＝1....n [Δt s,max ,Δt Cr,max ] (28)。

Images