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