Accurate quantum-corrected cubic equations of state for helium, neon, hydrogen, deuterium and their mixtures

ContentslistsavailableatScienceDirect

Fluid Phase Equilibria

journalhomepage:www.elsevier.com/locate/fluid

Accurate quantum-corrected cubic equations of state for helium, neon, hydrogen, deuterium and their mixtures

Ailo Aasen

^a

, Morten Hammer

^a

, Silvia Lasala

^b

, Jean-Noël Jaubert

^b

, Øivind Wilhelmsen

^a,c,∗

aSINTEF Energy Research, Trondheim NO-7465, Norway

bUniversité de Lorraine, Ecole Nationale Supérieure des Industries Chimiques, Laboratoire Ré et Génie des Procé dés (UMR CNRS 7274), 1 rue Grandville, Nancy 540 0 0, France

cNorwegian University of Science and Technology, Department of Energy and Process Engineering, Kolbjørn Hejes vei 1B, Trondheim NO-7491, Norway

a rt i c l e i n f o

Article history:

Received 1 June 2020 Revised 8 August 2020 Accepted 19 August 2020 Available online 20 August 2020 Keywords:

Cubic equation of state Quantum corrections Hydrogen

Helium Deuterium Neon

a b s t r a c t

Cubicequationsofstatehavethusfaryieldedpoorpredictionsofthethermodynamicpropertiesofquan- tumfluidssuchashydrogen,heliumanddeuteriumatlowtemperatures.Furthermore,theshapeofthe optimalα^{functionsofheliumand hydrogenhavebeenshowntonotdecay}monotonicallyasforother fluids.Inthiswork,wederivetemperature-dependentquantumcorrectionsforthecovolumeparameter ofcubicequationsofstatebymappingthemontotheexcludedvolumespredictedbyquantum-corrected Miepotentials.SubsequentregressionoftheTwuα ^{functionrecoversanearclassicalbehaviorwitha} monotonicdecayformostofthe temperaturerange.Thequantum correctionsresult inasignificantly betteraccuracy,especiallyforcaloricproperties.Whiletheaveragedeviationoftheisochoricheatcapac- ityofliquidhydrogenatsaturationexceeds80%withthepresentstate-of-the-art,theaveragedeviation is4%with quantumcorrections.Average deviationsfor thesaturationpressure arewellbelow 1%for allfourfluids.UsingPénelouxvolumeshiftsgivesaverageerrorsinsaturationdensitiesthatarebelow 2%forheliumandabout1%forhydrogen,deuteriumandneon.Parametersarepresentedfortwocubic equationsofstate:Peng–RobinsonandSoave–Redlich–Kwong.Thequantum-correctedcubicequationsof statearealsoabletoreproducethevapor-liquidequilibriumofbinarymixturesofquantumfluids,and theyarethefirstcubicequationsofstatethatareabletoaccuratelymodelthevapor-liquidequilibrium ofthehelium–neonmixture.Similartothequantum-correctedMiepotentialsthatwereusedtodevelop thecovolumecorrections,aninteractionparameterforthecovolumeisneededtorepresentthehelium–

hydrogenmixturetoahighaccuracy.Thequantum-correctedcubicequationofstatepavesthewayfor technologicalapplicationsofquantumfluidsthatrequiremodelswithbothhighaccuracyandcomputa- tionalspeed.

© 2020TheAuthors.PublishedbyElsevierB.V.

ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Since the first cubic equation of state (EoS) was introduced by van der Waals [1] more than 100 years ago, cubic EoS have remained the preferred choice when computational speed in combinationwithareasonableaccuracyisneeded.Theycanallbe formulatedintermsofthepressure

P

(

^T,

v )

= RT

v

−b− a

( v

−r1b

)( v

−r2b

)

^{, (1)}

∗ Corresponding author at: Norwegian University of Science and Technology, De- partment of Energy and Process Engineering, Kolbjørn Hejes vei 1B, NO-7491 Trond- heim, Norway.

E-mail address: oivind.wilhelmsen@ntnu.no (Ø. Wilhelmsen).

where T and v are the temperature andmolar volume, R is the gasconstant,andr₁andr₂areconstantsthatdefinetheparticular EoS.Thecovolume,b,andtheattractiveparameter,a,arefunctions oftemperatureandmixturecomposition.

Cubic EoS are usually preferred for computationally demand- ingprocesssimulations,optimizationstudies,andincomputational fluid dynamics [2,3]. Twoof the mostfrequently used cubic EoS are Soave–Redlich–Kwong (SRK) andPeng–Robinson (PR),named aftertheauthors ofthepaperswhere theseEoS wereintroduced [4,5]. Avariety of modifications havecontributedto improve the accuracyandwidespreaduseofcubicEoS,fornon-polarfluids[6], polarfluids[7],electrolytes[8]andmanyotherapplications[9,10]. Onetype offluidsthatsofar haseludedthe successfulmodeling by cubic EoS isfluids that exhibit vapor-liquid equilibrium(VLE) atverylowtemperatures,suchashydrogenandhelium.Sincehy- https://doi.org/10.1016/j.fluid.2020.112790

0378-3812/© 2020 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

drogenhasthepotentialtobecomeawidespreadenvironmentally friendlyenergycarrierinthefuture,theinterest insuch fluidsis increasing [11]. Forinstance, liquefaction is beinginvestigated as analternative forlarge-scale distributionofhydrogen acrosslong distances[12].Mixtures consistingofhelium, neonandhydrogen havebeen touted aspromising refrigerants withthe potential to significantlyimprovetheenergyefficiencyof thehydrogenlique- factionprocess[12,13].Unfortunately,even withadvanced mixing rules [14], cubic EoS have been unable to reproduce the VLE of helium–neon[12].The reasonisthatthesefluidsare stronglyin- fluenced by quantum-mechanical effects,as quantified by the de Brogliethermalwavelength

