🔬
Solver Purpose & Physical Scope
Compute flow velocity from pitot tube differential pressure. Includes Bernoulli incompressible and isentropic compressibility corrections.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 361 times
📄 Source File:
pitot_tube.f90
📁 calcul/FluidMechanics /
pitot_tube.f90
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 pitot_tube.f90 -o pitot_tube
2. Execution with input.txt redirection:
pitot_tube < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
Total (stagnation) pressure Pt\nStatic pressure Ps\nDensity Ï [kg/m³]\nGas constant R [J/(kg·K)]\nGas constant R [J/(kg·K)]\nStatic temperature T\nPipe diameter D\nViscosity μ [Pa·s]