Hydraulic design of restricted orifice surge tank

(25 March 2021)

This document describes design formulas and design example of restricted orifice surge tank.

Introduction

Function of surge tank

In general, a long pressure tunnel, the headrace / tailrace of pumped storage power plant needs a surge tank. The necessity of surge tank is judged whether length of pressure tunnel exceeds 500m or not. A surge tank has following functions.

To absorb abnormal pressure increase in the tunnel by water hammering when load rejection or input power rejection.
To absorb the abnormal pressure decrease in the tunnel by supplying water when rapid increase of discharge.

Surge tanks are categorized in following 4 types.

Simple surge tank
Differential surge tank
Restricted orifice surge tank (An orifice of a surge tank is called a ”port”)
Surge tank with chamber

These days Restricted Orifice Surge Tanks are usually adopted. This type can effectively attenuate amplitude of the water level in the tank and has comparatively simple design. In addition, in the case that a surge tank is constructed deep under the ground, such as a tailrace surge tank, the diameter of vertical shaft can be reduced by adding upper chamber.

Requirements for Hydraulic Design of Restricted Orifice Surge Tank

The following conditions are required for the design of Restricted Orifice Surge Tank.

The maximum amplitude of water level by surging is between the elevation of top and bottom of the tank.
To meet the required static and dynamic stability conditions of water level amplitude.
Maximum discharge does not exceed the following critical discharge.

Fundamental differential equations

Equation of motion

Euler's equation of fluid motion in one-dimensional problem is as follows.

$\begin{equation} \cfrac{dv}{dt}=-\cfrac{1}{\rho}\cfrac{dp}{dx}=-g\cfrac{dh}{dx} \end{equation}$

$\begin{align} &v & & \text{flow velocity} \\ &t & & \text{time} \\ &\rho & & \text{density of water} \\ &p & & \text{pressure} \\ &x & & \text{position} \\ &h & & \text{water head} \\ &g & & \text{gravity acceleration} \end{align}$

When the flow direction from reservoir to surge tank is positive, and the direction of water surface movement from up to dawn on the basis of reservoir water level in a surge tank is positive, positive gradient of water head is upward movement of water level in the surge tank.

$\begin{equation} \cfrac{dv}{dt}=-g\cfrac{dh}{dx}=g\cfrac{z}{L} \end{equation}$

$\begin{align} &z & & \text{variation of the water level in surge tank} \\ & & & \text{(positive is the downward from reservoir water level) } \\ &L & & \text{length of the pressure tunnel} \\ &v & & \text{flow velocity of the pressure tunnel} \\ \end{align}$

If head loss $\Delta h$ exists, the acceleration of the water in the pressure tunnel is reduced by $\Delta h$ .

$\begin{equation} \cfrac{dv}{dt}=g\cfrac{z}{L}-g\cfrac{\Delta h}{L}=\cfrac{z-\Delta h}{L/g} \end{equation}$

Because friction loss ([tex:h=cvv]) in the pressure tunnel and resistance loss of the port ( $k=v_p*v_p/2/g$ ) are caused in the case of Restricted Orifice Surge Tank, the equation of motion is given as follows taking into consideration the flow direction.

$\begin{equation} \cfrac{dv}{dt}=\cfrac{z-\Delta h}{L/g}=\cfrac{z-c\cdot |v|\cdot v-k}{L/g} \end{equation}$

Equation of continuity

When all flow direction of tunnel discharge, discharge in a surge tank and discharge of turbine is assumed positive, equation of continuity is given as follows.

$\begin{equation} Q=f \cdot v+F \cfrac{dz}{dt} \quad \rightarrow \quad \cfrac{dz}{dt}=\cfrac{Q-f \cdot v}{F} \end{equation}$

$\begin{align} &Q & & \text{discharge of turbine} \\ &f & & \text{section area of pressure tunnel} \\ &F & & \text{sectionarea of surge tank shaft} \end{align}$

Resistance of port

Because the port velocity ( $v_p$ ) is what discharge in a surge tank is divided by the effective area of the port, the resistance of the port ( $k$ ) can be given as follows, considering the flow direction.

$\begin{equation} v_p=\cfrac{Q-f \cdot v}{C_d \cdot F_p} \quad \rightarrow \quad k=\cfrac{|v_p| \cdot v_p}{2 g}=\cfrac{1}{2 g}\left|\cfrac{f \cdot v-Q}{C_d \cdot F_p}\right|\cdot \cfrac{f \cdot v-Q}{C_d \cdot F_p} \end{equation}$

Fundamental differential equations

$\begin{align} &\text{Equation of motion : } &\quad& \cfrac{dv}{dt}=\cfrac{z-c\cdot |v|\cdot v-k}{L/g}\\ &\text{Continuous equation : } &\quad& \cfrac{dz}{dt}=\cfrac{Q-f\cdot v}{F}\\ &\text{Resistance of the port : } &\quad& k=\cfrac{|v_p|\cdot v_p}{2g}=\cfrac{1}{2g}\cdot\left|\cfrac{f\cdot v-Q}{C_d\cdot F_p}\right|\cdot\cfrac{f\cdot v-Q}{C_d\cdot F_p} \end{align}$

$\begin{align} &v_p & & \text{flow velocity of port} \\ &F_p & & \text{section area of port} \\ &C_d & & \text{discharge coefficient of port} \\ &v & & \text{flow velocity of pressure tunnel (positive is from reservoir to surge tank)} \\ &z & & \text{variation of water in surge tank} \\ &g & & \text{gravity acceleration} \\ &c & & \text{head loss coefficient} \\ &L & & \text{length of pressure tunnel (from reservoir to port)} \\ &f & & \text{section area of pressure tunnel} \\ &F & & \text{section area of surge tank shaft} \\ &Q & & \text{discharge at interception} \end{align}$

In the surging analysis, above simultaneous differential equations are solved by Runge-Kutta method and the time series of the water surface profile in the surge tank can be obtained.

Equations for basic design of restricted orifice surge tank

Maximum water level rise in surge tank

For the estimation of the maximum water level rise from reservoir water level at the total load interception in the restricted orifice surge tank, Vogt-Forchheimer's equations shown below have been used.
These equations are one of the simple solution of the fundamental differential equation under the condition of instant interception of initial discharge.

