Axial Turbine Velocity Triangles & Efficiency (Soderberg)

program axial_turbine_velocity_triangles_stage
    implicit none
    integer :: iostat_val, turbine_type
    double precision :: m_dot_kgs, T01_K, P01_bar, r_mean_m, N_rpm, V_axial_ms
    double precision :: alpha2_deg, beta3_deg, tip_clearance_pct
    double precision :: cp_gas, gamma_k, R_gas
    double precision :: U_blade_ms, phi_flow, alpha2_rad, beta3_rad, V_theta2, W_theta2, beta2_rad, beta2_deg
    double precision :: W_theta3, V_theta3, alpha3_rad, alpha3_deg, V2_mag, W2_mag, W3_mag, V3_mag
    double precision :: delta_h0_Jkg, delta_h0_kJkg, psi_loading, degree_reaction, eta_tt, eta_ts, P_stage_MW
    character(len=32) :: reaction_type
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) turbine_type     ! 1=Gas Turbine HP Stage, 2=Steam Turbine Impulse Stage, 3=Turbine 50% Reaction, 4=Micro Gas Turbine
    read(*,*,iostat=iostat_val) m_dot_kgs        ! Mass Flow Rate [kg/s] (e.g. 45.0)
    read(*,*,iostat=iostat_val) T01_K            ! Total Inlet Temperature [K] (e.g. 1450.0)
    read(*,*,iostat=iostat_val) P01_bar          ! Total Inlet Pressure [bar a] (e.g. 18.0)
    read(*,*,iostat=iostat_val) r_mean_m         ! Mean Stage Radius [m] (e.g. 0.42)
    read(*,*,iostat=iostat_val) N_rpm            ! Shaft Rotational Speed [rpm] (e.g. 9500.0)
    read(*,*,iostat=iostat_val) V_axial_ms       ! Axial Velocity Va [m/s] (e.g. 180.0)
    read(*,*,iostat=iostat_val) alpha2_deg       ! Stator Exit Angle α2 [deg] (e.g. 68.0)
    read(*,*,iostat=iostat_val) beta3_deg        ! Rotor Exit Angle β3 [deg] (e.g. -62.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for axial turbine velocity triangles.'
        stop
    end if

    if (m_dot_kgs <= 0.0d0 .or. r_mean_m <= 0.0d0 .or. N_rpm <= 0.0d0 .or. V_axial_ms <= 0.0d0) then
        write(*,*) 'ERROR: Flow rate, radius, speed, and axial velocity must be positive.'
        stop
    end if

    ! Gas properties (Combustion gas vs Steam)
    if (turbine_type == 2) then ! Steam
        cp_gas = 2100.0d0; gamma_k = 1.30d0; R_gas = 461.5d0
    else ! Gas Turbine Combustion Products
        cp_gas = 1150.0d0; gamma_k = 1.33d0; R_gas = 287.0d0
    end if

    U_blade_ms = 2.0d0 * PI * (N_rpm / 60.0d0) * r_mean_m
    phi_flow = V_axial_ms / max(1.0d0, U_blade_ms)

    alpha2_rad = alpha2_deg * PI / 180.0d0
    beta3_rad = beta3_deg * PI / 180.0d0

    ! Velocity Triangles Computation
    V_theta2 = V_axial_ms * tan(alpha2_rad)
    W_theta2 = V_theta2 - U_blade_ms
    beta2_rad = atan2(W_theta2, V_axial_ms)
    beta2_deg = beta2_rad * 180.0d0 / PI

    W_theta3 = V_axial_ms * tan(beta3_rad)
    V_theta3 = W_theta3 + U_blade_ms
    alpha3_rad = atan2(V_theta3, V_axial_ms)
    alpha3_deg = alpha3_rad * 180.0d0 / PI

    V2_mag = sqrt(V_axial_ms**2 + V_theta2**2)
    W2_mag = sqrt(V_axial_ms**2 + W_theta2**2)
    W3_mag = sqrt(V_axial_ms**2 + W_theta3**2)
    V3_mag = sqrt(V_axial_ms**2 + V_theta3**2)

    ! Euler Specific Work & Loading
    delta_h0_Jkg = U_blade_ms * (V_theta2 - V_theta3)
    delta_h0_kJkg = delta_h0_Jkg / 1000.0d0
    psi_loading = delta_h0_Jkg / max(1.0d0, U_blade_ms**2)

    ! Degree of Reaction (R)
    degree_reaction = 1.0d0 - ((V_theta2 + V_theta3) / (2.0d0 * max(1.0d0, U_blade_ms)))

    if (abs(degree_reaction) < 0.15d0) then
        reaction_type = 'IMPULSE STAGE (R ≈ 0)'
    else if (abs(degree_reaction - 0.50d0) < 0.15d0) then
        reaction_type = '50% REACTION (PARSONS STAGE)'
    else
        reaction_type = 'HIGH REACTION STAGE'
    end if

    ! Soderberg / Smith Chart Efficiency Estimation
    eta_tt = 1.0d0 / (1.0d0 + (0.045d0 * (V2_mag**2) + 0.055d0 * (W3_mag**2)) / (2.0d0 * max(100.0d0, delta_h0_Jkg)))
    eta_tt = min(0.935d0, max(0.60d0, eta_tt))
    eta_ts = delta_h0_Jkg / ((delta_h0_Jkg / eta_tt) + 0.50d0 * (V3_mag**2))
    eta_ts = min(eta_tt, max(0.50d0, eta_ts))

    P_stage_MW = (m_dot_kgs * delta_h0_kJkg) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC -- AXIAL TURBINE STAGE VELOCITY TRIANGLES'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Blade Speed U / Axial Va     = ', U_blade_ms, ' m/s / ', V_axial_ms, ' m/s'
    write(*,'(A,F10.3,A,F10.3,A)') 'Flow Coeff (phi) / Loading   = ', phi_flow, ' / ', psi_loading, ''
    write(*,'(A,F10.2,A,A)')       'Degree of Reaction (R)       = ', degree_reaction, ' | ', trim(reaction_type)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'SPECIFIC WORK (dh0)          = ', delta_h0_kJkg, ' kJ/kg'
    write(*,'(A,F10.2,A)')  'STAGE POWER OUTPUT           = ', P_stage_MW, ' MW'
    write(*,'(A,F10.2,A)')  'Efficiency (Total-to-Total)  = ', eta_tt*100.0d0, ' %'
    write(*,'(A,F10.2,A)')  'Efficiency (Total-to-Static) = ', eta_ts*100.0d0, ' %'
    write(*,'(A,F10.1,A,F10.1,A)') 'Rotor Rel Flow Angles b2/b3  = ', beta2_deg, ' deg / ', beta3_deg, ' deg'
    write(*,'(A,F10.1,A)')  'Absolute Exit Swirl Angle a3 = ', alpha3_deg, ' deg'
    write(*,'(A)') '============================================================'

end program axial_turbine_velocity_triangles_stage