In [3]:

```
import math
from scipy.integrate import quad
#Given
H = 100; #heigth of dam
wb = 70; #width of base of dam
wt = 7; #width of top of dam
l = 1; #length of dam
hw = 98; #heigth of water in dam
hsu = 90; #heigth of slope on downstream side
s = 1/0.7; #slope on downstream side
gammad = 24; #unit weigth of dam
gammaw = 9.81; #unit weigth of water
E = 2.05e7; #modulus of elasticity
#(a) inertial forces and moments
alpha0 = 0.05; #from table 8.1
alphah = 2*alpha0;
#at 10m from top
def f0(y):
return 25.2-0.25*y
F10 = quad(f0,0,10)[0]
def f1(y):
return 25.2*(1-0.01*y)*(10-y)
M10 = quad(f1,0,10)[0]
#at 100m below top
def f2(y):
return 0.15*(1-0.01*y)*16.8*y
F100 = F10+ quad(f2,10,100)[0]
def f3(y):
return 0.15*(1-0.01*y)*16.8*y*(100-y)
M100 = M10+90*F10+ quad(f3,10,100)[0]
print "Inertial forces:At 10m from top: F = %.2f kn;M = %ikn-mAt 100m from top: F = %.2f kn;M = %ikn-m."%(F10,M10,F100,M100);
#(b) hydrodynamic pressure and moment
#at 10m from top
y = 8.;
W10 = 1680.;
alphah = F10/W10;
Cm = 0.735;
Cy = (Cm/2)*((y*(2-y/hw)/hw)+(y*(2-y/hw)/hw)**0.5);
p = Cy*alphah*gammaw*hw;
P10 = 0.726*p*y;
Mp10 = 0.299*p*y**2;
P10 = round(P10*100)/100;
Mp10 = round(Mp10*100)/100;
#at 100m from top
y = 98;
W100 = 84840;
alphah = F100/W100;
Cm = 0.735;
Cy = (Cm/2)*(y*(2-y/hw)/hw+(y*(2-y/hw)/hw)**0.5);
p = Cy*alphah*gammaw*hw;
P100 = 0.726*p*y;
Mp100 = 0.299*p*y**2;
print "Hydrodynamic forces:At 10m from top: F = %.2f kn;\
\nM = %.2fkn-mAt 100m from top: F = %i kn;\
\nM = %ikn-m."%(P10,Mp10,P100,Mp100);
# rounding off error.
```

In [5]:

```
import math
#Given
H = 100; #heigth of dam
wb = 70; #width of base of dam
wt = 7; #width of top of dam
l = 1; #length of dam
hw = 98; #heigth of water in dam
hsu = 90; #heigth of slope on downstream side
s = 1/0.7; #slope on downstream side
gammad = 24; #unit weigth of dam
gammaw = 9.81; #unit weigth of water
E = 2.05e7; #modulus of elasticity
beta = 1;
I = 2;
Fo = 0.25; #from table 8.2
#t = Sa/g;
t = 0.19; #from fig. 8.4
alphah = beta*I*Fo*t;
T = 5.55*H**2/wb*(gammad/(gammaw*E))**0.5;
#(a) Base shear
W = l*gammad*(wt*H+((hsu/s)*hsu)/2);
Fb = 0.6*W*alphah;
print "Base shear = %.2f KN."%(Fb);
#(b) Base moment
hbar = ((wt*H**2/2)+((hsu/s)*hsu**2/6))/((wt*H)+(hsu/s)*hsu/2);
Mb = 0.9*W*hbar*alphah;
print "Base moment = %.2f KN-m."%(Mb);
#(c) shear at 10m from top
Cv = 0.08;
F10 = Cv*Fb;
F10 = round(F10);
print "shear at 10m from top = %.2f KN."%(F10);
#(d) Moment at 10m from top
Cm = 0.02;
M10 = Cm*Mb;
M10 = round(M10);
print "moment at 10m from top = %.2f KN."%(M10);
#(e) Hydrodynamic pressure
#at 10m from top
y = 8;
W10 = 1680;
Cm = 0.735;
Cy = (Cm/2)*((y*(2-y/hw)/hw)+(y*(2-y/hw)/hw)**0.5);
p = Cy*alphah*gammaw*hw;
P10 = 0.726*p*y;
Mp10 = 0.299*p*y**2;
P10 = round(P10*100)/100;
Mp10 = round(Mp10*100)/100;
#at 100m from top
y = 98;
W100 = 84840;
Cm = 0.735;
Cy = (Cm/2)*(y*(2-y/hw)/hw+(y*(2-y/hw)/hw)**0.5);
p = Cy*alphah*gammaw*hw;
P100 = 0.726*p*y;
Mp100 = 0.299*p*y**2;
print "Hydrodynamic forces:At 10m from top: F = %.2f kn;\
\nM = %.2fkn-mAt 100m from top: F = %i \
kn;M = %ikn-m."%(P10,Mp10,P100,Mp100);
```

