damyarou

python, GMT などのプログラム

設計 ダムの安定計算

記事の最後に行く

はじめに

図に示すような単純形状のダムの安定計算プログラムです。 このような計算は擁壁などの場合も含めよく行うのですが、その都度プログラムを作っていました。 でも簡単な事例をどこかに載せておけば何かと便利なのでここに事例としてストックしておきます。

f:id:damyarou:20190516105639j:plain:w600

計算式


\begin{equation}
e=\cfrac{M}{V}-\cfrac{b}{2} \qquad
p_u=\cfrac{V}{b}\cdot\left(1-\cfrac{6\cdot e}{b}\right) \qquad
p_d=\cfrac{V}{b}\cdot\left(1+\cfrac{6\cdot e}{b}\right)
\end{equation}



\begin{equation}
SF_s=\cfrac{c\cdot b+V\cdot \tan{\phi}}{H} \qquad
SF_f=\cfrac{W}{U}
\end{equation}


 e Eccentricity from dam base center
 p_u Vertical stress at upstream bottom edge of dam
 p_d Vertical stress at downstream bottom edge of dam
 V Total vertical force
 H Total horizontal force
 M Total moment around upstream bottom edge of dam
 b dam base length in U/D direction
 c Cohesion of foundation rock
 \phiInternal friction angle of foundation rock
 W Sum of vertical force
 U Sum of uplift force
 SF_sSafety factor for sliding
 SF_fsafety factor for floatation



\begin{equation}
q_a=\cfrac{q_u}{SF_b} \qquad q_u=\cfrac{2\cdot c\cdot \cos{\phi}}{1-\sin{\phi}}
\end{equation}


 q_a Allowable bearing capacity
 q_u Ultimate bearing capacity
 SF_bSafety factor (=1.5: BS 8110-1)

Dynamic Water Pressure by Westergaard


\begin{align}
&\text{Pressure} & &p=\cfrac{7}{8}\cdot \gamma_w\cdot k_h\cdot \sqrt{H\cdot h} \\
&\text{Force} & &P=\cfrac{7}{12}\cdot \gamma_w\cdot k_h\cdot H^2 \qquad y=0.4\cdot H
\end{align}


 pdynamic water pressure at depth  h
 Phorizontal force due to dynamic water pressure
 \gamma_wunit weight of water
 k_hseismic coefficient
 Hdepth from water surface to foundation
 hdepth from water surface to acting point of dynamic water pressure
 yarm length from foundation to centroid of dynamic pressure force

Sedimantation Pressure Force


\begin{equation}
P=\cfrac{1}{2}\cdot \gamma_s\cdot C_e\cdot h_s^2 \qquad y=\cfrac{1}{3}\cdot hs
\end{equation}


 Phorizontal force due to sedimentation pressure
 \gamma_sunit weight of submerged sedimentation (=12 kN/m3)
 C_esedimentation pressure coefficient (=0.5)
 h_ssedimentation depth from foundation to surface of sedimentation
 yarm length from foundation to centroid of sedimentation pressure force

安定計算プログラム

このプログラムで扱っている荷重は以下の通りです。堆砂圧は含まれていません。

  • ダム自重および地震時慣性力
  • 上流側静水圧および地震時動水圧
  • 下流側静水圧および地震時動水圧
  • 下流側水重
  • 揚圧力

地震時動水圧については、上流側・下流側とも同じ方向にかけています。 例えば、地震時慣性力が上流側から下流側に作用する場合、上流側動水圧は上流から下流方向、下流側動水圧も上流から下流方向、すなわち堤体を引っ張る方向にかけています。

出力は、画面出力としていますが、計算はJupyter上で行っているので、テキストをコピーしてそのまま報告書に貼り付けられます。

プログラム

import numpy as np


