π¬
Solver Purpose & Physical Scope
Calculate atmospheric boundary layer wind velocity gradient using Log-Law and Power-Law, dynamic wind stagnation pressure on facades, and base overturning moments.
π Discipline: Cfd
β‘ Precision: IEEE-754 64-bit Real(`real(8)`)
π₯ Total Downloads: 413 times
π Source File:
atmospheric_boundary_layer_wind.f90
π calcul/CFD /
atmospheric_boundary_layer_wind.f90
program atmospheric_boundary_layer_wind
implicit none
integer :: iostat_val
double precision :: Uref_ms, zref_m, z0_rough_m, alpha_exp, Hbld_m, bld_width_m, rho_air
double precision :: u_star, U_top_log, U_top_pow, q_top_Pa, q_top_kPa, Iz_top_pct
double precision :: total_shear_kN, overturn_moment_kNm
double precision, parameter :: KAPPA = 0.40d0
! Read inputs
read(*,*,iostat=iostat_val) Uref_ms ! Reference Wind Speed at 10m [m/s] (e.g. 25.0)
read(*,*,iostat=iostat_val) zref_m ! Reference Height [m] (e.g. 10.0)
read(*,*,iostat=iostat_val) z0_rough_m ! Terrain Roughness Length z0 [m] (e.g. 0.30)
read(*,*,iostat=iostat_val) alpha_exp ! Power-Law Exponent alpha (e.g. 0.22)
read(*,*,iostat=iostat_val) Hbld_m ! Building Height [m] (e.g. 150.0)
read(*,*,iostat=iostat_val) bld_width_m ! Building Facade Width [m] (e.g. 35.0)
read(*,*,iostat=iostat_val) rho_air ! Air Density [kg/m3] (e.g. 1.225)
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid input data for atmospheric boundary layer calculation.'
stop
end if
if (Uref_ms <= 0.0d0 .or. z0_rough_m <= 0.0d0 .or. Hbld_m <= 0.0d0) then
write(*,*) 'ERROR: Wind speed, roughness, and building height must be positive.'
stop
end if
! Friction velocity u* [m/s]
u_star = (KAPPA * Uref_ms) / log(max(1.1d0, zref_m / z0_rough_m))
! Wind Velocity at Top of Building
U_top_log = (u_star / KAPPA) * log(max(1.05d0, Hbld_m / z0_rough_m))
U_top_pow = Uref_ms * ((Hbld_m / zref_m)**alpha_exp)
! Dynamic Stagnation Pressure at Top [Pa]
q_top_Pa = 0.5d0 * rho_air * (U_top_log**2)
q_top_kPa = q_top_Pa / 1000.0d0
! Turbulence Intensity at Building Top [%]
Iz_top_pct = (1.0d0 / log(max(1.05d0, Hbld_m / z0_rough_m))) * 100.0d0
! Total Integrated Base Shear Force [kN] and Overturning Moment [kN.m] (approximated)
total_shear_kN = (0.5d0 * rho_air * (Uref_ms**2) * bld_width_m * Hbld_m * (1.0d0 / (2.0d0 * alpha_exp + 1.0d0)) * &
((Hbld_m / zref_m)**(2.0d0 * alpha_exp))) / 1000.0d0
overturn_moment_kNm = total_shear_kN * (Hbld_m * ((2.0d0 * alpha_exp + 1.0d0) / (2.0d0 * alpha_exp + 2.0d0)))
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' THERMOFLUIDCALC β ATMOSPHERIC BOUNDARY LAYER WIND ENGINE'
write(*,'(A)') '============================================================'
write(*,'(A,F10.2,A,F10.1,A)') 'Reference Wind Uref / zref = ', Uref_ms, ' m/s / ', zref_m, ' m'
write(*,'(A,F10.4,A,F10.3)') 'Roughness z0 / Power exp a= ', z0_rough_m, ' m / ', alpha_exp
write(*,'(A,F10.1,A,F10.1,A)') 'Building Height / Width = ', Hbld_m, ' m / ', bld_width_m, ' m'
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,F10.2,A)') 'Friction Velocity (u*) = ', u_star, ' m/s'
write(*,'(A,F10.2,A)') 'TOP WIND SPEED (Log-Law) = ', U_top_log, ' m/s'
write(*,'(A,F10.2,A)') 'TOP WIND SPEED (Power-Law)= ', U_top_pow, ' m/s'
write(*,'(A,F10.2,A,F8.3,A)') 'Top Dynamic Wind Pressure = ', q_top_Pa, ' Pa (', q_top_kPa, ' kPa)'
write(*,'(A,F10.2,A)') 'Top Turbulence Intensity = ', Iz_top_pct, ' %'
write(*,'(A,F10.2,A)') 'Base Overturning Moment = ', overturn_moment_kNm, ' kN.m'
write(*,'(A,F10.2,A)') 'Total Wind Base Shear = ', total_shear_kN, ' kN'
write(*,'(A)') '============================================================'
end program atmospheric_boundary_layer_wind
π» How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 atmospheric_boundary_layer_wind.f90 -o atmospheric_boundary_layer_wind
2. Execution with input.txt redirection:
atmospheric_boundary_layer_wind < input.txt
π Sample input.txt File Structure
Sample Data:
26.0 10.0 1.20 0.33 200.0 40.0 1.225
Parameter Description:
Reference Wind Speed [m/s]\nReference Height [m]\nRoughness Length z0 [m]\nPower-Law Exponent alpha\nBuilding Height [m]\nBuilding Width [m]\nAir Density [kg/mΒ³]