$\begin{equation} \begin{cases} &m'\cdot k_0<1 & &(1+m'\cdot z_m)-\ln(1+m'\cdot z_m)=(1+m'\cdot h_0)-\ln(1-m'\cdot k_0) \\ &m'\cdot k_0>1 & &(m'\cdot |z_m|-1)+\ln(m'\cdot |z_m|-1)=\ln(m'\cdot k_0-1)-(m'\cdot h_0+1) \end{cases} \end{equation}$

$\begin{align} &h_0=c\cdot {v_0}^2 \\ &k_0=\cfrac{1}{2 g}\left(\cfrac{Q_0}{C_d F_p}\right)^2 \\ &m'=\cfrac{2 g F(h_0+k_0)}{L f {v_0}^2} \end{align}$

$\begin{align} &z_m & & \text{maximum water level rise in surge tank} \\ & & & \text{(positive is downward from reservoir water level)} \\ &h_0 & & \text{total head loss of pressure tunnel} \\ &k_0 & & \text{resistance of port} \\ &v_0 & & \text{flow velocity of pressure tunnel} \\ &c & & \text{head loss coefficient} \\ &Q_0 & & \text{maximum discharge} \\ &F_p & & \text{section area of port} \\ &C_d & & \text{discharge coefficient of port} \\ &L & & \text{length of pressure tunnel (from reservoir to port)} \\ &f & & \text{section area of pressure tunnel} \\ &F & & \text{section area of surge tank shaft} \\ &g & & \text{gravity acceleration} \end{align}$

To obtain the maximum water level rise of $z_m$ , Bisection method can be applied considering below equations.

$\begin{equation} f(z)=0 \end{equation}$

$\begin{equation} f(z)= \begin{cases} \{(1+m'\cdot z)-\ln(1+m'\cdot z)\}-\{(1+m'\cdot h_0)-\ln(1-m'\cdot k_0)\} &(m'\cdot k_0<1) \\ \{(m'\cdot |z|-1)+\ln(m'\cdot |z|-1)\}-\{\ln(m'\cdot k_0-1)-(m'\cdot h_0+1)\} &(m'\cdot k_0>1) \end{cases} \end{equation}$

In bove equations, the following conditions have to be satisfied everytime. These are for maintaining the positive values in the logarithm paragraph in the function of $f(z)$ .

$\begin{equation} \begin{cases} &0 < 1+m'\cdot z & (m'\cdot k_0 < 1) \\ &0 < m'\cdot |z|-1 & (m'\cdot k_0>1) \end{cases} \end{equation}$

Considering above conditions, the minimum limit of $z$ in Bisection method can be set as followings.

$\begin{equation} \begin{cases} z_1=-\cfrac{1}{m'}+0.0001 & (m'\cdot k_0<1) \\ |z_1|=\cfrac{1}{m'}+0.0001 & (m'\cdot k_0>1) \end{cases} \end{equation}$

And the maximum water level rise in the restricted orifice surge tank have to be smaller than the maximum water level rise of free surging. Therefore, the maximum water level rise of free surging shown below can be set as the maximum limit of $z$ in Bisection method.

$\begin{equation} z_2=v_0\cdot\sqrt{\cfrac{L\cdot f}{g\cdot F}} \end{equation}$

Using $z_1$ and $z_2$ as limit values of Bisection method, Vogt-Forchheimer's equations can be solved.

For the calculation of the water level rise, following issues shall be considered.

To adopt the most rigid condition of water level of the reservoir corresponding to the discharge conditions.
To set the small roughness in the case of load imterception or input interception, and to set the large roughness in the case of load rapid increase.

The actual conditions based on METI standard considered are shown below.

Case	water level and roughness	Headrace surge tank	Tailrace surge tank
Load interception	Reservoir water level Checked water level in shaft Variation of roughness	HWL of upper reservoir Up surge WL $-0.0015$	LWL of lower reservoir Down surge WL $-0.0015$
Load rapidly increase	Reservoir water level Checked water level in shaft Variation of roughness	LWL of upper reservoir Down surge WL $+0.0015$	HWL of lower reservoir Up surge WL $+0.0015$
Input interception	Reservoir water level Checked water level in shaft Variation of roughness	LWL of upper reservoir Down surge WL $-0.0015$	HWL of lower reservoir Up surge WL $-0.0015$

note

Above values of roughness variable are for concrete lining.
In case of steel lining, roughness variable is $0.001$ instead of $0.0015$ for concrete.
In case of no lining, roughness variable is $0.003$ instead of $0.0015$ for concrete.

Requirement for stability of water level vibration

For the stability against the water level vibration, Thoma-Schuller's stability conditions shown below shall be satisfied.

$\begin{align} &\text{Static stability conditions : } &\quad& h_0<\cfrac{H_g}{3}\sim\cfrac{H_g}{6} \\ &\text{Dynamic stability conditions : } &\quad& F>\cfrac{L f}{c(1+\eta)g H_g}\sim\cfrac{L f}{2 c g(H_g-z_m)} \end{align}$

where, $H_g$ means the gross head and $\eta=k_0 / h_0$ .

Critical discharge

The critical discharge by Calame-Garden is defined as follow.

$\begin{equation} Q_c=\cfrac{1}{c}\left(\cfrac{1}{2 g}\cdot\cfrac{L f^3}{F \eta}\right)^{1/2} \end{equation}$

If the critical discharge of $Q_c$ is smaller than the maximum discharge of $Q_0$ , the maximum water level rise for the critical discharge of $Q_c$ shall be checked. However, the surge tank shaft diameter is usually defined to satisfy the condition of $Q_0 \lt Q_c$ .

Required conditions for section area of surge tank shaft

Considering the conditions for the stability and the critical discharge, the section area of surge tank shaft of $F$ is defined to satisfy the following conditions.

$\begin{align} &F > F_1=\cfrac{h_0 L f}{c (h_0+k_0) g H_g}\\ &F > F_2=\cfrac{L f}{2 c g (H_g-z_m)}\\ &F < F_3=\cfrac{h_0 L f^3}{2 g c^2 k_0 {Q_0}^2} \end{align}$

Port diameter and shaft diameter

The resistance of the port $|z'_m|$ can be calculated as follow.

$\begin{equation} |z'_{rm}|=k_0-h_0 \end{equation}$

The relationship between port resistance $|z′_{rm}|$ and the maximum water Level rise $|z_m|$ are shown below.

$\begin{align} &|z_m| < |z′_{rm}| & & \text{port diameter is small because of large port resistance} \\ &|z_m| = |z′_{rm}| & & \text{port diameter is optimal} \\ &|z_m| > |z′_{rm}| & & \text{port diameter is large because of small port resistance} \end{align}$

where, $|z_m|$ can be obtained from Vogt-Forchheimer's equations.

Usually, the port diameter is defined as a value which is little bit larger than the optimal point.

f:id:damyarou:20210325030459p:plain

Design of chamber type surge tank

Since the Vogt-Forchheimer's equations give the maximum water level rise for a sinmle shaft type surge tank, ingenuity is needed to define the shape of chamber type surge tank. In case that a chamber type surge tank is required, the chamber size can be defined refering following figures. It shall be noted an air tunnel should be installed for a chamber type surge tank at the top of chamber.

Simple shaft type	Chamber type

Numerical solution of fundamental differential equations

Basic equations of Runge-Kutta method

Runge-Kutta method is a numerical integration method for the first order ordinary differential equation. The basic equatuins are shown below.