# global variants
elt=104.0  # dam top elevation
elb=70.5   # dam foundation elevation
n=0.65     # dpwnstream slope gradiant
tw=7.0     # top width of dam
bw=tw+(elt-elb)*n # bottom width of dam
wc=23.5    # unit weight of concrete (kN/m3)
ww=10.0    # unit weight of water (kN/m3)
cc=70.0*10 # cohesion of foundation rock (kN/m2)
phi=34.0   # internal friction angle of foundation rock (degrees)
qu=2*cc*np.cos(np.radians(phi))/(1-np.sin(np.radians(phi))) # ultimate bearing capacity
qa=qu/1.5  # allowable bearing capacity (SF=1.5)
print('* Basic Design Conditions')
print('Dam top level (EL)                            {0:10.3f} m'.format(elt))
print('Dam foundation level (EL)                     {0:10.3f} m'.format(elb))
print('Dam height                                    {0:10.3f} m'.format(elt-elb))
print('Dam upstream slope                                m={0:4.2f}'.format(0))
print('Dam downstream slope                              n={0:4.2f}'.format(n))
print('Dam top width                                 {0:10.3f} m'.format(tw))
print('Dam bottom width                              {0:10.3f} m'.format(bw))
print('Unit weight of concrete                       {0:10.3f} kN/m3'.format(wc))
print('Unit weight of water                          {0:10.3f} kN/m3'.format(ww))
print('Cohesion of foundation rock                   {0:10.3f} kN/m2'.format(cc))
print('Internal friction angle of foundation rock    {0:10.3f} degrees'.format(phi))
print('Allowable Bearing capacity of foundation rock {0:10.3f} kN/m2'.format(qa))
print()



def dambody():
    # dam self weight
    hh=elt-elb
    w1=wc*tw*hh        ; x1=tw/2         ; y1=hh/2
    w2=wc*(bw-tw)*hh/2 ; x2=tw+(bw-tw)/3 ; y2=hh/3
    wdb=w1+w2
    xx=(w1*x2+w2*x2)/wdb
    yy=(w1*y1+w2*y2)/wdb
    return wdb,xx,yy


def upwp(uwl):
    # upstream water pressure
    hh=uwl-elb
    uwp=ww*hh**2/2
    yy=hh/3
    return uwp,yy


def dnwp(dwl):
    # downstream water pressure
    hh=dwl-elb
    dwp=ww*hh**2/2
    yy=hh/3
    return dwp,yy


def dnww(dwl):
    # downstream water weight
    hh=dwl-elb
    dww=ww*n*hh**2/2
    xx=bw-n*hh/3
    return dww,xx


def updp(uwl,kh):
    # upstream dynamic water pressure
    hh=uwl-elb
    udp=ww*7/12*kh*hh**2
    yy=0.4*hh
    return udp,yy


def dndp(dwl,kh):
    # downstream dynamic water pressure
    hh=dwl-elb
    ddp=ww*7/12*kh*hh**2
    yy=0.4*hh
    return ddp,yy


def uplift(uwl,dwl):
    # uplift
    xd=8.0 # drain location
    cu=0.2 # uplift coefficient
    h1=uwl-elb
    h2=dwl-elb
    u1=ww*bw*h2               ; x1=bw/2
    u2=ww*xd*(h1-h2)*cu       ; x2=xd/2
    u3=ww*(bw-xd)*(h1-h2)*cu/2; x3=xd+(bw-xd)/3
    u4=ww*xd*(h1-h2)*(1-cu)/2 ; x4=xd/3
    upl=u1+u2+u3+u4
    xx=(u1*x1+u2*x2+u3*x3+u4*x4)/upl  
    return upl,xx


