Dynamics of monolayer physisorption in homogeneous mesoporous media

1

Dynamics of Monolayer Physisorption in Homogeneous

2

Mesoporous Media

3

Paul Papatzacos*

4Department of Mathematics and Physics, University of Stavanger, 4036 Stavanger, Norway

5 ABSTRACT: A model for monolayer physisorption of a one-

6 component gas on the pore surface of a homogeneous macroporous

7 or mesoporous porous medium is presented. It originates from an

8 averaging over many pores of a macroporous mediumfilled with a one-

9 componentfluid. The resulting model does not assume anything about

10 pore shape, but assumes that the pores are so large that capillary

11 condensation does not occur. Mathematically, the model gives coverage

12 as the solution of an ordinary,first-order, differential equation, where the

13 time derivative of coverage is proportional to the difference between the chemical potential of the adsorbate and the chemical

14 potential of the ambient gas. Coverage is determined by the ambient gas density, with temperature, adsorbate critical

15 temperature, and the Henry adsorption constant as parameters. The rest of this abstract describes what is deduced from the

16 equations of the model. Adsorbate phase transitions are built into the model by the use of van der Waals equations of state.

17 Equilibrium isotherms are derived from the equality of the chemical potentials. The differential equation for coverage makes it

18 possible to determine the mathematical stability of the equilibrium isotherms, and a number of properties of the isotherms are

19 derived, the most important being as follows: (i) an adsorbate phase transition is always accompanied by a well-defined

20 hysteresis loop, although“loop”is somewhat misleading as its vertical boundaries do not consist of equilibrium states; (ii) the

21 vertical boundaries are exactly located; (iii) the upper and lower boundaries consist of states that are mathematically stable,

22 while being either physically stable or metastable, and if physical metastability is the case, then the actual state of the adsorbate

23 (mono- or bi-phasic) will not be visible on the equilibrium isotherm. The shapes of the equilibrium isotherms are largely

24 determined by the value of the Henry constant, whether the isotherms are subcritical or supercritical. Expressions for the

25 location of an equilibrium isotherm’s region of fastest variation and for the locations of the vertical boundaries of its hysteresis

26 loop are found that also show the importance of Henry’s constant. Dynamical, that is, time-dependent isotherms are presented

27 for the case describing the variation of coverage resulting from forcing the ambient gas to undergo a compression−

28 decompression loop. Two subcases are considered: the subcritical and the supercritical adsorbate. It is shown that coverage in

29 terms of ambient pressure exhibits closed loops, even in supercritical isotherms. However, supercritical loops shrink when the

30 cycle time increases, reminiscent of rate-dependent hysteresis observed in piezoelectricity. The model is used to interpret two

31 experiments on the sorption of CO₂and CH₄on coal that showed hysteresis loops in isotherms of supercritical adsorbates and

32 that were originally interpreted as leading to different Henry constants for adsorption and for desorption. The interpretation set

33 forth here uses the inherent dynamics of the model and looks at the loop as just one isotherm evolving in time, thus leading to a

34 unique Henry constant.

1. INTRODUCTION

35In experiments on gas physisorption, one often observes a

36discontinuity in the equilibrium isotherms and a hysteresis

37loop. See Morishige and Shikimi¹ and references given there.

38The step and the loop occur at temperatures well below the

39critical temperature of the ambient gas, and at pressures well

40below its saturation pressure.

41 It has been shown by Hill, in an article published in 1947,²

42that hysteresis can be explained by assuming the existence of

43metastable adsorbed states in monolayer physisorption, no

44assumptions about the pores being necessary. Hill’s result is

45generalized in the present article, where monolayer phys-

46isorption is used to the exclusion of other processes. It must be

47mentioned that monolayer physisorption in a mesoporous or

48macroporous medium is, in a certain sense, in a class by itself,

49possibly together with multilayer physisorption if the number

50of layers is on the order of two or three. It has indeed been

shown³ that the size of the pore surface per unit volume of a 51

mesoporous or macroporous medium is such that the amount 52

adsorbed by monolayer physisorption is negligibly small when 53

compared to the amount thatflows in the pores. On the other 54

hand, physisorption by capillary condensation and/or chem- 55

isorption deal with adsorbed amounts that differ by orders of56

magnitude from those occurring in monolayer physisorption 57

and are essential to describe such processes as industrial 58

hydrocarbon recovery. Capillary condensation and chemisorp- 59

tion are not considered in any detail in this article. 60

The generalization of Hill’s result is done in the framework 61

of a sorption model, called M′ for convenience here. M′, a62

special case of a model M to be described presently, expresses 63

Received: September 11, 2019 Accepted: November 27, 2019

Article http://pubs.acs.org/journal/acsodf

© XXXX American Chemical Society A DOI:10.1021/acsomega.9b02956

ACS OmegaXXXX, XXX, XXX−XXX

64therate of change of coveragein terms of the coverage itself, of

65the ambient gas density, of temperature, of the critical

66temperatures of adsorbate and ambient gas, and of Henry’s

67adsorption constant. This is not the first appearance of an

68equation for the rate of change of coverage (see Alfé and

69Gillan⁴), but it is, to the author’s knowledge, thefirst time such

70an equation leads to an understanding of hysteresis loops in

71adsorption, and to a reinterpretation of experimental results.

72The three paragraphs below are short presentations of results

73that are described in more detail in Sections 2, 3, and4.

74 The first result concerns the placement of the vertical

75boundaries of the hysteresis loop (Section 2.5). Ball and

76Evans,⁵in their article on the mechanism for hysteresis, noted

77that the existence of physically metastable states will bring