In [6]:

```
import math
from scipy.integrate import quad
#Given
H = 100; #heigth of dam
wb = 70; #width of base of dam
wt = 7; #width of top of dam
l = 1; #length of dam
hw = 98; #heigth of water in dam
hsu = 90; #heigth of slope on downstream side
s = 1/0.7; #slope on downstream side
gammad = 24; #unit weigth of dam
gammaw = 9.81; #unit weigth of water
E = 2.05e7; #modulus of elasticity
#(a) Seismic coefficient method
alpha0 = 0.05; #from table 8.1
alphah = 2*alpha0;
alphav = 0.75*alphah;
#at 10m from top
def f4(y):
return alphav*168*(1-0.01*y)
F10 = quad(f4,0,10)[0]
#at 100m below top
def f5(y):
return alphav*(1-0.01*y)*16.8*y
F100 = F10+ quad(f5,10,100)[0]
print "Parta):At 10m from top: F = %.2f knAt 100m from top: F = %.2f kn."%(F10,F100);
#(b)Response spectrum method
beta = 1;
I = 2;
Fo = 0.25; #from table 8.2
#t = Sa/g;
t = 0.19; #from fig. 8.4
alphah = beta*I*Fo*t;
alphav = 0.75*alphah;
#at 10m from top
def f6(y):
return alphav*168*(1-0.01*y)
F10 = quad(f6,0,10)[0]
#at 100m below top
def f7(y):
return alphav*(1-0.01*y)*16.8*y
F100 = F10+ quad(f7,10,100)[0]
F100 = round(F100*100)/100;
print "Partb):At 10m from top: F = %.2f knAt 100m from top: F = %.2f kn."%(F10,F100);
```

In [8]:

```
import math
#Given
H = 100.; #heigth of dam
wb = 73.; #width of base of dam
wt = 7; #width of top of dam
l = 1.; #length of dam
hw = 98.; #heigth of water in dam
hsu = 90.; #heigth of slope on downstream side
s = 1/0.7; #slope on downstream side
gammad = 24.; #unit weigth of dam
gammaw = 9.81; #unit weigth of water
E = 2.05e7; #modulus of elasticity
#at 10m from top
y = 8;
alphah = 0.1;
Cm = 0.72;
Cy = (Cm/2)*((y*(2-y/hw)/hw)+(y*(2-y/hw)/hw)**0.5);
p10 = Cy*alphah*gammaw*hw;
F10 = 0.726*p10*y;
Mp10 = 0.299*p10*y**2;
#at 40m from top
y = 38;
alphah = 0.1;
Cm = 0.72;
Cy = (Cm/2)*((y*(2-y/hw)/hw)+(y*(2-y/hw)/hw)**0.5);
p40 = Cy*alphah*gammaw*hw;
F40 = 0.726*p40*y;
Mp40 = 0.299*p40*y**2;
#at 100m from top
y = 98;
alphah = 0.1;
Cm = 0.72;
Cy = (Cm/2)*((y*(2-y/hw)/hw)+(y*(2-y/hw)/hw)**0.5);
p100 = Cy*alphah*gammaw*hw;
F100 = 0.726*p100*y;
Mp100 = 0.299*p100*y**2;
p10 = round(p10*1000)/1000;
F10 = round(F10*1000)/1000;
Mp10 = round(Mp10*10)/10;
p40 = round(p40*1000)/1000;
F40 = round(F40*1000)/1000;
Mp40 = round(Mp40*10)/10;
p100 = round(p100*100)/100;
F100 = round(F100*1000)/1000;
Mp100 = round(Mp100*10)/10;
print "Hydrodynamic Forces:At 10m from top: P = %.2f KN/square m;\
\nF = %.2f KN;\
\nM = %.2f KN-m.At 40m from top: P = %.2f KN/square m.;\
\nF = %.2f KN;\
\nM = %.2f KN-m.At 100m from top: P = %.2f KN/square m;\
\nF = %.2f KN;\
\nM = %.2f KN-m."%(p10,F10,Mp10,p40,F40,Mp40,p100,F100,Mp100);
#vertical component of reservior water on horizontal section
s1 = 3./60;
Wh = (F100-F40)*s1;
Wh = round(Wh*100)/100;
print "vertical component of reservior water on horizontal section = %.2f kN/m."%(Wh);
```

In [9]:

