はじめに

ここで紹介しているのは、制水口式調圧水槽のサージングの基礎微分方程式（連立微分方程式）を Runge-Kutta 法により解き、水面振動の時刻歴を求めるプログラムです。以下のような出力が得られます。

f:id:damyarou:20190518111428j:plain

f:id:damyarou:20190518111441j:plain

計算理論

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}$

$v$	Flow velocity
$\rho$	Density of the water
$x$	Position
$g$	Acceleration of gravity
$t$	Time
$p$	Pressure
$h$	Water head

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}$

$z$	Variation of the water level in surge tank (positive is the downward as standard at reservoir water level)
$L$	Length of the pressure tunnel
$v$	Flow velocity of the pressure tunnel

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}$

$Q$	Discharge of the turbine
$f$	Area of the pressure tunnel
$F$	Area of the surge tank shaft

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

Following simultaneous differential equations can be solved by Runge-Kutta method.

$\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}$

$v_p$	Flow velocity of the port
$F_p$	Area of the port
$C_d$	Discharge coefficient of the port
$v$	Flow velocity of the pressure tunnel (positive is from reservoir to surge tank)
$z$	Variation of the water level in surge tank
$g$	Acceleration of gravity
$c$	Head loss coefficient
$L$	Length of the pressure tunnel (from reservoir to the port)
$f$	Area of the pressure tunnel
$F$	Area of the shaft
$Q$	Discharge (in the interception current of the time)

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{equation} \begin{cases} 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} \qquad \text{(resistance of port)}\\ \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} \qquad \text{(increment of pressure tunnel flow velocity)}\\ \Delta z=\Delta t\cdot f_z(q,v) & & f_z(q,v)=\cfrac{q-f\cdot v}{F} \qquad \text{(increment of surge tank water level)} \end{cases} \end{equation}$

After this, the calculation steps by Runge-Kutta method are shown.

Step-1	$\Delta z_0=\Delta t\cdot f_z(Q_0, V_0)$ (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)$
Step-2	$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)$
Step-3	$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)$
Step-4	$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)$
Step-5	$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)$

(Note)

Above equation and steps are only the basic parts of numerical solution. Please understand the method to treat the chamber room of the surge tank is included in the program code shiwn in this page.

入力データの準備

ここでは、導水路調圧水槽のみを有する一般水力の調圧水槽をモデルに入力データの準備を行います。導水路はコンクリートライニングされているものとします。

損失水頭計算は、通常、発電電力量を算出するために行われており、流量や粗度の条件が必ずしもサージング計算用の入力条件と一致しているとは限りません。そこで、サージング計算用の損失水頭を簡易に求めるため、「損失水頭は、流量の２乗および粗度係数の２乗に比例する」として、損失水頭の補正を行います。

ポート（制水口）の流量係数は、ポートの詳細形状により定まり、その形状によってはサージタンクへの流入時・流出時で値が異なります。しかしながら、設計上は、流入時・流出時共に、大きめの流量係数を入力しておくことにより、安全側、すなわち大きな水位上昇と水位低下が得られるため、ここでは $C_{d,in}=C_{d,out}=0.9$ を用いています。

圧力水路の損失係数 $c$ については、全負荷遮断、負荷急増ともに、最大流量に対する圧力水路損失水頭に対して計算するようにしています。これは、

損失水頭は（ほとんど）流量比の２乗に比例するため、いずれの流量を用いても同一の $c$ が得られること
通常、最大出力算定のため、最大流量に対する損失水頭が精度良く計算されていること

を考慮したものです。

事前準備

Case	Load interception	Load increase
* Conditions for Loss calculation
Discharge	49.0 m3/s	49.0 m3/s
Reservoir water Level	EL.282.5 (FSL)	EL.275.0 (MOL)
Roughness coefficient	n=0.014	n=0.014
Head loss of pressure tunnel	3.409 m	3.801 m
* Conditions for Surging Analysis
Discharge	49.0 m3/s	49.0 m3/s
Reservoir water Level	EL.282.5 (FSL)	EL.275.0 (MOL)
Modified roughness coefficient	n=0.0125 (0.014-0.0015)	n=0.0155 (0.014+0.0015)
Modified head loss of pressure tunnel $h_l$ *	2.718 m	4.659 m

$\begin{align} &\text{For load interception}& &h_l=3.409\times\left(\cfrac{49.0}{49.0}\right)^2\times\left(\cfrac{0.0125}{0.014}\right)^2=2.718 m \\ &\text{For load increase} & &h_l=3.801\times\left(\cfrac{49.0}{49.0}\right)^2\times\left(\cfrac{0.0155}{0.014}\right)^2=4.659 m \end{align}$

入力データの整理

Case	Load interception	Load increase
Gross head $H_g$	192.5 m	192.5 m
Discharge $Q$	49.0 => 0 m3/s	24.5 => 49.0 m3/s
Closing/Opening time $T$	8 sec	8 sec
Reservoir water level $RWL$	EL.282.5 (FSL)	EL.275.0 (MOL)
Length of pressure tunnel $L$	6576 m	6576 m
Section area of pressure tunnel $f$	26.162 m2	26.162 m2
Head loss of pressure tunnel $h_l$	2.718 m	4.659 m
Velocity head at bottom of surge tank $h_b$	0.179 m (v=1.873m/s)	0.179 m (v=1.873m/s)
Total head loss $h_0=h_l+h_b$	2.897 m	4.838 m
Average velocity of pressure tunnel $v_0$	1.873 m/s	1.873 m/s
Loss coefficient of pressure tunnel $c=h_0/v_0{}^2$	0.826	1.379
Section area of shaft $F$	63.585 m2 (Ds=9.0m)	63.585 m2 (Ds=9.0m)
Section area of port $F_p$	3.462 m2 (Dp=2.1m)	3.462 m2 (Dp=2.1m)
Discharge coefficient of port $C_d$	0.9	0.9

