program flow_cylinder implicit none integer :: i, iostat_val double precision :: D, Vinf, rho, mu, Re, St, f_vortex, CD, CL_osc double precision :: F_drag, F_lift_amp, omega_shed, V_lockin_lo, V_lockin_hi double precision :: fn_struct, zeta_struct, m_per_L double precision :: Re_i, CD_i, St_i, V_i, f_i double precision, parameter :: PI = 3.141592653589793d0 character(len=60) :: regime read(*,*,iostat=iostat_val) D if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid cylinder diameter input.' stop end if read(*,*,iostat=iostat_val) Vinf read(*,*,iostat=iostat_val) rho read(*,*,iostat=iostat_val) mu read(*,*,iostat=iostat_val) fn_struct read(*,*,iostat=iostat_val) zeta_struct read(*,*,iostat=iostat_val) m_per_L if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all cylinder flow inputs.' stop end if if (D <= 0.0d0 .or. Vinf <= 0.0d0) then write(*,*) 'ERROR: Diameter and velocity must be positive.' stop end if if (rho <= 0.0d0 .or. mu <= 0.0d0) then write(*,*) 'ERROR: Density and viscosity must be positive.' stop end if if (fn_struct <= 0.0d0) fn_struct = 1.0d0 if (zeta_struct < 0.0d0) zeta_struct = 0.01d0 if (m_per_L <= 0.0d0) m_per_L = 1.0d0 Re = rho * Vinf * D / mu ! Strouhal number correlation St = strouhal_number(Re) ! Vortex shedding frequency f_vortex = St * Vinf / D omega_shed = 2.0d0 * PI * f_vortex ! Drag coefficient correlation CD = drag_coeff(Re) ! Oscillating lift coefficient amplitude (approximate) CL_osc = lift_osc_coeff(Re) ! Forces per unit length F_drag = 0.5d0 * rho * Vinf**2 * D * CD F_lift_amp = 0.5d0 * rho * Vinf**2 * D * CL_osc ! Lock-in velocity range (f_vortex ≈ fn_struct) ! Lock-in when 0.8 fn < f_shed < 1.2 fn => V = f_n D / St if (St > 1.0d-8) then V_lockin_lo = 0.8d0 * fn_struct * D / St V_lockin_hi = 1.2d0 * fn_struct * D / St else V_lockin_lo = 0.0d0 V_lockin_hi = 0.0d0 end if ! Regime classification call classify_regime(Re, regime) write(*,'(A)') '============================================================' write(*,'(A)') ' FLOW AROUND CYLINDERS ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A)') '--- INPUTS --------------------------------------------------' write(*,'(A,ES12.4,A)') ' Cylinder Diameter = ', D, ' m' write(*,'(A,ES12.4,A)') ' Freestream Velocity = ', Vinf, ' m/s' write(*,'(A,ES12.4,A)') ' Fluid Density = ', rho, ' kg/m3' write(*,'(A,ES12.4,A)') ' Dynamic Viscosity = ', mu, ' Pa.s' write(*,'(A,ES12.4)') ' Reynolds Number = ', Re write(*,'(A,ES12.4,A)') ' Structural fn = ', fn_struct, ' Hz' write(*,'(A,ES12.4)') ' Structural Damping zeta = ', zeta_struct write(*,'(A,ES12.4,A)') ' Mass per Length = ', m_per_L, ' kg/m' write(*,*) write(*,'(A)') '--- VORTEX SHEDDING RESULTS ---------------------------------' write(*,'(A,ES12.4)') ' Strouhal Number = ', St write(*,'(A,ES12.4,A)') ' Shedding Frequency = ', f_vortex, ' Hz' write(*,'(A,ES12.4,A)') ' Shedding Angular Freq = ', omega_shed, ' rad/s' write(*,'(A,A)') ' Reynolds Regime = ', trim(regime) write(*,*) write(*,'(A)') '--- FORCE COEFFICIENTS --------------------------------------' write(*,'(A,ES12.4)') ' Drag Coefficient CD = ', CD write(*,'(A,ES12.4)') ' Osc Lift Coefficient CL = ', CL_osc write(*,'(A,ES12.4,A)') ' Drag Force per Length = ', F_drag, ' N/m' write(*,'(A,ES12.4,A)') ' Lift Amplitude per Length = ', F_lift_amp, ' N/m' write(*,*) write(*,'(A)') '--- LOCK-IN ANALYSIS ----------------------------------------' write(*,'(A,ES12.4,A)') ' Lock-in Velocity Low = ', V_lockin_lo, ' m/s' write(*,'(A,ES12.4,A)') ' Lock-in Velocity High = ', V_lockin_hi, ' m/s' if (Vinf >= V_lockin_lo .and. Vinf <= V_lockin_hi) then write(*,'(A)') ' *** WARNING: CURRENT VELOCITY IS IN LOCK-IN RANGE ***' else write(*,'(A)') ' Current velocity is outside lock-in range.' end if write(*,*) ! CD, St vs Re sweep write(*,'(A)') '--- CD AND ST VS RE SWEEP -----------------------------------' write(*,'(A)') ' Re CD St f_shed[Hz]' write(*,'(A)') ' -----------------------------------------------------------' do i = 1, 80 Re_i = 0.1d0 * (1.0d7/0.1d0)**(dble(i-1)/79.0d0) CD_i = drag_coeff(Re_i) St_i = strouhal_number(Re_i) V_i = Re_i * mu / (rho * D) f_i = St_i * V_i / D write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') Re_i, CD_i, St_i, f_i end do write(*,*) ! Shedding frequency vs velocity sweep write(*,'(A)') '--- FREQUENCY VS VELOCITY SWEEP -----------------------------' write(*,'(A)') ' V[m/s] f_shed[Hz] f/fn regime' write(*,'(A)') ' -----------------------------------------------------------' do i = 1, 50 V_i = Vinf * 0.1d0 * dble(i) Re_i = rho * V_i * D / mu St_i = strouhal_number(Re_i) f_i = St_i * V_i / D if (f_i/fn_struct > 0.8d0 .and. f_i/fn_struct < 1.2d0) then write(*,'(ES12.4,2X,ES12.4,2X,F10.4,2X,A)') V_i, f_i, f_i/fn_struct, 'LOCK-IN' else write(*,'(ES12.4,2X,ES12.4,2X,F10.4,2X,A)') V_i, f_i, f_i/fn_struct, 'safe' end if end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' St = f D / V (Strouhal number).' write(*,'(A)') ' St ~ 0.198(1 - 19.7/Re) for subcritical Re.' write(*,'(A)') ' CD from piecewise empirical fit vs Re.' write(*,'(A)') ' Lock-in: 0.8 fn < f_shed < 1.2 fn.' contains double precision function strouhal_number(Re_in) implicit none double precision, intent(in) :: Re_in if (Re_in < 40.0d0) then strouhal_number = 0.0d0 ! no periodic shedding else if (Re_in < 300.0d0) then strouhal_number = 0.198d0 * (1.0d0 - 19.7d0/Re_in) else if (Re_in < 3.0d5) then strouhal_number = 0.20d0 ! subcritical else if (Re_in < 3.5d6) then strouhal_number = 0.30d0 ! critical/supercritical else strouhal_number = 0.27d0 ! transcritical end if if (strouhal_number < 0.0d0) strouhal_number = 0.0d0 end function strouhal_number double precision function drag_coeff(Re_in) implicit none double precision, intent(in) :: Re_in if (Re_in < 1.0d0) then ! Oseen-Stokes regime drag_coeff = 24.0d0 / max(Re_in, 1.0d-10) * & (1.0d0 + 3.0d0/16.0d0*Re_in) if (drag_coeff > 200.0d0) drag_coeff = 200.0d0 else if (Re_in < 1.0d3) then drag_coeff = 10.0d0 / sqrt(Re_in) else if (Re_in < 3.0d5) then drag_coeff = 1.2d0 ! subcritical else if (Re_in < 5.0d5) then ! drag crisis drag_coeff = 1.2d0 - 0.9d0*(Re_in-3.0d5)/2.0d5 else if (Re_in < 3.5d6) then drag_coeff = 0.3d0 ! supercritical else drag_coeff = 0.6d0 ! transcritical recovery end if if (drag_coeff < 0.15d0) drag_coeff = 0.15d0 end function drag_coeff double precision function lift_osc_coeff(Re_in) implicit none double precision, intent(in) :: Re_in if (Re_in < 40.0d0) then lift_osc_coeff = 0.0d0 else if (Re_in < 300.0d0) then lift_osc_coeff = 0.3d0 else if (Re_in < 3.0d5) then lift_osc_coeff = 0.5d0 ! subcritical else if (Re_in < 5.0d5) then lift_osc_coeff = 0.1d0 ! crisis suppression else lift_osc_coeff = 0.2d0 ! supercritical end if end function lift_osc_coeff subroutine classify_regime(Re_in, rname) implicit none double precision, intent(in) :: Re_in character(len=60), intent(out) :: rname if (Re_in < 5.0d0) then rname = 'Creeping flow — no separation' else if (Re_in < 40.0d0) then rname = 'Steady twin vortices — no shedding' else if (Re_in < 200.0d0) then rname = 'Laminar vortex street (2D periodic)' else if (Re_in < 300.0d0) then rname = 'Transition to 3D vortex shedding' else if (Re_in < 3.0d5) then rname = 'Subcritical — turbulent wake, laminar BL' else if (Re_in < 3.5d6) then rname = 'Critical / supercritical — drag crisis' else rname = 'Transcritical — turbulent BL and wake' end if end subroutine classify_regime end program flow_cylinder