$\begin{equation} \cfrac{dy}{dt}=f(t,y) \qquad \rightarrow \qquad y_{i+1}=y_i+\cfrac{\Delta y_0+2\Delta y_1+2\Delta y_2+\Delta y_3}{6} \end{equation}$

$\begin{align} &\Delta y_0=\Delta t\cdot f(\:t_i,\: y_i\:) \\ &\Delta y_1=\Delta t\cdot f(\:t_i+\cfrac{\Delta t}{2},\: y_i+\cfrac{\Delta y_0}{2}\:) \\ &\Delta y_2=\Delta t\cdot f(\:t_i+\cfrac{\Delta t}{2},\: y_i+\cfrac{\Delta y_1}{2}\:) \\ &\Delta y_3=\Delta t\cdot f(\:t_i+\Delta t,\: y_i+\Delta y_2\:) \end{align}$

Solution by Runge-Kutta method

Subject

Using the known parameters of $Q_0$ , $V_0$ and $Z_0$ at time $t_0$ , to estimate the unknown parameters of $V$ and $Z$ at time $t=t_0+\Delta t$ , where the value of $Q$ at time $t=t_0+\Delta t$ is given (known) parameter.

Step of solving

Firstly, following functions are defined for the surging analysis.

$\begin{align} &q=Q(t) & & \text{(given value)} \\ &k=f_k(q,v) & & f_k(q,v)=\cfrac{1}{2g}\cdot \left|\cfrac{f\cdot v-q}{C_d\cdot F_p}\right|\cdot \cfrac{f\cdot v-q}{C_d\cdot F_p} \\ &\Delta v=\Delta t\cdot f_v(v,z,k) & & f_v(v,z,k)=\cfrac{z-c\cdot |v|\cdot v-k}{L/g} \\ &\Delta z=\Delta t\cdot f_z(q,v) & & f_z(q,v)=\cfrac{q-f\cdot v}{F} \end{align}$

Where,

$\begin{align} &q & & \text{discharge} \\ &k & & \text{resistance of port}\\ &\Delta v & & \text{increment of pressure tunnel flow velocity}\\ &\Delta z & & \text{increment of surge tank water level} \end{align}$

The calculation steps by Runge-Kutta method are shown.

Step-1

$\begin{align} &\Delta z_0=\Delta t\cdot f_z(Q_0, V_0) \\ &\text{(Judgement of in-flow or out-flow and set the discharge coefficient)}\\ &k_0=f_k(Q_0,V_0)\\ &\Delta v_0=\Delta t\cdot f_v(V_0, Z_0, k_0) \end{align}$

Step-2

$\begin{align} &v_1=V_0+1/2\cdot \Delta v_0 \\ &z_1=Z_0+1/2\cdot \Delta z_0 \\ &q_1=Q(t-1/2\cdot \Delta t) \\ &k_1=f_k(q_1,v_1) \\ &\Delta v_1=\Delta t\cdot f_v(v_1,z_1,k_1) \\ &\Delta z_1=\Delta t\cdot f_z(q_1,v_1) \end{align}$

Step-3

$\begin{align} &v_2=V_0+1/2\cdot \Delta v_2 \\ &z_2=Z_0+1/2\cdot \Delta z_1 \\ &k_2=f_k(q_1,v_2) \\ &\Delta v_2=\Delta t\cdot f_v(v_2,z_2,k_2) \\ &\Delta z_2=\Delta t\cdot f_z(q_1,v_2) \end{align}$

Step-4

$\begin{align} &v_3=V_0+\Delta v_2 \\ &z_3=Z_0+\Delta z_2 \\ &q_3=Q(t_0+\Delta t) \\ &k_3=f_k(q_3,v_3) \\ &\Delta v_3=\Delta t\cdot f_v(v_3,z_3,k_3) \\ &\Delta z_3=\Delta t\cdot f_z(q_3,v_3) \end{align}$

Step-5

$\begin{align} &Q=Q(t_0+\Delta t)=Q(t) \\ &V=V_0+1/6\cdot (\Delta v_0+2\Delta v_1+2\Delta v_2+\Delta v_3) \\ &Z=Z_0+1/6\cdot (\Delta z_0+2\Delta z_1+2\Delta z_2+\Delta z_3) \end{align}$

Design example

Headrace surge tank of PSPP

Preparation work

Since the head loss is calculated for the maximum plant discharge in many cases, the head loss for design of surge tank will be modified as shown below considering changes of the discharge and Manning's roughness coefficient also.

Item	Symbol	Unit	Load interception	Load increase	Input interception
Discharge	Q_org	m3/s	338	338	236.6
Reservoir water level	RWL	EL.m	1340 (FSL)	1315 (MOL)	1315 (MOL)
Roughness coefficient	n_org		0.013	0.013	0.013
Headloss of pressure tunnel	hl_org	m	13.065	13.065	13.065
Modified roughnee coefficient	n		0.0115 (-0.0015)	0.0145 (+0.0015)	0.0115 (-0.0015)
Modified head loss	hl	m	10.224	16.254	5.010

note

A midified head loss is calculated using following equation.

$\begin{equation} h_{l2}=h_{l1}\times\left(\cfrac{Q_2}{Q_1}\right)^2\times\left(\cfrac{n_2}{n_1}\right)^2 \end{equation}$

Where,

$h_{l1}$ is headloss at discharge of $Q_1$ and roughness coefficient of $n_1$ .
$h_{l2}$ is headloss at discharge of $Q_2$ and roughness coefficient of $n_2$ .

Actual calculations of modified head losses of pressure tunnel are shown below.

For load interception

$\begin{equation} h_l=13.065\times\left(\cfrac{338.0}{338.0}\right)^2\times\left(\cfrac{0.0115}{0.013}\right)^2=10.224 m \end{equation}$

For load increase

$\begin{equation} h_l=13.065\times\left(\cfrac{338.0}{338.0}\right)^2\times\left(\cfrac{0.0145}{0.013}\right)^2=16.254 m \end{equation}$

For input interception

$\begin{equation} h_l=13.065\times\left(\cfrac{236.6}{338.0}\right)^2\times\left(\cfrac{0.0115}{0.013}\right)^2=5.010 m \end{equation}$

Input data arrangement for surging analysis

Item	Symbol	Unit	Load interception	Load increase	Input interception
Gross head	Hg	m	677	677	677
Discharge	Q	m3/s	338 => 0	169 => 338	236.6 => 0
Closing time	Tc	sec	8	40	8
Reservoir water level	RWL	EL.m	1340 (FSL)	1315 (MOL)	1315 (MOL)
Diameter of pressure tunnel	d	m	8.2	8.2	8.2
Length of pressure tunnel	L	m	4800	4800	4800
Section area of pressure tunnel	f	m2	52.810	52.810	52.810
Average velocity of pressure tunnel	v0	m/s	6.400	6.400	4.478
Headloss of pressure tunnel	hl	m	10.224	16.254	5.010
Velocity head of pressure tunnel	hb	m	2.090	2.090	1.023
Total headloss	h0=hl+hb	m	12.314	18.344	6.033
Loss coefficient of pressure tunnel	c=h0/v0²		0.301	0.448	0.310
Diameter of shaft	DF	m	21.000	21.000	21.000
Diameter of port	Dp	m	4.500	4.500	4.500
Section area of shaft	F	m2	346.316	346.316	346.316
Section area of shaft	Fp	m2	15.904	15.904	15.904
Discharge coefficient of port	Cd		0.9	0.9	0.9

