🔬
Solver Purpose & Physical Scope
A body moving through a fluid experiences a drag force opposing its motion, defined by the drag equation:
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 489 times
📄 Source File:
drag_force_calculator.f90
📁 calcul/FluidMechanics /
drag_force_calculator.f90
program drag_force_calculator
implicit none
double precision :: U_inf, rho, mu, nu, A_ref, D_char
double precision :: Re, Cd, F_drag, q_dynamic
double precision :: pi
integer :: shape ! 1=sphere, 2=cylinder, 3=flat-plate-normal, 4=streamlined
character(len=30) :: shape_name
integer :: i
double precision :: Re_i, Cd_i
pi = 3.14159265358979d0
! Read inputs
read(*,*) shape ! Shape type
read(*,*) D_char ! Characteristic dimension [m] (diameter or length)
read(*,*) U_inf ! Free stream velocity [m/s]
read(*,*) rho ! Fluid density [kg/m3]
read(*,*) mu ! Dynamic viscosity [Pa.s]
nu = mu / rho
Re = rho * U_inf * D_char / mu
q_dynamic = 0.5d0 * rho * U_inf**2
! ============================================
! DRAG COEFFICIENT CORRELATIONS
! ============================================
select case(shape)
case(1)
shape_name = 'Sphere'
A_ref = pi * D_char**2 / 4.0d0
if (Re < 1.0d0) then
! Stokes flow
Cd = 24.0d0 / Re
else if (Re < 1000.0d0) then
! Schiller-Naumann
Cd = 24.0d0/Re * (1.0d0 + 0.15d0 * Re**0.687d0)
else if (Re < 2.0d5) then
! Standard drag
Cd = 0.44d0
else
! Post-critical (turbulent boundary layer)
Cd = 0.1d0
end if
case(2)
shape_name = 'Infinite Cylinder'
A_ref = D_char * 1.0d0 ! Per unit length
if (Re < 1.0d0) then
Cd = 8.0d0 * pi / (Re * (2.002d0 - log(Re)))
else if (Re < 1000.0d0) then
Cd = 1.0d0 + 10.0d0 * Re**(-2.0d0/3.0d0)
else if (Re < 2.0d5) then
Cd = 1.2d0
else
Cd = 0.3d0
end if
case(3)
shape_name = 'Flat Plate (normal)'
A_ref = D_char * D_char ! Square plate
Cd = 1.17d0 ! Approximately constant for Re > 1000
if (Re < 100.0d0) then
Cd = 2.0d0 ! Low Re approximation
end if
case(4)
shape_name = 'Streamlined Body'
A_ref = pi * (D_char/2.0d0)**2 ! Frontal area
Cd = 0.04d0 ! Well-streamlined body
if (Re < 1.0d5) then
Cd = 0.1d0
end if
case default
shape_name = 'Sphere'
A_ref = pi * D_char**2 / 4.0d0
Cd = 0.44d0
end select
! Drag force
F_drag = Cd * q_dynamic * A_ref
! ============================================
! OUTPUT
! ============================================
write(*,'(A)') '============================================================'
write(*,'(A)') ' EXTERNAL FLOW — DRAG FORCE CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- CONFIGURATION ---'
write(*,'(A,A)') ' Body Shape = ', trim(shape_name)
write(*,'(A,F12.4,A)') ' Characteristic Dim (D) = ', D_char, ' m'
write(*,'(A,ES12.4,A)') ' Reference Area = ', A_ref, ' m2'
write(*,*)
write(*,'(A)') '--- FLOW CONDITIONS ---'
write(*,'(A,F12.4,A)') ' Free Stream Velocity = ', U_inf, ' m/s'
write(*,'(A,F12.2,A)') ' Fluid Density = ', rho, ' kg/m3'
write(*,'(A,ES12.4,A)') ' Dynamic Viscosity = ', mu, ' Pa.s'
write(*,'(A,ES12.4)') ' Reynolds Number (Re) = ', Re
write(*,'(A,F12.2,A)') ' Dynamic Pressure (q) = ', q_dynamic, ' Pa'
write(*,*)
write(*,'(A)') '--- DRAG RESULTS ---'
write(*,'(A,F12.6)') ' Drag Coefficient (Cd) = ', Cd
write(*,'(A,F12.4,A)') ' Drag Force (F_D) = ', F_drag, ' N'
write(*,'(A,F12.4,A)') ' Drag Power (P_D) = ', F_drag * U_inf, ' W'
write(*,*)
! Terminal velocity for sphere
if (shape == 1) then
write(*,'(A)') '--- TERMINAL VELOCITY (sphere in gravity) ---'
! mg = F_drag -> (4/3)pi(D/2)3rho_s g = Cd x 0.5 x rho x V2 x pi(D/2)2
! V_t = sqrt(4rho_s g D / (3 Cd rho))
! Using rho_s = 2700 kg/m3 (aluminum) as example
write(*,'(A,F12.4,A)') ' V_terminal (steel) = ', &
sqrt(4.0d0 * 7800.0d0 * 9.81d0 * D_char / (3.0d0 * Cd * rho)), ' m/s'
write(*,'(A,F12.4,A)') ' V_terminal (aluminum) = ', &
sqrt(4.0d0 * 2700.0d0 * 9.81d0 * D_char / (3.0d0 * Cd * rho)), ' m/s'
write(*,'(A,F12.4,A)') ' V_terminal (wood) = ', &
sqrt(4.0d0 * 600.0d0 * 9.81d0 * D_char / (3.0d0 * Cd * rho)), ' m/s'
write(*,*)
end if
! Cd vs Re chart data
write(*,'(A)') '--- Cd vs Re PROFILE ---'
write(*,'(A)') ' Re Cd_sphere Cd_cylinder Cd_flat'
write(*,'(A)') ' ---'
do i = 1, 30
Re_i = 10.0d0**(-1.0d0 + dble(i-1) * 7.0d0 / 29.0d0)
! Sphere Cd
if (Re_i < 1.0d0) then
Cd_i = 24.0d0 / Re_i
else if (Re_i < 1000.0d0) then
Cd_i = 24.0d0/Re_i * (1.0d0 + 0.15d0 * Re_i**0.687d0)
else if (Re_i < 2.0d5) then
Cd_i = 0.44d0
else
Cd_i = 0.1d0
end if
write(*,'(ES12.4,2X,F12.6,2X,F12.6,2X,F12.6)') Re_i, &
Cd_i, &
merge(1.2d0, 1.0d0 + 10.0d0*Re_i**(-2.0d0/3.0d0), Re_i > 1000.0d0), &
merge(1.17d0, 2.0d0, Re_i > 100.0d0)
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---'
write(*,'(A)') ' F_D = Cd x (rhoU2/2) x A_ref'
write(*,'(A)') ' Sphere: Cd varies with Re (Stokes -> Schiller-Naumann)'
write(*,'(A)') ' Stokes (Re<1): Cd = 24/Re'
write(*,'(A)') ' Schiller-Naumann: Cd = 24/Re(1 + 0.15 Re^0.687)'
end program drag_force_calculator
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 drag_force_calculator.f90 -o drag_calc
2. Execution with input.txt redirection:
drag_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.1 10.0 1.225 1.81e-5
Parameter Description:
Object shape (1=Sphere, 2=Cylinder, 3=Flat Disk, 4=Streamlined body) Characteristic dimension (Diameter/Length) [m] Flow velocity [m/s] Fluid density [kg/m3] Fluid dynamic viscosity [Pa-s]