program pitot_tube implicit none integer :: i, iostat_val double precision :: P_total, P_static, rho, gamma, R_gas, T_static double precision :: D_pipe, A_pipe, mu double precision :: V, q_dynamic, M, mdot, Re, P_impact double precision :: V_incomp, V_comp, correction, a_sound double precision :: dP_i, V_i, M_i, mdot_i, q_i double precision, parameter :: PI = 3.141592653589793d0 character(len=60) :: comp_note read(*,*,iostat=iostat_val) P_total if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid total pressure input.' stop end if read(*,*,iostat=iostat_val) P_static read(*,*,iostat=iostat_val) rho read(*,*,iostat=iostat_val) gamma read(*,*,iostat=iostat_val) R_gas read(*,*,iostat=iostat_val) T_static read(*,*,iostat=iostat_val) D_pipe read(*,*,iostat=iostat_val) mu if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all pitot tube inputs.' stop end if if (P_total <= 0.0d0 .or. P_static <= 0.0d0) then write(*,*) 'ERROR: Pressures must be positive.' stop end if if (P_total < P_static) then write(*,*) 'ERROR: Total pressure must be >= static pressure.' stop end if if (rho <= 0.0d0) then write(*,*) 'ERROR: Density must be positive.' stop end if if (gamma <= 1.0d0) gamma = 1.4d0 if (R_gas <= 0.0d0) R_gas = 287.0d0 if (T_static <= 0.0d0) T_static = 293.15d0 if (D_pipe <= 0.0d0) D_pipe = 0.1d0 if (mu <= 0.0d0) mu = 1.81d-5 A_pipe = PI / 4.0d0 * D_pipe**2 q_dynamic = P_total - P_static P_impact = q_dynamic ! Speed of sound a_sound = sqrt(gamma * R_gas * T_static) ! Incompressible Bernoulli: V = sqrt(2 dP / rho) V_incomp = sqrt(2.0d0 * q_dynamic / rho) ! Mach number from isentropic relation: ! P_total/P_static = (1 + (gamma-1)/2 M^2)^(gamma/(gamma-1)) M = mach_from_pressure_ratio(P_total/P_static, gamma) ! Compressible velocity V_comp = M * a_sound ! Compressibility correction factor if (V_incomp > 1.0d-30) then correction = V_comp / V_incomp else correction = 1.0d0 end if ! Choose appropriate velocity if (M < 0.3d0) then V = V_incomp comp_note = 'Incompressible (M < 0.3) — Bernoulli valid' else if (M < 1.0d0) then V = V_comp comp_note = 'Subsonic compressible — isentropic correction applied' else V = V_comp comp_note = 'Supersonic — Rayleigh pitot formula recommended' end if ! Mass flow rate mdot = rho * V * A_pipe ! Reynolds number Re = rho * V * D_pipe / mu write(*,'(A)') '============================================================' write(*,'(A)') ' PITOT TUBE & PITOT-STATIC PROBE ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A)') '--- INPUTS --------------------------------------------------' write(*,'(A,ES12.4,A)') ' Total Pressure Pt = ', P_total, ' Pa' write(*,'(A,ES12.4,A)') ' Static Pressure Ps = ', P_static, ' Pa' write(*,'(A,ES12.4,A)') ' Impact Pressure = ', P_impact, ' Pa' write(*,'(A,ES12.4,A)') ' Fluid Density = ', rho, ' kg/m3' write(*,'(A,ES12.4)') ' Gamma = ', gamma write(*,'(A,ES12.4,A)') ' Gas Constant R = ', R_gas, ' J/(kg.K)' write(*,'(A,ES12.4,A)') ' Static Temperature = ', T_static, ' K' write(*,'(A,ES12.4,A)') ' Pipe Diameter = ', D_pipe, ' m' write(*,'(A,ES12.4,A)') ' Viscosity = ', mu, ' Pa.s' write(*,*) write(*,'(A)') '--- VELOCITY RESULTS ----------------------------------------' write(*,'(A,ES12.4,A)') ' Speed of Sound = ', a_sound, ' m/s' write(*,'(A,ES12.4,A)') ' Incompressible Velocity = ', V_incomp, ' m/s' write(*,'(A,ES12.4,A)') ' Compressible Velocity = ', V_comp, ' m/s' write(*,'(A,ES12.4)') ' Compressibility Factor = ', correction write(*,'(A,ES12.4,A)') ' Selected Velocity = ', V, ' m/s' write(*,'(A,ES12.4)') ' Mach Number = ', M write(*,'(A,A)') ' Flow Regime = ', trim(comp_note) write(*,*) write(*,'(A)') '--- FLOW PARAMETERS -----------------------------------------' write(*,'(A,ES12.4,A)') ' Dynamic Pressure q = ', q_dynamic, ' Pa' write(*,'(A,ES12.4,A)') ' Mass Flow Rate = ', mdot, ' kg/s' write(*,'(A,ES12.4)') ' Reynolds Number = ', Re write(*,'(A,ES12.4,A)') ' Pipe Area = ', A_pipe, ' m2' write(*,*) ! Velocity vs differential pressure sweep write(*,'(A)') '--- VELOCITY VS DELTA-P SWEEP -------------------------------' write(*,'(A)') ' dP[Pa] V_incomp[m/s] V_comp[m/s] Mach correction' write(*,'(A)') ' -----------------------------------------------------------------------' do i = 1, 60 dP_i = q_dynamic * 0.02d0 * dble(i) if (dP_i <= 0.0d0) cycle V_i = sqrt(2.0d0 * dP_i / rho) M_i = mach_from_pressure_ratio((P_static + dP_i)/P_static, gamma) if (M_i > 0.0d0) then write(*,'(ES12.4,2X,F12.4,2X,F12.4,2X,F10.5,2X,F10.5)') & dP_i, V_i, M_i*a_sound, M_i, & M_i*a_sound/max(V_i,1.0d-30) else write(*,'(ES12.4,2X,F12.4,2X,F12.4,2X,F10.5,2X,F10.5)') & dP_i, V_i, V_i, 0.0d0, 1.0d0 end if end do write(*,*) ! Mach number vs pressure ratio sweep write(*,'(A)') '--- MACH VS PRESSURE RATIO SWEEP ----------------------------' write(*,'(A)') ' Pt/Ps Mach V[m/s] T_total[K]' write(*,'(A)') ' -----------------------------------------------------------' do i = 1, 50 dP_i = 1.0d0 + 9.0d0 * dble(i-1) / 49.0d0 ! Pt/Ps from 1 to 10 M_i = mach_from_pressure_ratio(dP_i, gamma) V_i = M_i * a_sound write(*,'(F10.4,2X,F10.5,2X,F12.3,2X,F12.2)') & dP_i, M_i, V_i, T_static*(1.0d0+(gamma-1.0d0)/2.0d0*M_i**2) end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' Bernoulli: V = sqrt(2 dP / rho).' write(*,'(A)') ' Isentropic: Pt/Ps = (1 + (gamma-1)/2 M^2)^(gamma/(gamma-1)).' write(*,'(A)') ' Compressible V = M * a; a = sqrt(gamma R T).' write(*,'(A)') ' Correction factor = V_comp / V_incomp.' contains double precision function mach_from_pressure_ratio(pr, gam) implicit none double precision, intent(in) :: pr, gam double precision :: lo, hi, mid, pr_mid integer :: it if (pr <= 1.0d0) then mach_from_pressure_ratio = 0.0d0 return end if lo = 0.0d0 hi = 10.0d0 do it = 1, 200 mid = 0.5d0*(lo + hi) pr_mid = (1.0d0 + (gam-1.0d0)/2.0d0 * mid**2)**(gam/(gam-1.0d0)) if (pr_mid < pr) then lo = mid else hi = mid end if end do mach_from_pressure_ratio = 0.5d0*(lo + hi) end function mach_from_pressure_ratio end program pitot_tube