def calc(uwl,dwl,kh,icase):
    wdb,xdb,ydb=dambody()    # dam self weight
    uwp,yuwp=upwp(uwl)       # upstream water pressure
    dwp,ydwp=dnwp(dwl)       # downstream water pressure
    dww,xdww=dnww(dwl)       # downstream water weight
    udp,yudp=updp(uwl,kh)    # upstream dynamic water pressure
    ddp,yddp=dndp(dwl,kh)    # downstream dynamic water pressure
    upl,xupl=uplift(uwl,dwl) # uplift
    vv=np.zeros(7,dtype=np.float64)
    hh=np.zeros(7,dtype=np.float64)
    xx=np.zeros(7,dtype=np.float64)
    yy=np.zeros(7,dtype=np.float64)
    mm=np.zeros(7,dtype=np.float64)
    ss=['Dam self weight',
        'U/S water pressure',
        'D/S water pressure',
        'D/S water weight',
        'U/S dynamic pressure',
        'D/S dynamic pressure',
        'Uplift']
    i=0; vv[i]= wdb ; hh[i]= wdb*kh ; xx[i]=xdb  ; yy[i]=ydb   # Dam self weight
    i=1; vv[i]= 0   ; hh[i]= uwp    ; xx[i]=0    ; yy[i]=yuwp  # U/S water pressure
    i=2; vv[i]= 0   ; hh[i]=-dwp    ; xx[i]=0    ; yy[i]=ydwp  # D/S water pressure
    i=3; vv[i]= dww ; hh[i]=0       ; xx[i]=xdww ; yy[i]=0     # D/S water weight
    i=4; vv[i]= 0   ; hh[i]= udp    ; xx[i]=0    ; yy[i]=yudp  # U/S dynamic pressure
    i=5; vv[i]= 0   ; hh[i]= ddp    ; xx[i]=0    ; yy[i]=yddp  # D/S dynamic pressure
    i=6; vv[i]=-upl ; hh[i]=0       ; xx[i]=xupl ; yy[i]=0     # Uplift
    mm=vv*xx+hh*yy
    ee=np.sum(mm)/np.sum(vv)-bw/2
    pu=np.sum(vv)/bw*(1-6*ee/bw)
    pd=np.sum(vv)/bw*(1+6*ee/bw)
    fss=(cc*bw+np.sum(vv)*np.tan(np.radians(phi)))/np.sum(hh)
    fsf=(wdb+dww)/upl

    if icase==1: scase=' (usual operation)'
    if icase==2: scase=' (unusual operation)'
    if icase==3: scase=' (earthquake: SEE)'
    if icase==4: scase=' (earthquake: SEE)'
    if icase==5: scase=' (flood: PMF)'
    print('* Case-'+str(icase))
    s='-'*(21+12*5)
    print(s)
    print('Conditions: UWL=EL.{0:7.3f}m, DWL=EL.{1:7.3f}m, kh={2:6.3f}'.format(uwl,dwl,kh)+scase)
    print(s)
    print('{0:<21s}{1:>12s}{2:>12s}{3:>12s}{4:>12s}{5:>12s}'.format('Load Item','V (kN/m)','H (kN/m)','x (m)','y (m)','M (kN-m/m)'))
    print(s)
    for i in range(len(ss)):
        print('{0:<21s}{1:12.3f}{2:12.3f}{3:12.3f}{4:12.3f}{5:12.3f}'.format(ss[i],vv[i],hh[i],xx[i],yy[i],mm[i]))
    print('{0:<21s}{1:12.3f}{2:12.3f}{3:>12s}{4:>12s}{5:12.3f}'.format('Sum',np.sum(vv),np.sum(hh),'---','---',np.sum(mm)))
    print(s)
    print('{0:<21s}{1:>12s}{2:>12s}{3:>12s}{4:>12s}{5:>12s}'.format('Summary of Result','e (m)','pu (kN/m2)','pd (kN/m2)','FSs','FSf'))
    # Limit of acting point of resultant
    if icase==1: elim=bw/6 # middle third
    if icase==2: elim=bw/4 # middle half
    if icase==3: elim=bw/2 # within base
    if icase==4: elim=bw/2 # within base
    if icase==5: elim=bw/2 # within base
    se='e_limit={0:.3f}'.format(elim)
    print('{0:<21s}{1:12.3f}{2:12.3f}{3:12.3f}{4:12.3f}{5:12.3f}'.format(se,ee,pu,pd,fss,fsf))
    print(s)
    print()
    
    
