#pragma TextEncoding="UTF-8" #pragma rtGlobals=3 #pragma version=2.24 #pragma ModuleName=waterEOS // updater headers... maybe this will one day be cleaned up and released //static constant kProjectID=NaN // the project node on IgorExchange //static strconstant ksShortTitle="WaterEOS" // the project short title on IgorExchange // Note: you must choose sensible bracketing values for the property you want to calculate! // Some EOS behave more pathologigically than others and will require more careful choice of // bracketing values. // Bracketing values may be ignored in fugacity calculation. Sorry. // 2.20 March 2019 // fixed some bugs. // 2.10, 8/17/06 // added Zhang and Duan EOS // note: EOS_calculate_fugacity function must be updated when a new EOS is added to calculator. menu "Macros" "Water EOS calculator", EOScalculator() end function EOScalculator() EOS_init() EOS_makePanel() EOS_recalculate() end static function EOS_makePanel() wave EOSsettings=root:Packages:EOS:EOSsettings wave /T EOStextWave=root:Packages:EOS:EOStextWave DoWindow /K EOSpanel NewPanel /N=EOSpanel /W=(70,105,345,342) /K=1 as "High P EOS calculator for water" TabControl tabCalc,pos={0,5},size={275,20},tabLabel(0)="Calculator" TabControl tabCalc,tabLabel(1)="Units",tabLabel(2)="Limits",tabLabel(3)="Reference",value= 0 TabControl tabCalc,Proc=waterEOS#EOS_TabProc PopupMenu popupEOS,pos={40,40},size={200,21},title="Equation of state", Proc=waterEOS#EOSPopMenuProc PopupMenu popupEOS,mode=1,popvalue=EOStextWave[0][0],value=#"root:Packages:EOS:EOSlist", fsize=12 PopupMenu popupEOS, bodyWidth=125, disable=0 string unit switch(EOSsettings[%Pfactor]) case 1: unit="GPa" break case 10: unit="kbar" break case 1000: unit="MPa" break case 10000: unit="bar" break case 1e6: unit="kPa" break case 1e9: unit="Pa" break endswitch SetVariable setvarP,pos={75,75},size={150,16},title="Pressure ("+unit+")", Proc=waterEOS#EOSSetVarProc SetVariable setvarP,value= root:Packages:EOS:varP, noedit=0, limits={0,Inf,1 }, disable=0 switch(EOSsettings[%Toff]) case 0: unit="°C" break case 273: unit="K" break endswitch SetVariable setvarT,pos={75,100},size={150,16},title="Temperature ("+unit+")", Proc=waterEOS#EOSSetVarProc SetVariable setvarT,value= root:Packages:EOS:varT, noedit=0, limits={0,Inf,1 }, disable=0 switch(10*EOSsettings[%Vfactor]) case 10: unit="cc/mol" break case 1: unit="J/bar" break endswitch SetVariable setvarV,pos={75,125},size={150,16},title="Volume ("+unit+")", Proc=waterEOS#EOSSetVarProc SetVariable setvarV,value= root:Packages:EOS:varV, noedit=1, limits={0,Inf,0 }, disable=0 switch(EOSsettings[%Dfactor]) case 1: unit="g/cc" break case 1000: unit="g/l" break endswitch SetVariable setvarD,pos={75,150},size={150,16},title="Density ("+unit+") ", Proc=waterEOS#EOSSetVarProc SetVariable setvarD,value= root:Packages:EOS:varD, noedit=1, limits={0,Inf,0 }, disable=0 switch(EOSsettings[%Ffactor]) case 1: unit="GPa" break case 10: unit="kbar" break case 1000: unit="MPa" break case 10000: unit="bar" break case 1e6: unit="kPa" break case 1e9: unit="Pa" break endswitch SetVariable setvarF,pos={75,175},size={150,16},title="Fugacity ("+unit+")" SetVariable setvarF,value= root:Packages:EOS:varF, noedit=1, limits={0,50,0 }, bodyWidth=75 Button CalcFugButton, win=EOSpanel, disable=1, title="Calculate", Proc=waterEOS#EOS_calculate_fugacity, pos={150,175},size={75,16} CheckBox checkP,pos={50,75},size={16,14},title="",value= 0,mode=1, Proc=waterEOS#EOS_CheckProc, disable=0 CheckBox checkT,pos={50,100},size={16,14},title="",value= 0,mode=1, Proc=waterEOS#EOS_CheckProc, disable=0 CheckBox checkV,pos={50,125},size={16,14},title="",value= 1,mode=1, Proc=waterEOS#EOS_CheckProc, disable=0 GroupBox GroupUnits pos={30,50},size={215,170}, title="Select Units", disable=1 PopupMenu popupPunit pos={75,75 }, title="Pressure", value="Pa;kPa;bar;MPa;kbar;GPa", bodyWidth=60, size={120,16 }, disable=1, Proc=waterEOS#EOS_UnitsProc PopupMenu popupPunit mode=1+(EOSsettings[%Pfactor]<1e9)+(EOSsettings[%Pfactor]<1e6)+(EOSsettings[%Pfactor]<1e4)+(EOSsettings[%Pfactor]<1e3)+(EOSsettings[%Pfactor]<1e1) PopupMenu popupTunit pos={75,100 }, title="Temperature", value="°C;K", bodyWidth=60, size={120,16 }, disable=1, Proc=waterEOS#EOS_UnitsProc PopupMenu popupTunit mode =1+(EOSsettings[%Toff]>0) PopupMenu popupVunit pos={75,125 }, title="Volume", value="cc/mol;J/bar", bodyWidth=60, size={120,16 }, disable=1, Proc=waterEOS#EOS_UnitsProc PopupMenu popupVunit mode =1+(EOSsettings[%Vfactor]<1) PopupMenu popupDunit pos={75,150 }, title="Density", value="g/cc;g/l", bodyWidth=60, size={120,16 }, disable=1, Proc=waterEOS#EOS_UnitsProc PopupMenu popupDunit mode =1+(EOSsettings[%Dfactor]>1) PopupMenu popupFunit pos={75,175 }, title="Fugacity", value="Pa;kPa;bar;MPa;kbar;GPa;", bodyWidth=60, size={120,16 }, disable=1, Proc=waterEOS#EOS_UnitsProc PopupMenu popupFunit mode=1+(EOSsettings[%Pfactor]<1e9)+(EOSsettings[%Pfactor]<1e6)+(EOSsettings[%Pfactor]<1e4)+(EOSsettings[%Pfactor]<1e3)+(EOSsettings[%Pfactor]<1e1) GroupBox GroupLimits pos={30,50},size={215,170}, title="Settings for root finding", disable=1 PopupMenu popupLim pos={75,75 }, title="Root:", value="Volume;Pressure;Temperature", mode=1, bodyWidth=100, size={150,16 }, disable=1, Proc=waterEOS#EOS_PopLimProc SetVariable setvarlow pos={75,120 },size={150,1 }, title="High bracket (cc/mol)", value=root:Packages:EOS:EOSsettings[%Vhigh], disable=1, bodyWidth=50 SetVariable setvarhigh pos={75,145 },size={150,1 }, title="Low bracket (cc/mol)", value=root:Packages:EOS:EOSsettings[%Vlow], disable=1, bodyWidth=50 SetVariable setvartol pos={75,170 },size={150,1 }, title="Tolerance (cc/mol)", value=root:Packages:EOS:EOSsettings[%Vtol], disable=1, bodyWidth=50 TitleBox title0,pos={137.5,100},size={0,0},variable= root:Packages:EOS:ref, disable=1, anchor=mc ModifyPanel /W=EOSpanel, fixedSize= 1 end static function EOS_PopLimProc(s) : PopupMenuControl STRUCT WMPopupAction &s if(s.eventcode!=2) return 0 endif strswitch (s.popStr) case "Volume": SetVariable setvarlow win=EOSpanel, title="High bracket (cc/mol)", value=root:Packages:EOS:EOSsettings[%Vhigh] SetVariable setvarhigh win=EOSpanel, title="Low bracket (cc/mol)", value=root:Packages:EOS:EOSsettings[%Vlow] SetVariable setvartol win=EOSpanel, title="Tolerance (cc/mol)", value=root:Packages:EOS:EOSsettings[%Vtol] break case "Pressure": SetVariable setvarlow win=EOSpanel, title="High bracket (bar)", value=root:Packages:EOS:EOSsettings[%Phigh] SetVariable setvarhigh win=EOSpanel, title="Low bracket (bar)", value=root:Packages:EOS:EOSsettings[%Plow] SetVariable setvartol win=EOSpanel, title="Tolerance (bar)", value=root:Packages:EOS:EOSsettings[%Ptol] break case "Temperature": SetVariable setvarlow win=EOSpanel, title="High bracket (K)", value=root:Packages:EOS:EOSsettings[%Thigh] SetVariable setvarhigh win=EOSpanel, title="Low bracket (K)", value=root:Packages:EOS:EOSsettings[%Tlow] SetVariable setvartol win=EOSpanel, title="Tolerance (K)", value=root:Packages:EOS:EOSsettings[%Ttol] break endswitch end static function EOS_UnitsProc(STRUCT WMPopupAction &s) : PopupMenuControl if(s.eventcode!=2) return 0 endif wave EOSsettings=root:Packages:EOS:EOSsettings strswitch (s.ctrlName) case "popupPunit": strswitch (s.popStr) case "Pa": EOSsettings[%Pfactor]=1e9 break case "kPa": EOSsettings[%Pfactor]=1e6 break case "bar": EOSsettings[%Pfactor]=1e4 break case "MPa": EOSsettings[%Pfactor]=1e3 break case "kbar": EOSsettings[%Pfactor]=1e1 break case "GPa": EOSsettings[%Pfactor]=1 break endswitch break case "popupTunit": strswitch (s.popStr) case "°C": EOSsettings[%Toff]=0 break case "K": EOSsettings[%Toff]=273 break endswitch break case "popupVunit": strswitch (s.popStr) case "cc/mol": EOSsettings[%Vfactor]=1 break case "J/bar": EOSsettings[%Vfactor]=0.1 break endswitch break case "popupDunit": strswitch (s.popStr) case "g/cc": EOSsettings[%Dfactor]=1 break case "g/l": EOSsettings[%Dfactor]=1000 break endswitch break case "popupFunit": strswitch (s.popStr) case "GPa": EOSsettings[%Ffactor]=1 break case "kbar": EOSsettings[%Ffactor]=10 break case "MPa": EOSsettings[%Ffactor]=1000 break case "bar": EOSsettings[%Ffactor]=10000 break case "kPa": EOSsettings[%Ffactor]=1e6 break case "Pa": EOSsettings[%Ffactor]=1e9 break endswitch break endswitch end static function EOS_init() string savDF= GetDataFolder(1) // Save current DF for restore. NewDataFolder/O root:Packages NewDataFolder/O/S root:Packages:EOS // Make sure this exists. variable /G varP, varT, varV, varD, varF Make /O/T/N=(14,4) EOStextWave="" // EOStextWave[][2] contains names of functions that have the form P=f(V,T) // EOStextWave[][3] contains names of functions that are explicit in P: V=f(P,T) // order of rows must correspond to list in getVol function EOStextWave[0][0]="Brodholt and Wood" EOStextWave[0][1]="Brodholt, J.P. and Wood, B.J. (1993)\rSimulations of the structure and thermodynamic \r" EOStextWave[0][1]+="properties of water at high pressures and temperatures.\rJournal of Geophysical Research Solid Earth 98: 519-\r536." EOStextWave[0][2]="EOS_BrodholtWood" EOStextWave[0][3]="" EOStextWave[1][0]="Pitzer and Sterner" EOStextWave[1][1]="Pitzer, K.S. and Sterner, S.M. (1994).\rEquations of state valid continuously from zero to \r" EOStextWave[1][1]+="extreme pressures for H2O and CO2.\rJournal of Chemical Physics. 101: 3111-3116." EOStextWave[1][2]="EOS_PitzerSterner" EOStextWave[1][3]="" EOStextWave[2][0]="Holloway" EOStextWave[2][1]="Holloway, J.R. (1977)\rFugacity and activity of molecular species in super-\r" EOStextWave[2][1]+="critical fluids. \rIn:Thermodynamics in geology, pp 161-181." EOStextWave[2][2]="EOS_Holloway" EOStextWave[2][3]="" EOStextWave[3][0]="Saul and Wagner" EOStextWave[3][1]="Saul, A. and Wagner, W. (1989).\r" EOStextWave[3][1]+="A fundamental equation for water covering the range \r" EOStextWave[3][1]+="from the melting line to 1273 K at pressures up to \r" EOStextWave[3][1]+="25000 MPa.\r" EOStextWave[3][1]+="J. Phys. Chem. Ref. Data 18:1537-1564" EOStextWave[3][2]="EOS_SaulWagner" EOStextWave[3][3]="" EOStextWave[4][0]="Belonoshko and Saxena" EOStextWave[4][1]="Belonoshko, A. and Saxena, S.K. (1991).\rA molecular dynamics study of the pressure-volume- \r" EOStextWave[4][1]+="temperature properties of supercritical fluids: I. H2O.\rGeochimica et Cosmochimica Acta 55: 381-387." EOStextWave[4][2]="EOS_BelonoshkoSaxena" EOStextWave[4][3]="" EOStextWave[5][0]="Sakane" EOStextWave[5][1]="Sakane, S., Liu, W.B., Doren, D.J., Shock, E.L. and \rWood, R.H. (2001). \r" EOStextWave[5][1]+="Prediction of the Gibbs energies and an improved \r" EOStextWave[5][1]+="equation of state for water at extreme conditions \r" EOStextWave[5][1]+="from ab initio energies with classical simulations. \r" EOStextWave[5][1]+="Geochimica et Cosmochimica Acta, 65, 4067-4075." EOStextWave[5][2]="EOS_Sakane" EOStextWave[5][3]="" EOStextWave[6][0]="Zhang and Duan 05" EOStextWave[6][1]="Zhang and Duan (2005) Prediction of the PVT properties \r" EOStextWave[6][1]+="of water over wide range of temperatures and pressures \r" EOStextWave[6][1]+="from molecular dynamics simulation \r" EOStextWave[6][1]+="Physics of the Earth and Planetary Interiors, 149, 335-354" EOStextWave[6][2]="EOS_ZhangDuan05" EOStextWave[6][3]="" EOStextWave[7][0]="IAPWS95" EOStextWave[7][1]="The International Association for the \r" EOStextWave[7][1]+="Properties of Water and Steam \r" EOStextWave[7][1]+="Coefficients are from IAPWS R6-95(2018) " EOStextWave[7][2]="EOS_IAPWS95" EOStextWave[7][3]="" EOStextWave[8][0]="Zhang and Duan 09" EOStextWave[8][1]="Zhang and Duan (2009) A model for \r" EOStextWave[8][1]+="C-O-H fluid in the Earth's mantle\r" EOStextWave[8][1]+="\r" EOStextWave[8][1]+="" EOStextWave[8][2]="EOS_ZhangDuan09" EOStextWave[8][3]="" EOStextWave[9][0]="Holloway and Blank 94" EOStextWave[9][1]=" \r" EOStextWave[9][1]+="\r" EOStextWave[9][1]+="\r" EOStextWave[9][1]+="" EOStextWave[9][2]="EOS_Holloway94" EOStextWave[9][3]="" // edit edit calc fugacity function for cutoff row for pressure explicit eos EOStextWave[10][0]="Saxena and Fei" EOStextWave[10][1]="Saxena, S.K. and Fei, Y. (1997).\r" EOStextWave[10][1]+="High pressure and high temperature fluid \r" EOStextWave[10][1]+="fugacities. Geochimica et Cosmochimica \r" EOStextWave[10][1]+="Acta, 51, 793-791. " EOStextWave[10][2]="" EOStextWave[10][3]="EOS_V_SaxenaFei" EOStextWave[11][0]="Frost and Wood" EOStextWave[11][1]="Frost, D.J. and Wood, B.J. (1997). \r" EOStextWave[11][1]+="Experimental measurements of the properties of H2O- \r" EOStextWave[11][1]+="CO2 mixtures at high pressures and temperatures \r" EOStextWave[11][1]+="Geochimica et Cosmochimica Acta, 61, 3301-33010." EOStextWave[11][2]="" EOStextWave[11][3]="EOS_V_FrostWood" EOStextWave[12][0]="Holland and Powell 98" EOStextWave[12][1]="An internally consistent thermodynamic data\r" EOStextWave[12][1]+="set for phases of petrological interest\r" EOStextWave[12][1]+="J. metamorphic Geol., 16, 309–343\r" EOStextWave[12][1]+=" ** valid for supercritical conditions only! **" EOStextWave[12][2]="" EOStextWave[12][3]="EOS_V_HollandPowell98" EOStextWave[13][0]="Holland and Powell 91" EOStextWave[13][1]="A Compensated-Redlich-Kwong (CORK) equation for\r" EOStextWave[13][1]+="volumes and fugacities of CO2 and H20 in the \r" EOStextWave[13][1]+="range 1 bar to 50 kbar and 100-1600 °C\r" EOStextWave[13][1]+=" ** valid for supercritical conditions only! **" EOStextWave[13][2]="" EOStextWave[13][3]="EOS_V_HollandPowell91" if (!exists("EOSsettings")==1) Make /n=15 EOSsettings=1 SetDimLabel 0,0,Pfactor,EOSsettings SetDimLabel 0,1,Toff,EOSsettings SetDimLabel 0,2,Vfactor,EOSsettings SetDimLabel 0,3,Dfactor,EOSsettings SetDimLabel 0,4,Ffactor,EOSsettings SetDimLabel 0,5,Vhigh,EOSsettings SetDimLabel 0,6,Vlow,EOSsettings SetDimLabel 0,7,Vtol,EOSsettings SetDimLabel 0,8,Fres,EOSsettings SetDimLabel 0,9,Thigh,EOSsettings SetDimLabel 0,10,Tlow,EOSsettings SetDimLabel 0,11,Ttol,EOSsettings SetDimLabel 0,12,Phigh,EOSsettings SetDimLabel 0,13,Plow,EOSsettings SetDimLabel 0,14,Ptol,EOSsettings EOSsettings[%Toff]=0 EOSsettings[%Vhigh]=20 EOSsettings[%Vlow]=18.1 EOSsettings[%Vtol]=0.001 EOSsettings[%Fres]=1000 EOSsettings[%Thigh]=2000 EOSsettings[%Tlow]=500 EOSsettings[%Ttol]=0.01 EOSsettings[%Phigh]=100000 EOSsettings[%Plow]=1 EOSsettings[%Ptol]=0.1 endif string /G ref=EOStextWave[0][1] string /G EOSlist="" variable i for(i=0;i6) return 0 endif NVAR V=root:Packages:EOS:varV NVAR D=root:Packages:EOS:varD wave settings=root:Packages:EOS:EOSsettings if (stringmatch(s.ctrlName, "setvarV")) D=18.0153/(V*settings[%Vfactor])/settings[%Dfactor] endif if (stringmatch(s.ctrlName, "setvarD")) V=18.0153/(D*settings[%Dfactor])/settings[%Vfactor] endif EOS_recalculate() end static function EOS_recalculate() wave /T EOStextWave=root:Packages:EOS:EOStextWave wave settings=root:Packages:EOS:EOSsettings NVAR P=root:Packages:EOS:varP NVAR T=root:Packages:EOS:varT NVAR V=root:Packages:EOS:varV NVAR D=root:Packages:EOS:varD NVAR F=root:Packages:EOS:varF Make /free/n=1 w ControlInfo /W=EOSpanel checkV if (V_value) ControlInfo /W=EOSpanel popupEOS if (strlen(EOStextWave[V_Value-1][2])) FUNCREF EOStemplateFunction EOSfunc= $EOStextWave[V_Value-1][2]+"_V" w=T+273.15*(settings[%Toff]==0) FindRoots /H=(settings[%Vhigh])/L=(settings[%Vlow])/Q/T=(settings[%Vtol])/Z=(P/settings[%Pfactor]*10000) EOSfunc, w V = V_Flag ? NaN : V_Root else FUNCREF EOStemplateFunction2 ff= $EOStextWave[V_Value-1][3] V=ff(P/settings[%Pfactor]*10000,T+273.15*(settings[%Toff]==0)) endif D=18.0153/(V*settings[%Vfactor])/settings[%Dfactor] endif ControlInfo /W=EOSpanel checkT if (V_value) ControlInfo /W=EOSpanel popupEOS // function names are in different columns of EOStextWave depending on EOS format if (strlen(EOStextWave[V_Value-1][2])) // EOS calculates P(V,T), find T that satisfies P(V,T) FUNCREF EOStemplateFunction EOSfunc= $EOStextWave[V_Value-1][2]+"_T" w=V/settings[%Vfactor] FindRoots /H=(settings[%Thigh])/L=(settings[%Tlow])/Q/T=(settings[%Ttol])/Z=(P/settings[%Pfactor]*10000) EOSfunc, w T = V_Flag ? NaN : V_Root-273.15*(settings[%Toff]==0) else // EOS is explicit in P, calculates V(P,T) FUNCREF EOStemplateFunction EOSfunc= $EOStextWave[V_Value-1][3]+"_T" w=P/settings[%Pfactor]*10000 FindRoots /H=(settings[%Thigh])/L=(settings[%Tlow])/Q/T=(settings[%Ttol])/Z=(V/settings[%Vfactor]) EOSfunc, w T = V_Flag ? NaN : V_Root-273.15*(settings[%Toff]==0) endif endif ControlInfo /W=EOSpanel checkP if (V_value) ControlInfo /W=EOSpanel popupEOS // function names are in different columns of EOStextWave depending on EOS format if (strlen(EOStextWave[V_Value-1][2])) // EOS calculates P(V,T) ControlInfo /W=EOSpanel popupEOS FUNCREF EOStemplateFunction2 ff= $EOStextWave[V_Value-1][2] P=ff(T+273.15*(settings[%Toff]==0),V/settings[%Vfactor])/10000*settings[%Pfactor] else // EOS is explicit in P, so find P that satisfies V(P,T) FUNCREF EOStemplateFunction EOSfunc= $EOStextWave[V_Value-1][3]+"_P" w=T+273.15*(settings[%Toff]==0) FindRoots /H=(settings[%Phigh])/L=(settings[%Plow])/Q/T=(settings[%Ptol])/Z=(V/settings[%Vfactor]) EOSfunc, w P = V_Flag ? NaN : V_Root*(settings[%Pfactor])/10000 endif endif // (sometimes) recalculate fugacity ControlInfo /W=EOSpanel popupEOS // string strFugacityFunction = replacestring(" and ", S_value, "") // strFugacityFunction = replacestring(" ", strFugacityFunction, "") + "Fugacity" string strFugacityFunction = EOStextWave[V_Value-1][2] if (strlen(strFugacityFunction)==0) strFugacityFunction = EOStextWave[V_Value-1][3] endif strFugacityFunction = replacestring("EOS_V_", strFugacityFunction, "") strFugacityFunction = replacestring("EOS_", strFugacityFunction, "") strFugacityFunction += "Fugacity" if(exists(strFugacityFunction)==6) // must be threadsafe!!!! FUNCREF EOStemplateFunction2 ff= $strFugacityFunction F=ff(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] Button CalcFugButton, win=EOSpanel, disable=1 else Button CalcFugButton, win=EOSpanel, disable=0 endif // strswitch (S_value) // // // case "Holloway": // these EOS have an analytical solution (sort of) for fugacity // Button CalcFugButton, win=EOSpanel, disable=1 // F=HollowayFugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Belonoshko and Saxena": // Button CalcFugButton, win=EOSpanel, disable=1 // F=BelonoshkoSaxenaFugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Pitzer and Sterner": // Button CalcFugButton, win=EOSpanel, disable=1 // F=PitzerSternerFugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Sakane": // Button CalcFugButton, win=EOSpanel, disable=1 // F=SakaneFugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "IAPWS95": // Button CalcFugButton, win=EOSpanel, disable=1 // F=IAPWS95Fugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Saul and Wagner": // Button CalcFugButton, win=EOSpanel, disable=1 // F=SaulWagnerFugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Frost and Wood": // Button CalcFugButton, win=EOSpanel, disable=1 // F=FrostWoodFugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Zhang and Duan 09": // Button CalcFugButton, win=EOSpanel, disable=1 // F=ZhangDuan09Fugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Holland and Powell 91": // Button CalcFugButton, win=EOSpanel, disable=1 // F=HollandPowell91Fugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // case "Holland and Powell 98": // Button CalcFugButton, win=EOSpanel, disable=1 // F=HollandPowell98Fugacity(P/settings[%Pfactor], T-273.25*(settings[%Toff]>0))*settings[%Ffactor] // break // default: // Button CalcFugButton, win=EOSpanel, disable=0 // endswitch end static function EOS_calculate_fugacity(STRUCT WMButtonAction &s) : ButtonControl if(s.eventcode!=2) return 0 endif NVAR F=root:Packages:EOS:varF NVAR T=root:Packages:EOS:varT NVAR pressure=root:Packages:EOS:varP wave settings=root:Packages:EOS:EOSsettings wave /T EOStextWave=root:Packages:EOS:EOStextWave Button CalcFugButton, win=EOSpanel, title="Calculating...", disable=2 DoUpdate /W=$s.win // ---- do numerical integration variable TK, Pbar, Rgas, Vint, P0 TK=T+273.15*(settings[%Toff]==0) Pbar=pressure/settings[%Pfactor]*1e4 // bar Rgas=83.14510 // cc bar/K P0=1e-4 // GPa Make /free/n=2 pwave pwave[0]=TK ControlInfo /W=EOSpanel popupEOS // pwave[1] = V_Value>9 ? V_Value -1 + 1 : V_Value-1// EDIT 2 numbers HERE pwave[1] = V_Value-1 // 0-9, 10 - ? if (stringmatch(S_Value, "Saxena*")) // don't attempt to calculate Button CalcFugButton, win=EOSpanel, title="Calculate", disable=0 return 0 endif string strFugacityFunc=ReplaceString(" and ", S_Value, "") strFugacityFunc=ReplaceString(" ", strFugacityFunc, "")+"Fugacity" if(exists(strFugacityFunc)==6) FUNCREF EOStemplateFunction2 ff= $strFugacityFunc F=ff(Pbar/1e4, TK-273.15)*settings[%Ffactor] else // this is going to be slow. Very slow for some EOS! Vint=Integrate1D(waterEOS#EOS_GetVol,1, Pbar, 2, 0, pwave) F=P0*exp(Vint/Rgas/TK)*settings[%Ffactor] endif Button CalcFugButton, win=EOSpanel, title="Calculate", disable=1 end // wrapper function for fugacity calculation // this is a bad strategy - this function is called by Integrate1D. // The extra overhead of textwave and funcref really slows this down. static function EOS_GetVol(wave pwave, variable Pbars) // pwave[0] is temperature in K // pwave[1] tells us which EOS wave settings=root:Packages:EOS:EOSsettings Make /free/T/n=20 EOSfunctions EOSfunctions[0]={"EOS_BrodholtWood_V","EOS_PitzerSterner_V","EOS_Holloway_V","EOS_SaulWagner_V","EOS_BelonoshkoSaxena_V","EOS_Sakane_V","EOS_ZhangDuan05_V","EOS_IAPWS95_V","EOS_ZhangDuan09_V","EOS_Holloway94_V"} EOSfunctions[10]={"EOS_V_SaxenaFei","EOS_V_FrostWood","EOS_V_HollandPowell98","EOS_V_HollandPowell91"} if(pwave[1]>9) // EOS explicit FUNCREF EOStemplateFunction2 ff= $EOSfunctions[pwave[1]] return ff(Pbars,pwave[0]) // cc/mol else FUNCREF EOStemplateFunction EOSfunc=$EOSfunctions[pwave[1]] // find V that satisfies P(V,T) FindRoots /H=(settings[%Vhigh])/L=(settings[%Vlow])/Q/T=(settings[%Vtol])/Z=(Pbars) EOSfunc, pwave // EOSfunc will find T in pwave[0] return V_Flag ? NaN : V_Root // cc/mol endif end static function EOSPopMenuProc(STRUCT WMPopupAction &s) if(s.eventcode!=2) return 0 endif SVAR ref=root:Packages:EOS:ref wave /T EOStextWave=root:Packages:EOS:EOStextWave ref=EOStextWave[s.popNum-1][1] EOS_recalculate() // recalculate whatever needs calculating end static function EOS_CheckProc(STRUCT WMCheckboxAction &s) if(s.eventcode!=2) return 0 endif strswitch(s.ctrlName) case "checkP": CheckBox checkT,win=EOSpanel, value=0 CheckBox checkV,win=EOSpanel, value=0 SetVariable setvarP,win=EOSpanel, noedit=1, limits={0,Inf,0} SetVariable setvarT,win=EOSpanel, noedit=0, limits={0,Inf,1} SetVariable setvarV,win=EOSpanel, noedit=0, limits={0,Inf,0.1} SetVariable setvarD,win=EOSpanel, noedit=0, limits={0,Inf,0.1} break case "checkT": CheckBox checkP,win=EOSpanel, value=0 CheckBox checkV,win=EOSpanel, value=0 SetVariable setvarT,win=EOSpanel, noedit=1, limits={0,Inf,0} SetVariable setvarP,win=EOSpanel, noedit=0, limits={0,Inf,1} SetVariable setvarV,win=EOSpanel, noedit=0, limits={0,Inf,0.1} SetVariable setvarD,win=EOSpanel, noedit=0, limits={0,Inf,0.1} break case "checkV": CheckBox checkP,win=EOSpanel, value=0 CheckBox checkT,win=EOSpanel, value=0 SetVariable setvarV,win=EOSpanel, noedit=1, limits={0,Inf,0} SetVariable setvarD,win=EOSpanel, noedit=1, limits={0,Inf,0} SetVariable setvarP,win=EOSpanel, noedit=0, limits={0,Inf,1} SetVariable setvarT,win=EOSpanel, noedit=0, limits={0,Inf,1} break endswitch end // These templates allow selection of an EOS function to feed to findroots. // They describe the wrapper functions that have a format useable by findroots. // This template is for EOS functions that calculate P as a function of // V,T. The wrapper function can take either P or T in the parameter wave // w, and V as a variable. threadsafe function EOStemplateFunction(wave w, variable A) return 0 end // This template is for EOS functions that calculate V as a function of P,T threadsafe function EOStemplateFunction2(variable A, variable B) return 0 end static function EOS_TabProc(STRUCT WMTabControlAction &s) if(s.eventcode!=2) return 0 endif wave EOSsettings=root:Packages:EOS:EOSsettings PopupMenu popupEOS, win=EOSpanel, disable=(s.tab!=0) string unit switch(EOSsettings[%Pfactor]) case 1: unit="GPa" break case 10: unit="kbar" break case 1000: unit="MPa" break case 10000: unit="bar" break case 1e6: unit="kPa" break case 1e9: unit="Pa" break endswitch SetVariable setvarP,win=EOSpanel, title="Pressure ("+unit+")", disable=(s.tab!=0) switch(EOSsettings[%Toff]) case 0: unit="°C" break case 273: unit="K" break endswitch SetVariable setvarT,win=EOSpanel,title="Temperature ("+unit+")", disable=(s.tab!=0) switch(10*EOSsettings[%Vfactor]) case 10: unit="cc/mol" break case 1: unit="J/bar" break endswitch SetVariable setvarV,win=EOSpanel,title="Volume ("+unit+")", disable=(s.tab!=0) switch(EOSsettings[%Dfactor]) case 1: unit="g/cc" break case 1000: unit="g/l" break endswitch SetVariable setvarD,win=EOSpanel,title="Density ("+unit+") ", disable=(s.tab!=0) switch(EOSsettings[%Ffactor]) case 1: unit="GPa" break case 10: unit="kbar" break case 1000: unit="MPa" break case 10000: unit="bar" break case 1e6: unit="kPa" break case 1e9: unit="Pa" break endswitch ControlInfo /W=EOSpanel popupEOS SetVariable setvarF,win=EOSpanel,title="Fugacity ("+unit+")", disable=(s.tab!=0) Button CalcFugButton, win=EOSpanel, disable=(s.tab!=0) if (s.tab==0) if(stringmatch(S_Value, "Pitzer and Sterner")||stringmatch(S_Value, "Sakane")) Button CalcFugButton, win=EOSpanel, disable=1 else Button CalcFugButton, win=EOSpanel, disable=0 endif EOS_recalculate() endif CheckBox checkP, win=EOSpanel, disable=(s.tab!=0) CheckBox checkT, win=EOSpanel, disable=(s.tab!=0) CheckBox checkV, win=EOSpanel, disable=(s.tab!=0) GroupBox GroupUnits, win=EOSpanel, disable=(s.tab!=1) PopupMenu popupPunit, win=EOSpanel, disable=(s.tab!=1) PopupMenu popupTunit, win=EOSpanel, disable=(s.tab!=1) PopupMenu popupVunit, win=EOSpanel, disable=(s.tab!=1) PopupMenu popupDunit, win=EOSpanel, disable=(s.tab!=1) PopupMenu popupFunit, win=EOSpanel, disable=(s.tab!=1) GroupBox GroupLimits , win=EOSpanel, disable=(s.tab!=2) PopupMenu popupLim win=EOSpanel, disable=(s.tab!=2) SetVariable setvarlow win=EOSpanel, disable=(s.tab!=2) SetVariable setvarhigh win=EOSpanel, disable=(s.tab!=2) SetVariable setvartol win=EOSpanel, disable=(s.tab!=2) TitleBox title0, win=EOSpanel, disable=(s.tab!=3) return 0 end threadsafe function PitzerSternerFugacity(variable GPa, variable degC) variable Pa = GPa*1e9 variable T = degC+273.15 variable V = PitzerSternerVolume(GPa, degC) variable den = 1/V // mol/cc variable lnf, phi wave PS=PitzerSternerCoefficients() Make /D/free /n=10 c c=PS[p][0]*T^-4+PS[p][1]*T^-2+PS[p][2]*T^-1+PS[p][3]+PS[p][4]*T+PS[p][5]*T^2 // phi is dimensionless residual Helmholtz energy, A/(RT) phi = c[0]*den+(1/(c[1]+c[2]*den+c[3]*den^2+c[4]*den^3+c[5]*den^4)-1/c[1]) phi -= c[6]/c[7]*(exp(-c[7]*den)-1) phi -= c[8]/c[9]*(exp(-c[9]*den)-1) variable R = 8.314510 V /= 1e6 // m^3 lnf = (-ln(V) + phi + Pa*V/(R*T)) + ln(R*T) - 1 return exp(lnf)/1e9 // fugacity in GPa end threadsafe function SakaneFugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable Pa = GPa*1e9 variable T = degC+273.15 variable V = SakaneVolume(GPa, degC) variable den = 1/V wave Sak=SakaneCoefficients() Make /D/free/n=10 c c = Sak[p][0]*T^-4+Sak[p][1]*T^-2+Sak[p][2]*T^-1+Sak[p][3]+Sak[p][4]*T+Sak[p][5]*T^2 variable lnf, phi phi = c[0]*den+(1/(c[1]+c[2]*den+c[3]*den^2+c[4]*den^3+c[5]*den^4)-1/c[1]) phi -= c[6]/c[7]*(exp(-c[7]*den)-1) phi -= c[8]/c[9]*(exp(-c[9]*den)-1) variable R = 8.314510 V /= 1e6 // m^3 lnf = (-ln(V) + phi + Pa*V/(R*T)) + ln(R*T) - 1 return exp(lnf)/1e9 // fugacity in GPa end // there is an analytical form for this threadsafe function HollandPowell91FugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK = degC+273.15, Pbar = GPa*1e4, Rgas = 83.14510 // K, bar, cc bar/K/mol Make /free/n=1 w=TK variable Vint=Integrate1D(EOS_V_HollandPowell91_P,1, Pbar, 2, 0, w) return 1e-4*exp(Vint/Rgas/TK) // GPa end threadsafe function HollandPowell91Fugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK = degC+273.15, Pbar = GPa*1e4, Rgas = 83.14510 // K, bar, cc bar/K/mol variable V = HollandPowell98Volume(GPa, degC, MRK=1, excess=0) // cc/mol variable Z = Pbar*V/Rgas/TK variable a, b, c, d variable lnGammaExcess = 0, lnGammaMRK=0, P0=2 // calculate a and b for H + P 91 MRK make /free w_a={1113.4,-0.88517,4.53e-3,-1.3183e-5,-0.22291,-3.8022e-4,1.7791e-7,5.8487,-2.137e-2,6.8133e-5} if (TK<673) a = w_a[0] + w_a[1]*(673-TK) + w_a[2]*(673-TK)^2 + w_a[3]*(673-TK)^3 else a = w_a[0] + w_a[4]*(TK-673) + w_a[5]*(TK-673)^2 + w_a[6]*(TK-673)^3 endif // kJ^2/kbar K^(1/2) mol^-2 b=1.465 //kJ/kbar/mol a*=100000 b*=10 variable BB = b*Pbar/Rgas/TK variable AA = a/(b*Rgas*TK^1.5) lnGammaMRK = Z - 1 - ln(Z-BB) - AA*ln(1+BB/Z) // from appendix of H+P 91 if ((Pbar/1000)>P0) variable c0 = -3.025650e-2, c1 = -5.343144e-6 variable d0 = -3.2297554e-3, d1 = 2.2215221e-6 c = c0+c1*TK d = d0+d1*TK variable kbar = Pbar/1000 // V in kJ/kbar, P in kbar, so Rgas = 0.008314510 // kJ/K/mol lnGammaExcess = 1/Rgas/TK*( 2/3*c*(kbar-P0)^1.5 + d/2*(kbar-P0)^2 ) endif return GPa*exp(lnGammaMRK + lnGammaExcess) end threadsafe function HollandPowell98Volume(variable GPa, variable degC, [int MRK, int excess]) // pressure (GPa), temperature (°C) MRK = paramisdefault(MRK) ? 1 : MRK excess = paramisdefault(excess) ? 1 : excess variable Pbar = GPa*1e4, TK = degC + 273.15 // bar, K variable a = 1.9853e-3, b = -8.9090e-2, c = 8.0331e-2, P0 = 2.0 variable kbar = Pbar/1000 variable V_MRK=0, Vexcess=0 if (MRK) wave settings = root:Packages:EOS:EOSsettings Make /free/n=1 w=TK FindRoots /H=(settings[%Vhigh])/L=(settings[%Vlow])/Q/T=(settings[%Vtol])/Z=(Pbar) EOS_HollandPowell_MRK, w V_MRK = V_Flag ? NaN : V_Root endif if (excess && (Pbar/1000)>P0) Vexcess = a*(kbar-P0) + b*(kbar-P0)^0.5 + c*(kbar-P0)^0.25 Vexcess *= 10 // j/bar >> cc/mol endif return (V_MRK+Vexcess) // cc/mol end threadsafe function HollandPowell98Fugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK = degC + 273.15, Pbar = GPa*1e4, Rgas = 83.14510 // K, bar, cc bar/K/mol variable V = HollandPowell98Volume(GPa, degC, MRK=1, excess=0) // cc/mol variable Z = Pbar*V/Rgas/TK variable a, b, c, d variable lnGammaExcess = 0, lnGammaMRK=0, P0=2 // calculate a and b for H + P 91 MRK make /free w_a={1113.4,-0.88517,4.53e-3,-1.3183e-5,-0.22291,-3.8022e-4,1.7791e-7,5.8487,-2.137e-2,6.8133e-5} if (TK<673) a = w_a[0] + w_a[1]*(673-TK) + w_a[2]*(673-TK)^2 + w_a[3]*(673-TK)^3 else a = w_a[0] + w_a[4]*(TK-673) + w_a[5]*(TK-673)^2 + w_a[6]*(TK-673)^3 endif // kJ^2/kbar K^(1/2) mol^-2 b = 1.465 //kJ/kbar/mol a *= 100000 b *= 10 variable BB = b*Pbar/Rgas/TK variable AA = a/(b*Rgas*TK^1.5) lnGammaMRK = Z - 1 - ln(Z-BB) - AA*ln(1+BB/Z) // from appendix of H+P 91 if ((Pbar/1000)>P0) // switch to H+P 98 a, b, c a = 1.9853e-3 b = -8.9090e-2 c = 8.0331e-2 variable kbar=Pbar/1000 // Vexcess = a*(kbar-P0) + b*(kbar-P0)^0.5 + c*(kbar-P0)^0.25 // V in kJ/kbar, P in kbar, so Rgas = 0.008314510 // kJ/K/mol lnGammaExcess = 1/Rgas/TK*( a/2*(kbar-P0)^2 + 2/3*b*(kbar-P0)^1.5 + 4/5*c*(kbar-P0)^1.25 ) endif return GPa*exp(lnGammaMRK + lnGammaExcess) end // ------------------------ wrapper functions -------------------------------- // ------- wrapper functions for equations that calculate P(V,T) ------------ threadsafe function EOS_BrodholtWood_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_BrodholtWood(TK,w[0]) // bars end threadsafe function EOS_BrodholtWood_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_BrodholtWood(w[0],V) // bars end threadsafe function EOS_Holloway_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_Holloway(TK,w[0]) // bars end threadsafe function EOS_Holloway_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_Holloway(w[0],V) // bars end threadsafe function EOS_Holloway94_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_Holloway94(TK,w[0]) // bars end threadsafe function EOS_Holloway94_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_Holloway94(w[0],V) // bars end threadsafe function EOS_PitzerSterner_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_PitzerSterner(TK,w[0]) // bars end threadsafe function EOS_PitzerSterner_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_PitzerSterner(w[0],V) // bars end threadsafe function EOS_IAPWS95_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_IAPWS95(TK,w[0]) // bars end threadsafe function EOS_IAPWS95_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_IAPWS95(w[0],V) // bars end threadsafe function EOS_BelonoshkoSaxena_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_BelonoshkoSaxena(TK,w[0]) // bars end threadsafe function EOS_BelonoshkoSaxena_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_BelonoshkoSaxena(w[0],V) // bars end threadsafe function EOS_SaulWagner_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_SaulWagner(TK,w[0]) // bars end threadsafe function EOS_SaulWagner_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_SaulWagner(w[0],V) // bars end threadsafe function EOS_Sakane_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_Sakane(TK,w[0]) // bars end threadsafe function EOS_Sakane_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_Sakane(w[0],V) // bars end threadsafe function EOS_ZhangDuan05_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_ZhangDuan05(TK,w[0]) // bars end threadsafe function EOS_ZhangDuan05_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_ZhangDuan05(w[0],V) // bars end threadsafe function EOS_ZhangDuan09_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_ZhangDuan09(TK,w[0]) // bars end threadsafe function EOS_ZhangDuan09_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_ZhangDuan09(w[0],V) // bars end threadsafe function EOS_vanDerWaals_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_vanDerWaals(TK,w[0]) // bars end threadsafe function EOS_vanDerWaals_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_vanDerWaals(w[0],V) // bars end threadsafe function EOS_RedlichKwong_T(wave w, variable TK) // w[0]=V(cc/mol), TK=T(K) return EOS_RedlichKwong(TK,w[0]) // bars end threadsafe function EOS_RedlichKwong_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_RedlichKwong(w[0],V) // bars end // ------- wrapper functions for equations that calculate V(P,T) ------------ threadsafe function EOS_V_HollandPowell91_T(wave w, variable T) // w[0]=P (bar), T in K return EOS_V_HollandPowell91(w[0],T) // cc/mol end threadsafe function EOS_V_HollandPowell91_P(wave w, variable Pbar) // w[0]=T(K), Pbar = P (bar) return EOS_V_HollandPowell91(Pbar,w[0]) // cc/mol end threadsafe function EOS_V_HollandPowell98_T(wave w, variable T) // w[0]=P (bar), T in K return EOS_V_HollandPowell98(w[0],T) // cc/mol end threadsafe function EOS_V_HollandPowell98_P(wave w, variable Pbar) // w[0]=T(K), Pbar = P (bar) return EOS_V_HollandPowell98(Pbar,w[0]) // cc/mol end threadsafe function EOS_V_SaxenaFei_T(wave w, variable T) // w[0]=P (bar), T in K return EOS_V_SaxenaFei(w[0],T) // cc/mol end threadsafe function EOS_V_SaxenaFei_P(wave w, variable Pbar) // w[0]=T(K), Pbar = P (bar) return EOS_V_SaxenaFei(Pbar,w[0]) // cc/mol end threadsafe function EOS_V_FrostWood_T(wave w, variable T) // w[0]=P (bar), T in K return EOS_V_FrostWood(w[0], T) // cc/mol end threadsafe function EOS_V_FrostWood_P(wave w, variable Pbar) // w[0]=T(K), Pbar = P (bar) return EOS_V_FrostWood(Pbar,w[0]) // cc/mol end // -------------------------- equations of state for pure water -------------------------------------- // MRK type EOS for TIP4P water. For 10 kbars and over. // Calculates P in bars as a function of T(K) and V(cc/mol) threadsafe function EOS_BrodholtWood(variable T, variable V) variable a,b,Prep,Patr,R R = 83.144 // cc-bar/degree/mole a = -582468 -3038.79*T -9.24574e-3*T^2 + 3.02674e9/T^2 b = -3.90463E-02 - 0.991078*V Prep = R*T/(V-b) Patr = a/(sqrt(T)*V*(V+b)) return Prep - Patr + 36490.5/V -1.02451e7/V^2 -1.79681e8/V^3 + 2.18437e9/V^4 end // Holloway 1977, MRK type EOS // Calculates P in bars as a function of T(K) and V(cc/mol) threadsafe function EOS_Holloway(variable T, variable V) // temperature (K), volume (cc/mol) // ************* WHAT ABOUT THE FLOWERS CORRECTION??? ************ // i think the correction was for mixed fluids variable a, b, Prep, Patr, TempC, R TempC = T - 273.15 R = 83.1451 // cc-bar/degree/mole b = 14.6 a = 35e6 + 166.8e6 - 193080*TempC + 186.4*TempC^2 - 0.071288*TempC^3 Patr = a/(sqrt(T)*V*(V+b)) Prep = R*T/(V-b) return Prep - Patr end // Holloway 1994, MRK type EOS // Calculates P in bars as a function of T(K) and V(cc/mol) threadsafe function EOS_Holloway94(variable T, variable V) // temperature (K), volume (cc/mol) variable a, b, Prep, Patr, TempC, R R = 83.1451 // cc-bar/degree/mole b = 14.5 a = 115.98 - 0.0016295*T - 1.4984e-5*T^2 a *= 1e6 Patr = a/(sqrt(T)*V*(V+b)) Prep = R*T/(V-b) return Prep - Patr end // Pitzer, K.S. and Sterner, S.M., 1994. Equations of state valid // continuously from zero to extreme pressures for H2O and CO2. // Journal of Chemical Physics. 101: 3111-3116. threadsafe function EOS_PitzerSterner(variable T, variable V) // temperature (K), volume (cc/mol) wave PS=PitzerSternerCoefficients() Make /D/free/N=10 c c = PS[p][0]*T^-4+PS[p][1]*T^-2+PS[p][2]*T^-1+PS[p][3]+PS[p][4]*T+PS[p][5]*T^2 variable R, den, pressure, var_num, var_denom R = 8314510 // Pa.cc/K/mol den = 1/V // den in mol/cc var_num = c[2]+2*c[3]*den+3*c[4]*den^2+4*c[5]*den^3 var_denom = (c[1]+c[2]*den+c[3]*den^2+c[4]*den^3+c[5]*den^4)^2 pressure = den+c[0]*den^2-den^2*(var_num/var_denom) pressure += c[6]*den^2*exp(-c[7]*den)+c[8]*den^2*exp(-c[9]*den) pressure *= R * T return pressure/1e5 // return pressure in bar end threadsafe static function /WAVE PitzerSternerCoefficients() Make /free/D /N=(10,6) PS PS[0][0]= {0,0,0,0,0,0,3887865600000,0,-14182435000000,0} PS[0][1]= {0,0,0,0,0,0,-134948780,0,181653900,0} PS[0][2]= {246576.88,0.58638965,-6.278384,0,5665.4978,0,309165.64,-65537.898,-197690.68,92093.375} PS[0][3]= {51.359951,-0.0028646939,0.014791599,-0.42719875,-16.580167,0.10917883,7.5591105,188.10675,-23.530318,122.46777} PS[0][4]= {0,3.1375577e-05,0.00035779579,-1.6325155e-05,0.076560762,0,0,0,0,0} PS[0][5]= {0,0,1.5432925e-08,0,0,0,0,0,0,0} return PS end threadsafe function EOS_BelonoshkoSaxena(variable TK, variable V) // K, cc/mol variable a = 1.40203e5 - 41.2336*TK variable b = -7.72779e6 + 6.70124e3*TK variable c = 6.52012e9 - 0.45580e6*TK return (a*V^-1 + b*V^-2 + c*V^-4.586) // bars end // Saul, A. and Wagner, W. (1989). A fundamental equation for water covering the range from the // melting line to 1273 K at pressures up to 25000 MPa. J. Phys. Chem. Ref. Data 18:1537-1564 threadsafe function EOS_SaulWagner(variable Temperature, variable V) // K, cc/mol wave c=SaulWagnerCoefficients() variable density=18.01528/V*1e3 // kg/m^3 variable Tc=647.14, criticalDensity=322, R=0.46151805 // K,g/cc, J/gK variable tau=Tc/Temperature, delta=density/criticalDensity variable fi_RealDerivative=0, sum1=0, sum2=0, pressure variable Vgamma, d, t, a, i for(i=0;i<9;i+=1) a=C[i][%a]; d=C[i][%d]; t=C[i][%t] fi_RealDerivative+=a*d*delta^(d-1)*tau^t endfor for(i=9;i<54;i+=1) Vgamma=C[i][%g]; a=C[i][%a]; d=C[i][%d]; t=C[i][%t] fi_RealDerivative+=exp(-delta^Vgamma)*(a*delta^(d-1)*(d-Vgamma*delta^Vgamma)*tau^t) endfor for(i=55;i<58;i+=1) Vgamma=C[i][%g]; a=C[i][%a]; d=C[i][%d]; t=C[i][%t] sum1+=a*delta^(d+5)*tau^t sum2+=a*d*delta^(d-1)*tau^t endfor fi_RealDerivative+=(-2.4*exp(-0.4*delta^6)+12*exp(-2*delta^6))*sum1 fi_RealDerivative+=(exp(-0.4*delta^6)-exp(-2*delta^6))*sum2 pressure=density*R*Temperature*(1+delta*fi_RealDerivative) // kg/m^3, Pa m^3/K/g, K, so pressure in KPa return pressure/100 // bar end threadsafe static function /WAVE SaulWagnerCoefficients() Make /D/free/n=(58,4) SW SetDimLabel 1, 0, g, SW SetDimLabel 1, 1, d, SW SetDimLabel 1, 2, t, SW SetDimLabel 1, 3, a, SW SW[0][0]= {0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,0,0,0,0} SW[0][1]= {1,1,2,5,8,11,11,13,13,1,1,1,2,2,3,4,4,4,5,6,7,8,9,11,1,1,1,1,2,2,4,5,6,6,7,7,8,10,10,11,11,11,11,11,2,2,2,3,3,4,4,5,5,5,1,2,3,4} SW[0][2]= {0,2,0,9,0,0,12,7,13,0,1,3,1,5,5,2,3,5,6,4,1,8,0,1,0,9,10,11,0,8,5,4,2,12,3,10,3,2,8,0,1,3,4,6,13,14,15,14,16,13,26,15,23,25,50,40,32,26} SW[0][3]= {0.8216377478,-0.254389437,-0.0883086864,-8.90309724e-07,-1.241333357e-06,2.895590286e-09,1.403610309e-11,8.183943371e-13,-2.397905287e-13,-0.7519743341,-0.4151278588,-1.03051374,-1.648036888} SW[13][3]= {-0.4686350251,0.3560258142,-0.6364658294,0.2227482363,-0.08954849939,0.001557686788,0.001347719088,-0.001301353385,9.987368673e-07,0.0002263629476,2.89330495e-06,0.1995437169,-0.02707767662} SW[26][3]= {0.01849068216,-0.004402394357,-0.08546876737,0.1220538576,-0.2562237041,0.2555034636,-0.06323203907,3.351397575e-05,-0.06152834985,-0.0003533048208,0.03146309259,-0.002261795983} SW[38][3]= {0.00018689702,-0.001384614556,0.002713160073,-0.004866118539,0.003751789129,-0.0005692669373,-0.5876414555,0.56878383464,0.1642158198,0.5878635885,-0.2844301931,-0.2049198337,-0.004039233716} SW[51][3]= {0.05459049594,-0.008914260146,0.004974411254,-0.00709318338,0.01718796342,-0.01482653038,0.004517292884} return SW end threadsafe function EOS_IAPWS95(variable temperature, variable V) // temperature (K), volume (cc/mol) variable Tc=647.096 // critical temperature, K variable den_c=322 // critical density, kg/m^3 variable R=0.461518056e3 // J/kg/K - remove one significant figure to match output from C++ code variable mass=18.01528e-3 // kg/mol variable density=mass/V*1e6 // density, kg/m^3 variable delta=density/den_c variable tau=Tc/temperature wave w=IAPWS95Coefficients() // contains first 4 columns of the coefficient table variable phi=0 // derivative of residual Helmholtz free energy int i variable n, d, t, c for(i=1;i<8;i++) n=w[i][%n] d=w[i][%d] t=w[i][%t] phi+=n*d*delta^(d-1)*tau^t endfor for(i=8;i<52;i++) n=w[i][%n] d=w[i][%d] t=w[i][%t] c=w[i][%c] phi+=n*exp(-delta^c)*( delta^(d-1)*tau^t*(d-c*delta^c) ) endfor variable alpha_i, beta_i, gamma_i, epsilon_i for(i=52;i<55;i++) n=w[i][%n] d=w[i][%d] t=w[i][%t] alpha_i = 20 beta_i = 150 + 100*(i==54) gamma_i = 1.21 + 0.04*(i==54) epsilon_i = 1 phi+=n*delta^d*tau^t*exp( -alpha_i*(delta-epsilon_i)^2 - beta_i*(tau-gamma_i)^2 ) * (d/delta-2*alpha_i*(delta-epsilon_i) ) endfor // the final summation using coefficients in last two lines of table is written out in full // capitalised A, B, C, D represented as AA, BB, CC, DD. variable a, b, psi, bigDelta, theta, AA, BB, CC, DD, dpsi_ddelta, dbigDeltaeb_ddelta, dbigDelta_ddelta i = 55 a = 3.5; b = 0.85; BB = 0.2; CC = 28; DD = 700; AA = 0.32; beta_i = 0.3 n = w[i][%n] psi = exp(-CC*(delta-1)^2-DD*(tau-1)^2) theta = (1-tau)+AA*((delta-1)^2)^(1/(2*beta_i)) bigDelta = theta^2+BB*((delta-1)^2)^a dbigDelta_ddelta = (delta-1)*( AA*theta*2/beta_i*((delta-1)^2)^(1/(2*beta_i)-1) + 2*BB*a*((delta-1)^2)^(a-1) ) dbigDeltaeb_ddelta = b*bigDelta^(b-1)*dbigDelta_ddelta dpsi_ddelta = -2*CC*(delta-1)*psi phi += n*( bigDelta^b*(psi+delta*dpsi_ddelta) + dbigDeltaeb_ddelta*delta*psi ) i = 56 a = 3.5; b = 0.95; BB = 0.2; CC = 32; DD = 800; AA = 0.32; beta_i = 0.3 n = w[i][%n] psi = exp(-CC*(delta-1)^2-DD*(tau-1)^2) theta = (1-tau)+AA*((delta-1)^2)^(1/(2*beta_i)) bigDelta = theta^2+BB*((delta-1)^2)^a dbigDelta_ddelta = (delta-1)*( AA*theta*2/beta_i*((delta-1)^2)^(1/(2*beta_i)-1) + 2*BB*a*((delta-1)^2)^(a-1) ) dbigDeltaeb_ddelta = b*bigDelta^(b-1)*dbigDelta_ddelta dpsi_ddelta = -2*CC*(delta-1)*psi phi += n*( bigDelta^b*(psi+delta*dpsi_ddelta) + dbigDeltaeb_ddelta*delta*psi ) variable pressure = density*R*temperature*(1+delta*phi) return pressure/1e5 // Pa -> bar end threadsafe static function /WAVE IAPWS95Coefficients() Make /D/free/n=(57,4) w=nan SetDimLabel 1, 0, c, w SetDimLabel 1, 1, d, w SetDimLabel 1, 2, t, w SetDimLabel 1, 3, n, w w[8,22][0] = 1; w[23,42][0] = 2; w[43,46][0] = 3; w[47][0] = 4; w[48,51][0] = 6; w[55,56][0] = 3.5 w[1][1] = {1,1,1,2,2,3,4,1,1,1,2,2,3,4,4,5,7,9,10,11,13,15,1,2} w[25][1] = {2,2,3,4,4,4,5,6,6,7,9,9,9,9,9,10,10,12,3,4,4,5,14,3,6,6,6,3,3,3,0.85,0.95} w[1][2] = {-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4,2,13,9,3,4,11,4,13,1,7,1} w[25][2] = {9,10,10,3,7,10,10,6,10,10,1,2,3,4,8,6,9,8,16,22,23,23,10,50,44,46,50,0,1,4,0.2,0.2} w[1][3] = {0.12533547935523e-1,0.78957634722828e1,-0.87803203303561e1,0.31802509345418,-0.26145533859358} w[6][3] = {-0.78199751687981e-2,0.88089493102134e-2,-0.66856572307965,0.20433810950965,-0.66212605039687e-4} w[11][3] = {-0.19232721156002,-0.25709043003438,0.16074868486251,-0.40092828925807e-1,0.39343422603254e-6} w[16][3] = {-0.75941377088144e-5,0.56250979351888e-3,-0.15608652257135e-4,0.11537996422951e-8,0.36582165144204e-6} w[21][3] = {-0.13251180074668e-11,-0.62639586912454e-9,-0.10793600908932,0.17611491008752e-1,0.22132295167546} w[26][3] = {-0.40247669763528,0.58083399985759,0.49969146990806e-2,-0.31358700712549e-1,-0.74315929710341} w[31][3] = {0.47807329915480,0.20527940895948e-1,-0.13636435110343,0.14180634400617e-1,0.83326504880713e-2} w[36][3] = {-0.29052336009585e-1,0.38615085574206e-1,-0.20393486513704e-1,-0.16554050063734e-2,0.19955571979541e-2} w[41][3] = {0.15870308324157e-3,-0.16388568342530e-4,0.43613615723811e-1,0.34994005463765e-1,-0.76788197844621e-1} w[46][3] = {0.22446277332006e-1,-0.62689710414685e-4,-0.55711118565645e-9,-0.19905718354408,0.31777497330738} w[51][3] = {-0.11841182425981,-0.31306260323435e2,0.31546140237781e2,-0.25213154341695e4,-0.14874640856724,0.31806110878444} return w end threadsafe function EOS_Sakane(variable TK, variable V) // temperature (K), volume (cc/mol) wave Sak=SakaneCoefficients() Make /D/free/n=10 c c = Sak[p][0]*TK^-4+Sak[p][1]*TK^-2+Sak[p][2]*TK^-1+Sak[p][3]+Sak[p][4]*TK+Sak[p][5]*TK^2 variable var_num, var_denom, R, den, pressure R = 8314510 // Pa.cc/K/mol den = 1/V // mol/cc var_num = c[2]+2*c[3]*den+3*c[4]*den^2+4*c[5]*den^3 var_denom = (c[1]+c[2]*den+c[3]*den^2+c[4]*den^3+c[5]*den^4)^2 pressure = den+c[0]*den^2-den^2*(var_num/var_denom) pressure += c[6]*den^2*exp(-c[7]*den)+c[8]*den^2*exp(-c[9]*den) pressure *= R * TK return pressure/1e5 // return pressure in bar end threadsafe static function /WAVE SakaneCoefficients() Make /D/free/n=(10,6) Sak Sak[0][0]= {0,0,0,0,0,0,-12461703600000,0,154412618000000,0} Sak[0][1]= {0,0,0,0,0,0,-19579063.3,0,-771910139,0} Sak[0][2]= {274483.916,0.379124962,16.2544402,0,1732.81357,0,319957.478,-53213.8383,634276.326,197309.516} Sak[0][3]= {47.4274288,-0.00309382093,-0.0292684003,-0.686628542,-20.5333503,139.070662,-60.891445,198.958628,-98.990493,7.03943621} Sak[0][4]= {0,3.13284131e-05,0.00039342123,0.00029480696,0.0893992033,0,0,0,0,0} Sak[0][5]= {0,0,9.43483294e-09,0,0,0,0,0,0,0} return Sak end // for supercritical conditions only THIS BIT IS MRK ONLY - NOT // COMPENSATED Holland and Powell. A Compensated-Redlich-Kwong equation // for volumes and fugacities of CO2 and H2O in the range 1 bar to 50 // kbar and 100-1600 °C. Contrib. Mineral. petrol (1991) 109:265-273 // This function is used in EOS_V_HollandPowell91, EOS_V_HollandPowell98 threadsafe function EOS_HollandPowell_MRK(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) variable T=w[0] // K variable a, Patr, Prep make /free w_a={1113.4,-0.88517,4.53e-3,-1.3183e-5,-0.22291,-3.8022e-4,1.7791e-7,5.8487,-2.137e-2,6.8133e-5} variable b=1.465 if (T<673) a = w_a[0] + w_a[1]*(673-T) + w_a[2]*(673-T)^2 + w_a[3]*(673-T)^3 else a = w_a[0] + w_a[4]*(T-673) + w_a[5]*(T-673)^2 + w_a[6]*(T-673)^3 endif // a in kJ/kbar/mol // a*100000 >> bar cm^6 K^0.5 mol^-2 a *= 100000 b *= 10 variable R=83.1451 // cc-bar/K/mole Patr = a/(sqrt(T)*V*(V+b)) Prep = R*T/(V-b) return Prep - Patr // bars end // Zhang and Duan (2005) Prediction of the PVT properties of // water over wide range of temperatures and pressures from molecular dynamics simulation // Physics of the Earth and Planetary Interiors 149 (2005) 335–354. threadsafe function EOS_ZhangDuan05(variable TK, variable V) // temperature (K), volume (cc/mol) variable V_gamma Make /D/free/n=15 a a[0] = {NaN,3.49824207e-1,-2.91046273e0, 2.00914688e0,1.12819964e-1,7.48997714e-1,-8.73207040e-1,1.70609505e-2} a[8] = {-1.46355822e-2,5.79768283e-2,-8.41246372e-4,4.95186474e-3,-9.16248538e-3,-1.00358152e-1,-1.82674744e-3} V_gamma = 1.05999998e-2 variable R = 83.14467 // R in bar.cc/K/mol variable Tc = 647.25 //K variable Vc = 55.9480373 // cc/mol variable Tr = TK/Tc variable Vr = V/Vc variable B, C, D, E, F, G B = a[1] + a[2]/Tr^2 + a[3]/Tr^3 C = a[4] + a[5]/Tr^2 + a[6]/Tr^3 D = a[7] + a[8]/Tr^2 + a[9]/Tr^3 E = a[10] + a[11]/Tr^2 + a[12]/Tr^3 F = a[13] / Tr G = a[14] * Tr variable Z = 1 + B/Vr + C/Vr^2 + D/Vr^4 + E/Vr^5 + (F/Vr^2+G/Vr^4)*exp(-V_gamma/Vr^2) return Z * R * TK / V // return pressure in bar end threadsafe function EOS_ZhangDuan09(variable TK, variable V, [string phase]) // temperature (K), volume (cc/mol) // V /= 1e6 // m^3 //// phase = selectstring(paramisdefault(phase), phase, "H2O") // make /n=(8,2)/free wLenardJones // setdimlabel(0, 0, "null", wLenardJones) // setdimlabel(0, 1, "CH4", wLenardJones) // setdimlabel(0, 2, "H2O", wLenardJones) // setdimlabel(0, 3, "CO2", wLenardJones) // setdimlabel(0, 4, "H2", wLenardJones) // setdimlabel(0, 5, "CO", wLenardJones) // setdimlabel(0, 6, "O2", wLenardJones) // setdimlabel(0, 7, "C2H6", wLenardJones) // wLenardJones[][0] = {nan,154,510,235,31.2,105.6,124.5,246.1} // wLenardJones[][1] = {nan,3.691,2.88,3.79,2.93,3.66,3.36,4.35} // variable epsilon = 510, sigma = 2.88 // // if (paramisdefault(phase) == 0) // // // // strswitch (phase) // case "CH4": // epsilon = 154 // sigma = 3.691 // break // case "H2O": // epsilon = 510 // sigma = 2.88 // break // case "CO2": // epsilon = 235 // sigma = 3.79 // break // case "H2": // epsilon = 31.2 // sigma = 2.93 // break // case "CO": // epsilon = 105.6 // sigma = 3.66 // break // case "O2": // epsilon = 124.5 // sigma = 3.36 // break // case "C2H4": // epsilon = 246.1 // sigma = 4.35 // break // endswitch // endif phase = selectstring(paramisdefault(phase), phase, "H2O") make /n=(7,2)/free wLenardJones wLenardJones[][0] = {154,510,235,31.2,105.6,124.5,246.1} wLenardJones[][1] = {3.691,2.88,3.79,2.93,3.66,3.36,4.35} variable row = whichlistitem(phase, "CH4;H2O;CO2;H2;CO;O2;C2H6;", ";", 0, 0) variable epsilon = wLenardJones[row][0] variable sigma = wLenardJones[row][1] Make /D/free/n=16 a a[0]={NaN,2.95177298930e-2,-6.33756452413e3,-2.75265428882e5,1.29128089283e-3,-1.45797416153e2,7.65938947237e4,2.58661493537e-6} a[8]={0.52126532146e0,-1.39839523753e2,-2.36335007175e-8,5.35026383543e-3,-0.27110649951e0,2.50387836486e4,0.73226726041e0,1.54833359970e-2} // variable kB = 1.38064852e-23 // variable epsilon = 510*kB, sigma = 2.88e-10 // K, m variable Tm = 154*TK/epsilon variable Vm = V/1000*(3.691/sigma)^3 // does the factor of 1000 mean cc -> L? variable Z = 1 + (a[1] + a[2]/Tm^2 + a[3]/Tm^3)/Vm + (a[4] + a[5]/Tm^2 + a[6]/Tm^3)/Vm^2 + (a[7] + a[8]/Tm^2 + a[9]/Tm^3)/Vm^4 + (a[10] + a[11]/Tm^2 + a[12]/Tm^3)/Vm^5 + a[13]/Tm^3/Vm^2*(a[14]+a[15]/Vm^2)*exp(-a[15]/Vm^2) variable Rgas = 0.0831446261815324 // R in bar.L/K/mol variable Pm = Rgas*Tm*Z/Vm // bar return Pm*epsilon/3.0636/sigma^3 // return pressure in bar end // Liu et al (2019) Prediction of H2O PVT relations at high temperatures // by VHL equation of state // AIP Advances 9 (2019) 035334. doi: 10.1063/1.5090030 // this paper looks pretty dodgy function EOS_VHL(TK,V) // K, cc/mol variable TK, V variable pressure Make /free/n=11 coefC, coefD, coefE variable a, b, f coefC[0] = NaN coefC[1] = -0.96665 coefC[2] = -4.51039 coefC[3] = -2.58733 coefC[4] = -0.85487 coefC[5] = 2.08033 coefC[6] = -0.516311 coefC[7] = 2.02825 coefC[8] = -2.15543 coefC[9] = -0.12489 coefC[10] = -7.65161 coefD[0] = NaN coefD[1] = -2.31154 coefD[2] = -7.49608 coefD[3] = -2.56864 coefD[4] = -1.14947 coefD[5] = 2.09543 coefD[6] = -0.84757 coefD[7] = 2.81486 coefD[8] = -4.4448 coefD[9] = -0.30669 coefD[10] = -11.36945 coefE = 0 coefE[0] = NaN coefE[1] = 2.15701 coefE[2] = -2.52736 coefE[3] = 0.2708 coefE[4] = -1.0019 a = -0.002431 b = 0.0001776 f = 3.608 variable bCH4, R, Cstar, Dstar, Tstar, N, epsilon, d // what is value for d? bCH4 = 67.21 R = 8.3144e-5 // R in N = 6.022045e23 variable epsilon0,j,b0,k,c variable epsilonT, bT epsilon0 = 229.61778 j = 1188.7 b0 = 45.0936 k = -70.86645 c = 352.4585 epsilonT=epsilon0*(1+j/TK) bT=b0+k*exp(-TK/c) variable omega, Estar omega=bT^2/(bCH4*V) Tstar=TK/epsilonT variable i for(i=1;1<=10;i+=2) Estar+=coefE[i]*Tstar^coefE[i+1] Dstar+=coefD[i]*Tstar^coefD[i+1] Cstar+=coefC[i]*Tstar^coefC[i+1] endfor // can't figure out how this is supposed to be formulated // variable Bstar=2*pi*N*d^3/3 // variable Bstar=b variable Bstar variable Z = 1 + omega*Bstar + omega^2*Cstar + omega^3*Dstar + omega^4*Estar + (a+b*Tstar)*omega^f pressure=Z*R*TK/V return pressure // return pressure in bar end threadsafe function EOS_vanDerWaals(variable T, variable V) // temperature (K), volume (cc/mol) variable a, b, Tc, Pc, R R = 83.144 // cc-bar/degree/mole // use water critical point Tc = 647.096 Pc = 220.6 // bars a = 27/64*(R*Tc)^2/Pc // bar cc^2 b = R*Tc/(8*Pc) // cc/mol return R*T/(V-b)-a/V^2 //bars end threadsafe function EOS_RedlichKwong(variable T, variable V) // temperature (K), volume (cc/mol) variable a, b, Tc, Pc, R R = 83.144 // cc-bar/degree/mole // use water critical point Tc = 647.096 Pc = 220.6 // bars a = R^2*Tc^2.5/(9*(2^(1/3)-1)*Pc) b = (2^(1/3)-1)*R*Tc/(3*Pc) return R*T/(V-b)-a/(V*(V+b)*sqrt(T)) //bars end // --------------------- explicit equations of state ------------------ // the following equations are explicit in P // The names of these functions are prefixed by EOS_V_ // to indicate that they return volume, rather than pressure threadsafe function EOS_V_HollandPowell91(variable Pbar, variable TK) // for supercritical conditions only (needs some more code to work at lower PT) // Holland and Powell. A Compensated-Redlich-Kwong equation for volumes and fugacities // of CO2 and H2O in the range 1 bar to 50 kbar and 100-1600 °C. // Contrib. Mineral. Petrol. (1991) 109:265-273 variable c, d, V_MRK, Vexcess variable c0 = -3.025650e-2, c1 = -5.343144e-6 variable d0 = -3.2297554e-3, d1 = 2.2215221e-6 variable Po=2 wave settings=root:Packages:EOS:EOSsettings Make /free/n=1 w=TK FindRoots /H=(settings[%Vhigh])/L=(settings[%Vlow])/Q/T=(settings[%Vtol])/Z=(Pbar) EOS_HollandPowell_MRK, w V_MRK = V_Flag ? NaN : V_Root c = c0 + c1*TK d = d0 + d1*TK if ((Pbar/1000)>Po) Vexcess = (c*sqrt(Pbar/1000-Po)+d*(Pbar/1000-Po)) Vexcess *= 10 // j/bar >> cc/mol endif return (V_MRK+Vexcess) end threadsafe function EOS_V_HollandPowell98(variable Pbar, variable TK) variable a = 1.9853e-3, b = -8.9090e-2, c = 8.0331e-2, P0 = 2.0 variable kbar=Pbar/1000 variable V_MRK, Vexcess wave settings=root:Packages:EOS:EOSsettings Make /free/n=1 w=TK FindRoots /H=(settings[%Vhigh])/L=(settings[%Vlow])/Q/T=(settings[%Vtol])/Z=(Pbar) EOS_HollandPowell_MRK, w V_MRK = V_Flag ? NaN : V_Root if ((Pbar/1000)>P0) Vexcess = a*(kbar-P0) + b*(kbar-P0)^0.5 + c*(kbar-P0)^0.25 Vexcess *= 10 // j/bar >> cc/mol endif return (V_MRK+Vexcess) end threadsafe function EOS_V_FrostWood(variable Pbar, variable TK) // bar, K variable R = 83.14 variable a, b, c, d, V_MRK, Vexcess variable a0 = 5.3957e8, a1 = -6.362e5, a2 = 236.81 variable b0 = 27.732, b1 = -2.0179e-2, b2 = 9.2125e-6 variable c0 = -0.1244, c1 = 1.79e-4, c2 = -7.8587e-8 variable d0 = 2.186e-4, d1 = -3.6836e-7, d2 = 1.6127e-10 a = a0+a1*TK+a2*TK^2 b = b0+b1*TK+b2*TK^2 V_MRK = R*TK/Pbar+b-a*R*sqrt(TK)/((R*TK+b*Pbar)*(R*TK+2*b*Pbar)) c = c0+c1*TK+c2*TK^2 d = d0+d1*TK+d2*TK^2 Vexcess = (c*sqrt(Pbar)+d*Pbar) return (V_MRK+Vexcess) end threadsafe function EOS_V_SaxenaFei(variable Pbar, variable TK) // 2 kbar lower limit; be careful when calculating fugacities. variable Pa = Pbar*1e5 variable Tc = 647.14,Pc = 22064000,R = 8.31451 variable Tr = TK/Tc,Pr = Pa/Pc,Z variable A,B,C,D A = 1.4937-1.8626*Tr^-2+0.80003*Tr^-3-0.3941*ln(Tr) B = 4.241e-2*Tr^-1+2.4097e-2*Tr^-2-8.9634e-3*Tr^-3 C = -9.016e-7*Tr^-1-6.1345e-5*Tr^-2+2.2380e-5*Tr^-3+5.2335e-7*ln(Tr) D = -7.6707e-9*Tr^-1+4.1108e-8*Tr^-2-1.4798e-8*Tr^-3-6.3033e-21*Tr^3 Z = A+B*Pr+C*Pr^2+D*Pr^3 return R*TK*Z*1e6/Pa end // ------------------------- these are useful functions for manual use ----------------------------- threadsafe function BelonoshkoSaxenaVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.001)/Z=(GPa*10000) EOS_BelonoshkoSaxena_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function BelonoshkoSaxenaFugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) // V in cc/mol, P in bar, so R in bar cc/K/mol variable TK = degC + 273.15, Rgas = 83.14510, Pbar = GPa*1e4 variable a = 1.40203e5 - 41.2336*TK variable b = -7.72779e6 + 6.70124e3*TK variable c = 6.52012e9 - 0.45580e6*TK variable V = BelonoshkoSaxenaVolume(GPa, degC), V0 = BelonoshkoSaxenaVolume(1e-4, degC) // cc/mol variable intPdV = a*(ln(V) - ln(V0)) - b*(V^-1 - V0^-1) - c/3.586*(V^-3.586 - V0^-3.586) // bar.cc/mol return 1e-4*exp(1/Rgas/TK*(Pbar*V - V0 - intPdV)) //(a/V+b/V^2+c/V^4.586) // bars end threadsafe function BrodholtWoodVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.001)/Z=(GPa*10000) EOS_BrodholtWood_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function HollowayVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(15)/Q/T=(0.001)/Z=(GPa*10000) EOS_Holloway_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function Holloway94Volume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(15)/Q/T=(0.0001)/Z=(GPa*10000) EOS_Holloway94_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function SaulWagnerVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.0001)/Z=(GPa*10000) EOS_SaulWagner_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function PitzerSternerVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.0001)/Z=(GPa*10000) EOS_PitzerSterner_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function IdealVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable R=83.1451 // bar cc/K return R*(degC+273.15)/(GPa*10000) end threadsafe function SakaneVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.0001)/Z=(GPa*10000) EOS_Sakane_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function ZhangDuan05Volume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.0001)/Z=(GPa*10000) EOS_ZhangDuan05_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function ZhangDuan05FugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable Vint,Rgas,TK,p0,Pbar Rgas = 83.14467 // bar cc/K TK = degC+273.15 // K p0 = 1e-4 // GPa Pbar = GPa*1e4 // bars Make /free/n=1 w=degC Vint = Integrate1D(ZhangDuan05_wrapper,1,Pbar,2,0,w) return P0*exp(Vint/Rgas/TK) // GPa end // use DMW92a fugacity eq threadsafe function ZhangDuan05Fugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable Rgas, TK, Pbar Rgas = 83.14467 // bar cc/K TK = degC+273.15 // K Pbar = GPa*1e4 // bars variable volume = ZhangDuan05Volume(GPa, degC) // cc/mol variable Z = pBar * volume / (Rgas * TK) variable V_gamma Make/D/free/n=15 a a[0] = {NaN,3.49824207e-1,-2.91046273e0, 2.00914688e0,1.12819964e-1,7.48997714e-1,-8.73207040e-1,1.70609505e-2} a[8] = {-1.46355822e-2,5.79768283e-2,-8.41246372e-4,4.95186474e-3,-9.16248538e-3,-1.00358152e-1,-1.82674744e-3} V_gamma = 1.05999998e-2 variable Tc = 647.25 //K variable Vc = 55.9480373 // cc/mol variable Tr = TK/Tc variable B, C, D, E, F, v_beta // F and beta in DMW 92a variable F05, G05 // F and G in Zhang and Duan 2005 B = a[1] + a[2]/Tr^2 + a[3]/Tr^3 C = a[4] + a[5]/Tr^2 + a[6]/Tr^3 D = a[7] + a[8]/Tr^2 + a[9]/Tr^3 E = a[10] + a[11]/Tr^2 + a[12]/Tr^3 F05 = a[13] / Tr G05 = a[14] * Tr F = G05 / v_gamma v_beta = F05 / F variable dA, dB, dC, dD, dE, dF, dBeta, dGamma dB = 2 * B * Vc dC = 3 * C * Vc^2 dD = 5 * D * Vc^4 dE = 6 * E * Vc^5 dF = 2 * F * Vc^2 dBeta = v_beta dGamma = 3 * v_gamma * Vc^2 variable FVc2 = F * Vc^2 variable gammaVc2 = V_gamma * Vc^2 variable exptermV = exp(-gammaVc2 / volume^2) variable lnPhi lnPhi = -ln(Z) lnPhi += dB / volume lnPhi += dC / (2 * volume^2) lnPhi += dD / (4 * volume^4) lnPhi += dE / (5 * volume^5) lnPhi += (dF*v_beta + dBeta*FVc2) / (2*gammaVc2) * (1 - exptermV) lnPhi += (dF*gammaVc2 + dGamma*FVc2 - FVc2*v_beta*(dGamma - gammaVc2)) / (2*gammaVc2^2) * (1 - (gammaVc2/volume^2 + 1) * exptermV) lnPhi += - (dGamma - gammaVc2)*FVc2 / (2*gammaVc2^2) * (2 - (gammaVc2^2/volume^4 + 2*gammaVc2/volume^2 + 2) * exptermV) return exp(lnPhi) * GPa end threadsafe function ZhangDuan05_wrapper(wave w, variable Pbar) return ZhangDuan05Volume(Pbar/1e4, w[0]) end threadsafe function ZhangDuan09Volume(variable GPa, variable degC, [string phase]) // pressure (GPa), temperature (°C) Make /free/n=2 w w[0] = degC + 273.15 w[1] = 1 if (paramisdefault(phase) == 0 && strlen(phase)) w[1] = whichlistitem(phase, "CH4;H2O;CO2;H2;CO;O2;C2H6", ";", 0, 0) endif FindRoots /H=(50)/L=(10)/Q/T=(0.0001)/Z=(GPa*10000) EOS_ZhangDuan09_Vwrapper, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function EOS_ZhangDuan09_Vwrapper(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) string phase = stringfromlist(w[1], "CH4;H2O;CO2;H2;CO;O2;C2H6") return EOS_ZhangDuan09(w[0], V, phase=phase) // bars end // this wrapper has the right number of variables for the H2O EOS calculator to recognise it threadsafe function ZhangDuan09Fugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) return ZhangDuan09ComponentFugacity(GPa, degC) end threadsafe function ZhangDuan09ComponentFugacity(variable GPa, variable degC, [string phase]) // pressure (GPa), temperature (°C) phase = selectstring(paramisdefault(phase), phase, "H2O") variable TK = degC + 273.15 // K variable V = ZhangDuan09Volume(GPa, degC, phase=phase) // cc/mol Make /D/free/n=16 a a[0]={NaN,2.95177298930e-2,-6.33756452413e3,-2.75265428882e5,1.29128089283e-3,-1.45797416153e2,7.65938947237e4,2.58661493537e-6} a[8]={0.52126532146e0,-1.39839523753e2,-2.36335007175e-8,5.35026383543e-3,-0.27110649951e0,2.50387836486e4,0.73226726041e0,1.54833359970e-2} // variable epsilon = 510, sigma = 2.88 // if (paramisdefault(phase) == 0) make /n=(7,2)/free wLenardJones wLenardJones[][0] = {154,510,235,31.2,105.6,124.5,246.1} wLenardJones[][1] = {3.691,2.88,3.79,2.93,3.66,3.36,4.35} variable row = whichlistitem(phase, "CH4;H2O;CO2;H2;CO;O2;C2H6;", ";", 0, 0) variable epsilon = wLenardJones[row][0] variable sigma = wLenardJones[row][1] // endif variable Tm = 154*TK/epsilon variable Vm = V/1000*(3.691/sigma)^3 // does the factor of 1000 mean cc -> L? variable Z = 1 + (a[1] + a[2]/Tm^2 + a[3]/Tm^3)/Vm + (a[4] + a[5]/Tm^2 + a[6]/Tm^3)/Vm^2 + (a[7] + a[8]/Tm^2 + a[9]/Tm^3)/Vm^4 + (a[10] + a[11]/Tm^2 + a[12]/Tm^3)/Vm^5 + a[13]/Tm^3/Vm^2*(a[14]+a[15]/Vm^2)*exp(-a[15]/Vm^2) variable S1a = (a[1] + a[2]/Tm^2 + a[3]/Tm^3)/Vm + (a[4] + a[5]/Tm^2 + a[6]/Tm^3)/2/Vm^2 + (a[7] + a[8]/Tm^2 + a[9]/Tm^3)/4/Vm^4 + (a[10] + a[11]/Tm^2 + a[12]/Tm^3)/5/Vm^5 variable S1b = a[13]/2/a[15]/Tm^3*( a[14]+1-(a[14]+1+a[15]/Vm^2)*exp(-a[15]/Vm^2) ) variable lnGamma = Z - 1 - ln(Z) + S1a + S1b return GPa*exp(lnGamma) end // function ZhangDuan05FugacityTest(variable GPa, variable degC) // pressure (GPa), temperature (°C) // // variable TK = degC + 273.15 // K // variable V = ZhangDuan05Volume(GPa, degC) // cc/mol // // variable V_gamma // Make /D/free/n=15 a // a[0]={NaN,3.49824207e-1,-2.91046273e0, 2.00914688e0,1.12819964e-1,7.48997714e-1,-8.73207040e-1,1.70609505e-2} // a[8]={-1.46355822e-2,5.79768283e-2,-8.41246372e-4,4.95186474e-3,-9.16248538e-3,-1.00358152e-1,-1.82674744e-3} // V_gamma=1.05999998e-2 // // variable R=83.14467 // R in bar.cc/K/mol // variable Tc=647.25 //K // variable Vc=55.9480373 // cc/mol // // variable Tr=TK/Tc // variable Vr=V/Vc // //// variable epsilon = 510, sigma = 2.88 //// Tr = 154*TK/epsilon //// Vr = V/1000*(3.691/sigma)^3 // does the factor of 1000 mean cc -> L? // // variable B, C, D, E, F, G // B=a[1] + a[2]/Tr^2 + a[3]/Tr^3 // C=a[4] + a[5]/Tr^2 + a[6]/Tr^3 // D=a[7] + a[8]/Tr^2 + a[9]/Tr^3 // E=a[10] + a[11]/Tr^2 + a[12]/Tr^3 // F=a[13]/Tr // G=a[14]*Tr // // variable Z = 1 + B/Vr + C/Vr^2 + D/Vr^4 + E/Vr^5 + (F/Vr^2+G/Vr^4)*exp(-V_gamma/Vr^2) // // // // variable S1a = B/Vr + C/2/Vr^2 + D/4/Vr^4 + E/5/Vr^5 // variable S1b = 1/2/V_gamma*( F + 1 -(F + 1 + G/Vr^2)*exp(-V_gamma/Vr^2) ) // // a[13]/2/a[15]/Tm^3*( a[14]+1-(a[14]+1+a[15]/Vm^2)*exp(-a[15]/Vm^2) ) // // variable lnGamma = Z - 1 - ln(Z) + S1a + S1b // return GPa*exp(lnGamma) // // // // // // // //// Make /D/free/n=16 a //// a[0]={NaN,2.95177298930e-2,-6.33756452413e3,-2.75265428882e5,1.29128089283e-3,-1.45797416153e2,7.65938947237e4,2.58661493537e-6} //// a[8]={0.52126532146e0,-1.39839523753e2,-2.36335007175e-8,5.35026383543e-3,-0.27110649951e0,2.50387836486e4,0.73226726041e0,1.54833359970e-2} //// variable epsilon = 510, sigma = 2.88 //// variable Tm = 154*TK/epsilon //// variable Vm = V/1000*(3.691/sigma)^3 // does the factor of 1000 mean cc -> L? //// variable Z = 1 + (a[1] + a[2]/Tm^2 + a[3]/Tm^3)/Vm + (a[4] + a[5]/Tm^2 + a[6]/Tm^3)/Vm^2 + (a[7] + a[8]/Tm^2 + a[9]/Tm^3)/Vm^4 + (a[10] + a[11]/Tm^2 + a[12]/Tm^3)/Vm^5 + a[13]/Tm^3/Vm^2*(a[14]+a[15]/Vm^2)*exp(-a[15]/Vm^2) //// variable S1 = (a[1] + a[2]/Tm^2 + a[3]/Tm^3)/Vm + (a[4] + a[5]/Tm^2 + a[6]/Tm^3)/2/Vm^2 + (a[7] + a[8]/Tm^2 + a[9]/Tm^3)/4/Vm^4 + (a[10] + a[11]/Tm^2 + a[12]/Tm^3)/5/Vm^5 + a[13]/2/a[15]/Tm^3*( a[14]+1-(a[14]+1+a[15]/Vm^2)*exp(-a[15]/Vm^2) ) //// variable lnGamma = Z - 1 - ln(Z) + S1 //// return GPa*exp(lnGamma) //end threadsafe function vanDerWaalsVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(98)/L=(85)/Q/T=(0.01)/Z=(GPa*10000) EOS_vanDerWaals_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function RedlichKwongVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(100)/L=(23)/Q/T=(0.0001)/Z=(GPa*10000) EOS_RedlichKwong_V, w return V_Flag ? NaN : V_Root // cc/mol end threadsafe function DelanyHelgesonVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) return EOS_V_DelanyHelgeson(GPa*1e4,degC+273.15) // cc/mol end function VolumeFromEOS(string eos, variable GPa, variable degC) // use this for the "*_V" functions FUNCREF EOStemplateFunction EOSfunc=$eos Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(15)/Q/T=(0.001)/Z=(GPa*10000) EOSfunc, w return V_Flag ? NaN : V_Root end threadsafe function SaulWagnerFugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable Vint,Rgas,TK,p0,Pbar Rgas=83.1451 // bar cc/K TK=degC+273.15 // K p0=1e-4 // GPa Pbar=GPa*1e4 // bars Make /free/n=1 w=degC Vint=Integrate1D(SaulWagner_wrapper,1,Pbar,2,0,w) return P0*exp(Vint/Rgas/TK) // GPa end threadsafe function SaulWagner_wrapper(wave w, variable Pbar) return SaulWagnerVolume(Pbar/1e4, w[0]) end // JM Delaney and HC Helgeson (1978). Calculation of the thermodynamic consequences // of dehydration in subducting oceanic crust to 100 KB and < 800 °C. Am. J. Sci. 278:638-686 // takes pressure in bars, temperature in K, returns volume in cc/mol. threadsafe function EOS_V_DelanyHelgeson(variable Pbar, variable TK) Make /D/free/n=(5,5) coef coef[0][0]= {-56130.073,-15.285559,-0.026092451,1.7140501e-05,-6.0126987e-09} coef[0][1]= {0.38101798,0.0001375239,3.5988857e-08,-1.6860893e-11,NaN} coef[0][2]= {-2.1167697e-06,-1.5586868e-09,-2.791658e-14,NaN,NaN} coef[0][3]= {2.0266445e-11,6.6329577e-15,NaN,NaN,NaN} coef[0][4]= {-8.3225572e-17,NaN,NaN,NaN,NaN} variable degC=TK-273.15,volume=0 variable i,j for(i=0;i<5;i+=1) for(j=0;j<5-i;j+=1) volume+=j*coef[i][j]*(degC^i)*(Pbar^(j-1)) endfor endfor return volume/0.0239 // cal/mol/bar>>cc/mol end // 'compensated' Redlich Kwong type EOS threadsafe function FrostWoodVolume(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK = degC+273.15 variable P = GPa*1e4 variable R = 83.14 variable a, b, c, d, V_MRK, Vexcess variable a0 = 5.3957e8, a1 = -6.362e5, a2 = 236.81 variable b0 = 27.732, b1 = -2.0179e-2, b2 = 9.2125e-6 variable c0 = -0.1244, c1 = 1.79e-4, c2 = -7.8587e-8 variable d0 = 2.186e-4, d1 = -3.6836e-7, d2 = 1.6127e-10 a = a0+a1*TK+a2*TK^2 b = b0+b1*TK+b2*TK^2 V_MRK = R*TK/P+b-a*R*sqrt(TK)/((R*TK+b*P)*(R*TK+2*b*P)) c = c0+c1*TK+c2*TK^2 d = d0+d1*TK+d2*TK^2 Vexcess = (c*sqrt(P)+d*P) return (V_MRK+Vexcess) end threadsafe function FrostWoodFugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK = degC+273.15, Pa = GPa*1e9, Pbar = GPa*1e4 variable R = 83.14 variable a, b, c, d, V_MRK, Vexcess variable a0 = 5.3957e8, a1 = -6.362e5, a2 = 236.81 variable b0 = 27.732, b1 = -2.0179e-2, b2 = 9.2125e-6 variable c0 = -0.1244, c1 = 1.79e-4, c2 = -7.8587e-8 variable d0 = 2.186e-4, d1 = -3.6836e-7, d2 = 1.6127e-10 a = a0+a1*TK+a2*TK^2 b = b0+b1*TK+b2*TK^2 c = c0+c1*TK+c2*TK^2 d = d0+d1*TK+d2*TK^2 // analytical integration of the volume explicit approximation variable RTlnf = R*TK*ln(Pbar) + b*Pbar + a/(b*sqrt(TK))*(ln(R*TK + b*Pbar) - ln (R*TK + 2*b*Pbar)) + 2*c/3*(Pbar-1)^(3/2) + d/2*(Pbar-1)^2 return exp(RTlnf/R/TK)/1e4 // fugacity in GPa end threadsafe function FrostWoodFugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK, Pbar, Rgas, Vint, P0 TK=degC+273.15 // K Pbar=GPa*1e4 // bar Rgas=83.14510 // cc bar/K P0=1e-4 // GPa Make /free/n=1 w=degC Vint=Integrate1D(FrostWood_wrapper,1, Pbar, 2, 0, w) return P0*exp(Vint/Rgas/TK) // GPa end threadsafe function FrostWood_wrapper(wave w, variable bars) return FrostWoodVolume(bars*1e-4,w[0]) end threadsafe function HollowayFugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK=degC+273.15, Pa=GPa*1e9, Pbar=GPa*1e4 variable R=83.144621 variable V = HollowayVolume(GPa, degC) variable Z = Pbar*V/R/TK variable b = 14.6 variable a = 35e6 + 166.8e6 - 193080*degC + 186.4*degC^2 - 0.071288*degC^3 Variable AA = a/(b*R*TK^1.5) Variable BB = b*Pbar/R/TK variable lnGamma = Z - 1 - ln(z-BB) - AA* ln(1+BB/Z) // from appendix of H+P 91 return GPa*exp(lnGamma) // fugacity in GPa end threadsafe function Holloway94Fugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK=degC+273.15, Pa=GPa*1e9, Pbar=GPa*1e4 variable R=83.144621 variable V = Holloway94Volume(GPa, degC) variable Z = Pbar*V/R/TK variable b = 14.5 variable a = 115.98 - 0.0016295*TK - 1.4984e-5*TK^2 a *= 1e6 Variable AA = a/(b*R*TK^1.5) Variable BB = b*Pbar/R/TK variable lnGamma = Z - 1 - ln(z-BB) - AA* ln(1+BB/Z) return GPa*exp(lnGamma) // fugacity in GPa end threadsafe function Holloway_wrapper(wave w, variable bars) return HollowayVolume(bars*1e-4,w[0]) end threadsafe function HollowayFugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK, Pbar, Rgas, Vint, P0 TK = degC+273.15 // K Pbar = GPa*1e4 // bar Rgas = 83.14510 // cc bar/K P0 = 1e-4 // GPa Make /free/n=1 w=degC Vint = Integrate1D(Holloway_wrapper,1, Pbar, 2, 0, w) return P0*exp(Vint/Rgas/TK) // GPa end threadsafe function Holloway94_wrapper(wave w, variable bars) return Holloway94Volume(bars*1e-4,w[0]) end threadsafe function Holloway94FugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK, Pbar, Rgas, Vint, P0 TK=degC+273.15 // K Pbar=GPa*1e4 // bar Rgas=83.14510 // cc bar/K P0=1e-4 // GPa Make /free/n=1 w=degC Vint=Integrate1D(Holloway94_wrapper,1, Pbar, 2, 0, w) return P0*exp(Vint/Rgas/TK) // GPa end threadsafe function IAPWS95FugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable Vint,Rgas,TK,p0,Pbar Rgas=83.1451 // bar cc/K TK=degC+273.15 // K p0=1e-4 // GPa Pbar=GPa*1e4 // bars Make /free/n=1 w=degC Vint=Integrate1D(EOS_IAPWS95_wrapper,1,Pbar,2,0,w) return P0*exp(Vint/Rgas/TK) // GPa end threadsafe function IAPWS95Fugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable Rgas, TK, pressure Rgas=83.1451 // bar cc/K TK=degC+273.15 // K pressure=GPa*1e9 variable volume=IAPWS95Volume(GPa, degC) // cc/mol variable Tc=647.096 // critical temperature, K variable den_c=322 // critical density, kg/m^3 // variable R=0.461518056e3 // J/kg/K - remove one significant figure to match output from C++ code variable mass=18.01528e-3 // kg/mol variable density=mass/volume*1e6 // density, kg/m^3 variable delta=density/den_c variable tau=Tc/TK wave w=IAPWS95Coefficients() // contains first 4 columns of the coefficient table variable phi=0 // residual dimensionless Helmholtz free energy int i variable n, d, t, c for(i=1;i<8;i++) n=w[i][%n] d=w[i][%d] t=w[i][%t] phi+=n*delta^d*tau^t endfor for(i=8;i<52;i++) n=w[i][%n] d=w[i][%d] t=w[i][%t] c=w[i][%c] phi+=n*delta^d*tau^t*exp(-delta^c) endfor variable alpha_i, beta_i, gamma_i, epsilon_i for(i=52;i<55;i++) n=w[i][%n] d=w[i][%d] t=w[i][%t] alpha_i = 20 beta_i = 150 + 100*(i==54) gamma_i = 1.21 + 0.04*(i==54) epsilon_i = 1 phi+=n*delta^d*tau^t*exp( -alpha_i*(delta-epsilon_i)^2 - beta_i*(tau-gamma_i)^2 ) endfor // the final summation using coefficients in last two lines of table is written out in full // capitalised A, B, C, D represented as AA, BB, CC, DD. variable a, b, psi, bigDelta, theta, AA, BB, CC, DD, dpsi_ddelta, dbigDeltaeb_ddelta, dbigDelta_ddelta i = 55 a=3.5; b=0.85; BB=0.2; CC=28; DD=700; AA=0.32; beta_i=0.3 n=w[i][%n] psi=exp(-CC*(delta-1)^2-DD*(tau-1)^2) theta=(1-tau)+AA*((delta-1)^2)^(1/(2*beta_i)) bigDelta=theta^2+BB*((delta-1)^2)^a dbigDelta_ddelta=(delta-1)*( AA*theta*2/beta_i*((delta-1)^2)^(1/(2*beta_i)-1) + 2*BB*a*((delta-1)^2)^(a-1) ) dbigDeltaeb_ddelta=b*bigDelta^(b-1)*dbigDelta_ddelta dpsi_ddelta=-2*CC*(delta-1)*psi phi += n*bigDelta^b*delta*psi i = 56 a=3.5; b=0.95; BB=0.2; CC=32; DD=800; AA=0.32; beta_i=0.3 n=w[i][%n] psi=exp(-CC*(delta-1)^2-DD*(tau-1)^2) theta=(1-tau)+AA*((delta-1)^2)^(1/(2*beta_i)) bigDelta=theta^2+BB*((delta-1)^2)^a dbigDelta_ddelta=(delta-1)*( AA*theta*2/beta_i*((delta-1)^2)^(1/(2*beta_i)-1) + 2*BB*a*((delta-1)^2)^(a-1) ) dbigDeltaeb_ddelta=b*bigDelta^(b-1)*dbigDelta_ddelta dpsi_ddelta=-2*CC*(delta-1)*psi phi += n*bigDelta^b*delta*psi // phi is dimensionless residual Helmholtz energy, A/RT variable R=8.314510 // Pa m^3 /K /mol volume *= 1e-6 // volume in m^3, pressure in Pa variable lnf = ln(1/volume) + phi + pressure*volume/(R*TK) + ln(R*TK) - 1 return exp(lnf)/1e9 // fugacity in GPa end // THIS FUNCTION NOT TESTED!!!!!!!!!! threadsafe function RedlichKwongFugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK=degC+273.15, Pa=GPa*1e9, Pbar=GPa*1e4 variable V = RedlichKwongVolume(GPa, degC) // cc/mol variable V0 = RedlichKwongVolume(1e-4, degC) variable a, b, R, Tc, Pc R = 83.144 // cc-bar/degree/mole // use water critical point Tc = 647.096 Pc = 220.6 // bars // Tc = paramisdefault(Tc) ? 647.096 : Tc // Pc = paramisdefault(Pc) ? 220.6 : Pc // bars a = R^2*Tc^2.5/(9*(2^(1/3)-1)*Pc) b = (2^(1/3)-1)*R*Tc/(3*Pc) variable intPdV = R*TK*(ln(V-b)-ln(V0-b)) + a/(b*sqrt(TK)) * (ln((b+V)/V) - ln((b+V0)/V0)) variable intVdP = Pbar*V - 1*V0 - intPdV return exp(intVdP/R/TK)/1e4 // fugacity in ? GPa end threadsafe function SaulWagnerFugacity(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK = degC + 273.15, Pa = GPa * 1e9 variable V = SaulWagnerVolume(GPa, degC)*1e-6 // m^3/mol wave c = SaulWagnerCoefficients() variable density=18.01528/V/1e3 // kg/m^3 variable Tc=647.14, criticalDensity=322 //, R=0.46151805 // K,g/cc, J/gK variable tau=Tc/TK, delta=density/criticalDensity variable phi=0, sum1=0 variable Vgamma, d, t, a, i for(i=0;i<9;i+=1) a=C[i][%a]; d=C[i][%d]; t=C[i][%t] phi+=a*delta^d*tau^t endfor for(i=9;i<54;i+=1) Vgamma=C[i][%g]; a=C[i][%a]; d=C[i][%d]; t=C[i][%t] phi+=exp(-delta^Vgamma)*(a*delta^d*tau^t) endfor for(i=55;i<58;i+=1) Vgamma=C[i][%g]; a=C[i][%a]; d=C[i][%d]; t=C[i][%t] sum1+=a*delta^d*tau^t endfor phi+=(exp(-0.4*delta^6)-exp(-2*delta^6))*sum1 // phi is dimensionless residual Helmholtz energy, A/(RT) variable R=8.314510 // Pa m^3 /K /mol variable lnf = ln(1/V) + phi + Pa*V/(R*TK) + ln(R*TK) - 1 return exp(lnf)/1e9 // fugacity in GPa end threadsafe function EOS_IAPWS95_wrapper(w,Pbar) wave w; variable Pbar return IAPWS95Volume(Pbar/1e4, w[0]) end threadsafe function IAPWS95Volume(variable GPa, variable degC) // pressure (GPa), temperature (°C) Make /free/n=1 w=degC+273.15 FindRoots /H=(50)/L=(5)/Q/T=(0.0001)/Z=(GPa*10000) EOS_IAPWS95_V, w return V_Flag ? NaN : V_Root // cc/mol end // ------------------------------------------------------- // eq 29-30 of IAPWS R7-97(2012) // valid up to critical point, T = 647.096K threadsafe function SaturationPressure(variable Ts) // K Make /free/D /n=11 n n[0]=NaN n[1]=0.11670521452767e4 n[2]=-0.72421316703206e6 n[3]=-0.17073846940092e2 n[4]=0.12020824702470e5 n[5]=-0.32325550322333e7 n[6]=0.14915108613530e2 n[7]=-0.48232657361591e4 n[8]=0.40511340542057e6 n[9]=-0.23855557567849 n[10]=0.65017534844798e3 variable p0=10 // p0 in MPa. Use 10 to return answer in bars. variable s=Ts + n[9]/(Ts-n[10]) variable A,B,C A=s^2+n[1]*s+n[2] B=n[3]*s^2+n[4]*s+n[5] C=n[6]*s^2+n[7]*s+n[8] return p0*(2*C/(-B+sqrt(B^2-4*A*C)))^4 // bars end // ---------------------------------------------------------------------- // Duan and Zhang (2006) // Equation of state of the H2O, CO2, and H2O–CO2 systems up to 10 GPa // and 2573.15 K: Molecular dynamics simulations with ab initio potential // surface threadsafe function EOS_DuanZhang06(variable TK, variable V, variable xCO2, [int highP, wave /Z pwave, int globalParameters]) // temperature (K), volume (cc/mol) if (paramIsDefault(globalParameters)==0 && globalParameters) wave pwave = root:packages:EOS:wparam elseif (paramisdefault(pwave)) wave pwave = DuanZhang06Parameters(TK, XCO2) endif int lowP = paramisdefault(highP) ? 1 : highP == 0 highP = paramisdefault(highP) ? 1 : highP variable Rgas = 83.14467 // cc bar/(K mol) variable Pbar, Z if (lowP) Z = 1 Z += pwave[%BVc][0][0]/V Z += pwave[%CVc2][0][0]/V^2 Z += pwave[%DVc4][0][0]/V^4 Z += pwave[%EVc5][0][0]/V^5 Z += pwave[%FVc2][0][0]/V^2*(pwave[%vBeta][0][0]+pwave[%gammaVc2][0][0]/V^2)*exp(-pwave[%gammaVc2][0][0]/V^2) Pbar = Rgas*TK*Z/V if(Pbar<0.2e4 && Pbar>0) return Pbar endif endif if (highP) Z = 1 Z += pwave[%BVc][0][1]/V Z += pwave[%CVc2][0][1]/V^2 Z += pwave[%DVc4][0][1]/V^4 Z += pwave[%EVc5][0][1]/V^5 Z += pwave[%FVc2][0][1]/V^2*(pwave[%vBeta][0][1]+pwave[%gammaVc2][0][1]/V^2)*exp(-pwave[%gammaVc2][0][1]/V^2) Pbar = Rgas*TK*Z/V endif return Pbar end threadsafe static function croot(variable v) return sign(v)*(abs(v)^(1/3)) end // after executing this function, calling other DuanZhang06 functions with GlobalParameters=1 is valid function /wave DuanZhang06CreateGlobalParameterWave(variable TK, variable XCO2) wave wParam = DuanZhang06Parameters(TK, XCO2) NewDataFolder /O root:packages NewDataFolder /O root:packages:EOS DFREF dfr = root:packages:EOS KillWaves /Z root:packages:EOS:wParam MoveWave wParam, root:packages:EOS:wParam return wParam end function DuanZhang06Volume(variable GPa, variable degC, variable XCO2, [int Prange, int GlobalParameters, variable Vhigh, variable Vlow]) // pressure (GPa), temperature (°C) if (ParamIsDefault(GlobalParameters) || GlobalParameters==0) DuanZhang06CreateGlobalParameterWave(degC + 273.15, XCO2) endif Make/free/n=4 w w[0] = degC + 273.15 w[1] = XCO2 w[2] = paramisdefault(Prange) ? GPa > 0.2 : Prange w[3] = 1 // find suitable limits for findroots variable Pbar = GPa * 1e4 variable TK = degC+273.15 variable Rgas = 83.14467 if (paramisdefault(Vlow)) variable dev = Pbar - EOS_DuanZhang06_V(w, 6) Vlow = 6 variable devnew, V for(V=Vlow+2;V<=50;V+=2) devnew = Pbar - EOS_DuanZhang06_V(w, V) if (devnew <= dev) Vlow = V dev = devnew else break endif endfor endif if (paramisdefault(Vhigh)) Vhigh = 1000 variable Vmax = max(Rgas*TK/Pbar,1000) dev = Pbar - EOS_DuanZhang06_V(w, Vlow) for(V=Vlow+2;V<=Vmax;V+=2) devnew = Pbar - EOS_DuanZhang06_V(w, V) if (devnew > dev) Vhigh = V dev = devnew else break endif endfor endif FindRoots/H=(Vhigh)/L=(Vlow)/Q/T=(0.0001)/Z=(GPa*10000) EOS_DuanZhang06_V, w // FindRoots /H=(80)/L=(10+10*(GPa<1)+5*(XCO2>0.5))/Q/T=(0.0001)/Z=(GPa*10000) EOS_DuanZhang06_V, w return V_Flag ? NaN : V_Root // cc/mol end function DuanZhang06VolumeWrapper(wave w, variable Pbar) return DuanZhang06Volume(Pbar*1e-4, w[0], 0, GlobalParameters=w[1]) // for pure water end threadsafe function EOS_DuanZhang06_V(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) return EOS_DuanZhang06(w[0],V, w[1], highP=w[2], globalparameters = w[3]) // bars end // H2O fugacity function DuanZhang06FugacityNumerical(variable GPa, variable degC) // pressure (GPa), temperature (°C) variable TK, Pbar, Rgas, Vint, P0 TK=degC+273.15 // K Pbar=GPa*1e4 // bar Rgas=83.14467 // cc bar/K P0=1e-4 // GPa DuanZhang06CreateGlobalParameterWave(TK, 0) Make /free/n=2 w w[0]=degC w[1]=1 // use global parameters Vint = Integrate1D(DuanZhang06VolumeWrapper, 1, Pbar, 2, 0, w) return P0*exp(Vint/Rgas/TK) // GPa end function DuanZhang06Fugacity(variable GPa, variable degC, variable XCO2, int component) variable Rgas = 83.14467 variable TK = degC + 273.15 variable pBar = GPa * 10000 wave wparam = DuanZhang06CreateGlobalParameterWave(TK, XCO2) int lowP = 0 int highP = 1 // start with low P calculation variable InitialGPa = min(0.2, GPa) variable volume = DuanZhang06Volume(InitialGPa, degC, XCO2, Prange=0, GlobalParameters=1) // cc/mol variable Z = InitialGPa * 1e4 * volume / (Rgas * TK) // for initial calculation, Vref = inf, phiref = 1, Zref = 1 // use equation from Duan et al 1992 // equation 11 of Duan and Zhang 2006 is for EOS where F parameter is split into F and G! variable lnPhi variable dB = wparam[%dBVc][component][lowP] variable dC = wparam[%dCVc2][component][lowP] variable dD = wparam[%dDVc4][component][lowP] variable dE = wparam[%dEVc5][component][lowP] variable dF = wparam[%dFVc2][component][lowP] variable dBeta = wparam[%dvBeta][component][lowP] variable dGamma = wparam[%dgammaVc2][component][lowP] variable FVc2 = wparam[%FVc2][component][lowP] variable bet = wparam[%vBeta][component][lowP] variable gammaVc2 = wparam[%gammaVc2][component][lowP] variable exptermV = exp(-gammaVc2 / volume^2) lnPhi = -ln(Z) lnPhi += dB / volume lnPhi += dC / (2 * volume^2) lnPhi += dD / (4 * volume^4) lnPhi += dE / (5 * volume^5) lnPhi += (dF*bet + dBeta*FVc2) / (2*gammaVc2) * (1 - exptermV) lnPhi += (dF*gammaVc2 + dGamma*FVc2 - FVc2*bet*(dGamma - gammaVc2)) / (2*gammaVc2^2) * (1 - (gammaVc2/volume^2 + 1) * exptermV) lnPhi += - (dGamma - gammaVc2)*FVc2 / (2*gammaVc2^2) * (2 - (gammaVc2^2/volume^4 + 2*gammaVc2/volume^2 + 2) * exptermV) if (InitialGPa == GPa) return exp(lnPhi) * InitialGPa * (component ? XCO2 : 1-XCO2) endif variable lnPhiRef = lnPhi // !!!!!!! recalculate volume using high pressure coefficients variable Vref = DuanZhang06Volume(0.2, degC, XCO2, Prange=1, GlobalParameters=1) variable Zref = 0.2 * 1e4 * Vref / (Rgas * TK) dB = wparam[%dBVc][component][highP] dC = wparam[%dCVc2][component][highP] dD = wparam[%dDVc4][component][highP] dE = wparam[%dEVc5][component][highP] dF = wparam[%dFVc2][component][highP] dBeta = wparam[%dvBeta][component][highP] dGamma = wparam[%dgammaVc2][component][highP] FVc2 = wparam[%FVc2][component][highP] bet = wparam[%vBeta][component][highP] gammaVc2 = wparam[%gammaVc2][component][highP] volume = DuanZhang06Volume(GPa, degC, XCO2, Prange=1, GlobalParameters=1) // cc/mol Z = GPa * 1e4 * volume / (Rgas * TK) exptermV = exp(-gammaVc2 / volume^2) variable exptermVref = exp(-gammaVc2 / Vref^2) lnPhi = lnPhiRef -ln(Z) + ln(Zref) lnPhi += dB / volume - dB / Vref lnPhi += dC / (2 * volume^2) - dC / (2 * Vref^2) lnPhi += dD / (4 * volume^4) - dD / (4 * Vref^4) lnPhi += dE / (5 * volume^5) - dE / (5 * Vref^5) lnPhi += (dF*bet + dBeta*FVc2) / (2*gammaVc2) * (1 - exptermV) lnPhi -= (dF*bet + dBeta*FVc2) / (2*gammaVc2) * (1 - exptermVref) lnPhi += (dF*gammaVc2 + dGamma*FVc2 - FVc2*bet*(dGamma - gammaVc2)) / (2*gammaVc2^2) * (1 - (gammaVc2/volume^2 + 1) * exptermV) lnPhi -= (dF*gammaVc2 + dGamma*FVc2 - FVc2*bet*(dGamma - gammaVc2)) / (2*gammaVc2^2) * (1 - (gammaVc2/Vref^2 + 1) * exptermVref) lnPhi += - (dGamma - gammaVc2)*FVc2 / (2*gammaVc2^2) * (2 - (gammaVc2^2/volume^4 + 2*gammaVc2/volume^2 + 2) * exptermV) lnPhi -= - (dGamma - gammaVc2)*FVc2 / (2*gammaVc2^2) * (2 - (gammaVc2^2/Vref^4 + 2*gammaVc2/Vref^2 + 2) * exptermVref) return exp(lnPhi) * GPa * (component ? XCO2 : 1-XCO2) end // a function to fill all values, to be called once and then passed to all other functions threadsafe function /wave DuanZhang06Parameters(variable TK, variable XCO2) int lowP=0, highP=1, H2O=0, CO2=1 // create coefficients wave for EOS parameters from Table 4 of DZ06 Make/D/free/N=(16,2,2) wcoef=nan SetDimLabel 0, 13, alph, wcoef SetDimLabel 0, 14, bet, wcoef SetDimLabel 0, 15, gam, wcoef wcoef[1][H2O][lowP]={4.38269941E-02,-1.68244362E-01,-2.36923373E-01,1.13027462E-02,-7.67764181E-02} wcoef[6][H2O][lowP]={9.71820593E-02,6.62674916E-05,1.06637349E-03,-1.23265258E-03,-8.93953948E-06} wcoef[11][H2O][lowP]={-3.88124606E-05,5.61510206E-05,7.51274488E-03,2.51598931E+00,3.94000000E-02} wcoef[1][CO2][lowP]={1.14400435E-01,-9.38526684E-01,7.21857006E-01,8.81072902E-03,6.36473911E-02} wcoef[6][CO2][lowP]={-7.70822213E-02,9.01506064E-04,-6.81834166E-03,7.32364258E-03,-1.10288237E-04} wcoef[11][CO2][lowP]={1.26524193E-03,-1.49730823E-03,7.81940730E-03,-4.22918013E+00,1.58500000E-01} wcoef[1][H2O][highP]={4.68071541E-02,-2.81275941E-01,-2.43926365E-01,1.10016958E-02,-3.86603525E-02} wcoef[6][H2O][highP]={9.30095461E-02,-1.15747171E-05,4.19873848E-04,-5.82739501E-04,1.00936000E-06} wcoef[11][H2O][highP]={-1.01713593E-05,1.63934213E-05,-4.49505919E-02,-3.15028174E-01,1.25000000E-02} wcoef[1][CO2][highP]={5.72573440E-03,7.94836769E+00,-3.84236281E+01,3.71600369E-02,-1.92888994E+00} wcoef[6][CO2][highP]={6.64254770E+00,-7.02203950E-06,1.77093234E-02,-4.81892026E-02,3.88344869E-06} wcoef[11][CO2][highP]={-5.54833167E-04,1.70489748E-03,-4.13039220E-01,-8.47988634E+00,2.80000000E-02} // create a parameter wave for parameters in Table 5 of DZ06 make/D/free/N=(15,2,2) w // p is parameter, q is component, r is pressure range setdimlabel 0, 0, null, w setdimlabel 0, 1, BVc, w setdimlabel 0, 2, CVc2, w setdimlabel 0, 3, DVc4, w setdimlabel 0, 4, EVc5, w setdimlabel 0, 5, FVc2, w setdimlabel 0, 6, vBeta, w setdimlabel 0, 7, gammaVc2, w setdimlabel 0, 8, dBVc, w setdimlabel 0, 9, dCVc2, w setdimlabel 0, 10, dDVc4, w setdimlabel 0, 11, dEVc5, w setdimlabel 0, 12, dFVc2, w setdimlabel 0, 13, dvBeta, w setdimlabel 0, 14, dgammaVc2, w w[%null][][] = nan w[0][0][0] = TK variable Rgas = 83.14467 // cc bar/(K mol) make/D/free/N=2 wTc, wPc, wTr, wVc, wX make/D/free/N=(2,2) wB, wC, wD, wE, wF // p is component, q is prange wTc={647.25,304.1282} wPc={221.19,73.773} // bar wTr=TK/wTc wVc=Rgas*wTc/wPc wX={1-xCO2,xCO2} int i, j, k, l, m, n, nc nc=2 make/free/D/N=2 k1, k2, k3 make/free/D/N=(2,2) Vc, Coef // p is component, q is prange k1[lowP] = 3.131 - 5.0624e-3*TK + 1.8641e-6*TK^2 - 31.409/TK k2[lowP] = -46.646 + 4.2877e-2*TK - 1.0892e-5*TK^2 + 1.5782e4/TK k3[lowP] = 0.9 k1[highP] = 9.034 - 7.9212E-3*TK + 2.3285e-6*TK^2 - 2.4221e3/TK k2[highP] = -1.068 + 1.8756e-3*TK - 4.9371e-7*TK^2 + 6.6180e2/TK k3[highP] = 1.0 make/D/free/N=(2,2,2) kij make/D/free/N=(2,2,2,2) kijk wB[][] = wcoef[1][p][q] + wcoef[2][p][q]/wTr[p]^2 + wcoef[3][p][q]/wTr[p]^3 wC[][] = wcoef[4][p][q] + wcoef[5][p][q]/wTr[p]^2 + wcoef[6][p][q]/wTr[p]^3 wD[][] = wcoef[7][p][q] + wcoef[8][p][q]/wTr[p]^2 + wcoef[9][p][q]/wTr[p]^3 wE[][] = wcoef[10][p][q] + wcoef[11][p][q]/wTr[p]^2 + wcoef[12][p][q]/wTr[p]^3 wF[][] = wcoef[%alph][p][q]/wTr[p]^3 if (xCO2 == 0 || xCO2 == 1) int component = XCO2 // for volume calculation we always look for values for the mixture in the zeroth column w[%BVc][][] = wB[component][r]*wVc[component] // makes sense that BVc is BVc for the pure phase w[%CVc2][][] = wC[component][r]*wVc[component]^2 w[%DVc4][][] = wD[component][r]*wVc[component]^4 w[%EVc5][][] = wE[component][r]*wVc[component]^5 w[%FVc2][][] = wF[component][r]*wVc[component]^2 w[%vBeta][][] = wcoef[%bet][component][r] w[%gammaVc2][][] = wcoef[%gam][component][r]*wVc[component]^2 w[%dBVc][][] = 2 * wB[q][r]*wVc[q] w[%dCVc2][][] = 3 * wC[q][r]*wVc[q]^2 w[%dDVc4][][] = 5 * wD[q][r]*wVc[q]^4 w[%dEVc5][][] = 6 * wE[q][r]*wVc[q]^5 w[%dFVc2][][] = 2 * wF[q][r]*wVc[q]^2 w[%dvBeta][][] = wcoef[%bet][q][r] w[%dgammaVc2][][] = 3 * wcoef[%gam][q][r]*wVc[q]^2 else w[%BVc][][] = 0 kij[][][] = p==q ? 1 : k1[r] for(i=0;i 1 ? V_root2 : V_Root // cc/mol end //threadsafe function EOS_RedlichKwongCorrespondingStates(wave w, variable V) // w[0]=T(K), V=volume(cc/mol) // return EOS_RedlichKwong(w[0],V) // bars //end function EOS_RedlichKwongCorrespondingStates(wave w, variable V) // temperature (K), volume (cc/mol), Pc(bars), Tc (K) variable a, b, T, Tc, Pc, R R = 83.1446261815324 // cc-bar/degree/mole T = w[0] // K Pc = w[1] // bars Tc = w[2] // K a = R^2*Tc^2.5/(9*(2^(1/3)-1)*Pc) b = (2^(1/3)-1)*R*Tc/(3*Pc) // same as // a = 0.42748 * R^2 * Tc^2.5 / Pc // b = 0.08664 * R * Tc / Pc return R*T/(V-b)-a/(V*(V+b)*sqrt(T)) //bars end function BelonoshkoSaxenaN2() make /free/n=5 end function DMWvolume(wave w_x, variable GPa, variable degC, [int useK, int verbose]) useK = ParamIsDefault(useK) ? 1 : useK verbose = ParamIsDefault(verbose) ? 0 : verbose wave w_epsilon = DMWepsilonWave() wave w_sigma = DMWsigmaWave() wave mk1 = DMWk1Wave() // a symmetric matrix of interaction parameters wave mk2 = DMWk2Wave() setdimLabel 0, 0, null, w_epsilon, w_sigma, w_x setdimLabel 0, 1, H2O, w_epsilon, w_sigma, w_x setdimLabel 0, 2, CH4, w_epsilon, w_sigma, w_x setdimLabel 0, 3, CO2, w_epsilon, w_sigma, w_x setdimLabel 0, 4, CO, w_epsilon, w_sigma, w_x setdimLabel 0, 5, O2, w_epsilon, w_sigma, w_x setdimLabel 0, 6, N2, w_epsilon, w_sigma, w_x setdimLabel 0, 7, H2, w_epsilon, w_sigma, w_x setdimLabel 0, 8, Cl2, w_epsilon, w_sigma, w_x setdimLabel 0, 9, H2S, w_epsilon, w_sigma, w_x variable epsilon = 0, sigma = 0 int i, j for (i=1; i<=9; i++) for (j=1; j<=9; j++) epsilon += w_x[i] * w_x[j] * (useK ? mk1[i][j] : 1) * sqrt(w_epsilon[i] * w_epsilon[j]) sigma += w_x[i] * w_x[j] * (useK ? mk2[i][j] : 1) * (w_sigma[i] + w_sigma[j]) / 2 endfor endfor variable Pbar = GPa * 1e4 variable TK = degC + 273.15 variable Pm = 3.0626 * sigma^3 * Pbar / epsilon variable Tm = 154 * TK / epsilon Make/free/D w = {Tm} FindRoots /H=(50)/L=(0.000001)/Q/T=(0.0000001)/Z=(Pm) EOS_DMW, w variable Vm = V_numRoots == 2 ? v_root2 : v_root // Vm = max(v_root, v_root2) if (verbose) for(i=1;i<10;i++) if (w_x[i]) print GetDimLabel(w_x, 0, i), w_x[i] endif endfor print "sum", sum(w_x, 1, 9) print epsilon, sigma, Pm, Tm, Vm, Pbar, TK, "Use k =", useK print "roots: Vm", v_root, "V", 1000 * v_root * (sigma/3.691)^3, "Vm", v_root2, "V", 1000 * v_root2 * (sigma/3.691)^3 endif return V_Flag ? NaN : 1000 * Vm * (sigma/3.691)^3 end // Duan, Moller and Weare (1996). A general equation of state for // supercritical fluid mixtures and molecular dynamics simulation of // mixture PVTX properties. Geochimica et Cosmochimica Acta, 60, // 1209-1216. function EOS_DMW(wave w, variable Vm) wave w_a = DMWcoeffientWave() variable Tm = w[0] variable Rgas = 0.08314467 // L bar /K /mol variable Z = 1 Z += (w_a[1] + w_a[2]/Tm^2 + w_a[3]/Tm^3) / Vm Z += (w_a[4] + w_a[5]/Tm^2 + w_a[6]/Tm^3) / Vm^2 Z += (w_a[7] + w_a[8]/Tm^2 + w_a[9]/Tm^3) / Vm^4 Z += (w_a[10] + w_a[11]/Tm^2 + w_a[12]/Tm^3) / Vm^5 Z += w_a[13]/Tm^3/Vm^2 * (1 + w_a[14]/Vm^2) * exp(-w_a[14]/Vm^2) return Rgas * Tm * Z / Vm // Pm, volume of CH4 in l end function/wave DMWcoeffientWave() make/free/N=15/D w_a w_a[00] = NaN w_a[01] = 3.75504388e-02 w_a[02] = -1.08730273e+04 w_a[03] = 1.10964861e+06 w_a[04] = 5.41589372e-04 w_a[05] = 1.12094559e+02 w_a[06] = -5.92191393e+03 w_a[07] = 4.37200027e-06 w_a[08] = 4.95790731e-01 w_a[09] = -1.64902948e+02 w_a[10] = -7.07442825e-08 w_a[11] = 9.65727297e-03 w_a[12] = 4.87945175e-01 w_a[13] = 1.62257402e+04 w_a[14] = 8.99000000e-03 return w_a end function/wave DMWepsilonWave() make/free/D w = {NaN, 510, 154.0, 235.0, 98.0, 115.7, 101.0, 34.6, 348.7, 289.5} return w end function/wave DMWsigmaWave() make/free/D w = {NaN, 2.88, 3.691, 3.69, 3.66, 3.365, 3.63, 2.91, 3.692, 3.693} return w end function/wave DMWk1Wave() make/free/n=(10,10)/D w w = 1 w[1][3] = 0.840; w[3][1] = 0.840 // H2O-CO2 w[2][3] = 0.8563; w[3][2] = 0.8563 // CH4-CO2 w[2][6] = 0.9221; w[6][2] = 0.9221 // CH4-N2 return w end function/wave DMWk2Wave() make/free/n=(10,10)/D w w = 1 w[1][3] = 1.03 w[3][1] = 1.03 return w end function DMWexcessH(wave w_x, variable GPa, variable degC, [int useK, int verbose]) useK = ParamIsDefault(useK) ? 1 : useK verbose = ParamIsDefault(verbose) ? 0 : verbose wave w_a = DMWcoeffientWave() wave w_epsilon = DMWepsilonWave() wave w_sigma = DMWsigmaWave() wave mk1 = DMWk1Wave() // a symmetric matrix of interaction parameters wave mk2 = DMWk2Wave() setdimLabel 0, 0, null, w_epsilon, w_sigma, w_x setdimLabel 0, 1, H2O, w_epsilon, w_sigma, w_x setdimLabel 0, 2, CH4, w_epsilon, w_sigma, w_x setdimLabel 0, 3, CO2, w_epsilon, w_sigma, w_x setdimLabel 0, 4, CO, w_epsilon, w_sigma, w_x setdimLabel 0, 5, O2, w_epsilon, w_sigma, w_x setdimLabel 0, 6, N2, w_epsilon, w_sigma, w_x setdimLabel 0, 7, H2, w_epsilon, w_sigma, w_x setdimLabel 0, 8, Cl2, w_epsilon, w_sigma, w_x setdimLabel 0, 9, H2S, w_epsilon, w_sigma, w_x variable epsilon = 0, sigma = 0 int i, j for (i=1; i<=9; i++) for (j=1; j<=9; j++) epsilon += w_x[i] * w_x[j] * (useK ? mk1[i][j] : 1) * sqrt(w_epsilon[i] * w_epsilon[j]) sigma += w_x[i] * w_x[j] * (useK ? mk2[i][j] : 1) * (w_sigma[i] + w_sigma[j]) / 2 endfor endfor variable Pbar = GPa * 1e4 variable TK = degC + 273.15 variable Rgas = 0.08314467 // L bar /K /mol variable Pm = 3.0626 * sigma^3 * Pbar / epsilon variable Tm = 154 * TK / epsilon Make/free/D w = {Tm} FindRoots /H=(50)/L=(0.000001)/Q/T=(0.0000001)/Z=(Pm) EOS_DMW, w variable Vm = V_numRoots == 2 ? v_root2 : v_root // Vm = max(v_root, v_root2) if (verbose) for(i=1;i<10;i++) if (w_x[i]) print GetDimLabel(w_x, 0, i), w_x[i] endif endfor print "sum", sum(w_x, 1, 9) print epsilon, sigma, Pm, Tm, Vm, Pbar, TK, "Use k =", useK print "roots: Vm", v_root, "V", 1000 * v_root * (sigma/3.691)^3, "Vm", v_root2, "V", 1000 * v_root2 * (sigma/3.691)^3 endif variable ZZ = Pm * Vm / (Rgas * Tm) variable deltaH deltaH = (2*w_a[02]/Tm^3 + 3*w_a[03]/Tm^4) / Vm deltaH += (2*w_a[05]/Tm^3 + 3*w_a[06]/Tm^4) / (2 * Vm^2) deltaH += (2*w_a[08]/Tm^3 + 3*w_a[09]/Tm^4) / (4 * Vm^4) deltaH += (2*w_a[11]/Tm^3 + 3*w_a[12]/Tm^4) / (5 * Vm^5) deltaH += 3 * w_a[13] / (2*w_a[14]*Tm^4) * (2 - (2 + w_a[14]/Vm^2) * exp (-w_a[14]/Vm^2)) deltaH *= epsilon / 154 * Rgas * TK^2 deltaH += epsilon / 154 * Rgas * TK *(ZZ - 1) return deltaH end // H2O //Quantity Value Units Method Reference Comment //ΔfH°gas -241.826 ± 0.040 kJ/mol Review Cox, Wagman, et al., 1984 CODATA Review value //ΔfH°gas -241.83 kJ/mol Review Chase, 1998 Data last reviewed in March, 1979 //Quantity Value Units Method Reference Comment //S°gas,1 bar 188.835 ± 0.010 J/mol*K Review Cox, Wagman, et al., 1984 CODATA Review value //S°gas,1 bar 188.84 J/mol*K Review Chase, 1998 Data last reviewed in March, 1979 //Gas Phase Heat Capacity (Shomate Equation) // //Cp° = A + B*t + C*t2 + D*t3 + E/t2 //H° − H°298.15= A*t + B*t2/2 + C*t3/3 + D*t4/4 − E/t + F − H //S° = A*ln(t) + B*t + C*t2/2 + D*t3/3 − E/(2*t2) + G // Cp = heat capacity (J/mol*K) // H° = standard enthalpy (kJ/mol) // S° = standard entropy (J/mol*K) // t = temperature (K) / 1000. function EOS_ZD09(wave w, variable Vm) wave w_a = ZD09coefficientWave() variable Tm = w[0] variable Rgas = 0.08314467 // L bar /K /mol variable Z = 1 Z += (w_a[1] + w_a[2]/Tm^2 + w_a[3]/Tm^3) / Vm Z += (w_a[4] + w_a[5]/Tm^2 + w_a[6]/Tm^3) / Vm^2 Z += (w_a[7] + w_a[8]/Tm^2 + w_a[9]/Tm^3) / Vm^4 Z += (w_a[10] + w_a[11]/Tm^2 + w_a[12]/Tm^3) / Vm^5 Z += w_a[13]/Tm^3/Vm^2 * (w_a[14] + w_a[15]/Vm^2) * exp(-w_a[15]/Vm^2) return Rgas * Tm * Z / Vm // Pm, volume of CH4 in l end function/wave ZD09coefficientWave() make/free/N=16/D w_a w_a[00] = NaN w_a[01] = 2.95177298930e-2 w_a[02] = -6.33756452413e+3 w_a[03] = -2.75265428882e+5 w_a[04] = 1.29128089283e-3 w_a[05] = -1.45797416153e+2 w_a[06] = 7.65938947237e+4 w_a[07] = 2.58661493537e-6 w_a[08] = 0.52126532146e+0 w_a[09] = -1.39839523753e+2 w_a[10] = -2.36335007175e-8 w_a[11] = 5.35026383543e-3 w_a[12] = -0.27110649951e+0 w_a[13] = 2.50387836486e+4 w_a[14] = 0.73226726041e+0 w_a[15] = 1.54833359970e-2 return w_a end function/wave ZD09epsilonWave() make/free/D w = {NaN, 154.0, 510.0, 235.0, 31.2, 105.6, 124.5, 246.1} return w end // epsilon/kB in K function/wave ZD09sigmaWave() make/free/D w = {NaN, 3.691, 2.88, 3.79, 2.93, 3.66, 3.36, 4.35} return w end // 10^-10 m function/wave ZD09k1Wave() make/free/n=(8,8)/D w w = 1 w[2][3] = 0.85; w[3][2] = 0.85 // H2O-CO2 w[1][2] = 0.8; w[2][1] = 0.80// CH4-H2O return w end function/wave ZD09k2Wave() make/free/n=(8,8)/D w w = 1 w[2][3] = 1.02; w[3][2] = 1.02 // H2O-CO2 return w end function ZD09Volume(wave w_x, variable GPa, variable degC, [int useK, int verbose]) useK = ParamIsDefault(useK) ? 1 : useK verbose = ParamIsDefault(verbose) ? 0 : verbose wave w_epsilon = ZD09epsilonWave() wave w_sigma = ZD09sigmaWave() wave mk1 = ZD09k1Wave() // a symmetric matrix of interaction parameters wave mk2 = ZD09k2Wave() setdimLabel 0, 0, null, w_epsilon, w_sigma, w_x setdimLabel 0, 1, CH4, w_epsilon, w_sigma, w_x setdimLabel 0, 2, H2O, w_epsilon, w_sigma, w_x setdimLabel 0, 3, CO2, w_epsilon, w_sigma, w_x setdimLabel 0, 4, H2, w_epsilon, w_sigma, w_x setdimLabel 0, 5, CO, w_epsilon, w_sigma, w_x setdimLabel 0, 6, O2, w_epsilon, w_sigma, w_x setdimLabel 0, 7, C2H6, w_epsilon, w_sigma, w_x variable epsilon = 0, sigma = 0 int i, j for (i=1; i<=7; i++) for (j=1; j<=7; j++) epsilon += w_x[i] * w_x[j] * (useK ? mk1[i][j] : 1) * sqrt(w_epsilon[i] * w_epsilon[j]) sigma += w_x[i] * w_x[j] * (useK ? mk2[i][j] : 1) * (w_sigma[i] + w_sigma[j]) / 2 endfor endfor variable Pbar = GPa * 1e4 variable TK = degC + 273.15 variable Pm = 3.0626 * sigma^3 * Pbar / epsilon variable Tm = 154 * TK / epsilon Make/free/D w = {Tm} FindRoots /H=(50)/L=(0.0001)/Q/T=(0.0000001)/Z=(Pm) EOS_ZD09, w variable Vm = V_numRoots == 2 ? v_root2 : v_root // Vm = max(v_root, v_root2) if (verbose) for (i=1;i<=7;i++) if (w_x[i]) print GetDimLabel(w_x, 0, i), w_x[i] endif endfor print "sum", sum(w_x, 1, 7) print epsilon, sigma, Pm, Tm, Vm, Pbar, TK, "Use k =", useK print "roots: Vm", v_root, "V", 1000 * v_root * (sigma/3.691)^3, "Vm", v_root2, "V", 1000 * v_root2 * (sigma/3.691)^3 endif return V_Flag ? NaN : 1000 * Vm * (sigma/3.691)^3 end function/wave ZD09fugacity(wave w_x, variable GPa, variable degC, [int useK, int verbose]) useK = ParamIsDefault(useK) ? 1 : useK verbose = ParamIsDefault(verbose) ? 0 : verbose variable VL = ZD09volume(w_x, GPa, degC, useK=useK, verbose=verbose) / 1000 // litre variable TK = degC + 273.15 variable Pbar = GPa * 1e4 variable Rgas = 0.08314467 // L bar /K /mol make/N=10/free/D w_epsilon, w_sigma wave w_epsilon = ZD09epsilonWave() wave w_sigma = ZD09sigmaWave() wave w_a = ZD09coefficientWave() // lngamma = C1 - {C2} * (2*Tm / epsilon * (epsilon - epsilon1)) + 6/sigma * (1-Z) * (sigma - sigma1) int i, j variable Z = Pbar * VL / Rgas / TK variable epsilon, sigma, S1, S2 duplicate/free w_x, w_f, w_epsilon1, w_sigma1 w_f = 0 w_epsilon1 = 0 // w_sigma1 = 0 wave mk1 = ZD09k1Wave() wave mk2 = ZD09k2Wave() for (i=1; i<=7; i++) w_epsilon1 += (useK ? mk1[p][i] : 1) * w_x[i] * sqrt(w_epsilon[p] * w_epsilon[i]) w_sigma1 += (useK ? mk2[p][i] : 1) * w_x[i] * (w_sigma[p] + w_sigma[i]) / 2 for (j=1; j<=7; j++) epsilon += w_x[i] * w_x[j] * (useK ? mk1[i][j] : 1) * sqrt(w_epsilon[i] * w_epsilon[j]) sigma += w_x[i] * w_x[j] * (useK ? mk2[i][j] : 1) * (w_sigma[i] + w_sigma[j]) / 2 endfor endfor variable Pm = 3.0626 * sigma^3 * Pbar / epsilon variable Tm = 154 * TK / epsilon variable Vm = VL * (3.691/sigma)^3 S1 = (w_a[1] + w_a[2]/Tm^2 + w_a[3]/Tm^3) / Vm S1 += (w_a[4] + w_a[5]/Tm^2 + w_a[6]/Tm^3) / (2 * Vm^2) S1 += (w_a[7] + w_a[8]/Tm^2 + w_a[9]/Tm^3) / (4 * Vm^4) S1 += (w_a[10] + w_a[11]/Tm^2 + w_a[12]/Tm^3) / (5 * Vm^5) S1 += w_a[13] / (2*w_a[15]*Tm^3) * (w_a[14] + 1 - (w_a[14] + 1 + w_a[15]/Vm^2) * exp (-w_a[15]/Vm^2)) S2 = (2*w_a[02]/Tm^2 + 3 * w_a[03]/Tm^3) / Vm S2 += (2*w_a[05]/Tm^2 + 3 * w_a[06]/Tm^3) / (2 * Vm^2) S2 += (2*w_a[08]/Tm^2 + 3 * w_a[09]/Tm^3) / (4 * Vm^4) S2 += (2*w_a[11]/Tm^2 + 3 * w_a[12]/Tm^3) / (5 * Vm^5) S2 += 3 * w_a[13] / (2*w_a[15]*Tm^3) * (w_a[14] + 1 - (w_a[14] + 1 - w_a[15]/Vm^2) * exp (-w_a[15]/Vm^2)) // ln(gamma) NOTE: -2*S2*..., not +2*S2*...! w_f = Z - 1 - ln(Z) + S1 - 2*S2*(1 - w_epsilon1[p] / epsilon) + 6 * (1-Z) * (1 - w_sigma1[p] / sigma) // fugacity in bars w_f = Pbar * exp(w_f[p]) * w_x[p] if (verbose) print w_f endif if (verbose) for(i=1;i<=7;i++) if (w_x[i]) print GetDimLabel(w_x, 0, i), w_x[i], w_epsilon1[i], w_sigma1[i], w_f[i] endif endfor print epsilon, sigma, Pm, Tm, Vm, VL*1000, Pbar, TK print S1, S2, w_epsilon1[4], w_epsilon1[2], epsilon print (1 - w_epsilon1[4] / epsilon), 6 * (1-Z) * (1 - w_sigma1[4] / sigma) endif return w_f // partial fugacities in bar end function/wave ZD09lnGamma(wave w_x, variable GPa, variable degC, [int useK, int verbose]) useK = ParamIsDefault(useK) ? 1 : useK verbose = ParamIsDefault(verbose) ? 0 : verbose variable VL = ZD09volume(w_x, GPa, degC, useK=useK, verbose=verbose) / 1000 // litre variable TK = degC + 273.15 variable Pbar = GPa * 1e4 variable Rgas = 0.08314467 // L bar /K /mol make/N=10/free/D w_epsilon, w_sigma wave w_epsilon = ZD09epsilonWave() wave w_sigma = ZD09sigmaWave() wave w_a = ZD09coefficientWave() // lngamma = C1 - {C2} * (2*Tm / epsilon * (epsilon - epsilon1)) + 6/sigma * (1-Z) * (sigma - sigma1) int i, j variable Z = Pbar * VL / Rgas / TK variable epsilon, sigma, S1, S2 duplicate/free w_x, w_f, w_epsilon1, w_sigma1 w_f = 0 w_epsilon1 = 0 // w_sigma1 = 0 wave mk1 = ZD09k1Wave() wave mk2 = ZD09k2Wave() for (i=1; i<=7; i++) w_epsilon1 += (useK ? mk1[p][i] : 1) * w_x[i] * sqrt(w_epsilon[p] * w_epsilon[i]) w_sigma1 += (useK ? mk2[p][i] : 1) * w_x[i] * (w_sigma[p] + w_sigma[i]) / 2 for (j=1; j<=7; j++) epsilon += w_x[i] * w_x[j] * (useK ? mk1[i][j] : 1) * sqrt(w_epsilon[i] * w_epsilon[j]) sigma += w_x[i] * w_x[j] * (useK ? mk2[i][j] : 1) * (w_sigma[i] + w_sigma[j]) / 2 endfor endfor variable Pm = 3.0626 * sigma^3 * Pbar / epsilon variable Tm = 154 * TK / epsilon variable Vm = VL * (3.691/sigma)^3 S1 = (w_a[1] + w_a[2]/Tm^2 + w_a[3]/Tm^3) / Vm S1 += (w_a[4] + w_a[5]/Tm^2 + w_a[6]/Tm^3) / (2 * Vm^2) S1 += (w_a[7] + w_a[8]/Tm^2 + w_a[9]/Tm^3) / (4 * Vm^4) S1 += (w_a[10] + w_a[11]/Tm^2 + w_a[12]/Tm^3) / (5 * Vm^5) S1 += w_a[13] / (2*w_a[15]*Tm^3) * (w_a[14] + 1 - (w_a[14] + 1 + w_a[15]/Vm^2) * exp (-w_a[15]/Vm^2)) S2 = (2*w_a[02]/Tm^2 + 3 * w_a[03]/Tm^3) / Vm S2 += (2*w_a[05]/Tm^2 + 3 * w_a[06]/Tm^3) / (2 * Vm^2) S2 += (2*w_a[08]/Tm^2 + 3 * w_a[09]/Tm^3) / (4 * Vm^4) S2 += (2*w_a[11]/Tm^2 + 3 * w_a[12]/Tm^3) / (5 * Vm^5) S2 += 3 * w_a[13] / (2*w_a[15]*Tm^3) * (w_a[14] + 1 - (w_a[14] + 1 - w_a[15]/Vm^2) * exp (-w_a[15]/Vm^2)) // ln(gamma) NOTE: -2*S2*..., not +2*S2*...! w_f = Z - 1 - ln(Z) + S1 - 2*S2*(1 - w_epsilon1[p] / epsilon) + 6 * (1-Z) * (1 - w_sigma1[p] / sigma) return w_f end function/wave DMWfugacity(wave w_x, variable GPa, variable degC, [int useK, int verbose]) useK = ParamIsDefault(useK) ? 1 : useK verbose = ParamIsDefault(verbose) ? 0 : verbose variable VL = DMWvolume(w_x, GPa, degC, useK=useK, verbose=verbose) / 1000 // litre variable TK = degC + 273.15 variable Pbar = GPa * 1e4 variable Rgas = 0.08314467 // L bar /K /mol make/N=10/free/D w_epsilon, w_sigma wave w_epsilon = DMWepsilonWave() wave w_sigma = DMWsigmaWave() wave w_a = DMWcoeffientWave() // lngamma = Z - 1 - ln(Z) + C1 - {C2} * (2*Tm / epsilon * (epsilon - epsilon1)) + 6/sigma * (1-Z) * (sigma - sigma1) int i, j variable Z = Pbar * VL / Rgas / TK variable epsilon, sigma, C1, C2 duplicate/free w_x, w_f, w_epsilon1, w_sigma1 w_f = 0 w_epsilon1 = 0 // w_sigma1 = 0 wave mk1 = DMWk1Wave() wave mk2 = DMWk2Wave() for (i=1; i<=9; i++) w_epsilon1 += (useK ? mk1[p][i] : 1) * w_x[i] * sqrt(w_epsilon[p] * w_epsilon[i]) w_sigma1 += (useK ? mk2[p][i] : 1) * w_x[i] * (w_sigma[p] + w_sigma[i]) / 2 for (j=1; j<=9; j++) epsilon += w_x[i] * w_x[j] * (useK ? mk1[i][j] : 1) * sqrt(w_epsilon[i] * w_epsilon[j]) sigma += w_x[i] * w_x[j] * (useK ? mk2[i][j] : 1) * (w_sigma[i] + w_sigma[j]) / 2 endfor endfor variable Pm = 3.0626 * sigma^3 * Pbar / epsilon variable Tm = 154 * TK / epsilon variable Vm = VL * (3.691/sigma)^3 // C1 = Z - 1 - ln(Z) C1 = (w_a[1] + w_a[2]/Tm^2 + w_a[3]/Tm^3) / Vm C1 += (w_a[4] + w_a[5]/Tm^2 + w_a[6]/Tm^3) / (2 * Vm^2) C1 += (w_a[7] + w_a[8]/Tm^2 + w_a[9]/Tm^3) / (4 * Vm^4) C1 += (w_a[10] + w_a[11]/Tm^2 + w_a[12]/Tm^3) / (5 * Vm^5) // a14 or a4? a4 must be a typo C1 += w_a[13] / (2*w_a[14]*Tm^3) * (2 - (2 + w_a[14]/Vm^2) * exp (-w_a[14]/Vm^2)) C2 = (2*w_a[02]/Tm^2 + 3 * w_a[03]/Tm^3) / Vm C2 += (2*w_a[05]/Tm^2 + 3 * w_a[06]/Tm^3) / (2 * Vm^2) C2 += (2*w_a[08]/Tm^2 + 3 * w_a[09]/Tm^3) / (4 * Vm^4) C2 += (2*w_a[11]/Tm^2 + 3 * w_a[12]/Tm^3) / (5 * Vm^5) C2 += 3 * w_a[13] / (2*w_a[14]*Tm^3) * (2 - (2 + w_a[14]/Vm^2) * exp (-w_a[14]/Vm^2)) w_f = Z - 1 - ln(Z) + C1 - 2 * C2 * (1 / epsilon * (epsilon - w_epsilon1[p])) + 6/sigma * (1-Z) * (sigma - w_sigma1[p]) if (verbose) print w_f endif w_f = Pbar * exp(w_f[p]) * w_x[p] if (verbose) for(i=1;i<10;i++) if (w_x[i]) print GetDimLabel(w_x, 0, i), w_x[i], w_epsilon1[i], w_sigma1[i], w_f[i] endif endfor print epsilon, sigma, Pm, Tm, Vm, VL*1000, Pbar, TK print C1, C2, w_epsilon1[7], w_epsilon1[1], epsilon print (1 / epsilon * (epsilon - w_epsilon1[7])), 6/sigma * (1-Z) * (sigma - w_sigma1[7]) endif return w_f end