78about a transition to a physically stable state at some ambient

79pressure and thus produce a vertical boundary for the

80hysteresis loop at that pressure in adsorption as well as in

81desorption. They also remarked that determining the transition

82pressure is beyond the scope of an equilibrium theory, given

83that there are infinitely many physically metastable states. Now

84M′describes the evolution of isotherms with time, and it also

85determines the equilibrium isotherms. This implies that the

86mathematical stability of any point on an equilibrium isotherm

87can be found, and it turns out that physically metastable states

88are mathematically stable, except for just two such points, one

89for adsorption and one for desorption: these are mathemati-

90cally unstable and determine the transition pressures.

91 The second result concerns the new possibility implied by

92the ability of M′ to describe time-dependent isotherms. This

93has a direct relevance to measurements where adsorption and

94desorption follow different paths that join at a low and a high

95coverage, thus exhibiting a loop^7,6with no vertical boundaries.

96The explanation given by M′ is the one given by other

97workers: a genuine hysteresis loop must have two vertical

98boundaries, so their absence is explained by appealing to an

99insufficient equilibration time⁶ or waiting time.⁸ There is,

100however, a new possibility implied by the ability of M′ to

101describe time-dependent isotherms: that of considering the

102noncoinciding adsorption and desorption paths as being just

103one isotherm evolving in time under the application of a

104pressure-cycle consisting of compression followed by decom-

105pression of the ambient gas. The isotherms resulting from such

106a cycle are shown in Section 3, for the two cases of a

107supercritical and a subcritical isotherm.

108 The third result concerns the interpretation of experiments

109on sorption of methane (important as a source of energy) and

110on sorption of carbon dioxide (an important product to

111sequestrate). These two cases of sorption are exceptional in

112that they cannot be described in the framework of capillary

113condensation: the critical temperatures of the substances are

114low compared to storage temperatures, so that capillary

115condensation cannot occur. Wang et al.⁹ enumerate, and give

116references for, the hypotheses that have been made to explain

117the mechanism of methane and carbon dioxide sorption

118hysteresis: residual moisture in coal samples, surface geometry

119heterogeneity, chemical interaction, structural deformation,

120experimental inaccuracies, and insufficient waiting time. They

121conclude that the mechanism remains an open question. The

122most straightforward way to describe CH₄and CO₂ sorption

123has been to use the Langmuir model: see Jessen et al.⁷See also

124Wang et al.⁹ who look at two additional isotherms, one from

125the Dubinin−Radushkevich model and one from the dual

126sorption model, the latter allowing the inclusion of the effect of

coal swelling. Section 4 of the present article gives an 127

alternative description, based on the second result above, 128

that leads to a unique value for the Henry constant instead of 129

the two obtained byfitting separate equilibrium isotherms, one 130

for adsorption and another for desorption.⁷ 131

It is also worth mentioning that the mathematical expression132

of M′is simple enough to allow approximate expressions for a 133

number of useful quantities, such as the pressure at which the 134

isotherm is steepest and the width of the hysteresis loop. 135

A short description of model M′ and of the underlying 136

model M, follows. 137

M is a model for multiphaseflow in a porous medium, based 138

on the diffuse interface assumption.¹⁰ It is the result of an 139

averaging over many pores of the equations describing 140

Navier−Stokes flow in the pores. The averaging leads to a 141

new set of equations involving averaged quantities such as 142

density, velocity components, temperature, internal energy, 143

and entropy. M, and thereby M′, are based on the following 144

assumptions: (a1) the fluid-containing pores are connected; 145

(a2) the smallest pore-throat diameter is large when compared 146

to the average distance betweenfluid molecules, and also when 147

compared to their mean free path; (a3) adsorption occurs by 148

physisorption; (a4) adsorption is monolayer; (a5) the heat 149

generally released by adsorption does not appreciably change 150

the temperature; (a6) the averaged fluid quantities obey the151

same thermodynamical laws as the quantities of the original 152

fluid and, in particular, the averagedfluid has a well-defined 153

pressure obeying an equation of state that can be chosen 154

among the known ones. 155

M′contains three additional assumptions: (a7) the averaged 156

adsorbed fluid is assumed to have the thermodynamics of a 157

two-dimensionalfluid with, in particular, a spreading pressure 158

(the negative of the surface tension)¹¹obeying a van der Waals 159

equation of state; (a8) the ambientfluid is monophasic; (a9) 160

any externally applied changes, such as compressions or 161

decompressions, are done so slowly that the ambient fluid 162

velocity is negligibly small. 163

The consequences of these assumptions are discussed164

presently, after the introduction of the basic equations of M′. 165

These are obtained directly from M (see eq 17 in 166

Papatzacos¹⁰), with a slight modification in notation 167

c_Σ̇ = −Δμ (1) ₁₆₈

L( _f)

μ μ μ

Δ = _Σ− _{(2) 169}

The dot in thefirst equation denotes partial differentiation 170

with respect to time. A space dependence can also be included 171

for c_Σ, but is ignored here, as it is shown below that the 172

assumptions of model M′make it redundant. At the level of 173

model M, μf is the chemical potential of the ambient fluid, 174

modified by the addition of two terms: a term proportional to 175

the Laplacian of the ambient fluid density and a term176

proportional to the squared modulus of the ambient fluid 177

velocity. The Laplacian originates in the diffuse interface 178

framework, where large gradients of density exist in the 179

interfaces between phases, if two phases coexist. The squared 180

velocity accounts for the kinetic energy exchanged between 181

ambient and adsorbedfluids. The two equations above are the 182

core of model M′. They lead, after some preliminaries 183

presented as background material inSection 5, to a differential 184