def main():
    kh0=0.367                 # horizontal seismic coefficient
    upmf=101.35; dpmf=100.25  # PMF water level
    unwl= 98.00; dnwl= 91.37  # Normal water level
    ulwl= 95.00; dlwl= 87.40  # Dead water level (LWL)
    uwl=unwl; dwl=dnwl; kh= 0.000; calc(uwl,dwl,kh,1)
    uwl=unwl; dlw=dlwl; kh= 0.000; calc(uwl,dwl,kh,2)
    uwl=unwl; dwl=dnwl; kh= kh0  ; calc(uwl,dwl,kh,3)
    uwl=unwl; dwl=dnwl; kh=-kh0  ; calc(uwl,dwl,kh,4)
    uwl=upmf; dwl=dpmf; kh= 0.000; calc(uwl,dwl,kh,5)

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

画面出力

* Basic Design Conditions
Dam top level (EL)                               104.000 m
Dam foundation level (EL)                         70.500 m
Dam height                                        33.500 m
Dam upstream slope                                m=0.00
Dam downstream slope                              n=0.65
Dam top width                                      7.000 m
Dam bottom width                                  28.775 m
Unit weight of concrete                           23.500 kN/m3
Unit weight of water                              10.000 kN/m3
Cohesion of foundation rock                      700.000 kN/m2
Internal friction angle of foundation rock        34.000 degrees
Allowable Bearing capacity of foundation rock   1755.345 kN/m2

* Case-1
---------------------------------------------------------------------------------
Conditions: UWL=EL. 98.000m, DWL=EL. 91.370m, kh= 0.000 (usual operation)
---------------------------------------------------------------------------------
Load Item                V (kN/m)    H (kN/m)       x (m)       y (m)  M (kN-m/m)
---------------------------------------------------------------------------------
Dam self weight         14081.934       0.000      14.258      13.352  200784.914
U/S water pressure          0.000    3781.250       0.000       9.167   34661.458
D/S water pressure          0.000   -2177.785       0.000       6.957  -15150.121
D/S water weight         1415.560       0.000      24.253       0.000   34331.811
U/S dynamic pressure        0.000       0.000       0.000      11.000       0.000
D/S dynamic pressure        0.000       0.000       0.000       8.348       0.000
Uplift                  -6461.321       0.000      13.844       0.000  -89447.689
Sum                      9036.174    1603.465         ---         ---  165180.374
---------------------------------------------------------------------------------
Summary of Result           e (m)  pu (kN/m2)  pd (kN/m2)         FSs         FSf
e_limit=4.796               3.892      59.156     568.901      16.363       2.399
---------------------------------------------------------------------------------

* Case-2
---------------------------------------------------------------------------------
Conditions: UWL=EL. 98.000m, DWL=EL. 91.370m, kh= 0.000 (unusual operation)
---------------------------------------------------------------------------------
Load Item                V (kN/m)    H (kN/m)       x (m)       y (m)  M (kN-m/m)
---------------------------------------------------------------------------------
Dam self weight         14081.934       0.000      14.258      13.352  200784.914
U/S water pressure          0.000    3781.250       0.000       9.167   34661.458
D/S water pressure          0.000   -2177.785       0.000       6.957  -15150.121
D/S water weight         1415.560       0.000      24.253       0.000   34331.811
U/S dynamic pressure        0.000       0.000       0.000      11.000       0.000
D/S dynamic pressure        0.000       0.000       0.000       8.348       0.000
Uplift                  -6461.321       0.000      13.844       0.000  -89447.689
Sum                      9036.174    1603.465         ---         ---  165180.374
---------------------------------------------------------------------------------
Summary of Result           e (m)  pu (kN/m2)  pd (kN/m2)         FSs         FSf
e_limit=7.194               3.892      59.156     568.901      16.363       2.399
---------------------------------------------------------------------------------