λ

B[15,16]:

λ

^B=

^h

2

π

^mkBT, (2)

where h is Planck’s constant, m is the particle mass, and k_B is Boltzmann’sconstant. When

λ

B becomescomparable tothe typ- icalintermolecular separation, quantum effectscan influencethe fluidpropertiessignificantly.

Recently, the statistical associating fluid theory of quantum- corrected Mie potentials (SAFT-VRQ Mie) [17,18] was presented.

This EoS models quantum effects by adding the so-called Feynman–Hibbs-correctionstoaMiepotential,wheretheMiepo- tentialrepresentstheclassicalbehaviorofthefluid[19].Excellent agreementwithexperimentswasobtainedfortheVLEofhelium–

neon, except in the critical region due to the incomplete devel- opmentofthermodynamic perturbationtheory formixtures [20]. However, molecular simulations with the underlying interaction potentials of the EoS were in excellent agreement with the ex- periments,alsointhecriticalregion.Theinteractionpotentialsre- vealedthat a crucial effectto account foris “quantumswelling”, namelythattheeffectivesizeoftheparticlesincreasesatlowtem- peratures due to the wave-like nature of particles that becomes moreapparentas

λ

Bincreases.IncubicEoS,thesizeoftheparti- clesisrepresentedbythecovolumeb inEq.(1),whichisusually assumed to be temperature-independent and therefore does not accountforquantumswelling.Inthiswork,we shallderivequan- tumcorrectionsforthecovolumeparameterofcubicEoSbasedon thequantumswellingoftheFeynman–Hibbs-correctedMiepoten- tialsthatweresuccessfulinpreviouswork.

The shortcomings of cubic EoS for describing fluids that ex- hibitquantumeffects hasforlongbeen known.Inhis celebrated workthatformedthebasisfortheSRKEoS[4],Soavepointedout thatlessaccurate resultswereobtainedformixtures withhydro- gen.SomeeffortshavebeendirectedtowardsimprovingcubicEoS forsuch applications.Thomas etal.[21]evaluated theapplicabil- ityofseveralEoSformodelingpurehelium,demonstratingthatPR ingeneralgivesinaccuratepredictionsatlowtemperatures,except forpressuresbelow5bar.Kavianietal.[22] madeacubicEoSfor purehelium by fittingfive-parameter correlationsforboth the

α

functionand thecovolume, which were fittedto thermodynamic propertiesbelow15 K andup to 16 bar. Althoughhighaccuracy wasobtainedatlowtemperatures,theirEoSisnot applicablebe- yondtherangewhereitwasregressed.

Tkaczuk et al. [23] recently developed a Helmholtz energy mixture model for the helium–neon mixture. Their model relies onmultiparameterpure fluid EoScombinedwithmultiparameter mixturemodels.Tkaczuketal.usemorethan40binaryinteraction parameters to regress the properties of the binary helium–neon mixturemodel.Theirrepresentationoffluidpropertiesistherefore expectedto be much more computationally demanding than the descriptionprovidedinthiswork.

TheattractiveparameterinEq.(1),a,isusuallyassumedtode- pendontemperature,andiswrittenastheproductofitsvalueat thecriticalpoint,acandthe

α

^{function,as:a}(^T)=ac

α

(^T)^.Which functionalformtouseforthe

α

^{functionhasbeenwidelydebated.}

Fig. 1. Pair interaction potential for the Lennard-Jones (LJ) potential (black, solid line) and the LJ potential with second order Feynman–Hibbs (FH) quantum correc- tions (red, dashed). The temperature is given by k BT /= 1 , and the particle mass is that of helium. The FH correction increases the effective diameter σeff, and de- creases the well depth eff, compared to the classical LJ parameters, σand . The inset shows that the effective range of the potential increases slightly.

Recently,LeGuennecetal.[24]introducedconsistencycriteriafor the

α

^{functionsofcubicEoS.Theyshowedthatbyenforcingthese}

criteria,thepredictionsinthesupercriticaldomaincanbesignifi- cantlyimproved,withanegligiblelossofaccuracyinthesubcrit- ical domain [25].One of the consistency criteria wasthat the

α

function should exhibit a monotonic decayasa function oftem- perature.Following theserecommendations, cubicEoS havebeen reparametrized forthousands ofcomponents usingconsistent al- phafunctionsfittedtosaturationdata,includingforhelium,neon, hydrogenanddeuterium[6,26].LeGuennecetal.[24]showedthat the

α

^{function of hydrogen andhelium did not exhibit a mono-}

tonicdecaywithtemperature, anddidthusnotfollowtheconsis- tency criteriaderived forother fluids. The reasonforthisbehav- iorhasremainedunexplained,althoughithasbeenlinkedtotheir acentricfactorsbeingnegative[6].Wewillshowthatanearclas- sical behaviorwithamonotonic decaycan beobtainedforthe

α

functionbyincorporatingtemperature-dependentquantumcorrec- tionsforthecovolumeparameterthatexplicitlyaccountforquan- tumswelling.

