🔬
Solver Purpose & Physical Scope
Design Pelton, Francis, and Kaplan turbines. Compute specific speed, runner diameter, efficiency envelopes, shaft power, and cavitation safety.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 380 times
📄 Source File:
hydraulic_turbine.f90
📁 calcul/FluidMechanics /
hydraulic_turbine.f90
program hydraulic_turbine
implicit none
integer :: turbType, i, iostat_val
double precision :: H, Q, N_rpm, rho, g, eta_user
double precision :: P_hydro, P_shaft, eta, Ns, Nq, D_runner, U, C1
double precision :: omega, phi, psi, sigma_cavit, NPSH_req
double precision :: Q_i, H_i, P_i, eta_i, Ns_i
double precision, parameter :: PI = 3.141592653589793d0
character(len=40) :: turbName, regime
read(*,*,iostat=iostat_val) turbType
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid turbine type input.'
stop
end if
read(*,*,iostat=iostat_val) H
read(*,*,iostat=iostat_val) Q
read(*,*,iostat=iostat_val) N_rpm
read(*,*,iostat=iostat_val) rho
read(*,*,iostat=iostat_val) g
read(*,*,iostat=iostat_val) eta_user
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all turbine inputs.'
stop
end if
if (H <= 0.0d0 .or. Q <= 0.0d0 .or. N_rpm <= 0.0d0) then
write(*,*) 'ERROR: Head, flow and speed must be positive.'
stop
end if
if (rho <= 0.0d0) rho = 998.0d0
if (g <= 0.0d0) g = 9.81d0
if (eta_user <= 0.0d0 .or. eta_user > 1.0d0) eta_user = 0.0d0
omega = 2.0d0 * PI * N_rpm / 60.0d0
P_hydro = rho * g * Q * H
! Specific speed (metric) Ns = N sqrt(P) / H^(5/4) [P in kW]
Ns = N_rpm * sqrt(P_hydro/1000.0d0) / H**1.25d0
! Specific speed Nq = N sqrt(Q) / H^(3/4)
Nq = N_rpm * sqrt(Q) / H**0.75d0
select case(turbType)
case(1)
turbName = 'Pelton (impulse)'
if (eta_user > 0.01d0) then
eta = eta_user
else
eta = pelton_eta(Ns)
end if
! Pelton: bucket speed ~ 0.45-0.48 * jet velocity
C1 = sqrt(2.0d0*g*H) ! jet velocity
U = 0.46d0 * C1 ! peripheral speed
D_runner = 2.0d0 * U / omega
sigma_cavit = 0.0d0
NPSH_req = 0.0d0
case(2)
turbName = 'Francis (reaction, radial)'
if (eta_user > 0.01d0) then
eta = eta_user
else
eta = francis_eta(Ns)
end if
C1 = sqrt(2.0d0*g*H)
phi = 0.70d0 ! flow coefficient
U = phi * C1
D_runner = 2.0d0 * U / omega
sigma_cavit = 0.0432d0 * (Nq/100.0d0)**2
NPSH_req = sigma_cavit * H
case(3)
turbName = 'Kaplan (reaction, axial)'
if (eta_user > 0.01d0) then
eta = eta_user
else
eta = kaplan_eta(Ns)
end if
C1 = sqrt(2.0d0*g*H)
phi = 0.80d0
U = phi * C1
D_runner = 2.0d0 * U / omega
sigma_cavit = 0.28d0 + 0.0024d0 * Nq
NPSH_req = sigma_cavit * H
case default
write(*,*) 'ERROR: Turbine type must be 1 Pelton, 2 Francis, or 3 Kaplan.'
stop
end select
P_shaft = eta * P_hydro
psi = 2.0d0 * g * H / U**2 ! head coefficient
! Regime classification based on Ns
if (Ns < 60.0d0) then
regime = 'Low Ns — Pelton range'
else if (Ns < 300.0d0) then
regime = 'Medium Ns — Francis range'
else
regime = 'High Ns — Kaplan/propeller range'
end if
write(*,'(A)') '============================================================'
write(*,'(A)') ' HYDRAULIC TURBINE ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Turbine Type = ', trim(turbName)
write(*,'(A,ES12.4,A)') ' Net Head = ', H, ' m'
write(*,'(A,ES12.4,A)') ' Volume Flow = ', Q, ' m3/s'
write(*,'(A,F12.2,A)') ' Rotational Speed = ', N_rpm, ' rpm'
write(*,'(A,ES12.4,A)') ' Fluid Density = ', rho, ' kg/m3'
write(*,*)
write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
write(*,'(A,ES12.4)') ' Specific Speed Ns (metric)= ', Ns
write(*,'(A,ES12.4)') ' Specific Speed Nq = ', Nq
write(*,'(A,A)') ' Ns Regime = ', trim(regime)
write(*,'(A,ES12.4)') ' Efficiency = ', eta
write(*,'(A,ES12.4,A)') ' Hydraulic Power = ', P_hydro, ' W'
write(*,'(A,ES12.4,A)') ' Shaft Power = ', P_shaft, ' W'
write(*,'(A,ES12.4,A)') ' Shaft Power = ', P_shaft/1000.0d0, ' kW'
write(*,'(A,ES12.4,A)') ' Runner Diameter = ', D_runner, ' m'
write(*,'(A,ES12.4,A)') ' Peripheral Speed = ', U, ' m/s'
write(*,'(A,ES12.4,A)') ' Jet/Inlet Velocity = ', C1, ' m/s'
write(*,'(A,ES12.4)') ' Head Coefficient psi = ', psi
write(*,'(A,ES12.4)') ' Cavitation sigma = ', sigma_cavit
write(*,'(A,ES12.4,A)') ' NPSH Required = ', NPSH_req, ' m'
write(*,*)
! Efficiency vs Ns sweep
write(*,'(A)') '--- EFFICIENCY VS NS SWEEP ----------------------------------'
write(*,'(A)') ' Ns eta_Pelton eta_Francis eta_Kaplan'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, 60
Ns_i = 5.0d0 + 15.0d0*dble(i-1)
write(*,'(F10.2,2X,F10.5,2X,F10.5,2X,F10.5)') &
Ns_i, pelton_eta(Ns_i), francis_eta(Ns_i), kaplan_eta(Ns_i)
end do
write(*,*)
! Power vs Q sweep at given H and N
write(*,'(A)') '--- POWER VS FLOW SWEEP -------------------------------------'
write(*,'(A)') ' Q[m3/s] P_shaft[kW] eta'
write(*,'(A)') ' -------------------------------------------'
do i = 1, 40
Q_i = Q * 0.1d0 * dble(i)
P_i = rho * g * Q_i * H
Ns_i = N_rpm * sqrt(P_i/1000.0d0) / H**1.25d0
select case(turbType)
case(1); eta_i = pelton_eta(Ns_i)
case(2); eta_i = francis_eta(Ns_i)
case(3); eta_i = kaplan_eta(Ns_i)
end select
if (eta_user > 0.01d0) eta_i = eta_user
write(*,'(ES12.4,2X,ES12.4,2X,F10.5)') Q_i, eta_i*P_i/1000.0d0, eta_i
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' P_hydro = rho g Q H.'
write(*,'(A)') ' Ns = N sqrt(P_kW) / H^1.25 (metric specific speed).'
write(*,'(A)') ' Nq = N sqrt(Q) / H^0.75.'
write(*,'(A)') ' Euler: U = phi sqrt(2gH); D = 2U/omega.'
write(*,'(A)') ' Efficiency from empirical Ns-based envelope curves.'
contains
double precision function pelton_eta(Ns_in)
implicit none
double precision, intent(in) :: Ns_in
double precision :: x
! Bell-shaped efficiency curve peaking around Ns ~ 15-25
x = (Ns_in - 18.0d0) / 15.0d0
pelton_eta = 0.90d0 * exp(-0.5d0 * x**2)
if (pelton_eta < 0.30d0) pelton_eta = 0.30d0
if (pelton_eta > 0.92d0) pelton_eta = 0.92d0
end function pelton_eta
double precision function francis_eta(Ns_in)
implicit none
double precision, intent(in) :: Ns_in
double precision :: x
! Peak around Ns ~ 120-200
x = (Ns_in - 160.0d0) / 100.0d0
francis_eta = 0.93d0 * exp(-0.5d0 * x**2)
if (francis_eta < 0.30d0) francis_eta = 0.30d0
if (francis_eta > 0.94d0) francis_eta = 0.94d0
end function francis_eta
double precision function kaplan_eta(Ns_in)
implicit none
double precision, intent(in) :: Ns_in
double precision :: x
! Peak around Ns ~ 500-700
x = (Ns_in - 600.0d0) / 250.0d0
kaplan_eta = 0.92d0 * exp(-0.5d0 * x**2)
if (kaplan_eta < 0.30d0) kaplan_eta = 0.30d0
if (kaplan_eta > 0.93d0) kaplan_eta = 0.93d0
end function kaplan_eta
end program hydraulic_turbine
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 hydraulic_turbine.f90 -o hydraulic_turbine
2. Execution with input.txt redirection:
hydraulic_turbine < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
tt_init\nNet head H\nVolume flow Q [m³/s]\nRotational speed N [rpm]\nWater density Ï [kg/m³]\nGravity g [m/s²]\neu_init