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";