π» Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Water Hammer Calculator
Core Numerical Engine in Fortran 90 β’ 49 total downloads
water_hammer.f90
! =========================================================================
! Source File: water_hammer.f90
! =========================================================================
! ============================================================================
! ThermoFluidCalc β Water Hammer (Coup de bΓ©lier) Solver
! Reference: Wylie & Streeter, Fluid Transients in Systems; Joukowsky (1898)
! ============================================================================
program water_hammer
implicit none
! Inputs
double precision :: L ! Pipe length [m]
double precision :: D ! Pipe inner diameter [mm]
double precision :: e_thick ! Pipe wall thickness [mm]
double precision :: E_pipe ! Young's Modulus of pipe material [GPa]
double precision :: V0 ! Initial flow velocity [m/s]
double precision :: Kf ! Bulk modulus of fluid [GPa]
double precision :: rho ! Fluid density [kg/mΒ³]
double precision :: tc ! Valve closure time [s]
integer :: profile_option ! 1 = Instantaneous, 2 = Linear, 3 = Parabolic
double precision :: P0 ! Static pressure [kPa]
double precision :: Pv ! Vapor pressure [kPa]
! Constants
double precision, parameter :: g = 9.81d0
double precision, parameter :: pi = 3.141592653589793d0
! Outputs / Derived values
double precision :: a ! Wave speed [m/s]
double precision :: t_crit ! Critical time 2L/a [s]
double precision :: DP_inst ! Theoretical instantaneous pressure surge [kPa]
double precision :: DP_grad ! Theoretical gradual pressure surge (rigid column) [kPa]
double precision :: DP_max ! Absolute maximum pressure surge from simulation [kPa]
double precision :: P_max ! Maximum absolute pressure [kPa]
double precision :: P_min ! Minimum absolute pressure [kPa]
double precision :: F_support ! Force on pipe supports [kN]
double precision :: hoop_stress ! Stress in pipe wall [MPa]
logical :: cavitation_occurred
! MOC Solver variables
integer, parameter :: N = 20 ! Number of spatial reaches
double precision :: dx ! Spatial step [m]
double precision :: dt ! Time step [s]
integer :: N_steps ! Total time steps
double precision :: t_max ! Simulation duration [s]
double precision :: f_fric ! Darcy friction factor
double precision :: Re ! Reynolds number
double precision :: mu_visc ! Dynamic viscosity [Pa-s]
double precision :: A_area ! Pipe cross-sectional area [mΒ²]
double precision :: Q0 ! Initial volumetric flow rate [mΒ³/s]
double precision :: B ! Impedance parameter
double precision :: R_fric ! Friction resistance factor
! Grid arrays
double precision :: H(0:N) ! Hydraulic head at current time step [m]
double precision :: Q(0:N) ! Volumetric flow rate at current time step [mΒ³/s]
double precision :: H_new(0:N)
double precision :: Q_new(0:N)
! Time series tracking at the valve (node N)
integer, parameter :: MAX_TS_POINTS = 500
double precision :: ts_time(MAX_TS_POINTS)
double precision :: ts_pres(MAX_TS_POINTS)
integer :: ts_count, skip_step
! Temporary variables
integer :: i, nt
double precision :: t_val, tau, CP, CM, Cv, Q_N_temp, H_N_temp, Hv, H0_head, term, eps_over_D
character(len=20) :: closure_type
! Read inputs from stdin
read(*,*) L
read(*,*) D
read(*,*) e_thick
read(*,*) E_pipe
read(*,*) V0
read(*,*) Kf
read(*,*) rho
read(*,*) tc
read(*,*) profile_option
read(*,*) P0
read(*,*) Pv
! Defaults & safety bounds
if (L <= 0.0d0) L = 100.0d0
if (D <= 0.0d0) D = 100.0d0
if (e_thick <= 0.0d0) e_thick = 5.0d0
if (E_pipe <= 0.0d0) E_pipe = 200.0d0
if (Kf <= 0.0d0) Kf = 2.2d0
if (rho <= 0.0d0) rho = 1000.0d0
if (P0 <= 0.0d0) P0 = 200.0d0
if (Pv < 0.0d0) Pv = 2.34d0
if (tc < 0.0d0) tc = 0.0d0
! ββ 1. WAVE SPEED CALCULATIONS (Joukowsky) βββββββββββββββββ
! bulk modulus and Young's modulus converted to Pa (1e9)
! a = sqrt(Kf/rho) / sqrt(1 + (Kf * D)/(E * e))
a = sqrt((Kf * 1.0d9) / rho) / sqrt(1.0d0 + (Kf * D) / (E_pipe * e_thick))
! Critical time t_crit = 2L/a
t_crit = 2.0d0 * L / a
! Closure type
if (tc < t_crit) then
closure_type = "Instantaneous"
else
closure_type = "Gradual"
end if
! Theoretical surges
DP_inst = (rho * a * V0) / 1000.0d0 ! in kPa
if (tc > 0.0d0) then
DP_grad = (rho * L * V0 / tc) / 1000.0d0 ! rigid column in kPa
else
DP_grad = DP_inst
end if
! ββ 2. PREPARE METHOD OF CHARACTERISTICS (MOC) βββββββββββββ
dx = L / dble(N)
dt = dx / a
! Simulate for 5 wave cycles: T = 4L/a = 2 * t_crit
t_max = 5.0d0 * (4.0d0 * L / a)
N_steps = nint(t_max / dt)
if (N_steps < 100) N_steps = 100
! Pipe area
A_area = pi * (D / 1000.0d0)**2 / 4.0d0
Q0 = V0 * A_area
! Friction factor estimation (Haaland)
mu_visc = 0.001d0 ! water viscosity fallback
Re = rho * abs(V0) * (D / 1000.0d0) / mu_visc
if (Re < 2300.0d0) then
if (Re > 0.0d0) then
f_fric = 64.0d0 / Re
else
f_fric = 0.02d0
end if
else
eps_over_D = (0.05d0 / 1000.0d0) / (D / 1000.0d0)
term = (eps_over_D / 3.7d0)**1.11d0 + 6.9d0 / Re
f_fric = 1.0d0 / (-1.8d0 * log10(term))**2
end if
! MOC Parameters
B = a / (g * A_area)
R_fric = f_fric * dx / (2.0d0 * g * (D / 1000.0d0) * A_area**2)
! Initial Steady State Grade Line
! Reservoir at node 0 (fixed static pressure P0)
H0_head = P0 / (rho * g / 1000.0d0)
Hv = Pv / (rho * g / 1000.0d0)
do i = 0, N
H(i) = H0_head - dble(i) * R_fric * Q0**2
Q(i) = Q0
end do
! Initialize tracking
DP_max = 0.0d0
P_max = P0
P_min = P0
cavitation_occurred = .false.
ts_count = 0
skip_step = max(1, N_steps / MAX_TS_POINTS)
! Store initial point
ts_count = ts_count + 1
ts_time(ts_count) = 0.0d0
ts_pres(ts_count) = H(N) * (rho * g / 1000.0d0)
! ββ 3. MOC TRANSIENT LOOP ββββββββββββββββββββββββββββββββββ
do nt = 1, N_steps
t_val = dble(nt) * dt
! Determine valve opening tau
if (profile_option == 1) then
! Instantaneous
tau = 0.0d0
elseif (profile_option == 2) then
! Linear
if (t_val < tc) then
tau = 1.0d0 - (t_val / tc)
else
tau = 0.0d0
end if
else
! Parabolic
if (t_val < tc) then
tau = (1.0d0 - (t_val / tc))**2
else
tau = 0.0d0
end if
end if
! Interior nodes i = 1 ... N-1
do i = 1, N-1
CP = H(i-1) + B * Q(i-1) - R_fric * Q(i-1) * abs(Q(i-1))
CM = H(i+1) - B * Q(i+1) + R_fric * Q(i+1) * abs(Q(i+1))
H_new(i) = (CP + CM) / 2.0d0
Q_new(i) = (CP - CM) / (2.0d0 * B)
end do
! Reservoir Node i = 0 (Constant head)
H_new(0) = H0_head
CM = H(1) - B * Q(1) + R_fric * Q(1) * abs(Q(1))
Q_new(0) = (H_new(0) - CM) / B
! Valve Node i = N
CP = H(N-1) + B * Q(N-1) - R_fric * Q(N-1) * abs(Q(N-1))
Cv = (tau * Q0)**2 / H0_head
if (Cv > 0.0d0) then
! Quadratic formula for valve flow Q_N
Q_N_temp = (-Cv * B + sqrt((Cv * B)**2 + 4.0d0 * Cv * CP)) / 2.0d0
H_N_temp = CP - B * Q_N_temp
else
Q_N_temp = 0.0d0
H_N_temp = CP
end if
! Cavitation check at valve
if (H_N_temp < Hv) then
H_N_temp = Hv
if (tau > 0.0d0) then
Q_N_temp = (CP - Hv) / B
else
Q_N_temp = 0.0d0
end if
cavitation_occurred = .true.
end if
H_new(N) = H_N_temp
Q_new(N) = Q_N_temp
! Update grid values
H = H_new
Q = Q_new
! Track absolute pressure at the valve
H_N_temp = H(N)
H_N_temp = max(H_N_temp, Hv) ! Physical bound
! Convert to pressure
t_val = H_N_temp * (rho * g / 1000.0d0)
if (t_val > P_max) P_max = t_val
if (t_val < P_min) P_min = t_val
! Time series store
if (mod(nt, skip_step) == 0 .and. ts_count < MAX_TS_POINTS) then
ts_count = ts_count + 1
ts_time(ts_count) = dble(nt) * dt
ts_pres(ts_count) = t_val
end if
end do
! Ensure last point is captured
if (ts_count < MAX_TS_POINTS) then
ts_count = ts_count + 1
ts_time(ts_count) = dble(N_steps) * dt
ts_pres(ts_count) = H(N) * (rho * g / 1000.0d0)
end if
! Calculate pressure surge surge
DP_max = P_max - P0
if (DP_max < 0.0d0) DP_max = 0.0d0
! Force on pipe supports (F = A * DP_max)
! DP_max is in kPa, A_area in mΒ² -> F in kN
F_support = A_area * DP_max
! Stress in pipe wall (hoop stress)
! sigma = P * D / (2 * e_thick)
! P_max is in kPa -> P_max / 1000 is in MPa
hoop_stress = (P_max / 1000.0d0) * D / (2.0d0 * e_thick)
! ββ 4. OUTPUT RESULTS IN KEY-VALUE FORMAT βββββββββββββββββ
write(*, '(A, F14.2)') "Wave Speed = ", a
write(*, '(A, F14.4)') "Critical Time = ", t_crit
write(*, '(A, A)') "Closure Type = ", trim(closure_type)
write(*, '(A, F14.2)') "DP Inst = ", DP_inst
write(*, '(A, F14.2)') "DP Grad = ", DP_grad
write(*, '(A, F14.2)') "Pressure Surge = ", DP_max
write(*, '(A, F14.2)') "Max Pressure = ", P_max
write(*, '(A, F14.2)') "Min Pressure = ", P_min
write(*, '(A, A)') "Cavitation Risk = ", merge("Yes", "No ", cavitation_occurred)
write(*, '(A, F14.3)') "Support Force = ", F_support
write(*, '(A, F14.2)') "Hoop Stress = ", hoop_stress
! Timeline data points for PHP parsing
write(*, '(A)') "--- TIMELINE DATA ---"
do i = 1, ts_count
write(*, '(F10.5, A, F12.2)') ts_time(i), ",", ts_pres(i)
end do
end program water_hammer
Solver Description
Solves hydraulic transient pressure waves in pipes caused by valve closure using the 1D Method of Characteristics (MOC). Computes Joukowsky wave speeds, critical times, maximum surge pressures, hoop stresses, support forces, and flags cavitation risks.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
π οΈ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 water_hammer.f90 -o water_hammer_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
water_hammer_calc < input.txt
π₯ Downloads & Local Files
Preview of the required input file (input.txt):
! Pipe length [m]
150.0
! Pipe diameter [mm]
100.0
! Pipe wall thickness [mm]
5.0
! Pipe modulus [GPa]
200.0
! Flow velocity [m/s]
2.0
! Fluid bulk modulus [GPa]
2.2
! Fluid density [kg/m3]
1000.0
! Valve closure time [s]
0.5
! Valve closure profile (1=Inst, 2=Linear, 3=Parabolic)
2
! Static pressure [kPa]
300.0
! Vapor pressure [kPa]
2.34
150.0
! Pipe diameter [mm]
100.0
! Pipe wall thickness [mm]
5.0
! Pipe modulus [GPa]
200.0
! Flow velocity [m/s]
2.0
! Fluid bulk modulus [GPa]
2.2
! Fluid density [kg/m3]
1000.0
! Valve closure time [s]
0.5
! Valve closure profile (1=Inst, 2=Linear, 3=Parabolic)
2
! Static pressure [kPa]
300.0
! Vapor pressure [kPa]
2.34