Thephysicalvalidityofatemperature-dependentcovolumehas beenquestionedbyKalikhmanetal.[27],who demonstratedthat any repulsive term involving temperature-dependent covolumes resultsinanegativeinfinitevaluefortheisochoricheatcapacityat infinitepressure.We willshow that thequantum correctionsde- rivedinthisworkwillintroducenegativeheatcapacitiesinagree- mentwiththederivationsbyKalikhman,butonlyatpressureswell outsideregionsofindustrialapplicationsofthesubstance.Wealso verify that thermodynamic stability criteriaare fulfilled and that pressure–volumeisothermsdonotcrossclosetothecriticalpoint forthecovolumecorrectionsderived inthiswork.Parametersfor thequantumfluidswillbeprovided forthemostfrequentlyused cubic EoS, SRK and PR, butcan easily be extended to other cu- bicEoS.Thequantumcorrectionsyieldavastimprovementinthe accuracyofcubic EoS,bothforsingle-component fluidsandmix- tures.TheVLEofall examinedfluidmixturesare reproducedtoa highaccuracywithouttheneedforadvancedmixingrules[28].

2. Theory

We showed in previous work that Feynman–Hibbs-corrected Miepotentialsare capableofcapturingthethermodynamicprop- ertiesoffluidswithquantumeffectstoa relativelyhighaccuracy [17,18]. As shown in Fig. 1, the quantum corrections influence a classicalinteractionpotentialmainlythroughtwoeffects:

1. Theyincrease the distancewhere the classicalinteraction po- tential is zero,

σ

^{, to a larger,} temperature-dependent value,

σ

eff(^T)^.

2. Theyreduce the well-depth,

^{ofthe classical} interaction po- tentialtoasmaller,temperature-dependentvalue

eff(^T)^. Numerous studies have shown that cubic EoS are capable of representing the thermodynamic properties of classical fluids to a reasonable accuracy. The main hypothesesof thiswork are (1) thatthequantumswellingcanbeincorporatedintothecubicEoS through a temperature-dependent covolume correction, and (2) that the reduction in the well-depth can be captured by a suit- ablemodificationofthe

α

^{function.Inthefollowingwewillderive}

functional forms suitable to describe these quantum corrections based on the Feynman–Hibbs-corrected Mie potentials described inRefs.[17,18].

2.1. ThecovolumeofcubicEoS

Van derWaalsderived hisfamous equation ofstateunderthe assumption that the excluded volume per particle,known asthe covolume,isfourtimestheparticlevolume:

b=1 2

4

πσ

³

3

, (3)

where

σ

^{is the particle diameter. The covolume, b, used in the}

equation ofstate canbe found bymultiplying b withAvogadro’s number,

b=N_Ab. (4)

FromthereducedformofthevanderWaalsequationofstate,the covolumecan be relatedto thecritical propertiesofthe fluid by imposingthat (^dP/dV)Tc=0and(^d2P/dV²)Tc=0,whereTc isthe criticaltemperature.Thisgivesthefollowingexpressions:

bvdW =0.125T_critR

Pcrit , (5)

bSRK =0.08664T_critR

P_crit , (6)

bPR =0.07780T_critR

P_crit , (7)

forthecovolumesofthevanderWaals,theSRKandthePRequa- tionsofstaterespectively.

2.2. RelatingthecovolumeofcubicEoStointeractionpotentialsof classicalfluids

The intermolecular potential between two nearly spherical moleculescanbeapproximatedbyaMiepotential:

u_Mie

(

^r

)

^=C

σ

r

_λr

−

σ

r

_λa

, (8)

where

^{is thewell depth,}

σ

^{is the finitedistance atwhich the}

potential iszero,

λ

^{a and}

λ

^{r aretheattractive andrepulsive expo-}

nents,and C=

λ

^r

λ

^r−

λ

^a

λ

^r

λ

^a

_λr^λ₋^a_λa

. (9)

Forimpenetrablespheresasrepresentedbythehardspherepo- tential,theexcludedvolumeisgivenuniquelyasthevolumeofthe spheres.Thedefinitionof“excludedvolume” correspondingtothe repulsivepartoftheMiepotentialisnotunique,asthefiniteop- posingforce ofthe potentialatradiismaller than

σ

^{allows parti-}

clestobeatintermoleculardistancessmallerthan

σ

^.However,as

Fig. 2. The covolume b of SRK and PR divided by the excluded volume of the Mie potential models b Mie in Ref. [30] , plotted as a function of molecular weight for methane, argon, krypton and xenon.

therepulsivepartoftheinteractionpotentialisusuallyverysteep, theexcludedvolume canbeapproximatedby thepositivepartof thepotential.ForMiepotentials,theparticlediameteristhensim- ply

σ

^{inEq.(8),wherethe}interactionpotentialgoesfrompositive tonegative.Fig.2showsthattheexcludedvolumeoftheMiepo- tential,givenbyb_Mie=

πσ

³/6isproportionaltothecovolumepa- rameterbofcubicEoSforclassicalfluidsforbothSRKandPR.This observation,andthefactthattheproportionalityfactoriscloseto 0.5,hasbeenmadebefore[29].

2.3.Interactionpotentialsandquantumcorrections

AssumingthattheMiepotentialcanbeusedtodescribeclassi- calfluids,Feynman–Hibbscorrectionscanbeaddedtotheclassical potentialtorepresentquantumeffects,

u_MieFH

(

^r

)

^/

(

^C

)

^{= σ}r

λ^r

− ^σ_r

λ^a +D

T

Q1

( λ

^r

) σ

^λr

r^λr+2 −Q1

( λ

^a

) σ

^λa

r^λa+2

+D² T²

Q2

( λ

r

) σ

^λr

r^λr+4 −Q2

( λ

a

) σ

^λa

r^λa+4

,

(10)

where the term that has prefactor D is the first-order, and the termwithprefactorD² isthesecond-orderFeynman–Hibbsquan- tumcorrection,and

Q₁

( λ )

⁼

λ ( λ

⁻¹

)

^{, (11)}

Q2

( λ )

=1

2

( λ

⁺²

)( λ

⁺¹

) λ ( λ

⁻¹

)