note

Items shown by bold font in above table are defined by the analysis using Vogt-Forchheimer's equations.
The discharge coefficient of port for inflow ( $C_{d,in}$ ) is not always the same as that for outflow ( $C_{d,out}$ ), and the values depend on the detailed shape of the port. However, the safety side results can be obtained, if the larger discharge coefficients are used. Therefore, $C_{d,in}=C_{d,out}=0.9$ has been used for some studies.

Selection of shaft diameter and port diameter

For the headrace surge tank, The shaft diameter of 21.0m and the port diameter of 4.5m are selected considering the target water level rise of 35m.

f:id:damyarou:20210325030715p:plain

* Free surge water rise
Ds=20    z*=58.075
Ds=21    z*=55.310
Ds=22    z*=52.796
Ds=23    z*=50.500
Ds=24    z*=48.396
* Requirements for stability
Hg/6=112.833m  # for static stabity
Dsr1= 7.0m    # for dynamic stability (1)
Dsr2= 9.2m    # for dynamic stability (2)
Dsr3=43.9m    # for critical discharge

Result of surging analysis of headrace surge tank

f:id:damyarou:20210325030736p:plain

Tailrace surge tank of PSPP

Preparation work

Since the head loss is calculated for the maximum plant discharge in many case, the head loss for design of surge tank will be modified as shown below considering changes of the discharge and Manning's roughness coefficient also.

Item	Symbol	Unit	Load interception	Load increase	Input interception
Discharge	Q_org	m3/s	338	338	236.6
Reservoir water level	RWL	EL.m	630 (MOL)	670 (FSL)	670 (FSL)
Roughness coefficient	n_org		0.013	0.013	0.013
Headloss of pressure tunnel	hl_org	m	5.151	5.151	5.151
Modified roughnee coefficient	n		0.0115 (-0.0015)	0.0145 (+0.0015)	0.0115 (-0.0015)
Modified head loss	hl	m	4.031	6.408	1.975

note

A midified head loss is calculated using following equation.

$\begin{equation} h_{l2}=h_{l1}\times\left(\cfrac{Q_2}{Q_1}\right)^2\times\left(\cfrac{n_2}{n_1}\right)^2 \end{equation}$

Where, + [tex:h{l1}] is headloss at discharge of $Q_1$ and roughness coefficient of $n_1$ . + [tex:h{l2}] is headloss at discharge of $Q_2$ and roughness coefficient of $n_2$ .

Actual calculations of modified head losses of pressure tunnel are shown below.

For load interception

$\begin{equation} h_l=5.151\times\left(\cfrac{338.0}{338.0}\right)^2\times\left(\cfrac{0.0115}{0.013}\right)^2=4.031 m \end{equation}$

For load increase

$\begin{equation} h_l=5.151\times\left(\cfrac{338.0}{338.0}\right)^2\times\left(\cfrac{0.0145}{0.013}\right)^2=6.408 m \end{equation}$

For input interception

$\begin{equation} h_l=5.151\times\left(\cfrac{236.6}{338.0}\right)^2\times\left(\cfrac{0.0115}{0.013}\right)^2=1.975 m \end{equation}$

Input data arrangement for surging analysis

Item	Symbol	Unit	Load interception	Load increase	Input interception
Gross head	Hg	m	677	677	677
Discharge	Q	m3/s	338 => 0	169 => 338	236.6 => 0
Closing time	Tc	sec	8	40	8
Reservoir water level	RWL	EL.m	630 (MOL)	670 (FSL)	670 (FSL)
Diameter of pressure tunnel	d	m	8.2	8.2	8.2
Length of pressure tunnel	L	m	1749	1749	1749
Section area of pressure tunnel	f	m2	52.810	52.810	52.810
Average velocity of pressure tunnel	v0	m/s	6.400	6.400	4.480
Headloss of pressure tunnel	hl	m	4.031	6.408	1.975
Velocity head of pressure tunnel	hb	m	2.090	2.090	1.024
Total headloss	h0=hl+hb	m	6.121	8.498	2.999
Loss coefficient of pressure tunnel	c=h0/v0²		0.149	0.207	0.149
Diameter of shaft	DF	m	10.000	10.000	10.000
Diameter of port	Dp	m	4.500	4.500	4.500
Section area of shaft	F	m2	78.540	78.540	78.540
Section area of shaft	Fp	m2	15.904	15.904	15.904
Discharge coefficient of port	Cd		0.9	0.9	0.9

note

Items shown by bold font in above table are defined by the analysis using Vogt-Forchheimer's equations.
The discharge coefficient of port for inflow ( $C_{d,in}$ ) is not always the same as that for outflow ( $C_{d,out}$ ), and the values depend on the detailed shape of the port. However, the safety side results can be obtained, if the larger discharge coefficients are used. Therefore, $C_{d,in}=C_{d,out}=0.9$ has been used for some studies.

Selection of shaft diameter and port diameter

For the tailrace surge tank, The shaft diameter of 10.0m and the port diameter of 4.5m are selected considering the required minimum shaft diameter of 8.0m which is defined to satisfy the dynamic stability condition. For the tailrace surge tank, the target water level rise is not considered, because it can have a upper chamber.

f:id:damyarou:20210325030808p:plain

* Free surge water rise
Ds= 8    z*=87.641
Ds= 9    z*=77.903
Ds=10    z*=70.113
Ds=11    z*=63.739
Ds=12    z*=58.427
* Requirements for stability
Hg/6=112.833m  # for static stabity
Dsr1= 4.6m    # for dynamic stability (1)
Dsr2= 8.0m    # for dynamic stability (2)
Dsr3=37.6m    # for critical discharge

Dimensions of upper chamber

Since the water level rise to the upsurge water level from FSL of the lower reservoir is 52.363m which is estimated by Vogt-Forchheimer's equations, the corresponding water volume can be calculated as follow.

$\begin{equation} V_1=52.363 \times \cfrac{\pi\cdot 10.0^2}{4}=4113 m^3 \end{equation}$

In case that the invert level of upper chamber is set to 10m higher than the FSL of lower reservoir and the chamber dimensions are assumed as 15.0m wide x 7.0m high x 40.0m long, the water volume corresponding to these volumes can be calculated as follows.

$\begin{equation} V_a=10.0 \times \cfrac{\pi\cdot 10.0^2}{4}=785 m^3 \end{equation}$