```
import math
#no tension is permissible
#factor of safety against slidingis 1.5
#Given
wb = 3; #width of dam;
miu = 0.5; #coefficient of friction
Sg = 2.4; #specific gravity of masonary
gamma_w = 9.81; #unit weigth of water
c = 1;
#when uplift is considered
#when no tension is permissible then e = wb/6;
p1 = wb*Sg*gamma_w;
p2 = c*wb*gamma_w/2;
p3 = p1-p2;
p4 = p1*wb/2-p2*2;
p5 = gamma_w/6;
d1 = p4/p3; d2 = p5/p3;
d3 = 1.5-d1;
H = ((0.5-d3)/d2)**0.5;
H = round(H*100)/100;
print "when uplift is considered:"
print "Heigth of dam when no tension is permissible = %.2f m."%(H);
H = p3*0.5/(1.5*p5*3);
print "Heigth of dam when factor of safety against sliding is 1.5 = %.2f m."%(H);
#when uplift is not considered
p1 = wb*Sg*gamma_w;
p4 = p1*wb/2;
p5 = gamma_w/6;
d1 = p4/p1;
d2 = p5/p1;
H = (0.5/d2)**0.5;
H = round(H*100)/100;
print "when uplift is not considered:"
print "Heigth of dam when no tension is permissible = %.2f m."%(H);
H = p1*0.5/(1.5*p5*3);
print "Heigth of dam when factor of safety against sliding is 1.5 = %.2f m."%(H);
```

In [10]:

```
import math
#Given
c = 1;
hw = 6; #heigth of water in reservior
Bt = 1.5; #width of top of dam
H = 6; #heigth of the dam
wb = 4.5; #width of base of dam
Sg = 2.4; #specific gravity of masonary
gamma_w = 9.81; #weigth density of water
W1 = Bt*gamma_w*Sg*H;
W2 = gamma_w*Sg*H*(wb-Bt)/2;
L1 = (wb-Bt)+(Bt/2);
L2 = (2*(wb-Bt))/3
M1 = W1*L1
M2 = W2*L2
#Reaervior empty
SumW = W1+W2;
SumM = M1+M2;
x = SumM/SumW;
e = wb/2-x;
pnt = (SumW/wb)*(1+(6*e/wb));
pnh = (SumW/wb)*(1-(6*e/wb));
pnt = round(pnt*10)/10;
pnh = round(pnh*10)/10;
print "Reservior empty:";
print " Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
#Reservior full
W3 = gamma_w*H**2/2;
U = gamma_w*H*c*wb/2;
SumV = SumW-U;
L3 = hw/3;
L4 = 2*wb/3; #lever arm
M3 = W3*L3;
M4 = U*L4; #moment about toe
SumM1 = SumM-M4-M3;
x = SumM1/SumV;
e = wb/2-x;
pnt = (SumV/wb)*(1+(6*e/wb));
pnh = (SumV/wb)*(1-(6*e/wb));
pnt = round(pnt*10)/10;
pnh = round(pnh*10)/10;
print "Reservior full:";
print " Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
```

In [11]:

```
import math
from numpy import roots
#check the stability
#Given
c = 1;
hw = 6; #heigth of water in reservior
Bt = 1.5; #width of top of dam
H = 6; #heigth of the dam
gamma_m = 20; #unit weigth of masonary
gamma_w = 9.81; #weigth density of water
f = 1800; #compressive strength
miu = 0.6; #coefficient of friction
#to develop no tension e = b/6;x = b/3.
#hence on solving the relations we get
P = [1,2.944,-39.074]
wb = roots(P)[1]; #sign of coefficient is 2.944 is not taken correctly in book
#roots are 4.94 and -7.89
#math.since negative value cannot be taken
print "Neglecting the negative value.Width of base is = 4.94 m.";
W1 = Bt*gamma_m*H;
W2 = gamma_m*H*(wb-Bt)/2;
L1 = (wb-Bt)+(Bt/2);
L2 = (2*(wb-Bt))/3;
M1 = W1*L1,
M2 = W2*L2;
U = gamma_w*H*c*wb/2;
L4 = 2*wb/3;
M4 = U*L4;
W3 = gamma_w*H**2/2;
L3 = hw/3;
M3 = W3*L3;
SumW = W1+W2-U;
SumM = M1+M2-M4-M3;
pn = 2*SumW/wb;
pn = round(pn*10)/10;
print "Maximum stress = %.2f kN/square.m."%(pn);
print "Dam is safe against compression";
FOS = miu*SumW/W3;
FOS = round(FOS*100)/100;
print "Factor of safety against sliding = %.2f. <1"%(FOS);
print "Dam is unsafe against sliding.";
```

In [12]:

```
import math
from numpy import roots
#check the stability if uplift is neglected
#Given
c = 1;
hw = 6; #heigth of water in reservior
Bt = 1.5; #width of top of dam
H = 6; #heigth of the dam
gamma_m = 20; #unit weigth of masonary
gamma_w = 9.81; #weigth density of water
f = 1800; #compressive strength
miu = 0.6; #coefficient of friction
#to develop no tension e = b/6;x = b/3.
#hence on solving the relations we get
P = [1,1.5,-19.908]
wb = roots(P)[1];
#roots are 3.774 and -5.27
#math.since negative value cannot be taken
print "Neglecting the negative value.Width of base is = %.2f m."%wb;
W1 = Bt*gamma_m*H;
W2 = gamma_m*H*(wb-Bt)/2;
L1 = (wb-Bt)+(Bt/2);
L2 = (2*(wb-Bt))/3;
M1 = W1*L1,
M2 = W2*L2;
W3 = gamma_w*H**2/2;
L3 = hw/3;
M3 = W3*L3;
SumW = W1+W2;
SumM = M1+M2-M3;
pn = 2*SumW/wb;
pn = round(pn*10)/10;
print "Maximum stress = %.2f kN/square.m."%(pn);
print "Dam is safe against compression";
FOS = miu*SumW/W3;
FOS = round(FOS*1000)/1000;
print "Factor of safety against sliding = %.2f. > 1"%(FOS);
print "Dam is safe against sliding.";
```