, (12) D = ¯h²

12mkB. (13)

ForMiepotentialswithFeynman–Hibbs quantumcorrections, the repulsiveregionincreasesinsizeatlowertemperatures(cf.Fig.1).

Thiscanbequantifiedbytheeffectiveparticlediameter,

σ

eff,given by the separation where the quantum-corrected potential equals zero.Theeffectivediameterwasfoundbysolving

u_MieFH

(

^r;T

)

^{=0 (14)}

for r using a Newton–Raphson solver with analytical derivatives andthevalueof

σ

^{asinitialguess.}

2.4.Thequantum-correctedcovolumeparameter

Wewillrelate thequantumcovolumebq(T) totheclassicalco- volumebby

bq

(

^T

)

=b

β (

^T

)

,

β (

^T

)

=

σ

eff

(

^T

) σ

eff

(

^Tc

)

3

. (15)

Eq.(15)isbasedontheassumptionthat thecovolumeparameter bofaclassicalcubicEoSiscorrectatthecriticalpoint,asitisde- terminedfromthemeasuredcriticaltemperature,Tc,andpressure, P_c.Thequantumcorrection,

β

^{(T)capturesthequantumswellingof}

theparticlevolumewithrespecttoitsvolumeatthecriticalpoint.

Since

β

^{(T) is given fromthe} quantum-corrected Miepotential, it shouldthus be independent ofthe choice ofcubic EoS.This will bestudiedinfurtherdetailinSection4.

Clearly

σ

eff(^T−1=0)=

σ

^{;inother words,theclassicalvalue is}

recoveredathigh temperatures.Differentiating

σ

eff asdefinedby Eq.(14)withrespectto T⁻¹ givesthatthefollowingslopeshould besatisfiedinthelimitT⁻¹→0:

∂σ

eff

∂

^T−1

T⁻¹→0

=D

(

^Q1

( λ

^r

)

−Q₁

( λ

^a

) )

σ ( λ

^r−

λ

^a

)

^{. (16)}

Eq.(16) holds forboth the first andthe second order quantum- correctedMiepotentials.In thelimit T⁻¹→∞ however,thefirst order quantum-corrected Mie potential reaches the maximum value

σ

eff^max=

σ

Q1

( λ

r

)

Q1

( λ

a

)

1/(λ^r−λ^a)

. (17)

Thefollowingexpressioncansatisfyalltheseconstraints:

σ

eff=

σ

+ D T+T_adj

(

^Q1

( λ

r

)

−Q1

( λ

a

) )

σ ( λ

^r−

λ

^a

)

^{, (18)}

where

T_adj=c_FH D

(

^Q1

( λ

r

)

−Q1

( λ

a

) )

σ ( σ

eff^max−1

) ( λ

^r−

λ

^a

)

^{. (19)}

The reason for the prefactor c_FH=1.4 for first order-corrections andc_FH=0.5 for thesecond ordercorrections, isto compensate forthe inabilityofthelinear extrapolationinT⁻¹ athighertem- peratures.Theaboveexpressionhastherightderivativeinthehigh temperature limit, but underpredicts the maximum effective di- ameterslightly. However, since the limit 1/T → ∞ can never be reachedinpracticeduetosolidformation,thisisoflittlepractical relevance.

Eq. (18) is plotted together with the increase in the effective diameterfrom the first andthe second order corrected Mie po- tentials with parameters for hydrogen, deuterium, neon and he- liuminFig.3,wherethe curvesarefromEq.(18)andthe circles arefromthequantum-correctedMiepotentials.Thefiguredisplays excellentagreementbetweenthederivedexpressionandtheexact valuefor

σ

eff/

σ

^{obtainedbysolvingEq.(14),inparticularatsuper-}

critialtemperatures.

Thederived

β

^{functioncanbewrittenas}

β (

^T

)

=

1+ A T+B

3

/

1+ A Tc+B

3

, (20)

whereAareBarecomponent-specificparameters.SincecubicEoS arenotbasedonMiefluids,andsincetheFeynman–Hibbscorrec- tions become inaccurate at very low temperatures, we will also evaluatethe empiricalcorrelation givenby Eq. (20)whereA and Barefittedtoexperimentaldata.

2.5.The

α

^{functionandquantumeffects}

Inadditiontothehypothesisthatquantumswellingcanbecap- turedby Eq.(15), we also hypothesize that the reduction of the well-depthgivenby thequantum corrections (see Fig. 1) can be modeledbyproperlymodifyingthe

α

^{functionofcubicEoS.Since}

thequantum-corrected covolumeislarger atsubcriticaltempera- tures,theattractionandthusthe

α

^{functionmustbereduced,e.g.}

torecoverthesamesaturationpressures. Thisisbecausethetwo

Table 1

Triple point temperature T t, critical temperature T c and critical pressure P c for the four components considered.

For helium T trepresents the superfluid transition.

H 2[35] ⁴He [36] Ne [37] D 2[38]

T t(K) 13.957 2.17 24.556 18.724 T c(K) 33.19 5.1953 44.492 38.34 P c(bar) 12.964 2.276 26.79 16.796

termsofcubicEoS (cf.Eq.(1)) areregressedtogether toapproxi- matetheexperimentaldata.Forsupercriticaltemperatures,theco- volumeissmallerthantheclassicalcovolume,andthusthequan- tum

α

^{function should be smaller than the classical one. These}

effectswill changethe “bell-shaped”

α

function found previously for hydrogen andhelium forthe classical PR. However, it is not guaranteed that the change will make the

α

^{function consistent}

intheclassicalsense[24],andnotevenmonotonicallydecreasing.

Wewillexplorethisfurtherbycomparingtoexperimentaldatain Section4.1.

