🔬
Solver Purpose & Physical Scope
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:
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 502 times
📄 Source File:
boundary_layer.f90
📁 calcul/FluidMechanics /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 boundary_layer.f90 -o boundary_layer_calc
2. Execution with input.txt redirection:
boundary_layer_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
5.0 1.204 1.516e-5 1.0 500000.0
Parameter Description:
Free stream velocity U_inf [m/s] Fluid density [kg/m3] Kinematic viscosity nu [m2/s] Plate length L [m] Critical transition Reynolds number