In [13]:

```
import math
# calculate maximum permissible heigth of shutter so that no tension develops
#Given
Bt = 3; #width of top of dam
H = 12; #heigth of the dam
wb = 9; #width of base of dam
gamma_m = 21; #unit weigth of masonary
gamma_w = 9.81; #weigth density of water
W1 = Bt*gamma_m*H;
W2 = gamma_m*H*(wb-Bt)/2;
#taking moment about a point on base at 3m from toe
L1 = 3+Bt/2;
L2 = (2*(wb-Bt)/3)-3;
M1 = W1*L1
M2 = W2*L2;
M = M1+M2;
#net moment about this point should be zero for equilibrium
s = (M*6/gamma_w)**(1./3)-12;
s = round(s*100)/100;
# Results
print "maximum permissible heigth of shutter = %.2f m."%(s);
```

In [14]:

```
import math
#moment at 50m below water surface
#Given
c = 1;
H = 100; #heigth of dam
hw = 100; #heigth of water in reservior
FB = 1; #free board
s = 0.15; #slope of upstream face
gamma_w = 9.81; #unit weigth of water
alphah = 0.1;
theta = math.atan(s);
y = 50;
Cm = 0.735*(1-(theta*2/math.pi));
Cy = (Cm/2)*((y*(2-y/hw)/hw)+(y*(2-y/hw)/hw)**0.5);
pe = Cy*alphah*gamma_w*hw;
F = 0.726*pe*y;
M = 0.299*pe*y**2;
pe = round(pe*1000)/1000;
F = round(F*10)/10;
M = round(M*10)/10;
print "hydrodynamic earthquake pressure = %.2f kN/square.mshear = %.2f kN/m.Moment = %.2f kN-m/m."%(pe,F,M);
```

In [15]:

```
import math
from sympy import sec
#check stability
#Given
c = 1.;
H = 10.; #heigth of dam
hw = 10.; #heigth of water in reservior
wb = 8.25; #bottom width
Bt = 1.; #top width
Hs1 = 0.1; #slope on upstream side
gamma_w = 9.81; #unit weigth of water
gamma_m = 22.4; #unit weigth of masonary
f = 1400.; #permissible shear stress at joint
miu = 0.75; #coefficient of friction
fi = math.atan(0.625);
theta = math.atan(0.1);
W1 = Bt*H*gamma_m;
W2 = H*H*Hs1*gamma_m/2;
W3 = H*6.25*gamma_m/2;
W4 = hw*gamma_w*H*Hs1/2;
P = gamma_w*hw**2/2;
U = wb*gamma_w*hw*c/2;
SumV = W1+W2+W3+W4-U;
L3 = 2*(wb-(Hs1*H)-Bt)/3;
L1 = (wb-(Hs1*H)-Bt)+Bt/2;
L2 = (wb-(Hs1*H)-Bt)+Bt+(Hs1*H/3);
L4 = (wb-(Hs1*H)-Bt)+Bt+(2*Hs1*H/3);
L5 = 2*wb/3;L6 = hw/3;
M1 = W1*L1;
M2 = W2*L2;
M3 = W3*L3;
M4 = W4*L4;
M5 = U*L5;
M6 = P*L6;
SumM = M1+M2+M3+M4-M5-M6;
Mplus = M1+M2+M3+M4;
Mminus = M5+M6;
FOS = miu*SumV/P;
SFF = (miu*SumV+wb*1400)/P;
FOO = Mplus/Mminus;
FOS = round(FOS*100)/100;
SFF = round(SFF*10)/10;
FOO = round(FOO*100)/100;
print "Factor of safety against sliding = %.2f. >1 "%(FOS);
print "Shear friction factor = %.2f."%(SFF);
print "Factor of safety against overturning = %.2f. <1.5"%(FOO);
print "Dam is unsafe against overturning";
x = SumM/SumV;
e = wb/2-x;
p = hw*gamma_w;
pnt = (SumV/wb)*(1+(6*e/wb)); #calculation is done wrong in book;value of b is not taken correctly
pnh = (SumV/wb)*(1-(6*e/wb));
sigmat = pnt*sec(fi)**2;
sigmah = pnh*sec(theta)**2-p*math.tan(theta)**2;
taut = pnt*math.tan(fi);
tauh = -(pnh-p)*math.tan(theta);
pnt = round(pnt*10)/10;
pnh = round(pnh*10)/10;
sigmat = round(sigmat*10)/10;
sigmah = round(sigmah*10)/10;
taut = round(taut*10)/10;
tauh = round(tauh*10)/10;
print "Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
print "Principal stress at toe = %.2f kN/square.m."%(sigmat);
print "Principal stress at heel = %.2f kN/square.m."%(sigmah);
print "Shear stress at toe = %.2f kN/square.m."%(taut);
print "Shear stress at heel = %.2f kN/square.m."%(tauh);
```

