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