BEGIN "HAMMER: Water hammer analysis"
INTEGER Number!Of!Nodes;
REAL PROCEDURE Regula!Falsi(REAL X0,X1,Accuracy; REAL PROCEDURE F);
Begin"Regula Falsi According to Firouztale, Spencer & Wright"
Real X2;
While (abs(f(x2)-f(x1)) geq accuracy) do
x0:=((x2:=x1)-(x1:=x0))*f(x1)/(f(x1)-f(x2))+x1;
Return(x0);
End"Regula Falsi According to Firouztale, Spencer & Wright";
REAL PROCEDURE Redlich!Kwong(REAL Pressure!Reduced, Temperature!Reduced,
Accuracy);
BEGIN "Redlich!Kwong"
REAL h,k,z;
COMMENT : Z is Compressibility Factor, Z=PV/R/T, ie correction
factor leading to Redlich & Kwong Equation of State.
Smith & Van Ness, Introduction to Chemical
Engineering Thermodynamics;
z:=1;
While (abs(z-k) > Accuracy) do
z:=1/(1-(h:=.0867*Pressure!Reduced/Temperature!Reduced
/(k:=z)))-4.93*h*Temperature!Reduced^-1.5/(1+h);
Return(z);
END "Redlich!Kwong";
Real Procedure Bulk!Modulus (REAL Pressure, Critical!Pressure,
Temperature, Critical!Temperature,
Small!Number, R!Ideal!Constant);
BEGIN "Bulk Modulus from Thermodynamics"
REAL Z, dZ!dP;
COMMENT : k=rho(0)*dP/d.rho, P=Z*rho*r*T,rho dP/dZ=-Z(dZ/dP)=K;
Z:=Redlich!Kwong(Pressure/Critical!Pressure,
Temperature/Critical!Temperature,
Small!Number);
dZ!dP:=(Redlich!Kwong(Pressure/Critical!Pressure
*(1+Small!Number),
Temperature/Critical!Temperature,
Small/Number)
-Redlich!Kwong(Pressure/Critical!Pressure
*(1-Small!Number),
Temperature/Critical!Temperature,
Small/Number)
)/2/Small!Number*Critical!Pressure;
RETURN(-Z/dZ!dP);
END "Bulk Modulus from Thermodynamics";
Real Procedure Speed!Of!Wave(REAL Pressure, Critical!Pressure,
Temperature, Critical!Temperature,
Small!Number, R!Ideal!Constant,
Pipe!Diameter,Pipe!Thickness,Young!Modulus,
Pipe!Poisson!Ratio;
BOOLEAN Free!Pipe, Anchored!Pipe,
Expanding!Joints);
BEGIN "Speed of wave from thermodynamic and elastic pipe theories"
REAL Density,K,Poisson!Effect;
Density:=Pressure/R!Ideal!Constant/Temperature/
Redlich!Kwong(Pressure/Critical!Pressure,
Temperature/Critical!Temperature,
Small!Number);
K:=Bulk!Modulus( Pressure, Critical!Pressure,
Temperature, Critical!Temperature,
Small!Number,R!Ideal!Constant);
Poisson!Effect:= if Free!Pipe then (1.25-Pipe!Poisson!Ratio)
else if Anchored!Pipe then (1-Pipe!Poisson!Ratio^2)
else if Expanding!Joints then (1);
RETURN((Density*(1/K+Pipe!Diameter/Pipe!Thickness/Young!Modulus
*Poisson!Effect))^-.5);
END "Speed of wave from thermodynamic and elastic pipe theories";
Print(13&10&"How many nodes (number of sections +1)?>");
Number!Of!Nodes:=cvd(inchwl);
BEGIN "Water Hammer Internals"
REQUIRE "sai:ttyio.hdr" source!file; comment i/o procedures;
INTEGER I,J;
REAL Length,Diameter,Fanning!Friction,Wave!Speed,Initial!Velocity,
Initial!Head,Elevation!Head,End!Time,Valve!Delay,Discharge!Coefficient,
Delta!X,Delta!Time,Delta!Elevation,K!Loss,Time,g!c,Area!Factor;
REAL ARRAY Head,Velocity,X[1:Number!Of!Nodes];
BOOLEAN Non!Linear!Valve;
COMMENT : Central algorithm due to G.Z.Watters,Modern Analsis and
Control of Unsteady Flow in Pipelines and
Professor Richard Skalak CE4261X81 Applied Fluid Mechanics;
g!c:=TtFlt("Gravitational constant>","g sub c","32.174");
Length:=TtFlt("Length?>","homogeneous dimensions","10");
Diameter:=TtFlt("Pipe Diameter?>","homogenious dimension","1");
Fanning!Friction:=TtFlt("Fanning friction factor>","dim'less",".0001");
Initial!Velocity:=TtFlt("Initial velocity>");
Elevation!Head:=TtFlt("Elevation head>");
End!Time:=TtFlt("End time of simulation>","homogeneous dimensions","100");
Valve!Delay:=TtFlt("Time for valve closure>",
"homogeneous dimensions","50");
Non!Linear!Valve:=TtInt("Enter 0 if linear valve at end else discharge",
"Default is linear valve", "0");
If Non!Linear!Valve then
BEGIN "Input Discharge and Initial area fraction"
Discharge!Coefficient:=TtFlt(
"Valve discharge coefficient>","dim'less",".95");
Area!Factor:=TtFlt(
"Ratio of pipe to valve open areas>";
"dim'less",".75");
END "Input Discharge and Initial area fraction";
Wave!Speed:=TtFlt("Wave speed?>",
"default from thermodynamics, unusual # automatic", "517.1961");
If Wave!Speed=517.1961 then
BEGIN "Determine Wave Speed from physical properties"
REAL Pressure,Critical!Pressure,Temperature,
Critical!Temperature, R!Ideal!Constant,Pipe!Thickness,
Young!Modulus,Pipe!Poisson!Ration;
INTEGER Case!Number;
BOOLEAN Free!Pipe,Anchored!Pipe,Expanding!Joints;
REAL PROCEDURE Find!Pressure(REAL Pressure);
BEGIN "P"
Return(
Pressure-(Initial!Head*R!Ideal!Constant
*Temperature/g!c
*Redlich!Kwong(Pressure
/Critical!Pressure,
Temperature
/Critical!Temperature,
.001))^.5
);
END "P";
Temperature:=TtFlt("Temperature>",
"Room,Rankin is default","536.6")
Critical!Temperature:=TtFlt("Critical Temperature> ",
"Water,Rankin is default","1164.8");
R!Ideal!Constant:=TtFlt("Ideal Constant,R> ",
"English is default","10.73");
Pipe!Thickness:=TtFlt("Pipe Thickness>");
Young!Modulus:=TtFlt("Pipe Young's Modulus>",
"Steel is default","30000000");
Pipe!Poisson!Ratio:=TtFlt("Pipe Poisson Ratio>",
"Dimensionless, Steel is default",".3");
Case!Number:=TtInt(
"Free Pipe[1],Anchored[2],Expanding Joints[3]>",
"one of the three numbers,3 is default","3");
Case Case!Number Of
BEGIN "Choose Poisson Effect Pipe Type"
[1]Free!Pipe:=True;
[2]Anchored!Pipe:=True;
[3]Expanding!Joints:=True;
END "Choose Poisson Effect Pipe Type";
Critical!Pressure:=TtFlt("Critical Pressure>",
"Water,PSI is default","3198");
Pressure:=Regula!Falsi(
.8*(Initial!Head*R!Ideal!Constant
*Temperature/g!c)^.5,
1.2*(Initial!Head*R!Ideal!Constant
*Temperature/g!c)^.5,
.0001,Find!Pressure);
Wave!Speed:=Speed!Of!Wave(Pressure,Critical!Pressure,
Temperature,Critical!Temperature,
.0001,R!Ideal!Constant,
Diameter,Pipe!Thickness,Youg!Modulus,
Pipe!Poisson!Ratio,
Free!Pipe,Anchored!Pipe,
Expanding!Joints);
END "Determine Wave Speed from physical properties";
Delta!X:=Length/(Number!Of!Nodes-1);
Delta!Time:=Delta!X/Wave!Speed;
Delta!Elevation:=Elevation!Head/(Number!Of!Numes-1);
K!Loss:=24*Fanning!Friction*Delta!Time/Diameter;
For I:=1 step 1 until Number!Of!Nodes do
BEGIN "Set up arrays"
Velocity[i]:=Initial!Velocity;
Head[i]:=Initial!Head-(i-1)*48*Fanning!Friction*Delta!X
*Initial!Velocity^2/g!c/2/Diameter
X[i]:=(i-1)*Delta!X/Length;
END "Set up arrays";
COMMENT : Pseudo-condition to avoid zero initial difference
May be added by placing semicolon here ->
Velocity[2]:=Velocity[2]*Head[2]/Head[1];
for j:=1 step 1 until (End!Time/Delta!Time+1) do
BEGIN "Step ahead in time"
REAL Valve!loss,Valve!Area!Factor,Head!New,Head!Old!New,
Velocity!New,Old!Velocity!New;
Time:=Time+Delta!Time;
COMMENT : Data for gate valve, McCabe & Smith, Unit Operations
of Chemical Engineering: .5 open,k=56-full open k=.2
=> k:=-1.5*valve!area!factor+11
+Discharge term from Watters;
Valve!Area!Factor:=Area!factor*(1-time/valve!delay);
Valve!Loss:=-1.5*Valve!Area!Factor+11
+(Valve!Area!Factor/Discharge!Coefficient)^2;
For i:=1 step 1 until Number!Of!Nodes do
BEGIN "Each node"
Old!Head!New:=Head!New;
Old!Velocity!New:=Velocity!New;
Velocity!New:=
if i=1 then
Velocity[2]+g!c/Wave!Speed*(Head[1]-Head[2])
-K!Loss*Velocity[2]*abs(Velocity[2])
else if i=Number!Of!Nodes then
if Time>Valve!Delay then 0
else Initial!Velocity*(1-Time/Valve!Delay)
else .5*(Velocity[i-1]+Velocity[i+1]
+g!c/Wave!Speed*(Head[i-1]-Head[i+1]))
-K!Loss*(Velocity[i-1]*ABS(Velocity[i-1))
-Velocity[i+1]*ABS(Velocity[i+1]))
Head!New:=
if i=1 then Initial!Head
else if i=number!Of!Nodes then
(Head[Number!Of!Nodes-1]
+(Velocity[Number!Of!Nodes-1]
-Velocity[Number!Of!Nodes]
-(K!Loss+ Valve!Loss*
Non!Linear!Valve
/Non!Linear!Valve)
*Velocity[Number!Of!Nodes-1]
*abs(
Velocity[Number!Of!Nodes-1]
)
)*Wave!Speed/g!c)
else .5*(Head[i-1]+Head[i+1]
+Wave!Speed/g!c*(Velocity[i-1]
-Velocity[i+1])
-K!Loss*(Velocity[i-1]*ABS(Velocity[i-1])
-Velocity[i+1]*ABS(Velocity[i+1]))
*Wave!Speed/g!c
);
If i=Number!Of!Nodes then
BEGIN "Don't leave last one empty"
Head[i]:=Head!New;
Velocity[i]:=Velocity!New;
END "Don't leave last one empty";
if i=1 then
BEGIN "Don't go passed array bounds"
COMMENT : It's all taken care of;
END "Don't go passed array bounds"
else
BEGIN "Non-boundary indices"
Head[i-1]:=Old!Head!New;
Velocity[i-1]:=Old!Velocity!New;
END "Non-boundary indices";
END "Each node";
SetFormat(6,3);
Print (13&10&"For time=",time,",X,Velocity,Head are:");
For I:=1 step 1 until Number!Of!Nodes do
Print(13&10&null,X[i],Velocity[i],Head[i]);
END "Step ahead in time";
END "Water Hammer Internals";
END "HAMMER: Water hammer analysis";