equation for the coverage, presented inSection 5.4. 185

The consequences of assumptions (a1) to (a9) are as 186

follows. 187

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX B

188 As a result of the averaging process, the individual

189characteristics of the pores are lost, leaving only two

190parameters to characterize the medium as a whole, which are

191porosity and pore surface per unit volume.¹⁰

192 Assumptions (a1), (a2), and (a4) imply that adsorption

193does not affect the basic description of the ambientfluid by the

194equations of fluid mechanics expressing balance of mass,

195momentum, energy, and entropy. They also imply that

196adsorption induces negligible changes in the values of porosity

197and pore surface per unit volume. In fact, assumption (a4)

198implies, as stated above, that inside an arbitrary volume of

199porous medium, the total mass adsorbed is negligibly small

200when compared to the mass offluid that canflow freely in the

201pores.³ Assumption (a5) implies that the quantities character-

202izing the ambientfluid, such as its density and temperature and

203consequently its pressure and chemical potential, are not

204modified by sorption. (On the other hand, the quantities

205characterizing the adsorbed fluid are determined by the

206ambientfluid.) Assumption (a7) implies that phase transitions

207can occur in the adsorbate. Assumption (a8) implies that no

208interfaces exist in the ambientfluid, so that the Laplacian of the

209ambient density that modifies the chemical potential is

210negligible. Assumption (a9) implies that the other term

211modifying the chemical potential of the ambient fluid is

212negligibly small. With these last two assumptions, μf is the

213usual chemical potential of the ambient fluid. Another

214consequence of (a8) and (a9) is that eqs 1and 2 are space

215independent.

216 It is here emphasized that, as already mentioned, isotherm

217properties deduced in M′, like adsorbate phase transitions or

218hysteresis loops, do not depend on assumptions about the

219shape of the pores and, in particular, occur without the

220occurrence of capillary condensation.

221 The thermodynamical description of the fluids is given as

222background material in Section 5, where explicit expressions

223for the pressure, spreading pressure, and chemical potentials

224are given, and where Henry’s adsorption constant is

225introduced. As already mentioned, the adsorbate is a van der

226Waals fluid, whereas three alternatives are considered for the

227ambient fluid: ideal with zero volume particles (called ideal

228type 0), ideal with nonzero volume particles (called ideal type

2291), and van der Waals. TheΔμof eqs 1and2above is then

230derived as a function of coverage, ambient density, and some

231parameters that include temperature and Henry’s constant.

232The concepts of physical stability, metastability, and instability,

233known to occur in connection with the van der Waals equation

234of state, are illustrated infigures in the same section, for easy

235reference in the rest of the article.

236 The derivation of the basic differential equation defining M′

237is given as background material in Section 5.4. Other basic

238properties of M′, specifically the definitions of equilibrium

239isotherms and of mathematical stability are given inSection 2.

2. RESULT AND DISCUSSION 1: THEORY OF

240 EQUILIBRIUM ISOTHERMS

241Equilibrium isotherms are defined inSection 2.1. Mathematical

242stability is defined inSection 2.2, and is followed by a section

243on equilibrium isotherms given by analytical expressions, then

244by two sections on equilibrium isotherms given by numerical

245solutions.

246 2.1. Definition of an Equilibrium Isotherm. It is here

247referred to Section 5.4, where the differential equation giving

248the time rate of change of coverage is deduced, that is, eq 59.

The equilibrium isotherms follow from this equation: they are 249

the solutions that have zero rate of change. There are 250

additional requirements concerning stability and unicity: see 251

the last paragraph but one of this section. 252

An equilibrium solution is defined as follows. GivenT̃,T̃_Σc, 253

andψ, an equilibrium solution of eq 59is a set of points (r_e, 254

θe) in a Cartesian plane, satisfying 0 <r_e≤r_g, 0 <θe< 1, and 255

r T T

( ,_{e e}, , _c, ) 0

μ θ ψ

Δ ̃ ̃ _Σ̃ = _{(3) 256}

An equivalent form ofeq 3is found as follows. The equation 257

isfirst rewritten asμ̃Σ,red−T̃lnK̃_H=μ̃f,red, by usingeqs 56and 258

57. Dividing both sides byT̃, exponentiating, and referring to259

eq 45, onefinds 260

f̃ ( ,θ T̃, )τ =K f r T_{H f}̃ ̃( , ̃)

Σ (4) 261

This equation, with another definition for the proportion- 262

ality constant, is identical witheq 4in the article by Hoory and 263

Prausnitz.¹² 264

Equations 3and/or 4 give equilibrium isotherms provided 265

that (i) points that represent mathematically and/or physically 266

unstable states are discarded, and that (ii) non-unicity ofθefor 267

anyr_e agrees with observations of hysteresis. 268

Physical stability is considered inSection 5.1; mathematical 269

stability is defined inSection 2.2. 270

2.2. Mathematical Stability of Equilibrium Isotherms. 271

Mathematical stability, orm-stability, of equilibrium isotherms,272

is defined as follows.¹³ 273

For any given set {r_e,T̃,T̃_Σc,ψ}, a solutionθeofeq 3or4is274

said to be (asymptotically) m-stable if there is a neighborhood 275

ofθesuch that anyθinside this neighborhood approachesθeas 276

time increases. 277

It follows from the definition that a very simple criterion for 278

deciding whether any point on an equilibrium isotherm is m- 279

stable is as follows: 280

θe is an m-stable equilibrium solution if and only if Δμ̃ 281

changes from negative to positive values whenθ−θedoes so. 282 283 f1

Indeed, referring to Figure 1, and keeping eq 59in mind, one sees that, if the condition is satisfied, then any perturbation284