In [16]:

```
import math
#Check the stability and determine sliding factor and shear factor
#Given
c = 1;
miu = 0.75; #coefficient of friction
H = 90; #heigth of dam
wb = 73.1; #width of base
Bt = 7; #width of top of dam
hw = 89; #heigth of water in reservior
Hs1 = 28; #heigth of slope on upstream side
Hs2 = 83; #heigth of slope on downstream side
Cm = 0.735;
alphah = 0.1;
gamma_m = 23.5; #unit weigth of concrete
gamma_w = 9.81; #unit weigth of water
theta = math.atan(8./28);
fi = math.atan(0.7);
#self weigth of dam
W1 = (Hs1*8*gamma_m)/2
W2 = (Bt*H*gamma_m)
W3 = (Hs2**2*0.7*gamma_m)/2
#weigth of superimposed water
W4 = (Hs1*8*gamma_w)/2
W5 = (hw-Hs1)*8*gamma_w
U = hw*wb*2*gamma_w/6; #uplift force
wp = hw**2*gamma_w/2; #water pressure
hp = 0.726*Cm*alphah*gamma_w*hw**2; #hydrodynamic pressure
Mhp = 0.299*Cm*alphah*gamma_w*hw**3; #moment due to hydrodynamic pressure
#inertial load due to horizontal acceleration
I1 = W2/10;
I2 = W3/10;
I3 = W1/10;
SumV = W1+W2+W3+W4+W5-U;
SumH = wp+hp+I1+I2+I3;
L1 = (wb-8)+8/3
L2 = (0.7*Hs2)+(Bt/2)
L3 = (2*Hs2*0.7)/3.
L4 = (wb-8)+(2*8)/3.
L5 = (wb-8)+(8./2)
L6 = hw/3;
L7 = 2*wb/3;
M1 = W1*L1
M2 = W2*L2
M3 = W3*L3
M4 = W4*L4;
M5 = W5*L5;
M6 = wp*L6;
M7 = U*L7;
M8 = I1*45;
M9 = I2*83/3;
M10 = I3*28/3;
Mplus = M1+M2+M3+M4+M5;
Mminus = M6+M7+M8+M9+M10+Mhp;
SumM = Mplus-Mminus;
x = SumM/SumV;
e = wb/2-x;
pnt = (SumV/wb)*(1+(6*e/wb));
pnh = (SumV/wb)*(1-(6*e/wb));
sigmat = pnt*sec(fi)**2;
p = hw*gamma_w;
pe = Cm*alphah*gamma_w*hw;
sigmah = pnh*sec(theta)**2-(p+pe)*math.tan(theta)**2;
taut = pnt*math.tan(fi);
tauh = -(-pnh-(p+pe))*math.tan(theta);
print "Normal stress at toe = %i kN/square.m."%(pnt);
print "Normal stress at heel = %i kN/square.m."%(pnh);
print "Principal stress at toe = %i kN/square.m."%(sigmat);
print "Principal stress at heel = %i kN/square.m."%(sigmah);
print "Shear stress at toe = %i kN/square.m."%(taut);
print "Shear stress at heel = %i kN/square.m."%(tauh);
FOS = miu*SumV/SumH;
SFF = (miu*SumV+wb*1400)/SumH;
FOO = Mplus/Mminus;
Ffi = 1.2;Fc = 2.4;
F = (miu*SumV/Ffi+1400*wb/Fc)/SumH;
FOS = round(FOS*100)/100;
F = round(F*100)/100;
SFF = round(SFF*100)/100;
FOO = round(FOO*100)/100;
print "Factor of safety against sliding as per IS:6512-1972 = %.2f. <1.5"%(FOS);
print "Factor of safety against sliding as per IS:6512-1984 = %.2f. >1"%(F);
print "Shear friction factor = %.2f. <6"%(SFF);
print "Factor of safety against overturning = %.2f. <1.5"%(FOO);
print "Dam is unsafe for given loading conditions";
```

In [17]:

```
import math
#Check the stability and determine principal and shear stress at toe and heel
#Given
c = 1;
miu = 0.7; #coefficient of friction
H = 70; #heigth of dam
ht = 0; #heigth of tail water
Lf = 6.5; #location of foundation gallery from heel
wb = 52.5; #width of base
Bt = 7; #width of top of dam
hw = 70; #heigth of water in reservior
Hs1 = 35; #heigth of slope on upstream side
Hs2 = 60; #heigth of slope on downstream side
gamma_m = 24; #unit weigth of concrete
gamma_w = 9.81; #unit weigth of water
theta = math.atan(0.1);
fi = math.atan(0.7);
#self weigth of dam
W1 = (Hs1*3.5*gamma_m)/2
W2 = (Bt*H*gamma_m)
W3 = (Hs2**2*0.7*gamma_m)/2
#weigth of superimposed water
W4 = (Hs1*3.5*gamma_w)/2
W5 = (hw-Hs1)*3.5*gamma_w
wp = hw**2*gamma_w/2; #water pressure
Pt = gamma_w*ht
Ph = gamma_w*hw
Pg = (ht+(hw-ht)/3)*gamma_w
U = (Pt*(wb-Lf))+(Pg*Lf)+((Ph-Pg)*Lf/2)+((Pg-Pt)*(wb-Lf)/2)*c
l1 = (wb-Lf)/2
l2 = (2*(wb-Lf))/3
l3 = (wb-Lf)+(Lf/2)
l4 = (wb-Lf)+((2*Lf)/3)
L7 = (((Pt*(wb-Lf))*l1)+((Pg-Pt)*(wb-Lf)*l2/2)+((Pg*Lf)*l3)+((Ph-Pg)*Lf*l4/2))/U
L1 = (wb-3.5)+3.5/3
L2 = (0.7*Hs2)+(Bt/2)
L3 = (2*Hs2*0.7)/3
L4 = (wb-3.5)+(2*3.5)/3
L5 = (wb-3.5)+(3.5/2)
L6 = hw/3;
M1 = W1*L1
M2 = W2*L2
M3 = W3*L3
M4 = W4*L4;
M5 = W5*L5;
M6 = wp*L6;
M7 = U*L7;
SumV1 = W1+W2+W3;
SumM1 = M1+M2+M3;
SumV2 = SumV1+W4+W5;
SumM2 = SumM1+M4+M5-M6;
SumV3 = SumV2-U;
SumM3 = SumM2-M7;
Mplus = 1547377;
Mminus = 870421;
SumH = wp;
#case 1. Reservior empty
x = SumM1/SumV1;
e = wb/2-x;
pnt = (SumV1/wb)*(1+(6*e/wb));
pnh = (SumV1/wb)*(1-(6*e/wb));
sigmat = pnt*sec(fi)**2;
sigmah = pnh*sec(theta)**2;
taut = pnt*math.tan(fi);
tauh = pnh*math.tan(theta);
pnt = round(pnt*10)/10;
pnh = round(pnh*10)/10;
sigmat = round(sigmat*10)/10;
sigmah = round(sigmah*10)/10;
taut = round(taut*10)/10;
tauh = round(tauh*10)/10;
print "case 1. Reservior empty:";
print "Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
print "Principal stress at toe = %.2f kN/square.m."%(sigmat);
print "Principal stress at heel = %.2f kN/square.m."%(sigmah);
print "Shear stress at toe = %.2f kN/square.m."%(taut);
print "Shear stress at heel = %.2f kN/square.m."%(tauh);
#case2. reservior full without uplift
x = SumM2/SumV2;
e = wb/2-x;
p = hw*gamma_w;
pnt = (SumV2/wb)*(1+(6*e/wb));
pnh = (SumV2/wb)*(1-(6*e/wb));
sigmat = pnt*sec(fi)**2;
sigmah = pnh*sec(theta)**2-p*math.tan(theta)**2;
taut = pnt*math.tan(fi);
tauh = -(pnh-p)*math.tan(theta);
pnt = round(pnt*10)/10;
pnh = round(pnh*10)/10;
sigmat = round(sigmat*10)/10;
sigmah = round(sigmah*10)/10;
taut = round(taut*10)/10;
tauh = round(tauh*10)/10;
print "case 2. reservior full without uplift:";
print "Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
print "Principal stress at toe = %.2f kN/square.m."%(sigmat);
print "Principal stress at heel = %.2f kN/square.m."%(sigmah);
print "Shear stress at toe = %.2f kN/square.m."%(taut);
print "Shear stress at heel = %.2f kN/square.m."%(tauh);
#case3. reservior full with uplift
x = SumM3/SumV3;
e = wb/2-x;
p = hw*gamma_w;
pnt = (SumV3/wb)*(1+(6*e/wb));
pnh = (SumV3/wb)*(1-(6*e/wb));
sigmat = pnt*sec(fi)**2;
sigmah = pnh*sec(theta)**2-p*math.tan(theta)**2;
taut = pnt*math.tan(fi);
tauh = -(pnh-p)*math.tan(theta);
pnt = round(pnt);
pnh = round(pnh);
sigmat = round(sigmat);
sigmah = round(sigmah);
taut = round(taut);
tauh = round(tauh);
print "case 3. reservior full with uplift:";
print "Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
print "Principal stress at toe = %.2f kN/square.m."%(sigmat);
print "Principal stress at heel = %.2f kN/square.m."%(sigmah);
print "Shear stress at toe = %.2f kN/square.m."%(taut);
print "Shear stress at heel = %.2f kN/square.m."%(tauh);
FOS = miu*SumV3/SumH;
SFF = (miu*SumV3+wb*1400)/SumH;
FOO = Mplus/Mminus;
Ffi = 1.5;Fc = 3.6;
F = (miu*SumV3/Ffi+1400*wb/Fc)/SumH;
FOS = round(FOS*1000)/1000;
SFF = round(SFF*100)/100;
FOO = round(FOO*100)/100;
F = round(F*1000)/1000;
print "Factor of safety against sliding = %.2f."%(FOS);
print "Shear friction factor = %.2f."%(SFF);
print "Factor of safety against overturning = %.2f."%(FOO);
print "Factor of safety for load combination B = %.2f. > 1"%(F);
print "Dam is safe ";
#Case4.considering seismic forces
Cm = 0.712;
alphah = 0.1;
alphav = 0.08;
hp = 0.726*Cm*alphah*gamma_w*hw**2; #hydrodynamic pressure
Mhp = 0.299*Cm*alphah*gamma_w*hw**3; #moment due to hydrodynamic pressure
#inertial load due to horizontal acceleration
I1 = W2/10;
I2 = W3/10;
I3 = W1/10;
v = SumV1*alphav;
Mv = 116444;
SumV4 = SumV3-v;
SumH1 = SumH+I1+I2+I3+hp;
M8 = I1*35;
M9 = I2*20;
M10 = I3*35/3;
Mminus1 = 1161849;
SumM4 = SumM3-M8-M9-M10-Mhp-Mv;
x = SumM4/SumV4;
e = wb/2-x;
p = hw*gamma_w;
pe = Cm*alphah*gamma_w*hw;
pnt = (SumV4/wb)*(1+(6*e/wb));
pnh = (SumV4/wb)*(1-(6*e/wb));
sigmat = pnt*sec(fi)**2;
sigmah = pnh*sec(theta)**2-(p+pe)*math.tan(theta)**2;
taut = pnt*math.tan(fi);
tauh = (-pnh+(p+pe))*math.tan(theta);
pnt = round(pnt);
pnh = round(pnh);
sigmat = round(sigmat);
sigmah = round(sigmah);
taut = round(taut);
tauh = round(tauh);
print "case 4.considering seismic forces";
print "Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Normal stress at heel = %.2f kN/square.m."%(pnh);
print "Principal stress at toe = %.2f kN/square.m."%(sigmat);
print "Principal stress at heel = %.2f kN/square.m."%(sigmah);
print "Shear stress at toe = %.2f kN/square.m."%(taut);
print "Shear stress at heel = %.2f kN/square.m."%(tauh); #answer is wrong in book
FOS = miu*SumV4/SumH1;
SFF = (miu*SumV4+wb*1400)/SumH1;
FOO = Mplus/Mminus1;
Ffi = 1.2;Fc = 2.7;
F = (miu*SumV4/Ffi+1400*wb/Fc)/SumH1;
FOS = round(FOS*1000)/1000;
SFF = round(SFF*100)/100;
FOO = round(FOO*100)/100;
F = round(F*100)/100;
print "Factor of safety against sliding = %.2f."%(FOS);
print "Shear friction factor = %.2f."%(SFF);
print "Factor of safety against overturning = %.2f."%(FOO);
print "Factor of safety for load combination E = %.2f. > 1"%(F);
print "Dam is safe ";
```