Throughoutthiswork, thecubic

α

^{functioniscorrelatedusing}

thethree-parameterTwufunction[31]:

α (

^T

)

=T_r^N(M−1)exp[L

(

¹−T_r^MN

)

^]. (21) Eq.(21) is designedsuch that

α

(^Tc)=1 forall (L, M, N). Having three parameters,we assume that thisfunctional formisflexible enough tocapturethe quantumeffects onthe attractiveterm. In Section 4.1 we will investigate whether the optimal

α

^function

becomes a monotonic function when incorporating a quantum- correctedcovolume.

2.6. Mixingandcombiningrules

Foraweapplytheusualquadraticmixingrules:

a=

i

j

xixjai j, ai j=

aiiaj j

(

¹−ki j

)

, (22)

wherex_iandx_jarethemolefractionsofcomponentsiandj,and k_ij isabinaryinteractionparameter.

Themixturecovolumeiscalculatedas

b=

i

j

xixjbi j, bi j= bii+bj j

2

(

¹−li j

)

. (23)

Acommonchoiceistosetl_{i j}=0,whichconvertsthemixtureco- volume intoa linearcombinationofthesingle-component covol- umes.

Inaddition tothe

α

^and

β

^{functions, foreachcomponent we}

also fitted a (temperature-independent) Péneloux volume shift c [32]. Some properties such asenthalpies are affected by such a volume shift [33], and care should be taken to also shift these properties correctly. The impact of volume shifts on mixtures is morecomplicated[34],butonecanshowthatalinearmixingrule ofPénelouxparameters will notimpact thePxyphase envelopes.

Suchalinearmixingruleisadoptedinthiswork:

c=

i

x_ic_ii. (24)

3. Parameterestimationmethodology

In thiswork, we haveemployed thecommon methodof first fittingsingle-componentparameters,andinasubsequentstepfit- tingthebinaryinteractionparameters.Table1liststhetriplepoint temperature,criticaltemperature,andcriticalpressureforthecon- sideredcomponents.

Fig. 3. The relative increase in effective diameters for helium, hydrogen, neon and deuterium as computed exactly from the Mie-FH potentials (circles) and the analytic approximation Eq. (18) (curves).

Table 2

Weights in the objective function for the fitting pure fluid parameters. Square brackets indicate an interval within which the weight varied among the different components.

Saturation

Property p ^sat ρ^sat c ^sat_v c ^satp h ^sat Weight 1.0 0.5 [0,0.5] [0,0.5] 0.5 Single-phase

Property ρ^sup c ^supv c ^supp h ^sup w ^sup Weight 0.5 0.5 0.5 0.5 [0,0.1]

3.1. Objectivefunctions

3.1.1. Purefluids

The parameters (L, M, N) for the Twu

α

^{function (Eq. (21)),}

thePénelouxvolume shiftc,andoptionallythe parameters(A,B) in the covolume beta function (Eq. (20)), were fittedto pseudo- experimental data. These data were generated using the highly accurate referenceequations ofstate for helium [36], neon[37], hydrogen[35],anddeuterium[38],whichwere accessedthrough CoolProp[39].

The objective function isgiven by a weightedsumof relative deviations:

O^pure

(

^L,M,N;A,B;c

)

⁼

i∈{^states} ξ∈{^properties}

w_ξ

ξ

i^expt−

ξ

i^calc

ξ

i^expt

^{, (25)}

wherei indexesthe stateswherethedeviations arecalculated,

ξ

represents one of the ten thermodynamic properties considered, and w_ξ is the weight givento the property

ξ

^{. Table 2 lists the}

propertiesandweights usedintheobjectivefunction.We useda simplexalgorithm[40]fortheminimization.

Forthesaturationproperties,wechose20statesequispacedin temperature. The upper temperaturelimit was4.8 K forhelium, 30.0Kforhydrogen,34.5Kfordeuteriumand41.0Kforneon.The

lower temperature limits for neon and deuterium coincide with their triple point temperatures. Currently there are no technical applicationsofhydrogenattemperaturesbelowtheboiling point, 20.3K; to avoid sacrificing accuracy at higher temperatures, the lowertemperatureforhydrogenwaschosentobethesameasfor deuterium.The lower limit of helium waschosen as3 K, which isabove thesuperfluidtransitionat2.2Kbutstill wellbelowits boilingpointtemperatureof4.2K.

Thesingle-phaseproperties were evaluatedover a rectangular gridintemperature–pressurespace,withfivesupercriticaltemper- aturesand20pressuresfrom1barandup.Forneon,hydrogenand deuteriumthetemperatureswere50,100,150,200and300K,and themaximumpressurewas500bar. Forheliumthetemperatures were 20,30,40,100 and150 Kand themaximum pressurewas 300bar.

Residualpropertiesofpurefluids,suchastheresidualenthalpy h^R andtheresidual isochoricheat capacityc^R_v,are sometimesin- cludedinthe objectivefunctionwhen regressingEoSparameters.

We stronglyadvise against this practice:in our fits we haveen- countered meanaverage percentage errors (MAPEs) ofexceeding 200%forh^Rinthesupercriticaldomain,whiletheMAPEinh^supis lessthan2%.Thisisbecausethesepropertiesmaybecomezero,as showninFigs.5and6.

3.1.2. Binarymixtures

Using the optimalparameters forthe regression for pure flu- ids,binaryinteractionparameters (k_ij,l_ij) werefittedby minimiz- ingthetotalabsolutedeviationofmeasured(superscriptexpt)and calculated (superscript calc) molefractions, corresponding to the objectivefunction

O^binary

(

^ki j,l_{i j}

)

⁼

