🔬
Solver Purpose & Physical Scope
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$).
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 359 times
📄 Source File:
real_gas_eos.f90
📁 calcul/Thermodynamics /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 real_gas_eos.f90 -o real_gas_eos
2. Execution with input.txt redirection:
real_gas_eos < input.txt
📄 Sample input.txt File Structure
Sample Data:
3 350.0 8.0e6 304.13 7.377e6 0.2236 44.01
Parameter Description:
EOS model (1=van der Waals, 2=Redlich-Kwong, 3=Peng-Robinson) Temperature T [K] Pressure P [Pa] Critical temperature Tc [K] Critical pressure Pc [Pa] Acentric factor omega Molar mass [g/mol]