🔬
Solver Purpose & Physical Scope
Compute vortex shedding frequency, Strouhal number, drag/lift coefficients, and lock-in velocity range for circular cylinders in cross-flow.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 436 times
📄 Source File:
flow_cylinder.f90
📁 calcul/FluidMechanics /
flow_cylinder.f90
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 flow_cylinder.f90 -o flow_cylinder
2. Execution with input.txt redirection:
flow_cylinder < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
Cylinder diameter D\nFreestream velocity V∞ [m/s]\nFluid density Ï [kg/m³]\nDynamic viscosity μ [Pa·s]\nNatural frequency fn [Hz]\nStructural damping ratio ζ\nMass per unit length [kg/m]