💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Joule-Thomson Throttling
Core Numerical Engine in Fortran 90 • 22 total downloads
throttling_jt.f90
! =========================================================================
! Source File: throttling_jt.f90
! =========================================================================
program throttling_jt
implicit none
integer :: gas_type, iostat_val, i, n_sweep
double precision :: T_inlet, P_inlet, P_outlet
double precision :: Tc, Pc, omega, cp_gas, M_gas
double precision :: Tc_c, Pc_c, omega_c, cp_c, M_c
double precision :: R_gas, a_vdw, b_vdw
double precision :: mu_JT, T_outlet, dT, T_inv
double precision :: T_boil_est, P_H2O
double precision :: T_sw, mu_sw
character(len=40) :: gas_name, cooling_str, liq_str
R_gas = 8.314462d0 ! J/(mol.K)
! ── Read inputs ─────────────────────────────────────────────
read(*,*,iostat=iostat_val) gas_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid gas type input.'
stop
end if
read(*,*,iostat=iostat_val) T_inlet
read(*,*,iostat=iostat_val) P_inlet
read(*,*,iostat=iostat_val) P_outlet
read(*,*,iostat=iostat_val) Tc_c
read(*,*,iostat=iostat_val) Pc_c
read(*,*,iostat=iostat_val) omega_c
read(*,*,iostat=iostat_val) cp_c
read(*,*,iostat=iostat_val) M_c
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all inputs.'
stop
end if
! ── Input validation ────────────────────────────────────────
if (T_inlet <= 0.0d0) T_inlet = 300.0d0
if (P_inlet <= 0.0d0) P_inlet = 10.0d0
if (P_outlet <= 0.0d0) P_outlet = 0.1d0
if (P_outlet >= P_inlet) then
write(*,*) 'ERROR: P_outlet must be less than P_inlet.'
stop
end if
! ── Gas properties ──────────────────────────────────────────
if (gas_type == 1) then
gas_name = 'Nitrogen (N2)'
Tc = 126.2d0; Pc = 3.39d0; omega = 0.037d0
cp_gas = 29.1d0; M_gas = 28.014d0
else if (gas_type == 2) then
gas_name = 'Oxygen (O2)'
Tc = 154.6d0; Pc = 5.04d0; omega = 0.022d0
cp_gas = 29.4d0; M_gas = 32.0d0
else if (gas_type == 3) then
gas_name = 'Carbon Dioxide (CO2)'
Tc = 304.2d0; Pc = 7.38d0; omega = 0.224d0
cp_gas = 37.1d0; M_gas = 44.01d0
else if (gas_type == 4) then
gas_name = 'Methane (CH4)'
Tc = 190.6d0; Pc = 4.60d0; omega = 0.012d0
cp_gas = 35.7d0; M_gas = 16.04d0
else if (gas_type == 5) then
gas_name = 'Hydrogen (H2)'
Tc = 33.2d0; Pc = 1.30d0; omega = -0.217d0
cp_gas = 28.8d0; M_gas = 2.016d0
else if (gas_type == 6) then
gas_name = 'Air'
Tc = 132.5d0; Pc = 3.77d0; omega = 0.035d0
cp_gas = 29.1d0; M_gas = 28.97d0
else
gas_name = 'Custom Gas'
Tc = Tc_c; Pc = Pc_c; omega = omega_c
cp_gas = cp_c; M_gas = M_c
if (Tc <= 0.0d0) Tc = 200.0d0
if (Pc <= 0.0d0) Pc = 5.0d0
if (cp_gas <= 0.0d0) cp_gas = 30.0d0
if (M_gas <= 0.0d0) M_gas = 28.0d0
end if
! ── Van der Waals constants ─────────────────────────────────
! Pc in MPa -> Pa for calculation
a_vdw = 27.0d0 * R_gas**2 * Tc**2 / (64.0d0 * Pc * 1.0d6)
b_vdw = R_gas * Tc / (8.0d0 * Pc * 1.0d6)
! ── Joule-Thomson coefficient ───────────────────────────────
! mu_JT = (1/cp) * (2a/(RT) - b) in K/Pa
mu_JT = (1.0d0 / cp_gas) * (2.0d0 * a_vdw / (R_gas * T_inlet) - b_vdw)
! Convert to K/MPa
! mu_JT is already in K/Pa, multiply by 1e6 for K/MPa
! Actually: dT = mu_JT [K/Pa] * dP [Pa]
! P_inlet, P_outlet in MPa, so dP in Pa = (P_outlet - P_inlet)*1e6
! Temperature change
dT = mu_JT * (P_outlet - P_inlet) * 1.0d6
T_outlet = T_inlet + dT
! Inversion temperature (van der Waals)
T_inv = 2.0d0 * a_vdw / (R_gas * b_vdw)
! Boiling point estimate (rough: T_boil ~ 0.6 * Tc for most gases)
T_boil_est = 0.6d0 * Tc
! Cooling or heating?
if (mu_JT > 0.0d0) then
cooling_str = 'COOLING (mu_JT > 0)'
else if (mu_JT < 0.0d0) then
cooling_str = 'HEATING (mu_JT < 0)'
else
cooling_str = 'NO CHANGE (ideal gas)'
end if
! Liquefaction feasibility
if (T_outlet < T_boil_est .and. T_outlet > 0.0d0) then
liq_str = 'POSSIBLE — T_out < T_boil'
else
liq_str = 'NOT FEASIBLE at these conditions'
end if
! ── Output ──────────────────────────────────────────────────
write(*,'(A)') '============================================================'
write(*,'(A)') ' THROTTLING VALVE & JOULE-THOMSON EFFECT'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Gas = ', trim(gas_name)
write(*,'(A,F12.2,A)') ' Inlet Temperature T_in = ', T_inlet, ' K'
write(*,'(A,F12.4,A)') ' Inlet Pressure P_in = ', P_inlet, ' MPa'
write(*,'(A,F12.4,A)') ' Outlet Pressure P_out = ', P_outlet, ' MPa'
write(*,*)
write(*,'(A)') '--- GAS PROPERTIES ------------------------------------------'
write(*,'(A,F12.2,A)') ' Critical Temperature Tc = ', Tc, ' K'
write(*,'(A,F12.4,A)') ' Critical Pressure Pc = ', Pc, ' MPa'
write(*,'(A,F12.4)') ' Acentric Factor omega = ', omega
write(*,'(A,F12.4,A)') ' cp (molar) = ', cp_gas, ' J/(mol.K)'
write(*,'(A,F12.4,A)') ' Molar Mass M = ', M_gas, ' g/mol'
write(*,*)
write(*,'(A)') '--- VAN DER WAALS CONSTANTS ---------------------------------'
write(*,'(A,ES14.6,A)') ' a = ', a_vdw, ' Pa.m^6/mol^2'
write(*,'(A,ES14.6,A)') ' b = ', b_vdw, ' m^3/mol'
write(*,*)
write(*,'(A)') '--- JOULE-THOMSON COEFFICIENT -------------------------------'
write(*,'(A,ES14.6,A)') ' mu_JT = ', mu_JT, ' K/Pa'
write(*,'(A,F12.6,A)') ' mu_JT = ', mu_JT*1.0d6, ' K/MPa'
write(*,'(A,A)') ' Effect = ', trim(cooling_str)
write(*,*)
write(*,'(A)') '--- TEMPERATURE CHANGE --------------------------------------'
write(*,'(A,F12.4,A)') ' Delta_T = ', dT, ' K'
write(*,'(A,F12.2,A)') ' T_outlet = ', T_outlet, ' K'
write(*,'(A,F12.2,A)') ' T_outlet = ', T_outlet-273.15d0, ' C'
write(*,*)
write(*,'(A)') '--- INVERSION TEMPERATURE -----------------------------------'
write(*,'(A,F12.2,A)') ' T_inversion (vdW max) = ', T_inv, ' K'
write(*,'(A,F12.2,A)') ' T_inversion = ', T_inv-273.15d0, ' C'
if (T_inlet < T_inv) then
write(*,'(A)') ' T_in < T_inv => COOLING region'
else
write(*,'(A)') ' T_in > T_inv => HEATING region'
end if
write(*,*)
write(*,'(A)') '--- LIQUEFACTION FEASIBILITY --------------------------------'
write(*,'(A,F12.2,A)') ' T_boiling (est.) = ', T_boil_est, ' K'
write(*,'(A,A)') ' Liquefaction = ', trim(liq_str)
write(*,*)
write(*,'(A)') '--- IDEAL GAS COMPARISON ------------------------------------'
write(*,'(A)') ' Ideal gas: mu_JT = 0, no temperature change.'
write(*,'(A,F12.4,A)') ' Real gas Delta_T = ', dT, ' K'
write(*,*)
! ── Sensitivity sweep: T from 50 to 800 K ─────────────────
n_sweep = 40
write(*,'(A)') '--- SENSITIVITY: MU_JT VS INLET TEMPERATURE ----------------'
write(*,'(A)') ' T_in[K] mu_JT[K/MPa] T_out[K]'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, n_sweep
T_sw = 50.0d0 + dble(i-1) * (800.0d0 - 50.0d0) / dble(n_sweep - 1)
mu_sw = (1.0d0/cp_gas) * (2.0d0*a_vdw/(R_gas*T_sw) - b_vdw)
write(*,'(F10.2,4X,F14.6,4X,F12.2)') T_sw, mu_sw*1.0d6, &
T_sw + mu_sw*(P_outlet-P_inlet)*1.0d6
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Van der Waals: a=27R^2Tc^2/(64Pc), b=RTc/(8Pc)'
write(*,'(A)') ' mu_JT = (1/cp)*(2a/(RT) - b) [K/Pa]'
write(*,'(A)') ' T_inversion = 2a/(Rb) (maximum inversion temp)'
write(*,'(A)') ' Isenthalpic process: h_in = h_out'
write(*,'(A)') ' T_out = T_in + mu_JT * (P_out - P_in)'
write(*,'(A)') ' Ideal gas: mu_JT = 0 (no intermolecular forces)'
end program throttling_jt
Solver Description
Solves the Joule-Thomson throttling effect (isenthalpic expansion) for real gases using the van der Waals equation of state. Computes the Joule-Thomson coefficient ($\mu_{JT}$), outlet temperature change ($\Delta T$), maximum inversion temperature boundary, and determines whether the throttling process results in cooling, heating, or gas liquefaction.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 throttling_jt.f90 -o throttling_jt
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
throttling_jt < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Gas type (1=Nitrogen, 2=Oxygen, 3=CO2, 4=Methane, 5=Hydrogen, 6=Helium, 7=Custom)
1
! Inlet temperature Tin [K]
300.0
! Inlet pressure Pin [MPa]
20.0
! Outlet pressure Pout [MPa]
0.1
! Critical temperature Tc [K] for custom gas
200.0
! Critical pressure Pc [MPa] for custom gas
5.0
! Acentric factor omega for custom gas
0.04
! Molar heat capacity cp [J/(mol-K)] for custom gas
30.0
! Molar mass M [g/mol] for custom gas
28.0
1
! Inlet temperature Tin [K]
300.0
! Inlet pressure Pin [MPa]
20.0
! Outlet pressure Pout [MPa]
0.1
! Critical temperature Tc [K] for custom gas
200.0
! Critical pressure Pc [MPa] for custom gas
5.0
! Acentric factor omega for custom gas
0.04
! Molar heat capacity cp [J/(mol-K)] for custom gas
30.0
! Molar mass M [g/mol] for custom gas
28.0