πŸ’» Fortran Source Code Library

Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.

External Flow Drag Force Solver

Core Numerical Engine in Fortran 90 β€’ 0 total downloads

drag_force_calculator.f90
! =========================================================================
! Source File: 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


Solver Description

Calculates external drag force ($F_d = C_d \cdot \frac{1}{2}\rho V^2 A$). Computes Reynolds number and dynamically calculates drag coefficient ($C_d$) based on empirical flow correlations (such as Morrison's equation for spheres).

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 -ffree-line-length-none drag_force_calculator.f90 -o drag_calc

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

drag_calc < input.txt

πŸ“₯ Downloads & Local Files

Preview of the required input file (input.txt):

! Geometry Code (1=Sphere, 2=Cylinder, 3=Flat Plate)
1
! Characteristic Dimension (Diameter or Length) [m]
0.05
! Free-stream Velocity ($V$) [m/s]
10.0
! Fluid Density ($\rho$) [kg/mΒ³]
1.2
! Fluid Viscosity ($\mu$) [Pa-s]
1.8e-5