🔬
Solver Purpose & Physical Scope
Solve boundary layer thickness, displacement thickness, skin friction, and velocity profiles over a flat plate.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 461 times
📄 Source File:
blasius_boundary_layer.f90
📁 calcul/CFD /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 blasius_boundary_layer.f90 -o blasius_boundary_layer
2. Execution with input.txt redirection:
blasius_boundary_layer < input.txt
📄 Sample input.txt File Structure
Sample Data:
30 1.225 1.789e-5 0.5 350 300 0 1006 0.0262
Parameter Description:
U [m/s]\n[kg/m]\n[Pas]\nx_i\nTw_i\nTi_i\nPr\ncp [J/kgK]\nk_i