i

|

^xexpti −x^calc_i

|

⁺

i

|

^yexpti −y^calc_i

|

^{. (26)}

Here x_i and y_i are the mole fractions in the liquid and the va- por. The calculations were performed at the same temperature

Fig. 4. The optimal αand beta functions for four different cases of PR (defined in Section 4 ), all fitted to the same objective function. The lower temperature limits for each component are given by T tof Table 1 .

andpressureas themeasurements. Onlyregular VLE stateswere includedin the objective function, whereas liquid–liquid equilib- rium (LLE) and vapor-liquid–liquid equilibrium (VLLE) measure- mentswereexcluded.

4. Resultsanddiscussion

Toevaluateonaneutralbasistheinfluenceofthequantumcor- rectionsforthecovolumeparameter,weconsiderfivecases:

Classic: All calculationsareperformedwith thepresentstate- of-the-art forclassicalcubicEoSforquantumfluids,thetc- PREoS[6].

Classic-fit: We usethe classicalPR EoSbutfit theparameters ofthe

α

^{functionandthevolumeshift.Thiscaseallowsus}

toevaluate whethertheimprovementincomparisontothe tc-PREoScomesfromre-fittingagainstadifferentobjective function.

Case FH1:The quantumcorrection derived inSection 2based onfirst-orderFeynman–Hibbscorrections(FH1)withparam- etersfromRef.[18]areused.Newparametersareregressed forthe

α

^{functionandthevolumeshift.}

Case FH2:The quantumcorrection derived inSection 2based onthesecond-orderFeynman–Hibbscorrections(FH2)with

parametersfromRef.[18]areused.Newparametersarere- gressedforthe

α

^{functionandthevolumeshift.}

CaseEmpirical: We regress all parameters:(A, B) inEq. (20), the

α

^{function,andthevolumeshift.}

For the case Classic, we use the same values for the critical pressure andtemperature asRef. [6]. Forall other cases,we use thevaluesinTable1.

ThecasesClassicandClassic-fitarebenchmarksforgaugingthe improvement offered by quantum corrections. For Case FH1 and CaseFH2,thequantumcorrectionsforthecovolumeparameterare predicted. These cases have the samenumber of fittingparame- tersasclassicalcubicEoS.CaseEmpiricalallowsassessingwhether anythingcanbegainedbyalsofittingthequantumcorrectionsto thecovolume.

Table3liststheoptimalpure-componentparameters (L, M,N) for the Twu

α

^{function (Eq. (21)), the parameters (A, B) for the}

covolume correction (Eq. (20)), and the volume-shift parameter.

Wealsoobtainedtheoptimalbinaryinteractionparameters(k_ij,l_ij) presentedinTable4.

4.1. Single-componentsystems

Results for the pure components hydrogen, helium, neonand deuteriumaresummarizedforallthecaseslistedaboveinTable5.

Fig. 5. Thermodynamic properties of normal hydrogen (H 2): saturation diagrams (a and b) and supercritical isobars (c–f). The solid curves are calculated from the quantum PR EoS, and the dashed curves are calculated with the reference multiparameter EoS from Ref. [35] . The isobars are plotted between 20 and 300 K.

Evidently thequantum-corrected cubic EoSare vastlysuperior to theclassicalcubicEoS,evenforthecaseswithnonewfittingpa- rameters (CasesFH1 andFH2). The lowest value of theobjective function wasinall casesachievedforCase Empirical,sinceithas moreadjustableparameters.Interestingly,forhydrogen,deuterium andneon, thisdoesnot appreciablydecrease theMAPEsin com- parison toCaseFH1,andseemsto reflecttheparticularchoiceof objectivefunctioninsteadofconstituting animprovedrepresenta- tion of quantum fluids. For these fluids, we recommend the co- volume corrections basedon FH1 corrections, since thesehave a soundtheoretical basis.Theexception ishelium,whereCaseEm- piricalis clearlysuperiorandisthe recommendedmodel.Thisis expected, aswe showedinprevious work thatboth theFH1 and FH2 correctionswere unable to capture theproperties ofhelium below20K[18].

The temperature-dependence of the covolume corrections

β

and theregressed

α

^{functions are displayedin Fig. 4. The figure}

shows that the shape of the

α

^{function for the classic PR (Case}

Classic-fit) deviatesfromthe classicalmonotonically decayingbe- haviorforhelium andhydrogen.However,the

α

^{functionsforall}

cases that employ quantumcorrections forthe covolumeparam- eter recovera classical behavior, except forhelium atthe lowest temperatures.Using aquantum-correctedcovolumeforhelium,it isinfact possibletoachieve agoodrepresentationwithamono- tonically decreasing

α

^{function,asshownintherownamed“De-}

creasing

α

^{” in Table 5. Although this improves the accuracy of}

a few properties, most properties are less accurately predicted.

In particular,the MAPEinthe saturation pressureincreases from 0.67% to 1.18%. In a previous work [18], we showed that for he- lium,whichhasthestrongestquantumeffects,theFeynman–Hibbs corrections are only accurate down to about 20 K.We therefore hypothesize that the functional formof Eq.(20)must be further modifiedin orderto describe heliumto an even higheraccuracy atthelowesttemperatures.Suchamodificationhowever,fallsbe- yondthescopeofthepresentwork.

Theoptimalparameters(L,M,N)forthe

α

^{functionofhydrogen}

anddeuteriuminTable3arenotablydifferentfromtheothercom- ponents.Inparticular,thevalue ofLisone totwo ordersofmag- nitudelarger.Toexaminethismoreclosely,wehaveattemptedto refittheparameterswithfixedvaluesofL,whichrevealedthatthe objectivefunctioncanbereducedfurtherbyincreasingL.However, theachievablereductionoftheobjectivefunctionisinsignificantin comparisontotheobjectivefunctioncorresponding totheparam- etersinTable3.Wehypothesizethatthisstemsfromasuboptimal formoftheTwu