(Note-1) Loss coefficient of pressure tunnel

The recommendable condition of the calculation of the loss coefficient of $c$ is shown below.

A loss coefficient of $c_1$ can be obtain as follow equation using given discharge of $Q_1$ , pressure tunnel head loss of $h_{L1}$ and pressure tunnel section area of $A$ .

$\begin{equation} c_1=\cfrac{h_{L1}+\cfrac{Q_1{}^2}{2g A^2}}{\cfrac{Q_1{}^2}{A^2}}=h_{L1}\cdot \cfrac{A^2}{Q_1{}^2}+\cfrac{1}{2g} \end{equation}$

For another discharge of $Q_2$ , pressure head loss of $h_{L2}$ can be calculated as following.

$\begin{equation} h_{L2}=h_{L1}\cdot \left(\cfrac{Q_2}{Q_1}\right)^2 \end{equation}$

Using above, a loss coefficient of $c_2$ at discharge of $Q_2$ can be calculated.

$\begin{equation} c_2=\cfrac{h_{L2}+\cfrac{Q_2{}^2}{2g A^2}}{\cfrac{Q_2{}^2}{A^2}} =\cfrac{h_{L1}\left(\cfrac{Q_2}{Q_1}\right)^2+\cfrac{Q_2{}^2}{2g A^2}}{\cfrac{Q_2{}^2}{A^2}} =h_{L1}\cdot \cfrac{A^2}{Q_1{}^2}+\cfrac{1}{2g} \end{equation}$

From above, it can be understood $c_1=c_2$ .

Therefore, considering a matter that the pressure tunnel head loss is calculated for the maximum discharge in order to calculate the installed capacity in many cases, it is recommendable to apply the maximum discharge for the calculation of loss coefficient.

(Note-2) Discharge coefficient

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.

プログラム

プログラム内の入力データ用変数

strcom : comment (title of figure)
ICT    :  Calculation case (1: normal, 2: AFC)
AFCA   :  Discharge of half amplitude of AFC(m3/s)
AFCT   :  Period of AFC(s)
TMAX   :  End time of calculation
DT     :  Time increment for calculation
DTWR   :  Time increment for print out
PAA    :  Section area of Port
PCI    :  Discharge coefficient of Port for inflow
PCO    :  Discharge coefficient of Port for outflow
RWL    :  Reservoir Water Level
TNL    :  Length of pressure tunnel
TNA    :  Section area of pressure tunnel
TNC    :  Head loss coefficient of pressure tunnel
NST    :  Number of sections
SAA[]  :  Section area of specified level in the shaft
SEL[]  :  Level definition of sections of shaft (EL)
SLB[]  :  Label for section (text: for drawing)
NQT    :  Number of discharge input points
QTQ[]  :  Discharge at specified time
QTI[]  :  Specified time

入力データ・フォーマット

strcom
ICT, AFCA, AFCT
TMAX,dt,DTWR
PAA, PCI, PCO
RWL, TNL, TNA, TNC
NST
SAA[i], SEL[i], SLB[i]
......
NQT
QTQ[i], QTI[i]
......

入力データファイル作成

Mac なら bash のヒアドキュメントで作成するのですが、Windows で動かすため、Python の print 文でデータファイルを作成しています。

fnameW='inp_surge_11.csv'
f=open(fnameW,'w')
print('Load Interception (Cd_i=0.9,Cd_o=0.9)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('3.462,0.9,0.9',file=f)
print('282.5,6576,26.162,0.826',file=f)
print('2',file=f)
print('63.585,307.0,Top of Shaft (EL.307.0)',file=f)
print('63.585,257.0,Bottom of Shaft (EL.257.0)',file=f)
print('3',file=f)
print('48,0',file=f)
print('0,8',file=f)
print('0,9999',file=f)
f.close()

fnameW='inp_surge_12.csv'
f=open(fnameW,'w')
print('Load Increase (Cd_i=0.9,Cd_o=0.9)',file=f)
print('1,0,1',file=f)
print('600,0.01,0.1',file=f)
print('3.462,0.9,0.9',file=f)
print('275.0,6576,26.162,1.379',file=f)
print('2',file=f)
print('63.585,307.0,Top of Shaft (EL.307.0)',file=f)
print('63.585,257.0,Bottom of Shaft (EL.257.0)',file=f)
print('3',file=f)
print('24.5,0',file=f)
print('49.0,8',file=f)
print('49.0,9999',file=f)
f.close()

計算プログラム

#==========================
# 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

    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
    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
    ymin=250
    ymax=320
    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():
    lstr=['11','12']
    for str in lstr:
        start=time.time()
        fnameR='inp_surge_'+str+'.csv'
        fnameW='out_surge_'+str+'.csv'
        fnameF='fig_surge_'+str+'.jpg'
        #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()

Thank you.

記事の先頭に行く

damyarou

python, GMT などのプログラム

設計制水口式調圧水槽のサージング解析