π¬
Solver Purpose & Physical Scope
Evaluate 1D compressible flow with heat addition or combustion, computing exit Mach (M2), thermal choking limit (q_max), stagnation pressure loss, and entropy generation.
π Discipline: Cfd
β‘ Precision: IEEE-754 64-bit Real(`real(8)`)
π₯ Total Downloads: 241 times
π Source File:
compressible_rayleigh_flow_heat.f90
π calcul/CFD /
compressible_rayleigh_flow_heat.f90
program compressible_rayleigh_flow_heat
implicit none
integer :: iostat_val, iter
double precision :: M1_mach, T01_K, q_heat_kJkg, cp_kJkgK, P01_bar, gamma_ratio
double precision :: R_gas, T0_T0star_1, T0_star_K, q_max_kJkg, T02_K, T0_T0star_2
double precision :: M2_mach, M_low, M_high, M_mid, T0_mid
double precision :: T_Tstar_1, T_Tstar_2, P_Pstar_1, P_Pstar_2, P0_P0star_1, P0_P0star_2
double precision :: P1_bar, P2_bar, P02_bar, T1_K, T2_K, delta_s_JkgK
logical :: is_choked
! Read inputs
read(*,*,iostat=iostat_val) M1_mach ! Inlet Mach Number (e.g. 0.30)
read(*,*,iostat=iostat_val) T01_K ! Inlet Total Temp [K] (e.g. 300.0)
read(*,*,iostat=iostat_val) q_heat_kJkg ! Heat Added [kJ/kg] (e.g. 400.0)
read(*,*,iostat=iostat_val) cp_kJkgK ! Specific Heat [kJ/kg.K] (e.g. 1.005)
read(*,*,iostat=iostat_val) P01_bar ! Inlet Total Pressure [bar] (e.g. 5.0)
read(*,*,iostat=iostat_val) gamma_ratio ! Specific Heat Ratio (e.g. 1.40)
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid input data for Rayleigh flow calculation.'
stop
end if
if (M1_mach <= 0.01d0 .or. T01_K <= 0.0d0 .or. cp_kJkgK <= 0.0d0 .or. P01_bar <= 0.0d0) then
write(*,*) 'ERROR: Mach, temperature, cp, and pressure must be positive.'
stop
end if
R_gas = cp_kJkgK * (1.0d0 - 1.0d0 / gamma_ratio) * 1000.0d0 ! J/kg.K
! Inlet Total Temperature Ratio
T0_T0star_1 = (2.0d0 * (gamma_ratio + 1.0d0) * (M1_mach**2) * &
(1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M1_mach**2))) / &
((1.0d0 + gamma_ratio * (M1_mach**2))**2)
T0_star_K = T01_K / T0_T0star_1
q_max_kJkg = cp_kJkgK * (T0_star_K - T01_K)
if (q_heat_kJkg >= q_max_kJkg) then
is_choked = .true.
M2_mach = 1.0d0
T02_K = T0_star_K
T0_T0star_2 = 1.0d0
else
is_choked = .false.
T02_K = T01_K + (q_heat_kJkg / cp_kJkgK)
T0_T0star_2 = T02_K / T0_star_K
! Solve for M2 via Bisection
if (M1_mach < 1.0d0) then
M_low = M1_mach
M_high = 0.9999d0
else
M_low = 1.0001d0
M_high = M1_mach
end if
do iter = 1, 100
M_mid = 0.5d0 * (M_low + M_high)
T0_mid = (2.0d0 * (gamma_ratio + 1.0d0) * (M_mid**2) * &
(1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M_mid**2))) / &
((1.0d0 + gamma_ratio * (M_mid**2))**2)
if (M1_mach < 1.0d0) then
if (T0_mid < T0_T0star_2) then
M_low = M_mid
else
M_high = M_mid
end if
else
if (T0_mid < T0_T0star_2) then
M_high = M_mid
else
M_low = M_mid
end if
end if
end do
M2_mach = M_mid
end if
! Static and Stagnation Ratios
T_Tstar_1 = (M1_mach**2) * ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M1_mach**2)))**2
T_Tstar_2 = (M2_mach**2) * ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M2_mach**2)))**2
P_Pstar_1 = (gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M1_mach**2))
P_Pstar_2 = (gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M2_mach**2))
P0_P0star_1 = ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M1_mach**2))) * &
(((2.0d0 + (gamma_ratio - 1.0d0) * (M1_mach**2)) / (gamma_ratio + 1.0d0))**(gamma_ratio / (gamma_ratio - 1.0d0)))
P0_P0star_2 = ((gamma_ratio + 1.0d0) / (1.0d0 + gamma_ratio * (M2_mach**2))) * &
(((2.0d0 + (gamma_ratio - 1.0d0) * (M2_mach**2)) / (gamma_ratio + 1.0d0))**(gamma_ratio / (gamma_ratio - 1.0d0)))
P1_bar = P01_bar * (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M1_mach**2))**(-gamma_ratio / (gamma_ratio - 1.0d0))
P2_bar = P1_bar * (P_Pstar_2 / P_Pstar_1)
P02_bar = P01_bar * (P0_P0star_2 / P0_P0star_1)
T1_K = T01_K / (1.0d0 + 0.5d0 * (gamma_ratio - 1.0d0) * (M1_mach**2))
T2_K = T1_K * (T_Tstar_2 / T_Tstar_1)
delta_s_JkgK = (cp_kJkgK * 1000.0d0) * log(T2_K / T1_K) - R_gas * log(P2_bar / P1_bar)
! Output Results
write(*,'(A)') '============================================================'
write(*,'(A)') ' THERMOFLUIDCALC β COMPRESSIBLE RAYLEIGH HEAT FLOW ENGINE'
write(*,'(A)') '============================================================'
write(*,'(A,F10.2,A,F10.1,A)') 'Inlet Mach / Stagnation T01 = ', M1_mach, ' / ', T01_K, ' K'
write(*,'(A,F10.2,A,F10.2,A)') 'Heat Added / Max Choking q = ', q_heat_kJkg, ' kJ/kg / ', q_max_kJkg, ' kJ/kg'
write(*,'(A)') '------------------------------------------------------------'
write(*,'(A,F10.4)') 'EXIT MACH NUMBER (M2) = ', M2_mach
if (is_choked) then
write(*,'(A)') 'Flow Status = THERMALLY CHOKED (M2 = 1.0)'
else
write(*,'(A)') 'Flow Status = UNCHOKED HEAT ADDITION'
end if
write(*,'(A,F10.1,A,F10.1,A)') 'Total Temperature T01 -> T02= ', T01_K, ' K -> ', T02_K, ' K'
write(*,'(A,F10.3,A,F10.3,A)') 'Total Pressure P01 -> P02 = ', P01_bar, ' bar -> ', P02_bar, ' bar'
write(*,'(A,F10.1,A,F10.1,A)') 'Static Temp T1 -> T2 = ', T1_K, ' K -> ', T2_K, ' K'
write(*,'(A,F10.2,A)') 'Entropy Generation Delta s = ', delta_s_JkgK, ' J/(kg.K)'
write(*,'(A)') '============================================================'
end program compressible_rayleigh_flow_heat
π» How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 compressible_rayleigh_flow_heat.f90 -o compressible_rayleigh_flow_heat
2. Execution with input.txt redirection:
compressible_rayleigh_flow_heat < input.txt
π Sample input.txt File Structure
Sample Data:
0.22 650.0 850.0 1.15 15.0 1.33
Parameter Description:
Inlet Mach Number\nInlet Total Temp [K]\nHeat Added [kJ/kg]\nSpecific Heat cp [kJ/(kgΒ·K)]\nInlet Total Pressure [bar]\nSpecific Heat Ratio gamma