von Kármán Vortex Shedding (Strouhal Frequency)

program karman_vortex_shedding
    implicit none
    integer :: iostat_val
    double precision :: D_cyl_mm, U_inf_ms, rho_kgm3, mu_cP
    double precision :: D_m, mu_Pas, Re_D, St_num, fs_Hz, period_s
    double precision :: CL_rms, FL_rms_Nm, q_dyn_Pa, vortex_spacing_a_mm, vortex_width_h_mm
    character(len=36) :: wake_regime

    ! Read inputs
    read(*,*,iostat=iostat_val) D_cyl_mm      ! Cylinder Diameter [mm] (e.g. 50.0)
    read(*,*,iostat=iostat_val) U_inf_ms      ! Free-Stream Velocity [m/s] (e.g. 12.0)
    read(*,*,iostat=iostat_val) rho_kgm3      ! Fluid Density [kg/m3] (e.g. 1.204)
    read(*,*,iostat=iostat_val) mu_cP         ! Fluid Viscosity [cP] (e.g. 0.0182)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for von Karman vortex calculation.'
        stop
    end if

    if (D_cyl_mm <= 0.0d0 .or. U_inf_ms <= 0.0d0 .or. rho_kgm3 <= 0.0d0 .or. mu_cP <= 0.0d0) then
        write(*,*) 'ERROR: All input parameters must be positive.'
        stop
    end if

    D_m = D_cyl_mm * 1.0d-3
    mu_Pas = mu_cP * 1.0d-3

    ! Reynolds number Re_D
    Re_D = (rho_kgm3 * U_inf_ms * D_m) / mu_Pas

    ! Strouhal Number Correlation (Fey, Konig & Eckelmann / Roshko)
    if (Re_D < 47.0d0) then
        St_num = 0.0d0
        wake_regime = 'Steady Symmetric Separation Bubble'
        CL_rms = 0.0d0
    else if (Re_D <= 180.0d0) then
        St_num = 0.212d0 * (1.0d0 - 21.2d0 / Re_D)
        wake_regime = 'Laminar 2D von Karman Street'
        CL_rms = 0.15d0 * sqrt(Re_D / 180.0d0)
    else if (Re_D <= 300.0d0) then
        St_num = 0.212d0 * (1.0d0 - 12.7d0 / Re_D)
        wake_regime = '3D Mode-A & Mode-B Transition'
        CL_rms = 0.25d0
    else if (Re_D <= 2.0d5) then
        St_num = 0.198d0 * (1.0d0 - 19.7d0 / Re_D)
        wake_regime = 'Subcritical Turbulent Wake Street'
        CL_rms = 0.35d0
    else if (Re_D <= 3.5d6) then
        St_num = 0.28d0 ! Supercritical narrow wake
        wake_regime = 'Supercritical (Disrupted Shedding)'
        CL_rms = 0.10d0
    else
        St_num = 0.29d0 ! Transcritical
        wake_regime = 'Transcritical Turbulent Vortex Street'
        CL_rms = 0.28d0
    end if

    if (St_num > 0.0d0) then
        fs_Hz = (St_num * U_inf_ms) / D_m
        period_s = 1.0d0 / fs_Hz
    else
        fs_Hz = 0.0d0
        period_s = 0.0d0
    end if

    q_dyn_Pa = 0.5d0 * rho_kgm3 * (U_inf_ms**2)
    FL_rms_Nm = q_dyn_Pa * D_m * CL_rms ! N/m span

    ! von Karman ideal stable vortex street geometry (h/a = 0.281)
    if (fs_Hz > 0.0d0) then
        vortex_spacing_a_mm = ((0.85d0 * U_inf_ms) / fs_Hz) * 1000.0d0
        vortex_width_h_mm = 0.281d0 * vortex_spacing_a_mm
    else
        vortex_spacing_a_mm = 0.0d0
        vortex_width_h_mm = 0.0d0
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — VON KARMAN VORTEX SHEDDING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Cylinder Diameter / U_inf = ', D_cyl_mm, ' mm / ', U_inf_ms, ' m/s'
    write(*,'(A,F10.2,A,F10.4,A)') 'Fluid Density / Viscosity = ', rho_kgm3, ' kg/m3 / ', mu_cP, ' cP'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F12.1)')   'REYNOLDS NUMBER (Re_D)    = ', Re_D
    write(*,'(A,F10.4)')   'STROUHAL NUMBER (St)      = ', St_num
    write(*,'(A,F10.2,A,F8.4,A)') 'SHEDDING FREQUENCY (f_s)  = ', fs_Hz, ' Hz (Period = ', period_s, ' s)'
    write(*,'(A,A)')       'Wake Hydrodynamic Regime  = ', trim(wake_regime)
    write(*,'(A,F10.3,A,F10.2,A)')'Fluctuating Lift CL_rms   = ', CL_rms, ' (Force = ', FL_rms_Nm, ' N/m)'
    write(*,'(A,F10.1,A,F10.1,A)')'Vortex Spacing a / Width h= ', vortex_spacing_a_mm, ' mm / ', vortex_width_h_mm, ' mm'
    write(*,'(A)') '============================================================'

end program karman_vortex_shedding