💻 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.
Blasius Boundary Layer Solver
Core Numerical Engine in Fortran 90 • 28 total downloads
blasius_boundary_layer.f90
! =========================================================================
! Source File: blasius_boundary_layer.f90
! =========================================================================
program blasius_boundary_layer
implicit none
double precision :: U_inf, rho, mu, x_pos, T_wall, T_inf
double precision :: Pr_in, cp_val, k_cond
double precision :: nu, Re_x, Re_crit, x_trans, Pr_mol
! laminar
double precision :: d_lam, ds_lam, th_lam, shape_H_lam
double precision :: Cf_lam, Cfa_lam, tau_lam, ut_lam, FD_lam
double precision :: dT_lam, Nu_lam, Nua_lam, h_lam, q_lam, St_lam
! turbulent
double precision :: d_turb, ds_turb, th_turb, shape_H_turb
double precision :: Cf_turb, Cfa_turb, tau_turb, ut_turb, FD_turb
double precision :: Nu_turb, Nua_turb, h_turb, q_turb, St_turb
character(len=20) :: regime
integer :: i, n_pts, ios
double precision :: xc, rc, dl, dt2, dsl, dst, thl, tht, cfl, cft, nul2, nut2
read(*,*,iostat=ios) U_inf; if(ios/=0)then;write(*,*)'ERROR: Invalid U_inf.';stop;end if
read(*,*,iostat=ios) rho; if(ios/=0)then;write(*,*)'ERROR: Invalid rho.';stop;end if
read(*,*,iostat=ios) mu; if(ios/=0)then;write(*,*)'ERROR: Invalid mu.';stop;end if
read(*,*,iostat=ios) x_pos; if(ios/=0)then;write(*,*)'ERROR: Invalid x.';stop;end if
read(*,*,iostat=ios) T_wall; if(ios/=0)then;write(*,*)'ERROR: Invalid T_wall.';stop;end if
read(*,*,iostat=ios) T_inf; if(ios/=0)then;write(*,*)'ERROR: Invalid T_inf.';stop;end if
read(*,*,iostat=ios) Pr_in; if(ios/=0)then;write(*,*)'ERROR: Invalid Pr.';stop;end if
read(*,*,iostat=ios) cp_val; if(ios/=0)then;write(*,*)'ERROR: Invalid cp.';stop;end if
read(*,*,iostat=ios) k_cond; if(ios/=0)then;write(*,*)'ERROR: Invalid k_cond.';stop;end if
if(U_inf<=0d0)then;write(*,*)'ERROR: U_inf must be > 0.';stop;end if
if(rho<=0d0)then;write(*,*)'ERROR: rho must be > 0.';stop;end if
if(mu<=0d0)then;write(*,*)'ERROR: mu must be > 0.';stop;end if
if(x_pos<=0d0)then;write(*,*)'ERROR: x must be > 0.';stop;end if
if(cp_val<=0d0)then;write(*,*)'ERROR: cp must be > 0.';stop;end if
if(k_cond<=0d0)then;write(*,*)'ERROR: k_cond must be > 0.';stop;end if
nu = mu / rho
Re_x = U_inf * x_pos / nu
Re_crit = 5d5
x_trans = Re_crit * nu / U_inf
if (Pr_in > 0d0) then
Pr_mol = Pr_in
else
Pr_mol = mu * cp_val / k_cond
end if
if (Re_x < 3d5) then
regime = 'LAMINAR'
else if (Re_x > 1d6) then
regime = 'TURBULENT'
else
regime = 'TRANSITIONAL'
end if
! ---- LAMINAR (Blasius) -------------------------------------------
d_lam = 5d0 * x_pos / sqrt(Re_x)
ds_lam = 1.7208d0 * x_pos / sqrt(Re_x)
th_lam = 0.664d0 * x_pos / sqrt(Re_x)
shape_H_lam = ds_lam / th_lam
Cf_lam = 0.664d0 / sqrt(Re_x)
Cfa_lam = 1.328d0 / sqrt(Re_x)
tau_lam = 0.5d0 * Cf_lam * rho * U_inf**2
ut_lam = sqrt(tau_lam / rho)
FD_lam = 0.5d0 * Cfa_lam * rho * U_inf**2 * x_pos
dT_lam = d_lam / Pr_mol**(1d0/3d0)
Nu_lam = 0.332d0 * sqrt(Re_x) * Pr_mol**(1d0/3d0)
Nua_lam = 0.664d0 * sqrt(Re_x) * Pr_mol**(1d0/3d0)
h_lam = Nu_lam * k_cond / x_pos
q_lam = h_lam * (T_wall - T_inf)
if (rho*cp_val*U_inf > 0d0) then
St_lam = h_lam / (rho * cp_val * U_inf)
else
St_lam = 0d0
end if
! ---- TURBULENT (1/7 power law) -----------------------------------
d_turb = 0.37d0 * x_pos / Re_x**0.2d0
ds_turb = 0.0463d0 * x_pos / Re_x**0.2d0
th_turb = 0.0360d0 * x_pos / Re_x**0.2d0
shape_H_turb = ds_turb / th_turb
Cf_turb = 0.0592d0 / Re_x**0.2d0
Cfa_turb = 0.074d0 / Re_x**0.2d0
tau_turb = 0.5d0 * Cf_turb * rho * U_inf**2
ut_turb = sqrt(tau_turb / rho)
FD_turb = 0.5d0 * Cfa_turb * rho * U_inf**2 * x_pos
Nu_turb = 0.0296d0 * Re_x**0.8d0 * Pr_mol**(1d0/3d0)
Nua_turb = 0.037d0 * Re_x**0.8d0 * Pr_mol**(1d0/3d0)
h_turb = Nu_turb * k_cond / x_pos
q_turb = h_turb * (T_wall - T_inf)
if (rho*cp_val*U_inf > 0d0) then
St_turb = h_turb / (rho * cp_val * U_inf)
else
St_turb = 0d0
end if
! ---- Output ------------------------------------------------------
write(*,'(A)') '============================================================'
write(*,'(A)') ' BLASIUS BOUNDARY LAYER CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
write(*,'(A,F12.4,A)') ' U_inf = ', U_inf, ' m/s'
write(*,'(A,F12.4,A)') ' Density (rho) = ', rho, ' kg/m3'
write(*,'(A,ES14.6,A)') ' Dyn. Viscosity (mu) = ', mu, ' Pa.s'
write(*,'(A,ES14.6,A)') ' Kinematic Visc (nu) = ', nu, ' m2/s'
write(*,'(A,F12.6,A)') ' x (from leading edge) = ', x_pos, ' m'
write(*,'(A,F12.2,A)') ' T_wall = ', T_wall, ' K'
write(*,'(A,F12.2,A)') ' T_inf = ', T_inf, ' K'
write(*,'(A,F12.4)') ' Prandtl Number (Pr) = ', Pr_mol
write(*,'(A,F12.2,A)') ' cp = ', cp_val, ' J/kg.K'
write(*,'(A,F12.6,A)') ' k (conductivity) = ', k_cond, ' W/m.K'
write(*,*)
write(*,'(A)') '--- REYNOLDS NUMBER & REGIME --------------------------------'
write(*,'(A,ES14.6)') ' Re_x = ', Re_x
write(*,'(A,ES14.6)') ' Re_crit (assumed) = ', Re_crit
write(*,'(A,F12.6,A)') ' x_transition = ', x_trans, ' m'
write(*,'(A,A)') ' Flow Regime = ', trim(regime)
write(*,*)
write(*,'(A)') '--- BOUNDARY LAYER (Laminar - Blasius) ----------------------'
write(*,'(A,F12.6,A)') ' delta (BL thick) = ', d_lam*1d3, ' mm'
write(*,'(A,F12.6,A)') ' delta* (displacement) = ', ds_lam*1d3, ' mm'
write(*,'(A,F12.6,A)') ' theta (momentum) = ', th_lam*1d3, ' mm'
write(*,'(A,F12.6)') ' H = delta*/theta = ', shape_H_lam
write(*,'(A,F12.6,A)') ' delta_T (thermal) = ', dT_lam*1d3, ' mm'
write(*,*)
write(*,'(A)') '--- BOUNDARY LAYER (Turbulent - 1/7 law) -------------------'
write(*,'(A,F12.6,A)') ' delta (BL thick) = ', d_turb*1d3, ' mm'
write(*,'(A,F12.6,A)') ' delta* (displacement) = ', ds_turb*1d3, ' mm'
write(*,'(A,F12.6,A)') ' theta (momentum) = ', th_turb*1d3, ' mm'
write(*,'(A,F12.6)') ' H = delta*/theta = ', shape_H_turb
write(*,*)
write(*,'(A)') '--- SKIN FRICTION -------------------------------------------'
write(*,'(A,20X,A14,A14)') ' ','Laminar','Turbulent'
write(*,'(A,ES14.6,2X,ES14.6)') ' Cf (local) ', Cf_lam, Cf_turb
write(*,'(A,ES14.6,2X,ES14.6)') ' Cf (avg) ', Cfa_lam, Cfa_turb
write(*,'(A,ES14.6,2X,ES14.6,A)') ' tau_w ', tau_lam, tau_turb, ' Pa'
write(*,'(A,ES14.6,2X,ES14.6,A)') ' u_tau ', ut_lam, ut_turb, ' m/s'
write(*,'(A,ES14.6,2X,ES14.6,A)') ' F_D/width ', FD_lam, FD_turb, ' N/m'
write(*,*)
write(*,'(A)') '--- HEAT TRANSFER -------------------------------------------'
write(*,'(A,20X,A14,A14)') ' ','Laminar','Turbulent'
write(*,'(A,ES14.6,2X,ES14.6)') ' Nu_x (local)', Nu_lam, Nu_turb
write(*,'(A,ES14.6,2X,ES14.6)') ' Nu_x (avg) ', Nua_lam, Nua_turb
write(*,'(A,ES14.6,2X,ES14.6,A)') ' h (local) ', h_lam, h_turb, ' W/m2.K'
write(*,'(A,ES14.6,2X,ES14.6,A)') ' q (flux) ', q_lam, q_turb, ' W/m2'
write(*,'(A,ES14.6,2X,ES14.6)') ' St ', St_lam, St_turb
write(*,*)
! ---- Profile sweep -----------------------------------------------
write(*,'(A)') '--- PROFILE vs x -------------------------------------------'
write(*,'(A)') ' x[m] Re_x d_lam[mm] d_turb[mm] Cf_lam Cf_turb Nu_lam'
write(*,'(A)') ' ----------------------------------------------------------------------------------'
n_pts = 40
do i = 1, n_pts
xc = 0.01d0 + (max(x_pos*2d0, 1d0) - 0.01d0) * dble(i-1) / dble(n_pts-1)
rc = U_inf * xc / nu
if (rc < 1d0) rc = 1d0
dl = 5d0 * xc / sqrt(rc) * 1d3
dt2 = 0.37d0 * xc / rc**0.2d0 * 1d3
cfl = 0.664d0 / sqrt(rc)
cft = 0.0592d0 / rc**0.2d0
nul2 = 0.332d0 * sqrt(rc) * Pr_mol**(1d0/3d0)
write(*,'(F8.4,2X,ES12.4,2X,F10.4,2X,F10.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
xc, rc, dl, dt2, cfl, cft, nul2
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' Blasius: delta = 5x/sqrt(Re_x)'
write(*,'(A)') ' Blasius: Cf = 0.664/sqrt(Re_x)'
write(*,'(A)') ' Turb: delta = 0.37x/Re_x^(1/5)'
write(*,'(A)') ' Turb: Cf = 0.0592/Re_x^(1/5)'
write(*,'(A)') ' Nu_lam = 0.332*Re_x^0.5*Pr^(1/3)'
write(*,'(A)') ' Nu_turb = 0.0296*Re_x^0.8*Pr^(1/3)'
write(*,'(A)') '============================================================'
end program blasius_boundary_layer
Solver Description
Solve boundary layer thickness, displacement thickness, skin friction, and velocity profiles over a flat plate.
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 blasius_boundary_layer.f90 -o blasius_boundary_layer
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
blasius_boundary_layer < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! U [m/s]\n[kg/m]\n[Pas]\nx_i\nTw_i\nTi_i\nPr\ncp [J/kgK]\nk_i
30
! Parameter 2
1.225
! Parameter 3
1.789e-5
! Parameter 4
0.5
! Parameter 5
350
! Parameter 6
300
! Parameter 7
0
! Parameter 8
1006
! Parameter 9
0.0262
30
! Parameter 2
1.225
! Parameter 3
1.789e-5
! Parameter 4
0.5
! Parameter 5
350
! Parameter 6
300
! Parameter 7
0
! Parameter 8
1006
! Parameter 9
0.0262