求助：Sorption 计算吸附等温线 - 计算化学公社

################################################################################
#!perl #
# #
# PENG-ROBINSON EQUATION OF STATE #
# #
# Author: Reinier Akkermans (Biovia) #
# Version: 1.0 #
# Tested on: Materials Studio 7.0 #
# Required modules: Materials Visualizer #
# #
# Description: This script calculates the fugacity and pressure for a range of #
# volume using the Peng-Robinson Equation of state. The data is output in a #
# study table #
# #
################################################################################
use strict;
use Getopt::Long;
use MaterialsScript qw(:all);
use constant
{
GAS_CONSTANT => 8.314, # J/mol/K
# constants in the Peng-Robinson EoS
PR_a => 0.457235,
PR_b => 0.077796,
PR_kappa0 => 0.37464,
PR_kappa1 => 1.54226,
PR_kappa2 => -0.26992,
};
################################################################################
# BEGIN USER INPUT #
################################################################################
# example: CO2
my $Tc = 304.25; # temperature at critical point in K
my $Pc = 7380000; # pressure at critical point in Pa
my $omega = 0.228; # acentric factor
my $temperature = 298; # K
my $volume_start = 0.0004; # m^3/mol
my $volume_end = 0.1; # m^3/mol
my $volume_steps = 100;
################################################################################
# END USER INPUT #
################################################################################
my $Tr = $temperature/$Tc;
my $RTc = GAS_CONSTANT*$Tc; # J/mol
my $b = PR_b*$RTc/$Pc; # m^3/mol
my $kappa = PR_kappa0+PR_kappa1*$omega+PR_kappa2*$omega**2;
my $alpha = (1+$kappa*(1-sqrt($Tr)))**2;
my $aTc = PR_a*$RTc**2/$Pc; # Pa*(m^3/mol)^2
my $aT = $aTc*$alpha; # Pa*(m^3/mol)^2
my $RT = GAS_CONSTANT*$temperature;
my $std = Documents->New("PengRobinson.std");
$std->ColumnHeading(0) = "Volume (m^3/mol)";
$std->ColumnHeading(1) = "Pressure (kPa)";
$std->ColumnHeading(2) = "Fugacity (kPa)";
$std->ColumnHeading(3) = "A";
$std->ColumnHeading(4) = "B";
$std->ColumnHeading(5) = "Z";
$std->ColumnHeading(6) = "check";
my $scaleFactor = exp(log($volume_end/$volume_start)/($volume_steps-1));
my $volume = $volume_start;
for(my $i = 0; $i < $volume_steps; $i++)
{
$volume *= $scaleFactor if $i > 0;
my $pressure = $RT/($volume-$b)
-$aT/($volume*($b+$volume)+$b*($volume-$b));
my $A = $aT*$pressure/$RT**2;
my $B = $b*$pressure/$RT;
my $Z = $pressure*$volume/$RT;
# sanity check; should evaluate to zero
my $check = $Z**3+($B-1)*$Z**2+($A-2*$B-3*$B**2)*$Z+$B**3+$B**2-$A*$B;
my $lnfP = $Z-1
- log($Z-$B)
- $A/(sqrt(8)*$B)*log(($Z+(1+sqrt(2))*$B)/($Z-(1-sqrt(2))*$B));
my $fugacity = exp($lnfP)*$pressure;
$std->Cell($i,0) = $volume;
$std->Cell($i,1) = $pressure/1000;
$std->Cell($i,2) = $fugacity/1000;
$std->Cell($i,3) = $A;
$std->Cell($i,4) = $B;
$std->Cell($i,5) = $Z;
$std->Cell($i,6) = $check;
}

复制代码