* Case-3
---------------------------------------------------------------------------------
Conditions: UWL=EL. 98.000m, DWL=EL. 91.370m, kh= 0.367 (earthquake: SEE)
---------------------------------------------------------------------------------
Load Item                V (kN/m)    H (kN/m)       x (m)       y (m)  M (kN-m/m)
---------------------------------------------------------------------------------
Dam self weight         14081.934    5168.070      14.258      13.352  269787.014
U/S water pressure          0.000    3781.250       0.000       9.167   34661.458
D/S water pressure          0.000   -2177.785       0.000       6.957  -15150.121
D/S water weight         1415.560       0.000      24.253       0.000   34331.811
U/S dynamic pressure        0.000    1619.005       0.000      11.000   17809.057
D/S dynamic pressure        0.000     932.455       0.000       8.348    7784.132
Uplift                  -6461.321       0.000      13.844       0.000  -89447.689
Sum                      9036.174    9322.995         ---         ---  259775.663
---------------------------------------------------------------------------------
Summary of Result           e (m)  pu (kN/m2)  pd (kN/m2)         FSs         FSf
e_limit=14.388             14.361    -626.316    1254.374       2.814       2.399
---------------------------------------------------------------------------------

* Case-4
---------------------------------------------------------------------------------
Conditions: UWL=EL. 98.000m, DWL=EL. 91.370m, kh=-0.367 (earthquake: SEE)
---------------------------------------------------------------------------------
Load Item                V (kN/m)    H (kN/m)       x (m)       y (m)  M (kN-m/m)
---------------------------------------------------------------------------------
Dam self weight         14081.934   -5168.070      14.258      13.352  131782.814
U/S water pressure          0.000    3781.250       0.000       9.167   34661.458
D/S water pressure          0.000   -2177.785       0.000       6.957  -15150.121
D/S water weight         1415.560       0.000      24.253       0.000   34331.811
U/S dynamic pressure        0.000   -1619.005       0.000      11.000  -17809.057
D/S dynamic pressure        0.000    -932.455       0.000       8.348   -7784.132
Uplift                  -6461.321       0.000      13.844       0.000  -89447.689
Sum                      9036.174   -6116.064         ---         ---   70585.085
---------------------------------------------------------------------------------
Summary of Result           e (m)  pu (kN/m2)  pd (kN/m2)         FSs         FSf
e_limit=14.388             -6.576     744.629    -116.572      -4.290       2.399
---------------------------------------------------------------------------------

* Case-5
---------------------------------------------------------------------------------
Conditions: UWL=EL.101.350m, DWL=EL.100.250m, kh= 0.000 (flood: PMF)
---------------------------------------------------------------------------------
Load Item                V (kN/m)    H (kN/m)       x (m)       y (m)  M (kN-m/m)
---------------------------------------------------------------------------------
Dam self weight         14081.934       0.000      14.258      13.352  200784.914
U/S water pressure          0.000    4758.612       0.000      10.283   48934.399
D/S water pressure          0.000   -4425.312       0.000       9.917  -43884.349
D/S water weight         2876.453       0.000      22.329       0.000   64228.801
U/S dynamic pressure        0.000       0.000       0.000      12.340       0.000
D/S dynamic pressure        0.000       0.000       0.000      11.900       0.000
Uplift                  -8636.215       0.000      14.320       0.000 -123670.433
Sum                      8322.173     333.300         ---         ---  146393.332
---------------------------------------------------------------------------------
Summary of Result           e (m)  pu (kN/m2)  pd (kN/m2)         FSs         FSf
e_limit=14.388              3.203      96.041     482.390      77.275       1.964
---------------------------------------------------------------------------------

作図プログラム

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as tick


elt=104.0         # dam top elevation
elb=70.5          # dam foundation elevation
n=0.65            # dpwnstream slope gradiant
tw=7.0            # top width of dam
bw=tw+(elt-elb)*n # bottom width of dam
upmf=101.35; dpmf=100.25  # PMF water level
unwl= 98.00; dnwl= 91.37  # Normal water level
ulwl= 95.00; dlwl= 87.40  # Dead water level (LWL)


def barrow(x1,y1,x2,y2):
    col='#333333'
    arst='<->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)