$\begin{equation} V_b=15.0 \times 7.0 \times 40.0 =4200 m^3 \end{equation}$

From above,

$\begin{equation} V_1=4113 m^3 < V_2=V_a + V_b = 4985 m^3 \end{equation}$

Therefore, assumed chamber dimensions can be accepted. The wall height of upper chamber will be set as 8.0m with 1.0m margin to water surface.

Result of surging analysis of headrace surge tank

f:id:damyarou:20210325030828p:plain

Programs

Maximum water level rize

#==============================
# Vogt-Forchheimer's equation
# Bisection Method
#==============================
import numpy as np
import matplotlib.pyplot as plt


# global variables
g=9.8       # gravity acceleration


def inp_JH():
    q0=338       # discharge of pressure tunnel
    cd=0.9      # discharge coefficient of port
    h0=12.314    # head loss of pressure tunnel
    tl=4800 # length of pressure tunnel including untake tunnel
    ft=52.810   # section area of pressure tunnel
    hg=677    # gross head
    dds=21     # design shaft diameter
    ddp=4.5     # design port diameter
    shaft=np.arange(20,25,1).astype(np.float64)      # shaft diameter: parameter
    port=np.arange(0.1,10.05,0.05).astype(np.float64) # port diameter: parameter
    tstr='PSPP-J Headrace Surge Tank (h0={0:.3f})'.format(h0)
    xmin=0; xmax=10; dx=1
    ymin=0; ymax=50; dy=10
    fnameF='fig_JH_zmax.png'
    return q0,cd,h0,tl,ft,hg,dds,ddp,shaft,port,tstr,xmin,xmax,dx,ymin,ymax,dy,fnameF    
    
    
def inp_JT():
    q0=338       # discharge of pressure tunnel
    cd=0.9      # discharge coefficient of port
    h0=6.121    # head loss of pressure tunnel
    tl=1749 # length of pressure tunnel including untake tunnel
    ft=52.810   # section area of pressure tunnel
    hg=677    # gross head
    dds=10    # design shaft diameter
    ddp=4.5     # design port diameter
    shaft=np.arange(8,13,1).astype(np.float64)      # shaft diameter: parameter
    port=np.arange(0.1,10.05,0.05).astype(np.float64) # port diameter: parameter
    tstr='PSPP-J Tailrace Surge Tank (h0={0:.3f})'.format(h0)
    xmin=0; xmax=10; dx=1
    ymin=0; ymax=90; dy=10
    fnameF='fig_JT_zmax.png'
    return q0,cd,h0,tl,ft,hg,dds,ddp,shaft,port,tstr,xmin,xmax,dx,ymin,ymax,dy,fnameF