ofθethat takes place during time dt̃> 0, and that bringsθin 285

the interval (θe,θe+a), gives rise toΔμ̃> 0. It follows that dθ 286

= −Δμ̃ dt̃ < 0, so that θ is drawn back to θe. A similar 287

reasoning can be made for a perturbation that bringsθto the288

left ofθe, with the same conclusion thatθis drawn back toθe. 289

Necessity is also easily proven. 290

2.3. Analytical Equilibrium Isotherms and Their 291

Stability. Analytical solutions of eq 4 are considered here. 292

Figure 1.Figure used in definingm-stability inSection 2.2. IfΔμ̃vsθ behaves as shown in the vicinity ofθe, thenθeis m-stable.

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX C

293More precisely, one seeks to expressθeas a known function of

294r_e or vice versa. Obviously, this is only possible if one knows

295how to invert either f̃_Σ orf̃_f. No inverse off̃_Σ, given byeq 46

296(right), is known. Inverting f_f̃, however, may be possible as

297shown below. There are three acceptable approximations for

298the ambientfluid fugacity, and two of them are invertible: the

299ones obtained by assuming that the ambient fluid is ideal,

300either of type zero (zero volume particles) or of type 1

301(nonzero volume particles). The relevant expressions forf̃_fare

302given in eq 47(right) and 48(right). In fact, most analytical

303isotherms are expressions giving pressure as a function of

304coverage, see for example, Table I-1 in the book by Ross and

305Olivier.¹⁴This can easily be done in the present case by writing

306each invertible fugacity in terms of the corresponding pressure.

307Referring to the left and center equations of the set (36), it is

308easily seen that

i

kjjjjj y

