π» Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Flow Around Cylinders
Core Numerical Engine in Fortran 90 β’ 43 total downloads
flow_cylinder.f90
! =========================================================================
! Source File: 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
Solver Description
Compute vortex shedding frequency, Strouhal number, drag/lift coefficients, and lock-in velocity range for circular cylinders in cross-flow.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
π οΈ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 flow_cylinder.f90 -o flow_cylinder
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
flow_cylinder < input.txt
π₯ Downloads & Local Files
Preview of the required input file (input.txt):
! 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]
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0