๐Ÿ’ป 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.

Aerodynamic Lift & Drag Solver

Core Numerical Engine in Fortran 90 โ€ข 25 total downloads

lift_drag.f90
! =========================================================================
! Source File: lift_drag.f90
! =========================================================================

!==============================================================================
! ThermoFluidCalc โ€” Calculator #24 : Lift CL & Induced Drag CDi
!==============================================================================
! Physics :
!   CL  = CL_alpha * (alpha - alpha_0)       (linear lift model)
!   CDi = CL^2 / (pi * e * AR)               (induced drag)
!   CD  = CD0 + CDi                           (total drag, parabolic polar)
!   AR  = b^2 / S                             (aspect ratio)
!   L   = CL * q * S ,   D = CD * q * S
!   L/D = CL / CD
!
! Modes:
!   1 = Single point (one angle of attack)
!   2 = Alpha sweep  (range of angles)
!
! Reference : Gupta, ยง7.8, Eq. 7.87
!
! Build:
!   gfortran -O2 -o lift_drag lift_drag.f90
!==============================================================================
program lift_drag
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: PI = 3.141592653589793238_dp
  integer, parameter :: MAX_PTS = 5000

  integer  :: mode, npts, i
  real(dp) :: rho, V, S, b, e, CD0, CL_alpha, alpha0_deg, alpha_deg
  real(dp) :: alpha_min, alpha_max, dalpha
  real(dp) :: AR, q, CL, CDi, CD, L_force, D_force, LD
  real(dp) :: alpha_rad, alpha0_rad
  real(dp) :: best_LD, best_alpha, best_CL

  read(*,*) mode

  select case (mode)

  !==========================================================================
  ! MODE 1 : Single point
  !==========================================================================
  case (1)
    backspace(5)
    read(*,*) mode, rho, V, S, b, e, CD0, CL_alpha, alpha_deg, alpha0_deg

    call validate_inputs(rho, V, S, b, e)

    AR        = b**2 / S
    q         = 0.5_dp * rho * V**2
    alpha_rad = alpha_deg * PI / 180.0_dp
    alpha0_rad= alpha0_deg * PI / 180.0_dp
    CL        = CL_alpha * (alpha_rad - alpha0_rad)
    CDi       = CL**2 / (PI * e * AR)
    CD        = CD0 + CDi
    L_force   = CL * q * S
    D_force   = CD * q * S
    if (CD > 0.0_dp) then
      LD = CL / CD
    else
      LD = 0.0_dp
    end if

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Single Point'
    write(*,'(A,F12.6)')    'AR=', AR
    write(*,'(A,ES15.8)')   'Q=', q
    write(*,'(A,F12.6)')    'ALPHA_DEG=', alpha_deg
    write(*,'(A,F12.6)')    'ALPHA0_DEG=', alpha0_deg
    write(*,'(A,F12.6)')    'CL=', CL
    write(*,'(A,F12.6)')    'CDI=', CDi
    write(*,'(A,F12.6)')    'CD0=', CD0
    write(*,'(A,F12.6)')    'CD=', CD
    write(*,'(A,ES15.8)')   'LIFT=', L_force
    write(*,'(A,ES15.8)')   'DRAG=', D_force
    write(*,'(A,F12.4)')    'LD=', LD

  !==========================================================================
  ! MODE 2 : Alpha sweep
  !==========================================================================
  case (2)
    backspace(5)
    read(*,*) mode, rho, V, S, b, e, CD0, CL_alpha, alpha0_deg, &
              alpha_min, alpha_max, npts

    call validate_inputs(rho, V, S, b, e)
    if (npts < 2) npts = 2
    if (npts > MAX_PTS) npts = MAX_PTS

    AR         = b**2 / S
    q          = 0.5_dp * rho * V**2
    alpha0_rad = alpha0_deg * PI / 180.0_dp
    dalpha     = (alpha_max - alpha_min) / real(npts - 1, dp)

    best_LD    = -1.0e30_dp
    best_alpha = alpha_min
    best_CL    = 0.0_dp

    write(*,'(A,I1)')       'MODE=', mode
    write(*,'(A)')          'MODE_NAME=Alpha Sweep'
    write(*,'(A,F12.6)')    'AR=', AR
    write(*,'(A,ES15.8)')   'Q=', q
    write(*,'(A,F12.6)')    'CD0=', CD0
    write(*,'(A,I5)')       'NPTS=', npts

    write(*,'(A)') 'DATA_START'
    do i = 0, npts - 1
      alpha_deg = alpha_min + real(i, dp) * dalpha
      alpha_rad = alpha_deg * PI / 180.0_dp
      CL  = CL_alpha * (alpha_rad - alpha0_rad)
      CDi = CL**2 / (PI * e * AR)
      CD  = CD0 + CDi
      if (CD > 0.0_dp) then
        LD = CL / CD
      else
        LD = 0.0_dp
      end if

      if (LD > best_LD) then
        best_LD    = LD
        best_alpha = alpha_deg
        best_CL    = CL
      end if

      write(*,'(F8.3,4(A,F12.6))') alpha_deg, ',', CL, ',', CDi, ',', CD, ',', LD
    end do
    write(*,'(A)') 'DATA_END'

    write(*,'(A,F12.4)')   'MAX_LD=', best_LD
    write(*,'(A,F8.3)')    'ALPHA_MAXLD=', best_alpha
    write(*,'(A,F12.6)')   'CL_MAXLD=', best_CL

  case default
    write(*,'(A)') 'ERROR=Invalid mode (must be 1 or 2).'
    stop
  end select

contains

  subroutine validate_inputs(rho, V, S, b, e)
    real(dp), intent(in) :: rho, V, S, b, e
    if (rho <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Density must be positive.'; stop
    end if
    if (V <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Velocity must be positive.'; stop
    end if
    if (S <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Wing area must be positive.'; stop
    end if
    if (b <= 0.0_dp) then
      write(*,'(A)') 'ERROR=Wingspan must be positive.'; stop
    end if
    if (e <= 0.0_dp .or. e > 1.0_dp) then
      write(*,'(A)') 'ERROR=Oswald efficiency must be in (0, 1].'; stop
    end if
  end subroutine

end program lift_drag


Solver Description

Compute lift force, drag force, and respective coefficients for 2D airfoils and 3D wings.

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 lift_drag.f90 -o lift_drag

Execution Command:

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

lift_drag < input.txt

๐Ÿ“ฅ Downloads & Local Files

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

! (kg/mร‚ยณ)\nV (m/s)\nS wing area (mร‚ยฒ)\nb wingspan (m)\ne Oswald eff.\nCD0\nCL (per rad)\n0 (ร‚ยฐ)\n0 (ร‚ยฐ)
1.225
! Parameter 2
60.0
! Parameter 3
16.2
! Parameter 4
11.0
! Parameter 5
0.8
! Parameter 6
0.034
! Parameter 7
5.7
! Parameter 8
5.0
! Parameter 9
-2.0