π¬
Solver Purpose & Physical Scope
Calculate pipeline water hammer surge pressure (Joukowsky & Allievi): acoustic wave speed in elastic pipes (a), pipeline wave period (2L/a), maximum pressure surge (bar), and pipe wall hoop stress.
π Discipline: Fluid-mechanics
β‘ Precision: IEEE-754 64-bit Real(`real(8)`)
π₯ Total Downloads: 258 times
π Source File:
water_hammer_allievi_surge.f90
π calcul/FluidMechanics /
water_hammer_allievi_surge.f90
program water_hammer_allievi_surge
implicit none
integer :: iostat_val, pipe_material, fluid_type
double precision :: L_pipe_m, D_pipe_mm, wall_thick_mm, Q_flow_m3h, P_static_bar, Tc_valve_sec
double precision :: rho, K_bulk, E_pipe, nu_poisson, D_m, e_m, A_pipe, V0_ms
double precision :: a_wave_ms, Tp_pipe_sec, dP_joukowsky_bar, dH_joukowsky_m
double precision :: dP_actual_bar, P_peak_bar, hoop_stress_MPa, allievi_theta
character(len=32) :: closure_type, safety_status
double precision, parameter :: PI = 3.141592653589793d0
! Read inputs
read(*,*,iostat=iostat_val) pipe_material ! 1=Carbon Steel (E=205 GPa), 2=Ductile Iron (E=170 GPa), 3=HDPE (E=1.0 GPa), 4=PVC (E=3.0 GPa)
read(*,*,iostat=iostat_val) fluid_type ! 1=Water (K=2.19 GPa), 2=Crude Oil (K=1.5 GPa), 3=Ethanol (K=1.1 GPa)
read(*,*,iostat=iostat_val) L_pipe_m ! Pipeline Length L [m] (e.g. 1500.0)
read(*,*,iostat=iostat_val) D_pipe_mm ! Pipe Inner Diameter [mm] (e.g. 300.0)
read(*,*,iostat=iostat_val) wall_thick_mm ! Wall Thickness [mm] (e.g. 8.0)
read(*,*,iostat=iostat_val) Q_flow_m3h ! Flow Rate [m3/h] (e.g. 500.0)
read(*,*,iostat=iostat_val) P_static_bar ! Static Operating Pressure [bar] (e.g. 6.0)
read(*,*,iostat=iostat_val) Tc_valve_sec ! Valve Closure Time [s] (e.g. 2.5)
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid input data for water hammer calculation.'
stop
end if
if (L_pipe_m <= 0.0d0 .or. D_pipe_mm <= 0.0d0 .or. wall_thick_mm <= 0.0d0 .or. Q_flow_m3h <= 0.0d0) then
write(*,*) 'ERROR: Pipe geometry, length, and flow rate must be positive.'
stop
end if
! Fluid properties
if (fluid_type == 1) then ! Water
rho = 998.0d0; K_bulk = 2.19d9
else if (fluid_type == 2) then ! Oil
rho = 850.0d0; K_bulk = 1.50d9
else ! Ethanol
rho = 789.0d0; K_bulk = 1.10d9
end if
! Pipe material elasticity
if (pipe_material == 1) then ! Carbon Steel
E_pipe = 2.05d11; nu_poisson = 0.30d0
else if (pipe_material == 2) then ! Ductile Iron
E_pipe = 1.70d11; nu_poisson = 0.28d0
else if (pipe_material == 3) then ! HDPE
E_pipe = 1.0d9; nu_poisson = 0.45d0
else ! PVC
E_pipe = 3.0d9; nu_poisson = 0.38d0
end if
D_m = D_pipe_mm / 1000.0d0
e_m = wall_thick_mm / 1000.0d0
A_pipe = PI * (D_m**2) / 4.0d0
V0_ms = (Q_flow_m3h / 3600.0d0) / A_pipe
! Acoustic Wave Speed [m/s]
a_wave_ms = sqrt((K_bulk / rho) / (1.0d0 + (K_bulk / E_pipe) * (D_m / e_m) * (1.0d0 - nu_poisson**2)))
! Pipeline Characteristic Period [s]
Tp_pipe_sec = 2.0d0 * L_pipe_m / a_wave_ms
! Joukowsky Maximum Potential Surge (Rapid Closure)
dP_joukowsky_bar = (rho * a_wave_ms * V0_ms) / 1.0d5
dH_joukowsky_m = (a_wave_ms * V0_ms) / 9.80665d0
! Actual Surge based on Closure Time
if (Tc_valve_sec <= Tp_pipe_sec) then
closure_type = 'RAPID CLOSURE (Tc <= 2L/a)'
dP_actual_bar = dP_joukowsky_bar
else
closure_type = 'SLOW CLOSURE (Tc > 2L/a)'
dP_actual_bar = dP_joukowsky_bar * (Tp_pipe_sec / Tc_valve_sec)
end if
P_peak_bar = P_static_bar + dP_actual_bar
allievi_theta = a_wave_ms * Tc_valve_sec / (2.0d0 * L_pipe_m)
! Pipe Hoop Stress [MPa]
hoop_stress_MPa = ((P_peak_bar * 1.0d5) * D_m) / (2.0d0 * e_m * 1.0d6)
if (pipe_material == 1 .and. hoop_stress_MPa < 140.0d0) then
safety_status = 'SAFE (BELOW STEEL ALLOWABLE)'
else if (pipe_material == 3 .and. hoop_stress_MPa < 10.0d0) then
safety_status = 'SAFE (BELOW HDPE ALLOWABLE)'
else if (pipe_material == 4 .and. hoop_stress_MPa < 14.0d0) then
safety_status = 'SAFE (BELOW PVC ALLOWABLE)'
else
safety_status = 'WARNING: HIGH SURGE STRESS'
end if
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' THERMOFLUIDCALC β WATER HAMMER & SURGE PRESSURE ANALYSIS'
write(*,'(A)') '============================================================'
write(*,'(A,F10.1,A,F10.2,A)') 'Wave Speed (a) / Flow Velocity= ', a_wave_ms, ' m/s / ', V0_ms, ' m/s'
write(*,'(A,F10.2,A,F10.2,A)') 'Pipe Period (Tp) / Closure Tc = ', Tp_pipe_sec, ' s / ', Tc_valve_sec, ' s'
write(*,'(A,A)') 'Closure Regime Classification = ', trim(closure_type)
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,F10.2,A,F10.2,A)') 'MAX JOUKOWSKY SURGE (ΞP_max) = ', dP_joukowsky_bar, ' bar (', dH_joukowsky_m, ' m)'
write(*,'(A,F10.2,A)') 'ACTUAL SURGE PRESSURE (ΞP) = ', dP_actual_bar, ' bar'
write(*,'(A,F10.2,A)') 'PEAK PIPELINE PRESSURE (Pmax) = ', P_peak_bar, ' bar'
write(*,'(A,F10.1,A)') 'Pipe Wall Hoop Stress = ', hoop_stress_MPa, ' MPa'
write(*,'(A,A)') 'SAFETY INTEGRITY STATUS = ', trim(safety_status)
write(*,'(A)') '============================================================'
end program water_hammer_allievi_surge
π» How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 water_hammer_allievi_surge.f90 -o water_hammer_allievi_surge
2. Execution with input.txt redirection:
water_hammer_allievi_surge < input.txt
π Sample input.txt File Structure
Sample Data:
1 1 2000.0 400.0 8.0 800.0 6.0 3.0
Parameter Description:
Pipe Material (1=Carbon Steel, 2=Ductile Iron, 3=HDPE, 4=PVC)\nFluid Type (1=Water, 2=Crude Oil, 3=Ethanol)\nPipeline Length L [m]\nPipe Inner Diameter [mm]\nWall Thickness [mm]\nFlow Rate [mΒ³/h]\nStatic Operating Pressure [bar]\nValve Closure Time [s]