{zzzzz

f P f P P

, exp T

f

id0 id0

f

id1 id1 id1

̃ = ̃ ̃ = ̃ ̃ ̃

309 Assuming that the ambient gas can be approximated by an

310ideal gas of type 0,eq 4gives the known expression i

kjjjjj y

{zzzzz

P f

K T

K 1 exp 1

id0 27

H H

e e

e e

θ e

θ

θ θ

θ τ

̃ =

̃̃ = ̃̃ − − −

Σ

311 (5)

312 Assuming that the ambient gas can be approximated by an

313ideal gas of type 1,eq 4gives i

kjjjjj y

{zzzzz f

K T P

T P exp T

H

id1 id1

̃̃ ̃ = ̃ ̃ ̃ ̃

Σ

314 This equation is of the type Z exp Z = g and, g being

315positive, is equivalent toZ= Wp(g), where Wp is the principal

316part of the LambertW-function.¹⁵One thus obtains i

kjjjjj j

y {zzzzz z i

kjjjjj i

kjjjjj y

{zzzzzy {zzzzz

P T f

K T

T K

Wp ,

Wp 1

1 exp

1

27

id1

H

H e

e

e e

θ e

θ

θ θ

θ τ

̃ = ̃ ̃̃ ̃

= ̃ ̃ − − −

Σ

317 (6)

318 To the author’s knowledge, this equation has not previously

319appeared in the adsorption literature. A plot ofθeversusP̃^{id 0}is

320usually quite close to a plot of θe versus P̃^{id 1}, except for ψ-

321values less than 1 and forθe-values close to 1.

One now turns to the restrictions mentioned inSection 2.1, 322

so as to establish whether the expressions above are acceptable. 323

Plots of eq 6, say, show continuous curves, monotonic324 325 f2

increasing whenT̃ >T̃_Σc(Figure 2, left-hand plot), but with a region of three-valuedness when T̃ < T̃_Σc (Figure 2, right326

plot).^a 327

The m-stability of the supercritical isotherm (left plot of328

Figure 2) is established as follows. Lettingθmove on a vertical 329

line, say upward from a low value, thenθ−θeandΔμ̃both go 330

from negative to positive values when the isotherm is crossed. 331

The necessary and sufficient condition stated inSection 2.2 is332

satisfied, so that the supercritical isotherm is m-stable. It is also 333

physically stable (p-stable) because ambient gas and the 334

adsorbate behave nearly ideally. One thus recovers the known 335

result that the supercritical isotherms given by eq 5 are 336

physically correct. 337

Turning to the subcritical case (right plot ofFigure 2), one 338

sees that the isotherm has three parts: a central part, whose 339

states are p-unstable (orange line), connecting a lower to an 340

upper part. The lower is here called the adsorption branch, and 341

the upper is called the desorption branch. Each of these 342

branches contains a set of p-stable states (black line) and a set 343

of p-metastable states (green line). The m-stabilities of the 344

three parts are established by the same argument as used in the 345

supercritical case. One easilyfinds that the central part (orange346

line) consists of m-unstable states, so that there is no doubt in 347

discarding these points that already are p-unstable. However, 348

all points on the adsorption and desorption branches are, with 349

the exception of one point for each branch, m-stable, regardless 350

of the quality of their physical stability: the whole set of p- 351

stable points and almost the whole set of p-metastable points 352

are m-stable. The exception is, for each branch, the p- 353

metastable point having a vertical tangent, that is, the point on 354

the isotherm where θ = θm (desorption branch) or θ = θM355

(adsorption branch): indeed, repeating the m-stability argu- 356

ment, and letting θmove on a vertical tangent, thenΔμ̃does 357

not change sign when _m 358

M

θ−θ does. The m-instability of these

points is important in interpreting the region of two-valuedness 359

as a hysteresis loop with the following boundaries: its upper 360

and lower boundaries coincide with the desorption and 361

adsorption branches over a pressure interval [P̃₁, P̃₂]. P̃₂ is362

the abscissa of the point on the adsorption branch whose 363

ordinate is the left spinodal coverage, θM, whereas P̃₁ is the 364

abscissa of the point on the desorption branch whose ordinate 365

is the right spinodal coverage,θm. Note that the left and right 366

boundaries of the hysteresis loop can only be the vertical lines 367

Figure 2.Plots ofeq 6, for a supercritical (left) or subcritical (right) adsorbate for values ofT̃,T̃_Σc, andψ, as indicated. The right plot uses the notations ofFigure 14for the subscriptedθ′s, as well as for the meanings of the colors. Concerning the vertical dashed lines in the right plot, see Section 2.3.

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX D

368atP̃=P̃₁andP̃=P̃₂because (P̃₁,θm) and (P̃₂,θM) are the only

369m-unstable points on the desorption and the adsorption

370branches: the transition to a mathematically and physically

371stable state can only take place from (P̃₂, θM) during

372adsorption, and from (P̃₁, θm) during desorption. See Figure

3732, right plot. Note also that the vertical boundaries of the

374hysteresis loop are not places of equilibrium, as these are in

375regions where Δμ̃ ≠ 0. They are drawn as dashed lines in

376Figure 2to emphasize this fact.

377 A further remark on the p-metastable states follows: by

378definition, a p-metastable adsorbate state will, if perturbed,

379transit to a p-stable state of lower energy. Referring to Figure

38014(right), such a transition brings the adsorbate from a point

381on one of the green lines to the reconstructed dash-dotted line

382at the same value of coverage. This implies that it is not

383possible to say whether an equilibrium state represented by a

384point situated on one of the green lines ofFigure 2(right plot)

385is in a two-phase or in a one-phase adsorbate state: the two

386states have the same coverage and also the same ambient

387pressure. The equality of pressures is approximative and is due

388to the assumptions that characterize M′, implying that changes

389in the adsorbate are not “visible” in the ambient gas. As a

390consequence, p-stable and p-metastable adsorbate states are

391treated in the same way inSection 2.5below.

392 Thus, for equilibrium isotherms that can be obtained from

393analytical expressions, model M′ defines a subcritical

394equilibrium isotherm plotted against pressure as follows: it

395has an adsorption branch that spans all pressures up to a

396pressure P̃₂ at which the adsorbate is at its left spinodal

397coverage,θM. It has a desorption branch that spans pressures

398down to a pressure P̃₁ at which the adsorbate is at its right

399spinodal coverage θm. P̃₁ < P̃₂, hence, the interval (P̃₁, P̃₂)

400defines the pressure range of the hysteresis loop. The vertical

401boundaries of the loop at pressuresP̃₁andP̃₂do not consist of

402equilibrium points.

403 Sections 2.4 and2.5 look at supercritical and at subcritical

404isotherms in cases where no analytical solution is available, that

405is, whenμ̃f,redis given byeq 49(left).

406 2.4. Numerical Equilibrium Isotherms for Super-

407critical Adsorbates and Their Stability. Model M′ is

408now considered for the general case of numerically solvingeq 3

409forT̃ ≥T̃_Σc.

410 According toeq 56, theΔμ̃ versus θ curve, at givenr_e, T̃,

411T̃_Σc, andψ, is theμ̃Σ,redversusθcurve with a vertical translation

412induced by μ̃f,red(r_e,T̃) and byψ(T̃). Theμ̃Σ,redversusθcurve

413has the shape of the monotonically increasing curve shown in

414Figure 15 (right). One can immediately conclude that there

415will always be a solution (as the curve spans the vertical axis

416from −∞to +∞), and that it is unique. Using mathematical

417stability, as defined inSection 2.2, together with the plots of

418Δμ̃ versusθ, one concludes thatθe(r_e) is m-stable in addition

419to being p-stable. For later reference, this equilibrium solution

420is written as r T T

T T

( , , , ) unique solution of eq 3,

( )

e e c

c

θ ̃ ̃ ψ =

̃ ≥ ̃

Σ

421 Σ (7)

422 Its general shape is given inFigure 2(left).

423 The interpretation ofeq 3as the intersection with theθ-axis

424of a vertically translatedμ̃Σ,redversusθ curve has some useful

425consequences concerning the general shape of an equilibrium

426supercritical isotherm, especially the location of its point of

427inflection.

Knowing the set {T̃,T̃_Σc,ψ}, it is possible to roughly predict 428

the position of the region where the isotherm is steepest, 429

assuming that the parameters in the set above are such that the 430

isotherm flattens out. This follows from the observation that431

the μ̃Σ,red versus θ curve can be translated upward by an432

arbitrary amount by lettingr_ego arbitrarily close to 0. Letting433

r_e increase from such a value, the μ̃Σ,red versus θ curve is 434

translated downward, and its intersection with theθaxis gives 435

values ofθethat increase, the fastest increase taking place when436

the inflection point of theμ̃Σ,redversusθcurve is close to theθ- 437

axis. This can only happen if ψ (the other quantity that is 438

subtracted fromμ̃Σ,redin the expression ofΔμ̃) is large enough. 439

It is then useful to define a functionΨu(T̃,T̃_Σc) as follows: ifψ440

= Ψu, then the inflection point of the Δμ̃ versus θ curve, at 441

givenT̃,T̃_Σc, and atr=r_g(T̃), is on theθ-axis. The inflection 442

point, easily found by equating to zero the second derivative of 443

μ̃Σ,redwith respect toθ, occurs atθ= 1/3, independently ofT̃ 444

and of T̃_Σc, so that 445

T T T T r T T

( , ) (1/3, , ) ( ( ), )

u c μ ,red c μf,red g

Ψ ̃ _Σ̃ = ̃ ̃ ̃ − ̃ ̃ ̃

Σ Σ

(8) 446

SeeFigure 4 for the general behavior ofΨu. Note that Ψu 447

does not depend on the values in the set+. 448

449 f3

Figure 3, shows how the shape of the equilibrium isotherm changes for values ofψin the neighborhood ofΨu. Note that 450

the tendency towardflattening occurs whenψ>Ψu. 451

On an equilibrium isotherm plotted asθeversusr_e, one can 452

find the approximate position of its point of inflection. It is the 453

value of r_e, here denoted r_i, that corresponds to θe = 1/3. 454

According to eq 4, it is given by 455

f r T_f̃( ,_i ̃ =) e^−ψ/T̃f̃ (1/3,T̃, )τ

Σ (9) 456

One can obtain a good approximation for thisr_i if one can 457

assume that it is sufficiently close to 0 thatf̃_f≈T̃r_i(seeeq 47 458

(right)). Then, usingeq 46(right), one gets 459

r K

e

2 e e

i 1/2 9/4 T 2

/ 1/2 9/4

H

= =

̃

τ ψ

− τ

− ̃ −

(10) 460

InFigure 3, the values ofT̃r_i/P̃₀given by this expression for461

ψ= 1.5Ψuand 2.0Ψuare indicated by the vertical segments. 462

Figure 3.Equilibrium supercritical isotherms,θevsP̃/P̃₀, whereT̃_Σc= 0.5 and T̃ = 0.6. Curves originate from eq 7. The value of ψ is indicated for each curve (Ψu= 1.2188). The horizontal dashed line is atθe= 1/3. Concerning the vertical segments, see the end ofSection 2.4.

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX E

463 2.5. Numerical Equilibrium Isotherms for Subcritical

464Adsorbates and Their Stability. Returning to the

465interpretation of eq 3 as the intersection with the θ-axis of a

466vertically translatedμ̃Σ,redversusθcurve, one now assumesT̃<

467T̃_Σc, so that theμ̃Σ,redversusθcurve has a local maximum and a

468local minimum as shown inFigure 15(right). There can now

469be one, two, or three values ofθe, depending on the position of

470the Δμ̃ versusθcurve in relation to theθ-axis. This position

471depends, in turn, on the values ofμ̃f,redandψ.

472 The use of mathematical stability in this case leads to the

473following. If there is just one value ofθe, then it occurs as the

474intersection with theθ-axis of that part of the translatedμ̃Σ,red 475versusθcurve that is either to the left of the local maximum or

476to the right of the local minimum and is therefore m-stable. If

477there are two distinct values of θe, then one of them is m-

478stable, the other (being the abscissa of the local maximum or

479minimum) is m-unstable and discarded. If there are three

480values of θe, then: (i) the smallest and largest are m-stable,

481whereas the middle one is m-unstable and discarded; (ii) the

482smallest is in the interval (0, θM), and the largest is in the

483interval (θm, 1).

484 For later reference, the m-stable equilibrium solutions are

485written r T T T T

( , , , ) solution of eq 3 in (0, ),

( )

ea e c M

c

θ ̃ ̃ ψ = θ

̃ < ̃

Σ

486 Σ (11)

r T T T T

( , , , ) solution of eq 3 in ( , 1),

( )

ed e c m

c

θ ̃ ̃ ψ = θ

̃ < ̃

Σ

487 Σ (12)

488 Because θM < θm, there are two separate branches as in

489Section 2.3, an adsorption branch, θea, and a desorption

490branch,θed, creating a double-valuedness for certain pressures.

491As in the right plot ofFigure 2, this is interpreted to mean that

492there is a hysteresis loop and a phase transition for the

493adsorbate.

494 The sizes ofψand ofr_gbeing important, one is led to define

495the following two functions of T̃ andT̃_Σc

T T T T r T

( , ) ( , , ) ( , )

m c μ ,red θm c μf,red g

Ψ ̃ _Σ̃ = ̃ ̃ ̃ − ̃ ̃

Σ Σ

496 (13)

T T T T r T

( , ) ( , , ) ( , )

M c μ ,red θM c μf,red g

Ψ ̃ _Σ̃ = ̃ ̃ ̃ − ̃ ̃

Σ Σ

497 (14)

498whereθmandθMare as defined inFigures 14and15, and byeq

49937. It is easily seen thatΨm<ΨM. Both functions depend onT̃

f4 500andT̃_Σcand not on the set+.Figure 4shows the behavior of

501Ψm,ΨM, andΨuversusT̃ for T̃_Σc= 0.5.

502 It is now clear that the shape of an isotherm critically

503depends on the value of ψ. To make a more detailed

504description, plots ofΔμ̃(θ,r,T̃,T̃_Σc,ψ) versusθare shown in the

f5 505third row of Figure 5, where the parameters are chosen as

506follows:T̃= 0.4,T̃_Σc= 0.5;ris given its maximum value,r_g(T̃),

507consistent with the assumption that the ambientfluid is a gas;

508finally,ψ<Ψm(left-hand plot),Ψm<ψ<ΨM(center plot),

509and ψ > ΨM (right-hand plot). Obviously, the value of ψ

510determines the position of the curve relative to theθ-axis, given

511that r=r_g.

512 One can now drawθe versusr_efor the three ψ-cases given

513above by gradually reducing r_e from r_g to 0 and getting the

514corresponding value(s) of θe. Graphically, this means trans-

515lating the curves shown in the third row vertically upward, and

516noting the values ofθeat which the black lines cross theθ-axis:

517the intersection of the black line to the left givesθea, whereas

the black line to the right givesθed. Numerically, by usingeqs 518

11 and 12. Functions θea and θed plotted against r_e or, 519

equivalently against the pressure, are shown as solid lines in the 520

fourth row ofFigure 5. It is graphically obvious thatθeddoes 521

not exist for the case 0 < ψ<Ψm, and that both θeaand θed 522

exist for larger ψ-values, albeit for a limited range for ror P̃ 523

values. The dotted lines that connect the desorption and the 524

adsorption branches are the isotherms given byeq 6. See, for 525

comparison, the right plot in Figure 2, and note, in particular 526

that P̃(r₁,T̃) =P̃₁and thatP̃(r₂,T̃) =P̃₂. 527

One sees that the adsorption branches, shown in the fourth 528

row, only reach the valuer=r₂<r_g, ifψis large enough that 529

the local maximum of theΔμ̃versusθcurve is below theθ-axis 530

whenr=r_g, so that lifting the curve by reducingrbrings the 531

local maximum on theθ-axis whenr=r₂: see the third column 532

of the fourth row. Note also that, when this occurs, the 533

adsorption branch stops being defined and that it seems to534

have a vertical tangent atr=r₂. Similar statements hold for the 535

desorption branch. 536

The statements that the adsorption branch stops at r₂, 537

whereas the desorption branch stops at r₁, both with vertical 538

tangents, are correct because an equilibrium isotherm is a curve 539

satisfying μ̃Σ,red(θe,T̃,T̃_Σc) − μ̃f,red(r_e,T̃) − ψ = 0, whose 540

tangents are given by 541

r d r d

/ /

e e

f,red e ,red e

θ μ

μ θ

= ∂ ̃ ∂

∂ ̃_Σ ∂

Thus, the local extrema of theμ̃Σ,redversusθcurve produce 542

vertical tangents on the isotherms, and there is agreement with 543

the case of the analytical equilibrium isotherms, Section 2.3. 544

The agreement goes, in fact, deeper because it is possible to 545

repeat the stability analysis of Section 2.3 for the isotherms, 546

and also for the points with vertical tangents, with the same 547

conclusions. 548

It is now possible to interpret the multivaluedness in the 549

equilibrium isotherms, by introducing a hysteresis loop. 550

Referringfirst to the plots of the third column, fourth and551

fifth rows of Figure 5, one interprets the two-valuedness of 552

coverage as in Section 2.3: there is an adsorption curve 553

abcdβα, and a desorption curve αβγδba. As pointed out in 554

Figure 4.Specialψ′s versus dimensionless temperature,T̃. The upper thick line isΨMdefined byeq 14, the lower thick line isΨmdefined by eq 13, and the thin line isΨu, defined by eq 8. The first two are defined for 0 <T̃ < T̃_Σc, whereas the last is defined forT̃ > 0. All curves are drawn withT̃_Σc= 0.5.

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX F

555Section 2.3, no solution ofeq 11or ofeq 12will be situated on

556one of the vertical boundaries of the hysteresis loop bcdβγδ.

557 The ambient densitiesr₁andr₂, identifying the left and right

558boundaries of the hysteresis loop (see the fourth row inFigure

5595) are functions of T̃, T̃_Σc, and ψ, easily calculable by the

560following expressions

r x (0, ) such thatr ( , ,x T T, , ) 0,

( )

1 g m c

m

θ ψ

ψ

= ∈ Δμ ̃ ̃ =

> Ψ

Σ

561 (15)

r x (0, ) such thatr ( , ,x T T, , ) 0,

( )

2 g M c

M

θ ψ

ψ

= ∈ Δμ ̃ ̃ =

> Ψ

Σ

562 (16)

Turning now to the plots of the second column, fourth and 563

fifth rows ofFigure 5, it is seen that the upper solid line of the 564

plot in the fourth row must be discarded as a possible 565

equilibrium isotherm, at least if one assumes that desorption 566

occurs after adsorption, because adsorption stops before the 567

mathematically unstable point (whereθ=θM) is reached, and 568

the adsorbate has not transited to the desorption branch. 569

Adsorption occurs along abc, stops at c because of theP̃ <P̃₀ 570

condition, and then desorption follows the adsorption branch 571

but in the reverse direction, that is, along cba. 572

An important conclusion follows: there are two main classes 573

for the values of ψ: the class ψ≤ ΨMin which the isotherms 574

show low coverage and are structureless (fifth row, left-hand 575

and center plots inFigure 5) and the classψ>ΨMin which the 576

Figure 5.Illustrating the numerical calculation of equilibrium isotherms for subcritical temperatures. All plots are done withT̃= 0.4 andT̃_Σc= 0.5.

The values ofψare, column-wise:ψ<Ψm(left),Ψm<ψ<ΨM(center), andψ>ΨM(right). For the plots of the third row, it is reminded that +

T T r T T

( , , ) ( , ) ( , )

,red c f,red

μ μ θ μ ψ

Δ ̃ = ̃ ̃ ̃ − ̃ ̃ − ̃

Σ Σ . See text inSection 2.5, aftereq 14.

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX G

577isotherms show an adsorbate phase transition and a hysteresis

578loop (fifth row, right-hand plot).

579 When ψ > ΨM, it is of some interest to predict the

580approximate pressure of the center of the loop, and also the

581pressures of its left and right vertical boundaries. The pressure

582in the approximate middle of the loop isP̃(r_i,T̃) wherer_iis the

583previously found density at which the ambient chemical

584potential has an inflection point, see eq 10. The pressures

585approximating the left and right vertical boundaries areP̃(r₁,T̃)

586andP̃(r₂,T̃) where approximate expressions can be found for r₁

587and r₂ in a manner similar to the one that led to eq 10.

588Equations 15, and16are equivalent to

f T K f r T

f T K f r T

( , , ) ( , ),

( , , ) ( , )

m H f 1

M H f 2

θ τ

θ τ

̃ ̃ = ̃ ̃ ̃

̃ ̃ = ̃ ̃ ̃

Σ Σ

589and assuming thatr₁andr₂are small, onefinds

r T T E E

K

r T T E E

K

( , , ) ( ) e ( )

,

( , , ) ( ) e ( )

T

T

1 c m

/ m

H

2 c M / M

H

ψ τ τ

ψ ψ τ

̃ ̃ = =

̃

̃ ̃ = =

̃

ψ

ψ Σ

− ̃

Σ − ̃

590 (17)

591where

i

kjjjjj y

{zzzzz

E f T

( ) ( T, , )

1 exp

1

27

m 4

m m

m

m m

τ θ τ θ m

θ

θ θ

θ

= τ

̃ ̃̃ =

− − −

Σ

592 (18)

593(E_Mbeing obtained from above by substituting M to m) and

594where θmand θMare the functions ofτgiven byeq 37.

595 In the same approximation for which the above expressions

596are valid (low values of r), one can write the width of the

597hysteresis loop, in units of pressure, asw= 8P_cT̃(r₂−r₁) (see

598eq 33, right), that is

w P T E E

P T E E

K

8 ( ) ( ) e

8 ( ) ( )

T

c M m

/

c M m

H

τ τ

τ τ

= ̃[ − ]

= ̃[ − ]

̃

ψ

− ̃

599 (19)

f6 600 E_M−E_mversusτis shown inFigure 6. The three pressures

601P̃(r_i,T̃),P̃(r₁,T̃), andP̃(r₂,T̃) are indicated by vertical lines on

602the plot of thefifth row and third column inFigure 5.

3. RESULT AND DISCUSSION 2: THEORY OF TIME-DEPENDENT ISOTHERMS 603

Time-dependent solutions of eq 59 are now considered, that 604

simulate the thought experiment described in the next 605

paragraph below, leading to a description of dynamical 606

isotherms, and to their behavior in relation to equilibrium 607

isotherms and to hysteresis loops. Two subcases are 608

considered: the supercritical and the subcritical adsorbate. It 609

is believed that the actual values ofT̃ andT̃_Σcare of secondary 610

importance as compared to their relative sizes. T̃_Σc= 0.5 has 611

been used repeatedly above, and is also used in what follows, 612

whereasT̃ = 0.6 andT̃ = 0.4 are used for the two cases. 613

The thought-experiment considered is as follows. An 614

amount of mesoporous or macroporous medium is placed in 615

a containerfilled with a one-component gas at low pressureP 616

and at uniform constant temperatureT, which is less that the 617

critical temperature. After the gas in the pores has settled into a 618

state of zero velocity and uniform pressure P_i, which is less 619

than the saturation pressureP₀, the pressure in the gas is slowly 620

increased to a pressureP_f≤P₀, then slowly decreased back to 621

P_i. Recordings of the amounts adsorbed, and of the 622

corresponding ambient densities or pressures are assumed to 623

occur continuously. The duration of the cycleP_i→P_f→P_iis624

assumed to be large enough for the ambient density to remain 625

nearly uniform, and for the ambient velocity to be nearly zero, 626

in the macroporous medium at all times. (See the description 627

of model M′ in theIntroduction.) 628

The cycle of applied pressure forces the ambient density to629

follow a similar cycle. In M′, r(t̃) is needed as extra input to 630

solveeq 59. Mathematically, one can introduce a functionP̃(t)̃ 631

andfind the resultingr(t̃) by solving the equation of state. A632

simplification is described in what follows, that avoids time- 633

consuming calculations by starting directly with a functionr(t)̃ 634

that goes through a cycle, starting from a low-value r_ϵ, 635

increasing to a high-value r_t≤ r_g, then decreasing back to r_ϵ. 636 637 f7

Such a function is shown inFigure 7, where the increasing and

decreasing parts are linear. The minimumr_ϵis introduced so as 638

to avoid the singularity of the chemical potential atr= 0, and a 639

value r_ϵ = r_t/10³ is in general sufficient. The figure defines 640

7

r( ; )t̃, where7is the set of three parameters {r_ϵ,r_t,t_d̃}. The time parameter,t_d̃ will be called the cycle duration. 641

Equation 59 is solved^b with an initial condition, 642

r 7 T

(0) ( ; 0) exp( / )

θ = ψ ̃ . Concerning the arguments ofΔμ̃, r( ; )7 t̃ is substituted to r and, as already mentioned,T̃_Σc= 0.5 and two values are considered for T̃, that is 0.4 and 0.6. The 643

central purpose of the calculations is to determine the 644

influence of the values of ψ and of the cycle durationt_d̃ on 645

the shapes of the time-dependent isotherms. 646

Figure 6. For τ= T̃/T̃_Σc < 1, the width of the hysteresis loop is proportional toE_M−E_m(seeeq 19), a function ofτthat vanishes atτ

= 0 andτ= 1. Its maximum, 0.03281 is attained whenτ= 0.5917.

Figure 7.Functionr( ; )7 t̃, where7is the set of three parameters {r_ϵ, rt,t̃_d}. The function simulates a compression−decompression cycle.

DOI:10.1021/acsomega.9b02956 ACS OmegaXXXX, XXX, XXX−XXX H