Fortran Source Code & Offline Testing

Boundary Layer Theory

Core Numerical Engine in Fortran 90 • 63 total downloads

boundary_layer.f90

! =========================================================================
! Source File: boundary_layer.f90
! =========================================================================

program boundary_layer
    implicit none
    
    double precision :: U_inf, rho, nu, L, Re_cr
    double precision :: Re_L, x_tr, Re_x
    double precision :: delta, delta_star, theta, Cf, tau_w
    double precision :: q_inf, CD_total, F_drag, A_corr
    double precision :: pi, x, dx
    integer :: i, n_stations
    character(len=20) :: regime
    
    ! Turbulent virtual origin
    double precision :: x_virt
    double precision :: delta_lam_tr, delta_turb
    
    pi = 3.14159265358979d0
    n_stations = 40
    
    ! Read inputs
    read(*,*) U_inf      ! Free-stream velocity [m/s]
    read(*,*) rho        ! Fluid density [kg/m3]
    read(*,*) nu         ! Kinematic viscosity [m2/s]
    read(*,*) L          ! Plate length [m]
    read(*,*) Re_cr      ! Critical Reynolds number for transition
    
    ! Global calculations
    Re_L = U_inf * L / nu
    x_tr = Re_cr * nu / U_inf
    q_inf = 0.5d0 * rho * U_inf**2
    
    ! ============================================
    ! HEADER
    ! ============================================
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   BOUNDARY LAYER ANALYSIS — FLAT PLATE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUT PARAMETERS ---'
    write(*,'(A,F12.4,A)')  '  Free-stream velocity U∞ = ', U_inf, ' m/s'
    write(*,'(A,F12.4,A)')  '  Fluid density ρ         = ', rho, ' kg/m³'
    write(*,'(A,ES12.4,A)') '  Kinematic viscosity ν   = ', nu, ' m²/s'
    write(*,'(A,F12.4,A)')  '  Plate length L          = ', L, ' m'
    write(*,'(A,ES12.4)')   '  Critical Re_x           = ', Re_cr
    write(*,*)
    write(*,'(A)') '--- GLOBAL RESULTS ---'
    write(*,'(A,ES12.4)')   '  Reynolds at plate end Re_L = ', Re_L
    
    if (x_tr < L) then
        write(*,'(A,F12.6,A)') '  Transition location x_tr   = ', x_tr, ' m'
        write(*,'(A,F8.2,A)')  '  Transition at              = ', (x_tr/L)*100.0d0, ' % of plate'
        write(*,'(A)')         '  Flow type                  = Mixed (laminar + turbulent)'
    else
        write(*,'(A,F12.6,A)') '  Transition location x_tr   = ', x_tr, ' m (beyond plate)'
        write(*,'(A)')         '  Flow type                  = Fully laminar'
    end if
    write(*,*)
    
    ! ============================================
    ! TOTAL DRAG COEFFICIENT
    ! ============================================
    write(*,'(A)') '--- TOTAL DRAG COEFFICIENT ---'
    
    if (Re_L < Re_cr) then
        ! Fully laminar
        CD_total = 1.328d0 / sqrt(Re_L)
        write(*,'(A,ES12.6)')   '  C_D (laminar)           = ', CD_total
        write(*,'(A)')          '  Formula: C_D = 1.328 / √Re_L'
    else
        ! Mixed boundary layer
        ! A depends on Re_cr: A = Re_cr × (0.074/Re_cr^0.2 − 1.328/√Re_cr)
        A_corr = Re_cr * (0.074d0 / Re_cr**0.2d0 - 1.328d0 / sqrt(Re_cr))
        CD_total = 0.074d0 / Re_L**0.2d0 - A_corr / Re_L
        write(*,'(A,ES12.6)')   '  C_D (mixed BL)          = ', CD_total
        write(*,'(A,F12.4)')    '  Correction A            = ', A_corr
        write(*,'(A)')          '  Formula: C_D = 0.074/Re_L^0.2 − A/Re_L'
    end if
    
    F_drag = CD_total * q_inf * L * 1.0d0  ! per unit width
    write(*,'(A,F12.6,A)')  '  Drag force (per m width)= ', F_drag, ' N/m'
    write(*,'(A,F12.4,A)')  '  Dynamic pressure q∞     = ', q_inf, ' Pa'
    write(*,*)
    
    ! ============================================
    ! BOUNDARY LAYER AT PLATE END
    ! ============================================
    write(*,'(A)') '--- BOUNDARY LAYER AT x = L ---'
    Re_x = U_inf * L / nu
    
    if (Re_L < Re_cr) then
        delta = 5.0d0 * L / sqrt(Re_x)
        delta_star = 1.7208d0 * L / sqrt(Re_x)
        theta = 0.664d0 * L / sqrt(Re_x)
        Cf = 0.664d0 / sqrt(Re_x)
        regime = 'Laminar'
    else
        delta = 0.37d0 * L / Re_x**0.2d0
        delta_star = delta / 8.0d0
        theta = 7.0d0 * delta / 72.0d0
        Cf = 0.0592d0 / Re_x**0.2d0
        regime = 'Turbulent'
    end if
    
    tau_w = Cf * q_inf
    
    write(*,'(A,A)')        '  Regime                  = ', trim(regime)
    write(*,'(A,F12.6,A)')  '  δ  (BL thickness)       = ', delta*1000, ' mm'
    write(*,'(A,F12.6,A)')  '  δ* (displacement)       = ', delta_star*1000, ' mm'
    write(*,'(A,F12.6,A)')  '  θ  (momentum)           = ', theta*1000, ' mm'
    write(*,'(A,ES12.6)')   '  C_f (skin friction)     = ', Cf
    write(*,'(A,F12.4,A)')  '  τ_w (wall shear)        = ', tau_w, ' Pa'
    write(*,'(A,F12.6)')    '  Shape factor H = δ*/θ   = ', delta_star / theta
    write(*,*)
    
    ! ============================================
    ! PROFILE TABLE
    ! ============================================
    write(*,'(A)') '--- BOUNDARY LAYER PROFILE ALONG PLATE ---'
    write(*,'(A)') '  x [m]      Re_x         δ [mm]     δ* [mm]    θ [mm]     Cf          τ_w [Pa]   Regime'
    write(*,'(A)') '  -------    ----------   --------   --------   --------   ----------  --------   ------'
    
    do i = 1, n_stations
        x = L * dble(i) / dble(n_stations)
        Re_x = U_inf * x / nu
        
        if (Re_x < Re_cr) then
            ! Laminar (Blasius)
            delta = 5.0d0 * x / sqrt(Re_x)
            delta_star = 1.7208d0 * x / sqrt(Re_x)
            theta = 0.664d0 * x / sqrt(Re_x)
            Cf = 0.664d0 / sqrt(Re_x)
            regime = 'Laminar'
        else
            ! Turbulent (1/7 power law)
            delta = 0.37d0 * x / Re_x**0.2d0
            delta_star = delta / 8.0d0
            theta = 7.0d0 * delta / 72.0d0
            Cf = 0.0592d0 / Re_x**0.2d0
            regime = 'Turbulent'
        end if
        
        tau_w = Cf * q_inf
        
        write(*,'(F8.4,2X,ES12.4,2X,F8.4,3X,F8.4,3X,F8.4,3X,ES10.4,2X,F8.4,3X,A)') &
            x, Re_x, delta*1000, delta_star*1000, theta*1000, Cf, tau_w, trim(regime)
    end do
    
    write(*,*)
    
    ! ============================================
    ! CHART DATA (machine-readable)
    ! ============================================
    write(*,'(A)') '--- CHART_DATA ---'
    do i = 1, n_stations
        x = L * dble(i) / dble(n_stations)
        Re_x = U_inf * x / nu
        
        if (Re_x < Re_cr) then
            delta = 5.0d0 * x / sqrt(Re_x)
            delta_star = 1.7208d0 * x / sqrt(Re_x)
            theta = 0.664d0 * x / sqrt(Re_x)
            Cf = 0.664d0 / sqrt(Re_x)
            write(*,'(F10.6,A,F10.6,A,F10.6,A,F10.6,A,ES12.6,A,I1)') &
                x, ',', delta*1000, ',', delta_star*1000, ',', theta*1000, ',', Cf, ',', 0
        else
            delta = 0.37d0 * x / Re_x**0.2d0
            delta_star = delta / 8.0d0
            theta = 7.0d0 * delta / 72.0d0
            Cf = 0.0592d0 / Re_x**0.2d0
            write(*,'(F10.6,A,F10.6,A,F10.6,A,F10.6,A,ES12.6,A,I1)') &
                x, ',', delta*1000, ',', delta_star*1000, ',', theta*1000, ',', Cf, ',', 1
        end if
    end do
    write(*,'(A)') '--- END_CHART_DATA ---'
    
    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ---'
    write(*,'(A)') '  LAMINAR (Blasius):'
    write(*,'(A)') '    δ = 5.0 x / √Re_x'
    write(*,'(A)') '    δ* = 1.7208 x / √Re_x'
    write(*,'(A)') '    θ = 0.664 x / √Re_x'
    write(*,'(A)') '    C_f = 0.664 / √Re_x'
    write(*,'(A)') ''
    write(*,'(A)') '  TURBULENT (1/7 power law):'
    write(*,'(A)') '    δ = 0.37 x / Re_x^0.2'
    write(*,'(A)') '    δ* = δ / 8'
    write(*,'(A)') '    θ = 7δ / 72'
    write(*,'(A)') '    C_f = 0.0592 / Re_x^0.2'
    write(*,'(A)') ''
    write(*,'(A)') '  MIXED DRAG:'
    write(*,'(A)') '    C_D = 0.074/Re_L^0.2 − A/Re_L'
    write(*,'(A)') '    where A = Re_cr(0.074/Re_cr^0.2 − 1.328/√Re_cr)'
    
end program boundary_layer

Solver Description

Fluid flow over a flat plate creates a boundary layer due to shear stresses. The boundary layer thickness ($\delta$) and skin friction coefficient ($C_f$) are modeled by Blasius (laminar) and empirical turbulent equations:

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 boundary_layer.f90 -o boundary_layer_calc

Execution Command:

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

boundary_layer_calc < input.txt

📥 Downloads & Local Files

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

! Free stream velocity U_inf [m/s]
5.0
! Fluid density [kg/m3]
1.204
! Kinematic viscosity nu [m2/s]
1.516e-5
! Plate length L [m]
1.0
! Critical transition Reynolds number
500000.0

📦 Download Offline Bundle (ZIP) 20 dl 📄 Download input.txt 23 dl 📥 Download boundary_layer.f90 20 dl

💻 Fortran Source Code Library

Boundary Layer Theory

Solver Description

Key Numerical Methods & Architecture

🛠️ Local Compilation

📥 Downloads & Local Files