💻 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.
Real Gas Equations of State
Core Numerical Engine in Fortran 90 • 27 total downloads
real_gas_eos.f90
! =========================================================================
! Source File: real_gas_eos.f90
! =========================================================================
program real_gas_eos
implicit none
integer :: eos_type, i, iostat_val
double precision :: T, P, Tc, Pc, omega, Rgas, Mmol
double precision :: a_eos, b_eos, kappa, alpha_pr
double precision :: v, Z, v_ideal, Z_ideal
double precision :: dh_dep, ds_dep, dg_dep, fugacity, phi_fug
double precision :: P_i, v_i, Z_i, T_i
double precision, parameter :: Ru = 8.314462d0 ! J/(mol K)
character(len=40) :: eos_name
read(*,*,iostat=iostat_val) eos_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid EOS type input.'
stop
end if
read(*,*,iostat=iostat_val) T
read(*,*,iostat=iostat_val) P
read(*,*,iostat=iostat_val) Tc
read(*,*,iostat=iostat_val) Pc
read(*,*,iostat=iostat_val) omega
read(*,*,iostat=iostat_val) Mmol
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all EOS inputs.'
stop
end if
if (T <= 0.0d0 .or. P <= 0.0d0 .or. Tc <= 0.0d0 .or. Pc <= 0.0d0) then
write(*,*) 'ERROR: T, P, Tc, Pc must be positive.'
stop
end if
if (Mmol <= 0.0d0) Mmol = 28.97d0
Rgas = Ru / (Mmol * 1.0d-3) ! J/(kg K) if Mmol in g/mol
! Ideal gas
v_ideal = Ru * T / P ! m3/mol
Z_ideal = 1.0d0
select case(eos_type)
case(1)
eos_name = 'van der Waals'
a_eos = 27.0d0 * Ru**2 * Tc**2 / (64.0d0 * Pc)
b_eos = Ru * Tc / (8.0d0 * Pc)
call solve_cubic_vdw(T, P, a_eos, b_eos, Ru, v)
case(2)
eos_name = 'Redlich-Kwong'
a_eos = 0.42748d0 * Ru**2 * Tc**2.5d0 / Pc
b_eos = 0.08664d0 * Ru * Tc / Pc
call solve_cubic_rk(T, P, a_eos, b_eos, Ru, v)
case(3)
eos_name = 'Peng-Robinson'
kappa = 0.37464d0 + 1.54226d0*omega - 0.26992d0*omega**2
alpha_pr = (1.0d0 + kappa*(1.0d0 - sqrt(T/Tc)))**2
a_eos = 0.45724d0 * Ru**2 * Tc**2 / Pc * alpha_pr
b_eos = 0.07780d0 * Ru * Tc / Pc
call solve_cubic_pr(T, P, a_eos, b_eos, Ru, v)
case default
write(*,*) 'ERROR: EOS type must be 1 vdW, 2 RK, or 3 PR.'
stop
end select
Z = P * v / (Ru * T)
! Departure functions (molar basis, J/mol)
call departure_functions(eos_type, T, P, v, a_eos, b_eos, Ru, Tc, &
dh_dep, ds_dep, dg_dep, fugacity, phi_fug)
write(*,'(A)') '============================================================'
write(*,'(A)') ' REAL GAS — EQUATIONS OF STATE ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Equation of State = ', trim(eos_name)
write(*,'(A,ES12.4,A)') ' Temperature T = ', T, ' K'
write(*,'(A,ES12.4,A)') ' Pressure P = ', P, ' Pa'
write(*,'(A,ES12.4,A)') ' Critical Temperature Tc = ', Tc, ' K'
write(*,'(A,ES12.4,A)') ' Critical Pressure Pc = ', Pc, ' Pa'
write(*,'(A,ES12.4)') ' Acentric Factor omega = ', omega
write(*,'(A,ES12.4,A)') ' Molar Mass = ', Mmol, ' g/mol'
write(*,'(A,ES12.4)') ' Reduced Temperature Tr = ', T/Tc
write(*,'(A,ES12.4)') ' Reduced Pressure Pr = ', P/Pc
write(*,*)
write(*,'(A)') '--- EOS PARAMETERS ------------------------------------------'
write(*,'(A,ES12.4)') ' a (attraction) = ', a_eos
write(*,'(A,ES12.4)') ' b (co-volume) = ', b_eos
write(*,*)
write(*,'(A)') '--- STATE RESULTS -------------------------------------------'
write(*,'(A,ES12.4,A)') ' Molar Volume v = ', v, ' m3/mol'
write(*,'(A,ES12.4,A)') ' Ideal Gas Volume = ', v_ideal, ' m3/mol'
write(*,'(A,ES12.6)') ' Compressibility Factor Z = ', Z
write(*,'(A,ES12.4,A)') ' Specific Volume = ', v/(Mmol*1.0d-3), ' m3/kg'
write(*,'(A,ES12.4,A)') ' Density = ', (Mmol*1.0d-3)/v, ' kg/m3'
write(*,*)
write(*,'(A)') '--- DEPARTURE FUNCTIONS (molar) -----------------------------'
write(*,'(A,ES12.4,A)') ' Enthalpy Departure dh = ', dh_dep, ' J/mol'
write(*,'(A,ES12.4,A)') ' Entropy Departure ds = ', ds_dep, ' J/(mol.K)'
write(*,'(A,ES12.4,A)') ' Gibbs Departure dg = ', dg_dep, ' J/mol'
write(*,'(A,ES12.4,A)') ' Fugacity = ', fugacity, ' Pa'
write(*,'(A,ES12.6)') ' Fugacity Coefficient phi = ', phi_fug
write(*,*)
! Z vs P sweep (isothermal)
write(*,'(A)') '--- Z VS PRESSURE SWEEP (isothermal) ------------------------'
write(*,'(A)') ' P[Pa] Z v[m3/mol] rho[kg/m3]'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, 60
P_i = P * 0.01d0 * (200.0d0)**(dble(i-1)/59.0d0)
select case(eos_type)
case(1); call solve_cubic_vdw(T, P_i, a_eos, b_eos, Ru, v_i)
case(2); call solve_cubic_rk(T, P_i, a_eos, b_eos, Ru, v_i)
case(3); call solve_cubic_pr(T, P_i, a_eos, b_eos, Ru, v_i)
end select
Z_i = P_i * v_i / (Ru * T)
write(*,'(ES12.4,2X,F12.6,2X,ES12.4,2X,ES12.4)') &
P_i, Z_i, v_i, (Mmol*1.0d-3)/v_i
end do
write(*,*)
! Z vs T sweep (isobaric)
write(*,'(A)') '--- Z VS TEMPERATURE SWEEP (isobaric) -----------------------'
write(*,'(A)') ' T[K] Z v[m3/mol] Tr'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, 50
T_i = Tc * 0.5d0 + (Tc*3.0d0 - Tc*0.5d0)*dble(i-1)/49.0d0
! Recompute a for PR (temperature-dependent alpha)
if (eos_type == 3) then
alpha_pr = (1.0d0 + kappa*(1.0d0 - sqrt(T_i/Tc)))**2
a_eos = 0.45724d0 * Ru**2 * Tc**2 / Pc * alpha_pr
end if
if (eos_type == 2) then
a_eos = 0.42748d0 * Ru**2 * Tc**2.5d0 / Pc
end if
select case(eos_type)
case(1); call solve_cubic_vdw(T_i, P, a_eos, b_eos, Ru, v_i)
case(2); call solve_cubic_rk(T_i, P, a_eos, b_eos, Ru, v_i)
case(3); call solve_cubic_pr(T_i, P, a_eos, b_eos, Ru, v_i)
end select
Z_i = P * v_i / (Ru * T_i)
write(*,'(F12.2,2X,F12.6,2X,ES12.4,2X,F10.4)') T_i, Z_i, v_i, T_i/Tc
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' vdW: P = RT/(v-b) - a/v^2.'
write(*,'(A)') ' RK: P = RT/(v-b) - a/(T^0.5 v(v+b)).'
write(*,'(A)') ' PR: P = RT/(v-b) - a alpha/(v(v+b)+b(v-b)).'
write(*,'(A)') ' Z = Pv/(RT).'
contains
subroutine solve_cubic_vdw(Tin, Pin, ain, bin, R, vout)
implicit none
double precision, intent(in) :: Tin, Pin, ain, bin, R
double precision, intent(out) :: vout
double precision :: v_guess, Pv, dPdv, dv
integer :: it
v_guess = R*Tin/Pin
do it = 1, 200
Pv = R*Tin/(v_guess-bin) - ain/v_guess**2
dPdv = -R*Tin/(v_guess-bin)**2 + 2.0d0*ain/v_guess**3
dv = -(Pv - Pin)/dPdv
v_guess = v_guess + dv
if (v_guess < bin*1.01d0) v_guess = bin*1.01d0
if (abs(dv) < 1.0d-15*abs(v_guess)) exit
end do
vout = v_guess
end subroutine solve_cubic_vdw
subroutine solve_cubic_rk(Tin, Pin, ain, bin, R, vout)
implicit none
double precision, intent(in) :: Tin, Pin, ain, bin, R
double precision, intent(out) :: vout
double precision :: v_guess, Pv, dPdv, dv, sqT
integer :: it
sqT = sqrt(Tin)
v_guess = R*Tin/Pin
do it = 1, 200
Pv = R*Tin/(v_guess-bin) - ain/(sqT*v_guess*(v_guess+bin))
dPdv = -R*Tin/(v_guess-bin)**2 + ain/(sqT)* &
(2.0d0*v_guess+bin)/(v_guess*(v_guess+bin))**2
dv = -(Pv - Pin)/dPdv
v_guess = v_guess + dv
if (v_guess < bin*1.01d0) v_guess = bin*1.01d0
if (abs(dv) < 1.0d-15*abs(v_guess)) exit
end do
vout = v_guess
end subroutine solve_cubic_rk
subroutine solve_cubic_pr(Tin, Pin, ain, bin, R, vout)
implicit none
double precision, intent(in) :: Tin, Pin, ain, bin, R
double precision, intent(out) :: vout
double precision :: v_guess, Pv, dPdv, dv, denom
integer :: it
v_guess = R*Tin/Pin
do it = 1, 200
denom = v_guess*(v_guess+bin)+bin*(v_guess-bin)
Pv = R*Tin/(v_guess-bin) - ain/denom
dPdv = -R*Tin/(v_guess-bin)**2 + ain*(2.0d0*v_guess+2.0d0*bin-bin)/ &
denom**2 * (2.0d0*v_guess+2.0d0*bin)
! Numerical derivative fallback
dPdv = -R*Tin/(v_guess-bin)**2 + &
ain*(2.0d0*v_guess)/(denom**2)
dv = -(Pv - Pin)/dPdv
if (abs(dv) > 0.5d0*v_guess) dv = sign(0.5d0*v_guess, dv)
v_guess = v_guess + dv
if (v_guess < bin*1.01d0) v_guess = bin*1.01d0
if (abs(dv) < 1.0d-15*abs(v_guess)) exit
end do
vout = v_guess
end subroutine solve_cubic_pr
subroutine departure_functions(etype, Tin, Pin, vin, ain, bin, R, Tcin, &
dh, ds, dg, fug, phi)
implicit none
integer, intent(in) :: etype
double precision, intent(in) :: Tin, Pin, vin, ain, bin, R, Tcin
double precision, intent(out) :: dh, ds, dg, fug, phi
double precision :: Zval, lnphi
Zval = Pin*vin/(R*Tin)
select case(etype)
case(1)
! van der Waals departure
dh = R*Tin*(Zval - 1.0d0) - 2.0d0*ain/vin + Pin*vin - R*Tin
dh = Pin*vin - R*Tin - ain/vin
ds = R*log(max((vin-bin)*Pin/(R*Tin), 1.0d-30))
lnphi = bin/(vin-bin) - 2.0d0*ain/(R*Tin*vin) - log(Zval)
case(2)
! Redlich-Kwong departure
dh = R*Tin*(Zval-1.0d0) - 1.5d0*ain/(bin*sqrt(Tin))* &
log((vin+bin)/vin)
ds = R*log(max(Zval*(1.0d0-bin/vin),1.0d-30)) - &
0.5d0*ain/(bin*Tin*sqrt(Tin))*log((vin+bin)/vin)
lnphi = Zval - 1.0d0 - log(max(Zval*(1.0d0-bin/vin),1.0d-30)) - &
ain/(bin*R*sqrt(Tin)*Tin)*log((vin+bin)/vin)
case(3)
! Peng-Robinson departure
dh = R*Tin*(Zval-1.0d0) + (Tin*0.0d0 - ain)/(2.0d0*sqrt(2.0d0)*bin)* &
log((vin+(1.0d0+sqrt(2.0d0))*bin)/(vin+(1.0d0-sqrt(2.0d0))*bin))
ds = R*log(max(Zval*(vin-bin)/vin, 1.0d-30))
lnphi = Zval - 1.0d0 - log(max(Zval-bin*Pin/(R*Tin),1.0d-30)) - &
ain/(2.0d0*sqrt(2.0d0)*bin*R*Tin)* &
log((vin+(1.0d0+sqrt(2.0d0))*bin)/ &
(vin+(1.0d0-sqrt(2.0d0))*bin))
end select
phi = exp(lnphi)
fug = phi * Pin
dg = R*Tin*lnphi
end subroutine departure_functions
end program real_gas_eos
Solver Description
Evaluates real gas thermodynamic properties using three cubic equations of state (van der Waals, Redlich-Kwong, and Peng-Robinson). Computes compressibility factor Z (solving the cubic cardano equation for vapor and liquid roots), molar volume, density, enthalpy departure ($\Delta h$), entropy departure ($\Delta s$), fugacity, and fugacity coefficient ($\phi$).
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 real_gas_eos.f90 -o real_gas_eos
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
real_gas_eos < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! EOS model (1=van der Waals, 2=Redlich-Kwong, 3=Peng-Robinson)
3
! Temperature T [K]
350.0
! Pressure P [Pa]
8.0e6
! Critical temperature Tc [K]
304.13
! Critical pressure Pc [Pa]
7.377e6
! Acentric factor omega
0.2236
! Molar mass [g/mol]
44.01
3
! Temperature T [K]
350.0
! Pressure P [Pa]
8.0e6
! Critical temperature Tc [K]
304.13
! Critical pressure Pc [Pa]
7.377e6
! Acentric factor omega
0.2236
! Molar mass [g/mol]
44.01