In [18]:

```
import math
#design practical profile of gravity dam
#Given
c = 1;
rlb = 1450; #R.L of base of dam
rlw = 1480.5; #R.L of water level
Sg = 2.4; #specific gravity of masonary
gamma_w = 9.81; #unit weigth of water
w = 1; #heigth of waves
f = 1200; #safe compressive stress for masonary
FB = 1.5*w;
rlt = FB+rlw; #R.L of top of dam
H = rlt-rlb; #heigth of dam
LH = f/(gamma_w*(Sg+1))
LH = round(LH*100)/100;
print "Heigth of dam = %.2f m."%(H);
print "limiting heigth of dam = %.2f m."%(LH);
print "Dam is low gravity dam";
hw = rlw-rlb;
#keep top width,a = 4.5.
a = 4.5;
P = hw/(Sg**0.5);
P = round(P*10)/10;
print "Base width of elementary profile = %.2f m."%(P);
uo = a/16;
wb = uo+P;
wb = round(wb);
print "Base width = %.2f m."%(wb);
D = 2*a*(Sg**0.5);
D = round(D);
print "Dismath.tance upto which u/s slope is vertical from water level = %.2f m."%(D);
```

In [24]:

```
import math
#determine if dam is safe against sliding
#Given
hw = 97.; #heigth of water in reservior
Bt = 7.; #width of top of dam
H = 100.; #heigth of the dam
Hs2 = 90.; #heigth of slope on downstream side
wb = 75.; #width of base of dam
miu = 0.75; #coefficient of friction
gamma_d = 2.4; #weigth density of concrete
gamma_w = 1000.; #weigth density of water
# Calculations
P = gamma_w*hw**2/(2*1000);
W1 = Bt*gamma_d*H;
W2 = (wb-Bt)*Hs2*gamma_d/2;
W = W1+W2;
FOS = miu*W/P;
FOS = round(FOS*1000)/1000;
# Results
print "Factor of safety against sliding = %.2f."%(FOS);
print "Dam is safe against sliding";
```