α

^{function,butitseemsunlikelythatamoreflex-}

ibleformwouldyieldasignificantlybetterfit.

Thethermodynamicpropertiesofdeuteriumareknownlessac- curately than for the other components considered in thiswork, andtheMAPEsinTable5shouldbe interpretedwithcaution.For example,theuncertaintiesinthe vaporpressures are 2% andthe uncertaintiesinthesaturatedliquiddensitiesare3%[38].Fordeu- terium,speedsofsoundwere not includedinthe objectivefunc- tionasthereislargescatterintheexperimentaldata[38].

We have also regressed parameters forthe SRK EoS, and the correspondingparametersandMAPEsareprovidedinthesupple- mentaryinformation.The main conclusion fromthisis that, also for SRK, the quantum corrections yield greatly improved agree- mentwithexperimentaldatawithnonewfittingparameters.This stronglyindicatesthatthehypothesisstatedinSection2.4isvalid;

namelythatthequantumcorrectionofthecovolumeparameteris independentofthefunctionalformofthecubicEoS.Wefindthat the performance of SRK is worse than forPR, especially for the liquid-phasedensities.

Figs. 5 and 6 illustrate the high accuracy of the quantum- correctedcubic EoSfor heliumandhydrogen, whichare thetwo componentsexhibitingthestrongestquantumeffects.

4.2.Thermodynamicconsistency

Temperature-dependent covolume corrections are usually avoided in the development of cubic EoS since they can lead

Fig. 6. Thermodynamic properties of helium ( ⁴He): saturation diagrams (a and b) and supercritical isobars (c–f). The solid curves are calculated from the quantum PR EoS, and the dashed curves are calculated with the reference multiparameter EoS from Ref. [36] . The isobars are plotted between 20 and 150 K.

Fig. 7. Contour plots of the molar isochoric heat capacity c v(Figs. (a)–(d)), and the molar isobaric heat capacity c p(Figs. (e)–(h)) in the temperature–density plane. The saturation curve is given by the black curve, and negative values are illustrated by gray regions. The lower temperature is the triple point temperature, except for ⁴He where it is 2.5 K. Isobars for 100 bar, 250 bar and 500 bar are also indicated.

to unphysical results. In particular, they lead to diverging and negative isochoric heat capacities at sufficiently high pressures [27].InFig.7,we showthatthisisalsothecaseforthequantum corrected cubic EoS. However, thisonly occurs atpressures that are well beyond the domain of industrial interest, and into the solid-formationregime forsomeofthecomponents.Wehavealso includedcontourplotsoftheisobaricheatcapacityc_p.Anegative cp indicates thermodynamic instability of the pure fluid [12,41],

and such unstable regions should not occur in the regions of known fluid stability; Fig. 7 shows that no spurious instabilities arepredicted.

Interestingly,Fig.7f andhillustrate spurious regionsof stabil- ity for hydrogenand deuterium. Since the isobaric heat capacity becomes positive within the spinodal region, it follows that the equation of state predicts the existence of a pseudostable phase within the spinodalregion foradiabatic andisobaric systems(PS

Fig. 8. Pxy diagrams for the binary mixtures of the fluids considered in this work. The crosses are experimental measurements. The full lines are VLE, the dashed lines are LLE, and the dotted lines indicate VLLE, all calculated with the optimal parameters given in Tables 3 and 4 . The colors signify different temperatures.

ensemble) [12,41]. Such a pseudostable phase isin fact alsopre- dictedbytheclassicalSRKandPRequationsofstate,andisthere- fore not an artefact of quantum corrections. Aursand et al have also found a region ofpositive Cp within the spinodal region for methane, even for the simpler van der Waals cubic equation of state[12].

Another consistency criterion is that pressure–density isotherms should not cross for pure components close to the criticalpoint.Althoughthisisnotarigorousconsistencycriterion, such a crossing of near-criticalisotherms hasnever been experi- mentallyobserved.Forthefluidsconsideredinthiswork,wehave indeedcheckedthatthiscriterionissatisfied.

4.3.Theperformanceformixtures

Fig. 8 compares calculations using the quantum-corrected PR with available experiments [42–53] for the six binary mixtures:

hydrogen-neon, deuterium-neon, hydrogen-deuterium, helium- neon, helium-deuterium and helium-hydrogen. The figure shows that quantum-correctedcubic EoSis inexcellent agreement with availableexperimentalresults,evenforpressuresupto200bar.

Forthehydrogen-heliummixture,wefoundthatitwasneces- sarytoincludeaninteractionparameter,l_ij,thataccountsfornon- idealmixingof thecovolumeparameter. Thiswasexpected, asa similarparameterisneededforthecross-interactionparameterof

Table 3

Optimal parameters ( L, M, N ) for the Twu αfunction, ( A, B ) for the beta function, and c for the Péneloux volume shift. The parameters for the recommended quantum-corrected PR is shaded gray.

Model L ( −) M ( −) N ( −) A (K) B (K) c (cm ³/mol) Normal hydrogen (H 2)

Classic-fit 2.8994 −0 . 61791 −0 . 42846 0 0 −4 . 1101 Classic 1.5147 −3 . 7959 −0 . 1377 0 0 −5 . 3386 FH1 156.21 −0 . 0062072 5.047 3.0696 12.682 −3 . 8139 FH2 347.52 −0 . 0027936 8.2946 5.8821 14.791 −2 . 9125 Empirical 158.54 −0 . 0061196 5.2105 3.477 15 −3 . 8140 Helium ( ⁴He)

