💻 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 Compressibility Factor Z
Core Numerical Engine in Fortran 90 • 27 total downloads
compressibility_z.f90
! =========================================================================
! Source File: compressibility_z.f90
! =========================================================================
program compressibility_z
implicit none
integer::gas,iostat_val,i,ns,j
double precision::T,P,Tc,Pc,omega,Tc_c,Pc_c,om_c
double precision::Tr,Pr,Z_pit,Z_vdw,Z_pr,B0,B1
double precision::a_vdw,b_vdw,Ap,Bp,Z_iter
double precision::a_pr,b_pr,kappa,alpha_pr,A_pr_coeff,B_pr_coeff
double precision::P_sw,Z1,Z2,Z3
double precision,parameter::R=8.314462d0
character(len=30)::gname
read(*,*,iostat=iostat_val) gas; read(*,*,iostat=iostat_val) T
read(*,*,iostat=iostat_val) P; read(*,*,iostat=iostat_val) Tc_c
read(*,*,iostat=iostat_val) Pc_c; read(*,*,iostat=iostat_val) om_c
if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
select case(gas)
case(1);gname='Nitrogen (N2)';Tc=126.2d0;Pc=3.39d0;omega=0.037d0
case(2);gname='Oxygen (O2)';Tc=154.6d0;Pc=5.04d0;omega=0.022d0
case(3);gname='CO2';Tc=304.2d0;Pc=7.38d0;omega=0.224d0
case(4);gname='Methane (CH4)';Tc=190.6d0;Pc=4.60d0;omega=0.012d0
case(5);gname='Hydrogen (H2)';Tc=33.2d0;Pc=1.30d0;omega=-0.217d0
case(6);gname='Air';Tc=132.5d0;Pc=3.77d0;omega=0.035d0
case default;gname='Custom';Tc=Tc_c;Pc=Pc_c;omega=om_c;gas=7
end select
if(Tc<=0)Tc=300.0d0; if(Pc<=0)Pc=5.0d0; if(T<=0)T=300.0d0; if(P<=0)P=1.0d0
Tr=T/Tc; Pr=P/Pc
! Pitzer
B0=0.083d0-0.422d0/Tr**1.6d0
B1=0.139d0-0.172d0/Tr**4.2d0
Z_pit=1.0d0+(B0+omega*B1)*Pr/Tr
! Van der Waals
a_vdw=27.0d0/64.0d0*(R*Tc)**2/(Pc*1.0d6)
b_vdw=R*Tc/(8.0d0*Pc*1.0d6)
Ap=a_vdw*P*1.0d6/(R*T)**2
Bp=b_vdw*P*1.0d6/(R*T)
Z_vdw=1.0d0
do j=1,50; Z_vdw=1.0d0+Bp-Ap/(Z_vdw+Bp/Z_vdw); end do
if(Z_vdw<0.1d0)Z_vdw=0.1d0
! Peng-Robinson
kappa=0.37464d0+1.54226d0*omega-0.26992d0*omega**2
alpha_pr=(1.0d0+kappa*(1.0d0-sqrt(Tr)))**2
a_pr=0.45724d0*(R*Tc)**2/(Pc*1.0d6)*alpha_pr
b_pr=0.07780d0*R*Tc/(Pc*1.0d6)
A_pr_coeff=a_pr*P*1.0d6/(R*T)**2
B_pr_coeff=b_pr*P*1.0d6/(R*T)
Z_pr=1.0d0
do j=1,80
Z_pr=1.0d0+B_pr_coeff-(A_pr_coeff-B_pr_coeff-B_pr_coeff**2)/(Z_pr*(Z_pr+B_pr_coeff)+B_pr_coeff*(Z_pr-B_pr_coeff))
if(Z_pr<0.05d0)Z_pr=0.05d0
end do
write(*,'(A)')'============================================================'
write(*,'(A)')' COMPRESSIBILITY FACTOR (Z) CALCULATOR'
write(*,'(A)')'============================================================'
write(*,*)
write(*,'(A)')'--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' Gas = ',trim(gname)
write(*,'(A,F10.2,A)') ' Temperature T = ',T,' K'
write(*,'(A,F10.4,A)') ' Pressure P = ',P,' MPa'
write(*,'(A,F10.2,A)') ' Tc = ',Tc,' K'
write(*,'(A,F10.4,A)') ' Pc = ',Pc,' MPa'
write(*,'(A,F10.4)') ' Acentric Factor omega = ',omega
write(*,*)
write(*,'(A)')'--- REDUCED PROPERTIES --------------------------------------'
write(*,'(A,F10.4)') ' Tr = T/Tc = ',Tr
write(*,'(A,F10.4)') ' Pr = P/Pc = ',Pr
write(*,*)
write(*,'(A)')'--- COMPRESSIBILITY FACTOR Z --------------------------------'
write(*,'(A,F10.6)') ' Z (Ideal Gas) = ',1.0d0
write(*,'(A,F10.6)') ' Z (Pitzer Correlation) = ',Z_pit
write(*,'(A,F10.6)') ' Z (Van der Waals) = ',Z_vdw
write(*,'(A,F10.6)') ' Z (Peng-Robinson) = ',Z_pr
write(*,*)
ns=40
write(*,'(A)')'--- Z VS PRESSURE SWEEP -------------------------------------'
write(*,'(A)')' P[MPa] Z_Pitzer Z_VdW Z_PR'
write(*,'(A)')' -----------------------------------------------------------'
do i=1,ns
P_sw=0.1d0+dble(i-1)*(50.0d0-0.1d0)/dble(ns-1)
Pr=P_sw/Pc
B0=0.083d0-0.422d0/Tr**1.6d0;B1=0.139d0-0.172d0/Tr**4.2d0
Z1=1.0d0+(B0+omega*B1)*Pr/Tr
Ap=a_vdw*P_sw*1.0d6/(R*T)**2;Bp=b_vdw*P_sw*1.0d6/(R*T)
Z2=1.0d0;do j=1,50;Z2=1.0d0+Bp-Ap/(Z2+Bp/Z2);end do;if(Z2<0.05d0)Z2=0.05d0
A_PR=a_pr*P_sw*1.0d6/(R*T)**2;B_PR=b_pr*P_sw*1.0d6/(R*T)
Z3=1.0d0;do j=1,80;Z3=1.0d0+B_PR-(A_PR-B_PR-B_PR**2)/(Z3*(Z3+B_PR)+B_PR*(Z3-B_PR));if(Z3<0.05d0)Z3=0.05d0;end do
write(*,'(F8.3,4X,F10.6,2X,F10.6,2X,F10.6)') P_sw,Z1,Z2,Z3
end do
write(*,*)
write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)')' Pitzer: Z=1+(B0+omega*B1)*Pr/Tr'
write(*,'(A)')' Van der Waals: iterative cubic solution'
write(*,'(A)')' Peng-Robinson: iterative cubic with kappa correlation'
end program compressibility_z
Solver Description
Calculates real gas compressibility factor Z using Pitzer generalized correlations, Van der Waals, and Peng-Robinson cubic equations of state.
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 compressibility_z.f90 -o compressibility_z
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
compressibility_z < 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=Air, 7=Custom)
1
! Temperature T [K]
300.0
! Pressure P [MPa]
10.0
! Custom critical temperature Tc [K]
0.0
! Custom critical pressure Pc [MPa]
0.0
! Custom acentric factor omega
0.0
1
! Temperature T [K]
300.0
! Pressure P [MPa]
10.0
! Custom critical temperature Tc [K]
0.0
! Custom critical pressure Pc [MPa]
0.0
! Custom acentric factor omega
0.0