In [23]:

```
import math
#Factor of safety against overturning
#Factor of safety against sliding
#Shear friction factor
#Given
c = 1.;
H = 10.; #heigth of dam
hw = 10.; #heigth of water in reservior
wb = 8.25; #bottom width
Bt = 1.; #top width
Hs1 = 0.1; #slope on upstream side
gamma_w = 9.81; #unit weigth of water
gamma_m = 22.4; #unit weigth of masonary
f = 1400.; #permissible shear stress at joint
miu = 0.75; #coefficient of friction
fi = math.atan(0.625);
theta = math.atan(0.1);
# Calculations
W1 = Bt*H*gamma_m;
W2 = H*H*Hs1*gamma_m/2;
W3 = H*6.25*gamma_m/2;
W4 = hw*gamma_w*H*Hs1/2;
P = gamma_w*hw**2/2;
U = wb*gamma_w*hw*c/2;
SumV = W1+W2+W3+W4-U;
L3 = 2*(wb-(Hs1*H)-Bt)/3;
L1 = (wb-(Hs1*H)-Bt)+Bt/2;
L2 = (wb-(Hs1*H)-Bt)+Bt+(Hs1*H/3);
L4 = (wb-(Hs1*H)-Bt)+Bt+(2*Hs1*H/3);
L5 = 2*wb/3;L6 = hw/3;
M1 = W1*L1;M2 = W2*L2;M3 = W3*L3;M4 = W4*L4;
M5 = U*L5;M6 = P*L6;
SumM = M1+M2+M3+M4-M5-M6;
Mplus = M1+M2+M3+M4;
Mminus = M5+M6;
FOS = miu*SumV/P;
SFF = (miu*SumV+wb*1400)/P;
FOO = Mplus/Mminus;
FOS = round(FOS*100)/100;
SFF = round(SFF*10)/10;
FOO = round(FOO*100)/100;
# Results
print "Factor of safety against sliding = %.2f."%(FOS);
print "Shear friction factor = %.2f."%(SFF);
print "Factor of safety against overturning = %.2f."%(FOO);
print "Dam is unsafe against overturning";
```

In [22]:

```
import math
#Given
c = 1.;
hw = 80.; #heigth of water in reservior
Bt = 6.; #width of top of dam
H = 84.; #heigth of the dam
Hs2 = 75.; #heigth of slope on downstream side
wb = 56.; #width of base of dam
Lf = 8.; #dismath.tance of foundation gallery from heel
gamma_d = 23.5; #weigth density of concrete
gamma_w = 9.81; #weigth density of water
ht = 6.; #heigth of tail water
# Calculations
W1 = Bt*gamma_d*H;
W2 = gamma_d*Hs2*(wb-Bt)/2;
W3 = gamma_w*ht*4/2;
W4 = gamma_w*hw**2/2;
W5 = gamma_w*ht**2/2;
Pt = gamma_w*ht
Ph = gamma_w*hw
Pg = (ht+(hw-ht)/3)*gamma_w
U = (Pt*(wb-Lf))+(Pg*Lf)+((Ph-Pg)*Lf/2)+((Pg-Pt)*(wb-Lf)/2)*c
l1 = (wb-Lf)/2
l2 = (2*(wb-Lf))/3
l3 = (wb-Lf)+(Lf/2)
l4 = (wb-Lf)+((2*Lf)/3)
L6 = (((Pt*(wb-Lf))*l1)+((Pg-Pt)*(wb-Lf)*l2/2)+((Pg*Lf)*l3)+((Ph-Pg)*Lf*l4/2))/U
L1 = (wb-Bt)+(Bt/2)
L2 = (2*(wb-Bt))/3
L3 = 4./3;
L4 = hw/3;
L5 = ht/3;
M1 = W1*L1
M2 = W2*L2
M3 = W3*L3
M4 = W4*L4
M5 = W5*L5
M6 = U*L6;
SumV = W1+W2+W3-U;
SumH = W4-W5;
SumM = M1+M2+M3-M4+M5-M6;
x = SumM/SumV;
e = wb/2-x;
pnt = (SumV/wb)*(1+(6*e/wb));
pnh = (SumV/wb)*(1-(6*e/wb));
pnt = round(pnt*10)/10;
pnh = round(pnh*10)/10;
# Results
print "Maximum Normal stress at toe = %.2f kN/square.m."%(pnt);
print "Maximum Normal stress at heel = %.2f kN/square.m."%(pnh);
```