def drawfig(shaft,port,zz,zd,dzz,ddp,xmin,xmax,dx,ymin,ymax,dy,tstr,fnameF):
    fsz=16
    xsize=8; ysize=6 # size of figure
    fig=plt.figure(figsize=(xsize,ysize),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Port diameter (m)')
    plt.ylabel('Maximum water level rise (m)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)
    for i,sh in enumerate(shaft):
        # Since zm is monotone increasing function for port diameter, 
        # non-converged values of zm are omitted from plot data
        n=np.argmin(zz[i,:])
        plt.plot(port[n:],zz[i,n:],color='#000000',lw=1)
        xs=xmax
        ys=zz[i,-1]
        ss='Ds={0:.0f}m'.format(sh)
        plt.text(xs,ys,ss,color='#000000',ha='right',va='bottom',fontsize=fsz-2)
    plt.plot([xmin],[ymin],'-',color='#000000',lw=1,label="Water level rise")
    plt.plot(port,zd,'--',color='#000000',lw=1,label="$|z_m'|=k_0-h_0$")
    plt.plot([ddp],[dzz],'o',color='#000000',ms=8,label='Design point')   
    plt.plot([xmin,ddp],[dzz,dzz],'-',color='#000000',lw=0.5)
    plt.plot([ddp,ddp],[ymin,dzz],'-',color='#000000',lw=0.5)
    ss1='Dp={0:.1f}m'.format(ddp)
    ss2='Zmax={0:.3f}m'.format(dzz)
    plt.text(ddp,ymin,ss1,color='#000000',ha='left',va='bottom',rotation=90,fontsize=fsz-2)
    plt.text(xmin,dzz,ss2,color='#000000',ha='left',va='bottom',rotation=0,fontsize=fsz-2)
#    plt.legend(loc='lower right',fontsize=fsz-2,shadow=True)
    plt.legend(loc='lower right',fontsize=fsz-2,shadow=True,numpoints=1,title='Ds: shaft diameter').get_title().set_fontsize(fsz-2)
    plt.tight_layout()
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def func(z,param):
    # h0 is global variable
    h0=param[0]
    md=param[1]
    k0=param[2]
    if md*k0<1: f=((1+md*z)-np.log(1+md*z))-((1+md*h0)-np.log(1-md*k0))
    if md*k0>1: f=((md*z-1)+np.log(md*z-1))-(np.log(md*k0-1)-(md*h0+1))
    return f


def bisection(x01,x02,param):
    imax=100    # maximum number of iterations
    err=1.0e-10 # allowable error of the zero value
    x1=x01
    x2=x02
    for i in range(0,imax):
        xi=0.5*(x1+x2)
        f1=func(x1,param)
        f2=func(x2,param)
        fi=func(xi,param)
        if abs(x2-x1) < err: break # converged
        if 0<f1*f2: break          # no solution
        if f1*fi < 0.0: x2=xi
        if fi*f2 < 0.0: x1=xi
    return xi

           
def cal_z(param,zs):
    h0=param[0]
    md=param[1]
    k0=param[2]
    eps1=1.0e-10
    zmax=0
    z2=zs # maximum limit of assumed solution
    if md*k0<1: z1=-1/md+eps1 # minimum limit of assumed solution
    if md*k0>1: z1= 1/md+eps1 # minimum limit of assumed solution
    zmax=bisection(z1,z2,param)
    return zmax
    

def main():
    for iii in [1,2]:
        if iii==1:
            q0,cd,h0,tl,ft,hg,dds,ddp,shaft,port,tstr,xmin,xmax,dx,ymin,ymax,dy,fnameF=inp_JH()
        if iii==2:
            q0,cd,h0,tl,ft,hg,dds,ddp,shaft,port,tstr,xmin,xmax,dx,ymin,ymax,dy,fnameF=inp_JT()
        
        v0=q0/ft    # flow velocity of pressure tunnel
        fs=np.pi*shaft**2/4 # shaft diameter
        fp=np.pi*port**2/4  # port diameter
        zz=np.zeros((len(shaft),len(port)),dtype=np.float64) # array for |zm|
        zd=np.zeros(len(port),dtype=np.float64)              # array for |zrm'|
        eps=0.0001
        dzz=0
        for i in range(len(fs)):
            for j in range(len(port)):
                k0=1/2/g*(q0/cd/fp[j])**2        # k0
                md=2*g*fs[i]*(h0+k0)/tl/ft/v0**2 # m'
                zs=v0*np.sqrt(tl*ft/g/fs[i])     # free surge water rise
                param=np.array([h0,md,k0])
                zz[i,j]=np.abs(cal_z(param,zs))     # calculation of |zm|
                # set design water level rize at design point
                #     dzz : design water level rise
                #     dds : design shaft diameter 
                #     ddp : design port diameter
                if dds-eps<shaft[i] and shaft[i]<dds+eps:
                    if ddp-eps<port[j] and port[j]<ddp+eps:
                        dzz=zz[i,j]
        zd=1/2/g*(q0/cd/fp)**2-h0 # |zrm'|=k0-h0
        drawfig(shaft,port,zz,zd,dzz,ddp,xmin,xmax,dx,ymin,ymax,dy,tstr,fnameF)

        # Free surge water rise
        zs=v0*np.sqrt(tl*ft/g/fs) # free surge water rise
        print('* Free surge water rise')
        for i in range(len(fs)):
            print('Ds={0:2.0f}    z*={1:6.3f}'.format(shaft[i],zs[i]))

        # required diameter
        fp=np.pi*ddp**2/4
        cc=h0/v0**2
        k0=1/2/g*(q0/cd/fp)**2
        f1d=h0*tl*ft/cc/(h0+k0)/g/hg;       dr1=np.sqrt(4*f1d/np.pi) # Thoma-Schuller's dynamic stability (1)
        f2d=tl*ft/2/cc/g/(hg-dzz);          dr2=np.sqrt(4*f2d/np.pi) # Thoma-Schuller's dynamic stability (2)
        f3d=h0*tl*ft**3/2/g/cc**2/k0/q0**2; dr3=np.sqrt(4*f3d/np.pi) # Calame-Garden's critical discharge
        hh=hg/6
        print('* Requirements for stability')
        print('Hg/6={0:6.3f}m  # for static stabity'.format(hh))
        print('Dsr1={0:4.1f}m    # for dynamic stability (1)'.format(dr1))
        print('Dsr2={0:4.1f}m    # for dynamic stability (2)'.format(dr2))
        print('Dsr3={0:4.1f}m    # for critical discharge'.format(dr3))


#==============
# Execution
#==============
if __name__ == '__main__': main()

Creation of Input data for surging analysis

The positive direction of discharge is defined to surgetank from reservoir. Therefore, followings shall be noted.

The discharge in the case of input intersection for headrace surge tank has negative sign.
The discharge in the case of load interception and load increase for tailrace surge tank have negative signs.

Input data format is shown below.

Variables	Description
strcom	comment (title of figure)
ICT, AFCA, AFCT	ICT=1(normal), AFCA=0 (AFC discharge), AFCT=1 (period of AFC)
TMAX, dt, DTWR	calculation time, time increment for calc, time increment for print
PAA, PCI, PCO	section area of port, C for port inflow, C for port outflow
RWL, TNL, TNA, TNC	Reservoir WL, tunnel length, tunnel section area, c of tunnel
NST	number of section
SAA[i], SEL[i], SLB[I] ......	section area, elevation, label
NQT	number of discharge input points
QTQ[i], QTI[I] ......	discharge at specified time, time
ymin, ymax	range of y-axis of drawing

The case of each created file is shown below.

Case	File name
Headrace surge tank: load interception	inp_surge_JH1.csv
Headrace surge tank: load increase	inp_surge_JH2.csv
Headrace surge tank: input interception	inp_surge_JH3.csv
Tailrace surge tank: load interception	inp_surge_JT1.csv
Tailrace surge tank: load increase	inp_surge_JT2.csv
Tailrace surge tank: input interception	inp_surge_JT3.csv

# strcom                   comment (title of figure)
# ICT, AFCA, AFCT          ICT=1(normal), AFCA=0 (AFC discharge), AFCT=1 (period of AFC)
# TMAX,dt,DTWR             calculation time, time increment for calc, time increment for print
# PAA, PCI, PCO            section area of port, C for port inflow, C for port outflow
# RWL, TNL, TNA, TNC       Reservoir WL, tunnel length, tunnel section area, c of tunnel
# NST                      number of section
# SAA[i], SEL[i], SLB[i]   section area, elevation, label
# ......
# NQT                      number of discharge input points
# QTQ[i], QTI[i]           discharge at specified time, time
# ......
# ymin, ymax               for drawing


# PSPP-J headrace ST
fnameW='inp_surge_JH1.csv'
f=open(fnameW,'w')
print('Load Interception (PSPP-J headrace ST)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('15.904,0.9,0.9',file=f)
print('1340,4800,52.810,0.301',file=f)
print('2',file=f)
print('346.313,1379.0,Top of Shaft (EL.1379.0)',file=f)
print('346.313,1275.0,Bottom of Shaft (EL.1275.0)',file=f)
print('3',file=f)
print('338,0',file=f)
print('0,8',file=f)
print('0,9999',file=f)
print('1270,1390',file=f)
f.close()

fnameW='inp_surge_JH2.csv'
f=open(fnameW,'w')
print('Load Increase (PSPP-J headrace ST)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('15.904,0.9,0.9',file=f)
print('1315,4800,52.810,0.448',file=f)
print('2',file=f)
print('346.313,1379.0,Top of Shaft (EL.1379.0)',file=f)
print('346.313,1275.0,Bottom of Shaft (EL.1275.0)',file=f)
print('3',file=f)
print('169,0',file=f)
print('338,40',file=f)
print('338,9999',file=f)
print('1270,1390',file=f)
f.close()

fnameW='inp_surge_JH3.csv'
f=open(fnameW,'w')
print('Input Interception (PSPP-J headrace ST)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('15.904,0.9,0.9',file=f)
print('1315,4800,52.810,0.301',file=f)
print('2',file=f)
print('346.313,1379.0,Top of Shaft (EL.1379.0)',file=f)
print('346.313,1275.0,Bottom of Shaft (EL.1275.0)',file=f)
print('3',file=f)
print('-236.6,0',file=f)
print('0,8',file=f)
print('0,9999',file=f)
print('1270,1390',file=f)
f.close()

# PSPP-J tailrace ST
fnameW='inp_surge_JT1.csv'
f=open(fnameW,'w')
print('Load Interception (PSPP-J headrace ST)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('15.904,0.9,0.9',file=f)
print('630,1749,52.810,0.149',file=f)
print('3',file=f)
print('600.000,688.0,Top of chamber wall (EL.688.0)',file=f)
print('600.000,680.0,Bottom of chamber (EL.680.0)',file=f)
print('78.540,556.3,Bottom of shaft (EL.556.3)',file=f)
print('3',file=f)
print('-338,0',file=f)
print('0,8',file=f)
print('0,9999',file=f)
print('550,700',file=f)
f.close()

fnameW='inp_surge_JT2.csv'
f=open(fnameW,'w')
print('Load Interception (PSPP-J headrace ST)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('15.904,0.9,0.9',file=f)
print('670,1749,52.810,0.207',file=f)
print('3',file=f)
print('600.000,688.0,Top of chamber wall (EL.688.0)',file=f)
print('600.000,680.0,Bottom of chamber (EL.680.0)',file=f)
print('78.540,556.3,Bottom of shaft (EL.556.3)',file=f)
print('3',file=f)
print('-169,0',file=f)
print('-338,40',file=f)
print('-338,9999',file=f)
print('550,700',file=f)
f.close()

fnameW='inp_surge_JT3.csv'
f=open(fnameW,'w')
print('Load Interception (PSPP-J headrace ST)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('15.904,0.9,0.9',file=f)
print('670,1749,52.810,0.149',file=f)
print('3',file=f)
print('600.000,688.0,Top of chamber wall (EL.688.0)',file=f)
print('600.000,680.0,Bottom of chamber (EL.680.0)',file=f)
print('78.540,556.3,Bottom of shaft (EL.556.3)',file=f)
print('3',file=f)
print('236.6,0',file=f)
print('0,8',file=f)
print('0,9999',file=f)
print('550,700',file=f)
f.close()

Surging analysis

#==========================
# Surging Analysis
#==========================
import numpy as np
import sys
import matplotlib.pyplot as plt
import time


def wout(dts,zw,vw,qn,fout):
    #print of Time, WL of Surge tank, Velocity of Tunnel, Discharge, k
    cf=0.5*(PCI+PCO)*PAA
    k=fk(vw,qn,cf)
    print('{0:15e},{1:15e},{2:15e},{3:15e},{4:15e}'.format(dts,zw,vw,qn,k),file=fout)


def fz(v,q,fa):
    return (q-TNA*v)/fa


def fk(v,q,cf):
    return np.abs((TNA*v-q)/cf)*(TNA*v-q)/cf/2.0/g


def fv(z,v,wrk):
    return (z-TNC*np.abs(v)*v-wrk)/TNL*g


def qncal(tt):
    qn=QI+AFCA*np.sin(2.0*np.pi/AFCT*tt)
    if tt>=AFCS:
        for i in range(1,NQT):
            if QTI[i-1]<tt and tt<=QTI[i]:
                qn=QTQ[i-1]+(QTQ[i]-QTQ[i-1])/(QTI[i]-QTI[i-1])*(tt-QTI[i-1])
                break
    return qn


def runge(qo,vo,zo,fa,dtt,dts):
    wrk=0.0
    dz0=dtt*fz(vo,qo,fa)
    cf=PCI*PAA              # Water level rising (in-flow from port)
    if dz0>0.0: cf=PCO*PAA  # Water level drop (out-flow from port)
    wrk=fk(vo,qo,cf)
    dv0=dtt*fv(zo,vo,wrk)

    v1=vo+0.5*dv0
    z1=zo+0.5*dz0
    dtw=dts-dtt*0.5
    q1=qncal(dtw)
    wrk=fk(v1,q1,cf)
    dv1=dtt*fv(z1,v1,wrk)
    dz1=dtt*fz(v1,q1,fa)

    v1=vo+0.5*dv1
    z1=zo+0.5*dz1
    wrk=fk(v1,q1,cf)
    dv2=dtt*fv(z1,v1,wrk)
    dz2=dtt*fz(v1,q1,fa)

    v1=vo+dv2
    z1=zo+dz2
    dtw=dts
    q2=qncal(dtw)
    wrk=fk(v1,q2,cf)
    dv3=dtt*fv(z1,v1,wrk)
    dz3=dtt*fz(v1,q2,fa)

    dv=(dv0+2.0*dv1+2.0*dv2+dv3)/6.0
    dz=(dz0+2.0*dz1+2.0*dz2+dz3)/6.0
    return dv,dz,q2


def datainp(fnameR):
    # data input
    global g
    global strcom
    global ICT,AFCA,AFCT
    global TMAX,DT,DTWR
    global PAA,PCI,PCO
    global RWL,TNL,TNA,TNC
    global NST,SEL,SAA,SLB
    global NQT,QTQ,QTI
    global AFCS, QI
    global ymin,ymax

    g=9.8

    fin=open(fnameR,'r')
    text=fin.readline()
    text=text.strip()
    strcom=text         # Comment

    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    ICT =int(text[0])   # Calculation case (1: normal, 2: AFC)
    AFCA=float(text[1]) # Discharge of half amplitude of AFC(m3/s)
    AFCT=float(text[2]) # Period of AFC(s)

    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    TMAX=float(text[0]) # End time of calculation
    DT  =float(text[1]) # Time increment for calculation
    DTWR=float(text[2]) # Time increment for print out

    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    PAA=float(text[0])  # Section area of Port
    PCI=float(text[1])  # Discharge coefficient of Port for inflow
    PCO=float(text[2])  # Discharge coefficient of Port for outflow

    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    RWL=float(text[0])  # Reservoir Water Level
    TNL=float(text[1])  # Length of pressure tunnel
    TNA=float(text[2])  # Section area of pressure tunnel
    TNC=float(text[3])  # Head loss coefficient of pressure tunnel

    SAA=[]
    SEL=[]
    SLB=[]
    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    NST=int(text[0])             # Number of sections
    for i in range(0,NST):
        text=fin.readline()
        text=text.strip()
        text=text.split(',')
        SAA=SAA+[float(text[0])] # Section area of specified level in the shaft
        SEL=SEL+[float(text[1])] # Level definition of sections of shaft (EL)
        SLB=SLB+[text[2]]        # Label for section (text: for drawing)

    QTQ=[]
    QTI=[]
    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    NQT=int(text[0])             # #Number of discharge input points
    for i in range(0,NQT):
        text=fin.readline()
        text=text.strip()
        text=text.split(',')
        QTQ=QTQ+[float(text[0])] # Discharge at specified time point
        QTI=QTI+[float(text[1])] # Time at specified time point

    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    ymin=float(text[0])
    ymax=float(text[1])
    fin.close()

    # Preprocessing
    QI=QTQ[0]
    if ICT==1: AFCS=-1.0
    if ICT==2: AFCS=QTI[1]
    _list = [[-1]*3 for i in range(0,NST)]
    for i in range(0,NST):
        _list[i][0]=SEL[i]
        _list[i][1]=SAA[i]
        _list[i][2]=SLB[i]
    _list.sort()         # Sort section numbers in small order of elevation
    for i in range(0,NST):
        SEL[i]=_list[i][0]
        SAA[i]=_list[i][1]
        SLB[i]=_list[i][2]
    del _list
    SAA=SAA+[SAA[NST-1]] # Dummy for judgement of top level (section area)
    SEL=SEL+[SEL[NST-1]] # Dummy for judgement of top level (EL)


def calc(fnameW):
    # Main calculation
    fout=open(fnameW,'w')
    print('Time,WL of Surge tank,Velocity of Tunnel,Discharge,k',file=fout)
    nnp=int(DTWR/DT+0.00001) # Step number of print out
    dto=DT                   # Time increment
    qn=QI                    # Initial discharge
    vi=qn/TNA                # Initial velocity of tunnel discharge
    zi=TNC*vi*vi             # Initial head loss due to initial velocity
    if vi<0.0: zi=-TNC*vi*vi # Definition of sign of initial head loss (direction of flow)
    zo=RWL-zi                # Initial water level in Surge Tank
    # Print out initial values (time=0)
    dts=0.0
    zw=zo
    vw=vi
    zi=zi
    qn=qn
    wout(dts,zw,vw,qn,fout)
    #To search the section number with initial water level
    for j in range(NST-1,-1,-1):
        if SEL[j]<zo: break
    npe=j      #Number of section with initial water level
    npw=0      #Initialization of counter for print out
    dts=0.0    #Initialization of elapsed time
    iret=0
    while dts<=TMAX:
        dts=dts+dto  #The time found solution
        npw=npw+1
        fa=SAA[npe]
        dv,dz,q1=runge(qn,vi,zi,fa,dto,dts)
        zw=RWL-zi-dz
        iw=npe+1
        if iw>NST-1: iw=NST-1
        if zw<SEL[npe] or SEL[iw]<zw:
            dti=0.1*dto
            for ic in range(1,11):
                dtw=dts-dto+dti*float(ic)
                zw=RWL-zi
                if SEL[npe+1]<zw: npe=npe+1
                if zw<SEL[npe]: npe=npe-1
                if NST-1<npe: iret=1 # water level greater than USWL
                if npe<0:     iret=2 # water level less than DSWL
                fa=SAA[npe]
                dv,dz,q2=runge(q1,vi,zi,fa,dti,dtw)
                vi=vi+dv
                zi=zi+dz
                q1=q2
        else:
            vi=vi+dv
            zi=zi+dz
            qn=q1
        # Print out results
        if npw>=nnp:
            npw=0
            zw=RWL-zi
            vw=vi
            qn=q1
            wout(dts,zw,vw,qn,fout)
    fout.close()
    return iret


def drawfig(fnameW,fnameF):
    data=np.loadtxt(fnameW,delimiter=',', skiprows=1, usecols=(0,1))
    ti=data[:,0]
    wl=data[:,1]
    xmin=0.0
    xmax=TMAX
    fsz=14
    fig=plt.figure(figsize=(10, 5),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlabel('Time (second)')
    plt.ylabel('Elevation (EL.m)')
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.grid(linestyle='--')
    plt.plot(ti,wl,color='black',lw=2.0,label='WL in the shaft')
    plt.hlines(RWL,xmin,xmax,alpha=0.3,color='black',lw=1.5,linestyles='-.')
    plt.text(xmax-10,RWL,'RWL (EL.{0:.1f})'.format(RWL),fontsize=fsz,color='black',ha='right',va='bottom')
    for i in range(0,NST):
        if SLB[i].find('xxx')<0:
            plt.hlines(SEL[i],xmin,xmax,alpha=0.3,color='black',lw=1.5,linestyles='-')
            plt.text(xmax-10,SEL[i],SLB[i],fontsize=fsz,color='black',ha='right',va='bottom')
    n1=np.argmax(wl)
    n2=np.argmin(wl)
    plt.text(ti[n1],wl[n1],'max:{0:.3f}({1:.1f}sec)'.format(wl[n1],ti[n1]),fontsize=fsz,color='black',ha='left',va='bottom')
    plt.text(ti[n2],wl[n2],'min:{0:.3f}({1:.1f}sec)'.format(wl[n2],ti[n2]),fontsize=fsz,color='black',ha='left',va='top')
    plt.title(strcom,loc='left',fontsize=fsz)
    plt.savefig(fnameF,dpi=200)
    plt.show()
    plt.close()


def main():
    lstr1=['JH1','JH2','JH3']
    lstr2=['JT1','JT2','JT3']
    lstr=lstr1+lstr2
    for str in lstr:
        start=time.time()
        fnameR='inp_surge_'+str+'.csv'
        fnameW='out_surge_'+str+'.csv'
        fnameF='_fig_surge_'+str+'.png'
        #param=sys.argv
        param=['surge', fnameR,fnameW,fnameF]
        fnameR=param[1] # input data file
        fnameW=param[2] # output data file
        fnameF=param[3] # output image data file

        datainp(fnameR)
        iret=calc(fnameW)
        time1=time.time()-start    
        drawfig(fnameW,fnameF)
        time2=time.time()-start-time1

        print('iret={0:d}'.format(iret))
        if iret==0: print('succsessful termination')
        if iret==1: print('water level is greater than Top of the shaft')
        if iret==2: print('water level is less than Bottom of the shaft')

        print(time1)
        print(time2)
    

#==============
# Execution
#==============
if __name__ == '__main__': main()

Marge images

# 画像結合
# https://note.nkmk.me/python-pillow-concat-images/

from PIL import Image
import os
import glob


def comb_h(im1, im2):
    dst = Image.new('RGB', (im1.width + im2.width, im1.height))
    dst.paste(im1, (0, 0))
    dst.paste(im2, (im1.width, 0))
    return dst

def comb_v(im1, im2):
    dst = Image.new('RGB', (im1.width, im1.height + im2.height))
    dst.paste(im1, (0, 0))
    dst.paste(im2, (0, im1.height))
    return dst

def multi_comb_h(im_list):
    _im=im_list.pop(0)
    for im in im_list:
        _im=comb_h(_im, im)
    return _im

def multi_comb_v(im_list):
    _im=im_list.pop(0)
    for im in im_list:
        _im=comb_v(_im, im)
    return _im


    
fig_01='_fig_surge_JH1.png'; im_01 = Image.open(fig_01)
fig_02='_fig_surge_JH2.png'; im_02 = Image.open(fig_02)
fig_03='_fig_surge_JH3.png'; im_03 = Image.open(fig_03)
fnameF='fig_surge_JH.png'
im_list_1=[im_01, im_02, im_03]
multi_comb_v(im_list_1).save(fnameF)
    

fig_01='_fig_surge_JT1.png'; im_01 = Image.open(fig_01)
fig_02='_fig_surge_JT2.png'; im_02 = Image.open(fig_02)
fig_03='_fig_surge_JT3.png'; im_03 = Image.open(fig_03)
fnameF='fig_surge_JT.png'
im_list_1=[im_01, im_02, im_03]
multi_comb_v(im_list_1).save(fnameF)
file_list = glob.glob('_*.png')


for file in file_list:
    print("remove：{0}".format(file))
    os.remove(file)

That's all. Thank you.