def main():
    fsz=12
    xmin=-20; xmax=40; dx=5
    ymin=65; ymax=110; dy=5
    plt.figure(figsize=(8,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Basic Shape of Concrete Gravity Dam (unit: m)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Distance (m)')
    plt.ylabel('Elevation (EL.m)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.gca().spines['right'].set_visible(False)
    plt.gca().spines['top'].set_visible(False)
#    plt.gca().yaxis.set_ticks_position('left')
#    plt.gca().xaxis.set_ticks_position('bottom')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz+2)

    damx=np.array([0,0,tw,bw])
    damy=np.array([elb,elt,elt,elb])
    plt.fill(damx,damy,facecolor='#eeeeee',edgecolor='#000000',lw=2) # dam body
    plt.fill([-5,bw+5,bw+5,-5],[elb,elb,elb-2,elb-2],fill=False, hatch='//',lw=0) # hatching of foundation
    plt.plot([-5,bw+5],[elb,elb],'-',lw=2) # foundation line
    
    plt.plot([xmin+5,0],[elt,elt],'-',lw=1)
    plt.plot([xmin+5,0],[elb,elb],'-',lw=1)
    plt.plot([xmin+5,0],[upmf,upmf],'-',lw=1)
    plt.plot([xmin+5,0],[unwl,unwl],'-',lw=1)
    plt.plot([xmin+5,0],[ulwl,ulwl],'-',lw=1)
    plt.text(xmin+5.0,elt,'EL.{0:.3f}'.format(elt),va='bottom',ha='left',fontsize=fsz)
    plt.text(xmin+5.0,elb,'EL.{0:.3f}'.format(elb),va='bottom',ha='left',fontsize=fsz)
    plt.text(xmin+5.0,upmf,'PMF: EL.{0:.3f}'.format(upmf),va='bottom',ha='left',fontsize=fsz)
    plt.text(xmin+5.0,unwl,'NWL: EL.{0:.3f}'.format(unwl),va='bottom',ha='left',fontsize=fsz)
    plt.text(xmin+5.0,ulwl,'LWL: EL.{0:.3f}'.format(ulwl),va='bottom',ha='left',fontsize=fsz)

    plt.plot([xmax-5,tw+(elt-dpmf)*n],[dpmf,dpmf],'-',lw=1)
    plt.plot([xmax-5,tw+(elt-dnwl)*n],[dnwl,dnwl],'-',lw=1)
    plt.plot([xmax-5,tw+(elt-dlwl)*n],[dlwl,dlwl],'-',lw=1)
    plt.text(xmax-5.0,dpmf,'PMF: EL.{0:.3f}'.format(dpmf),va='bottom',ha='right',fontsize=fsz)
    plt.text(xmax-5.0,dnwl,'NWL: EL.{0:.3f}'.format(dnwl),va='bottom',ha='right',fontsize=fsz)
    plt.text(xmax-5.0,dlwl,'LWL: EL.{0:.3f}'.format(dlwl),va='bottom',ha='right',fontsize=fsz)

    plt.plot([0,0],[elt,elt+3],'-',lw=1,color='#333333')
    plt.plot([tw,tw],[elt,elt+3],'-',lw=1,color='#333333')
    x1=0; x2=tw; y1=elt+2.5; y2=y1; barrow(x1,y1,x2,y2)
    xs=0.5*(x1+x2); ys=y1
    plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz)

    plt.plot([0,0],[elb,elb-4],'-',lw=1,color='#333333')
    plt.plot([bw,bw],[elb,elb-4],'-',lw=1,color='#333333')
    x1=0; x2=bw; y1=elb-3.5; y2=y1; barrow(x1,y1,x2,y2)
    xs=0.5*(x1+x2); ys=y1
    plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='top',ha='center',fontsize=fsz)

    x1=2; x2=x1; y1=elb; y2=elt; barrow(x1,y1,x2,y2)
    xs=x1; ys=0.5*(y1+y2)
    plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='left',fontsize=fsz,rotation=90)
    
    ang=np.degrees(np.arctan(1/n))
    ys=80.0; xs=bw-(ys-elb)*n+1
    plt.text(xs,ys,'1:{0:.2f}'.format(n),va='center',ha='center',fontsize=fsz,rotation=-ang)
    
    fnameF='fig_dam.jpg'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)    
    plt.show()    
        
#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く