Francis & Kaplan Turbine Cavitation Sizer (Thoma σ)

Standalone, self-contained numerical routine. Verify algorithms, inspect boundary condition equations, or compile locally for batch parametric runs.

📦 Download Complete ZIP 📥 Download .f90 Code 📄 Download input.txt

🔬

Solver Purpose & Physical Scope

Calculate Francis and Kaplan reaction turbine cavitation limits: Thoma critical cavitation coefficient sigma, metric specific speed nq, and maximum setting height above tailwater.

📂 Discipline: Turbomachinery ⚡ Precision: IEEE-754 64-bit Real(`real(8)`) 📥 Total Downloads: 356 times 📄 Source File: francis_kaplan_turbine_cavitation_thoma.f90

📁 calcul/Turbomachinery / francis_kaplan_turbine_cavitation_thoma.f90

program francis_kaplan_turbine_cavitation_thoma
    implicit none
    integer :: iostat_val, turbine_select
    double precision :: H_net_m, Q_m3s, N_rpm, alt_m, Tw_C
    double precision :: P_atm_Pa, Pv_Pa, nq_metric, Ns_power, P_hyd_MW, P_shaft_MW
    double precision :: sigma_c, H_atm_m, H_vap_m, zs_max_m, cavitation_margin_m
    character(len=32) :: turbine_rec, cav_status
    double precision, parameter :: G_ACC = 9.80665d0
    double precision, parameter :: RHO_W = 1000.0d0

    ! Read inputs
    read(*,*,iostat=iostat_val) turbine_select  ! 1=Francis Turbine, 2=Kaplan / Bulb Turbine
    read(*,*,iostat=iostat_val) H_net_m         ! Net Available Water Head [m] (e.g. 85.0)
    read(*,*,iostat=iostat_val) Q_m3s           ! Water Flow Rate Q [m3/s] (e.g. 24.0)
    read(*,*,iostat=iostat_val) N_rpm           ! Rotational Speed N [rpm] (e.g. 375.0)
    read(*,*,iostat=iostat_val) alt_m           ! Powerhouse Altitude above sea level [m] (e.g. 450.0)
    read(*,*,iostat=iostat_val) Tw_C            ! River Water Temperature [deg C] (e.g. 15.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for Francis/Kaplan cavitation calculation.'
        stop
    end if

    if (H_net_m <= 0.0d0 .or. Q_m3s <= 0.0d0 .or. N_rpm <= 0.0d0) then
        write(*,*) 'ERROR: Head, flow rate, and speed must be positive.'
        stop
    end if

    ! Barometric Atmospheric Pressure per standard atmosphere
    P_atm_Pa = 101325.0d0 * ((1.0d0 - 2.25577d-5 * alt_m)**5.25588d0)
    
    ! Vapor pressure of water
    Pv_Pa = 611.2d0 * exp((17.67d0 * Tw_C) / (Tw_C + 243.5d0))

    H_atm_m = P_atm_Pa / (RHO_W * G_ACC)
    H_vap_m = Pv_Pa / (RHO_W * G_ACC)

    ! Specific speed Metric nq = N * sqrt(Q) / H^0.75
    nq_metric = N_rpm * sqrt(Q_m3s) / (H_net_m**0.75d0)

    P_hyd_MW = (RHO_W * G_ACC * Q_m3s * H_net_m) / 1.0d6
    P_shaft_MW = P_hyd_MW * 0.93d0 ! 93% average reaction turbine efficiency
    Ns_power = N_rpm * sqrt(P_shaft_MW * 1000.0d0) / (H_net_m**1.25d0)

    if (nq_metric < 35.0d0) then
        turbine_rec = 'PELTON IMPULSE RECOMMENDED'
    else if (nq_metric <= 110.0d0) then
        turbine_rec = 'HIGH-HEAD FRANCIS RUNNER'
    else if (nq_metric <= 260.0d0) then
        turbine_rec = 'MEDIUM/LOW-HEAD FRANCIS'
    else
        turbine_rec = 'KAPLAN / BULB AXIAL RUNNER'
    end if

    ! Thoma Critical Cavitation Coefficient (USBR / IEC standards)
    if (turbine_select == 1) then ! Francis
        sigma_c = 0.0437d0 * ((nq_metric / 100.0d0)**1.64d0) + 0.015d0
    else ! Kaplan
        sigma_c = 0.105d0 * ((nq_metric / 100.0d0)**1.75d0) + 0.080d0
    end if

    ! Maximum permissible suction setting height zs [m]
    zs_max_m = (H_atm_m - H_vap_m) - sigma_c * H_net_m - 0.50d0 ! 0.5m safety margin

    if (zs_max_m >= 0.0d0) then
        cav_status = 'ABOVE TAILWATER (POSITIVE ZS)'
    else
        cav_status = 'SUBMERGED RUNNER (NEGATIVE ZS)'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — HYDRAULIC TURBINE CAVITATION (THOMA σ)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Specific Speed (nq / Ns)     = ', nq_metric, ' rpm / ', Ns_power, ''
    write(*,'(A,A)')               'Turbine Selection Topology   = ', trim(turbine_rec)
    write(*,'(A,F10.4,A)')         'Thoma Critical Cavitation σc = ', sigma_c, ''
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'MAX RUNNER SETTING HEIGHT zs = ', zs_max_m, ' m'
    write(*,'(A,A)')        'Powerhouse Excavation Regime = ', trim(cav_status)
    write(*,'(A,F10.2,A)')  'Barometric Pressure Head     = ', H_atm_m, ' m w.c.'
    write(*,'(A,F10.2,A)')  'Cavitation Head Drop (σc·H)  = ', sigma_c * H_net_m, ' m'
    write(*,'(A,F10.2,A)')  'Shaft Power Generated        = ', P_shaft_MW, ' MW'
    write(*,'(A)') '============================================================'

end program francis_kaplan_turbine_cavitation_thoma

💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 francis_kaplan_turbine_cavitation_thoma.f90 -o francis_kaplan_turbine_cavitation_thoma

2. Execution with input.txt redirection:

francis_kaplan_turbine_cavitation_thoma < input.txt

📄 Sample input.txt File Structure

Sample Data:

Parameter Description:

Turbine Type (1=Francis, 2=Kaplan)
Net Head Hnet [m]
Water Flow Q [m³/s]
Rotational Speed N [rpm]
Altitude [m]
Water Temp [°C]