Classic-fit −0 . 046019 1.2618 0.69755 0 0 −3 . 4875

Classic 0.0063 1.2175 1.0909 0 0 −4 . 8915

FH1 0.18976 1.3964 0.58143 1.8774 7.7564 −2 . 9291 FH2 1.1393 93.272 0.0044747 2.7979 5.2677 −3 . 9406 Empirical 0.48558 1.7173 0.30271 1.4912 3.2634 −3 . 1791 Neon (Ne)

Classic-fit 0.40805 0.98441 0.78674 0 0 −2 . 6039

Classic 0.1887 0.947 1.4698 0 0 −2 . 3573

FH1 0.40453 0.95861 0.8396 0.4673 2.4634 −2 . 4665 FH2 0.38356 0.94695 0.87127 0.4679 0.88094 −2 . 4556 Empirical 0.3981 0.96535 0.82696 0.22069 −0 . 65243 −2 . 5676 Normal deuterium (D 2)

Classic-fit 0.3089 1.0716 0.6551 0 0 −4 . 4250

Classic 0.1486 0.9968 1.0587 0 0 −4 . 5555

FH1 55.007 −0 . 016981 3.1621 1.6501 7.309 −3 . 8718 FH2 63.647 −0 . 014525 3.283 1.9086 3.4071 −3 . 6319 Empirical 52.586 −0 . 017779 3.2179 2.2117 12.768 −3 . 8717

Table 4

Optimal interaction parameters ( k ij, l ij) for the quantum PR EoS.

D 2 H 2 4He

H 2 (0,0) − −

4He (0.45,0) (0 . 17 , −0 . 16) − Ne (0.18,0) (0.18,0) (−0 . 17 , 0)

hydrogen-helium as described by quantum-corrected Mie poten- tials[17].Thesemiclassicalforcefieldforthehelium–hydrogenin- teractionusesanon-additivityfortheinteractiondiameterof1.05 [17],i.e.

σ

12=1.05(

σ

1+

σ

2)/2.Thistranslatesintoanon-additivity forthecovolumeof1.05³≈1.16,whichcorrespondstol_{i j}=−0.16. Usingl_{i j}=−0.16andrefittingk_ij yieldeda significantlyimproved fitfor thehelium–hydrogen mixture.We foundthat little can be gained by further tuning the l_ij parameter. For the helium–neon mixture,weverifiedthatincludingl_ijinthefittingprocessyieldsa negligibleimprovement.Thisisonceagaininagreementwiththe

Table 5

Mean average percentage error (MAPE) relative to the reference equations of state [35–38] for each com- ponent. The properties are saturated liquid density ρ^sat, saturation pressure p ^sat, saturated liquid isobaric heat capacity c ^satp, supercritical density ρ^sup, supercritical isobaric heat capacity c ^supp , and supercritical speed of sound w ^sup. Results for the model “Classic” were generated with parameters from Ref. [6] . The MAPEs for the recommended quantum-corrected PR is shaded gray.

Model p ^sat ρ^{sat c sat}v c ^sat_p h ^sat ρ^{sup c sup}v c ^supp h ^sup w ^sup Normal hydrogen (H 2)

Classic-fit 0.52 3.59 40.72 14.93 1.26 3.18 4.12 1.39 0.52 5.14 Classic 2.02 1.98 83.88 19.12 2.12 6.64 4.75 3.37 2.34 6.76 FH1 0.33 1.10 4.11 11.16 0.93 0.71 1.04 1.05 0.60 3.29 FH2 0.60 1.00 6.58 14.77 1.99 3.43 1.53 1.79 1.01 6.22 Empirical 0.36 1.11 3.94 10.90 0.90 0.78 0.97 1.07 0.61 3.38 Helium ( ⁴He)

Classic-fit 1.94 7.61 2.21 65.61 5.18 6.38 4.52 3.39 3.78 15.08 Classic 2.52 4.30 18.06 51.47 3.79 10.50 8.77 2.76 4.74 23.24 FH1 0.81 4.67 4.79 33.59 2.87 1.55 2.58 1.16 0.49 4.62 FH2 0.29 4.88 5.15 10.46 2.84 1.12 1.64 1.50 0.39 1.46 Empirical 0.67 1.70 2.17 12.26 1.76 0.45 1.64 0.74 0.48 2.57 Decreasing α 1.18 2.18 1.74 12.19 1.45 0.51 1.77 0.70 0.54 2.59 Neon (Ne)

Classic-fit 0.57 1.51 1.69 8.66 0.96 0.38 1.66 0.53 0.24 1.53 Classic 0.86 1.78 13.18 6.13 0.82 1.07 5.66 2.22 0.98 1.74 FH1 0.25 1.18 1.99 8.16 0.59 0.57 2.25 0.65 0.23 2.01 FH2 0.26 1.17 2.17 8.46 0.61 0.59 2.30 0.67 0.24 2.04 Empirical 0.29 1.29 1.74 7.83 0.65 0.40 1.92 0.56 0.21 1.71 Normal deuterium (D 2)

Classic-fit 2.04 1.89 6.27 20.14 2.31 2.07 3.09 0.62 0.67 15.79 Classic 1.89 1.43 7.67 18.64 1.88 2.82 3.65 0.88 0.78 17.10 FH1 0.61 0.83 6.55 14.23 0.90 0.60 0.90 0.84 0.44 10.47 FH2 0.78 1.17 11.45 20.52 1.22 1.07 1.01 0.99 0.35 9.07 Empirical 0.60 0.81 5.66 12.80 0.85 0.70 0.82 0.86 0.44 10.05