Robust Three-Phase Vapor-Liquid-Asphaltene Equilibrium

Subscriber access provided by Lancaster University Library

Thermodynamics, Transport, and Fluid Mechanics

Robust Three-Phase Vapor-Liquid-Asphaltene Equilibrium Calculation Algorithm for Isothermal CO2 Flooding Applications Ruixue Li, and Huazhou Andy Li Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.9b02740 • Publication Date (Web): 05 Aug 2019 Downloaded from pubs.acs.org on August 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

Industrial & Engineering Chemistry Research

Robust Three-Phase Vapor-Liquid-Asphaltene Equilibrium Calculation Algorithm for Isothermal CO2 Flooding Applications

Ruixue Li1, 2 and Huazhou Li1* 1School of Mining and Petroleum Engineering, Faculty of Engineering, University of Alberta, Edmonton, Canada T6G 1H9 2 College of Energy, Chengdu University of Technology, Chengdu, 610059, PR China

*Corresponding Author: Dr. Huazhou Andy Li Assistant Professor, Petroleum Engineering University of Alberta Phone: 1-780-492-1738 Email: [email protected] 1 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 55

45

Abstract

46

CO2 flooding is an effective enhanced oil recovery process for light oil reservoirs.

47

Asphaltenes can easily precipitate during CO2 flooding, leading to the appearance of

48

three-phase vapor-liquid-asphaltene (VLS) equilibria. A prerequisite for accurately

49

simulating the CO2 flooding process is developing a robust three-phase VLS equilibrium

50

calculation algorithm. In this study, we develop a robust and efficient three-phase VLS

51

equilibrium calculation algorithm with the use of asphaltene-precipitation model proposed by

52

Nghiem et al. [1]. To develop this algorithm, a two-phase flash calculation algorithm is first

53

developed to split the mixture into an asphaltene phase and a non-asphaltene phase.

54

Moreover, two different three-phase VLS flash calculation algorithms are developed and

55

incorporated into our three-phase equilibrium calculation algorithm. New initialization

56

approaches for both stability test and flash calculations are proposed. The performance of this

57

three-phase

58

pressure-composition (P-X) diagrams for several reservoir fluids mixed with pure or impure

59

CO2.

VLS

equilibrium

calculation

algorithm

60

2 ACS Paragon Plus Environment

is

tested

by

generating

Page 3 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


61

1. Introduction

62

Carbon dioxide (CO2) flooding becomes one of the most effective enhanced oil recovery

63

processes for light oil reservoirs [2]. During CO2 flooding, asphaltenes can easily precipitate

64

due to the change in the composition of reservoir fluids or the change in pressure/temperature

65

conditions [3, 4]. Therefore, three-phase vapor-liquid-asphaltene (VLS) equilibria can

66

frequently appear in light oil reservoirs where CO2 flooding is implemented. To have a better

67

understanding of the CO2 flooding process in light oil reservoirs, it is necessary to have a

68

compositional simulator that is capable of handling three-phase VLS equilibria. In

69

compositional simulations, phase equilibrium calculation will be conducted in each grid at

70

each time step, implying that phase equilibrium calculations will be conducted at varying

71

conditions. In a CO2 flooding where reservoir temperature can be considered as constant, the

72

varying conditions correspond to varying pressures and gas concentrations. Therefore, a

73

prerequisite for building a robust compositional simulator is having an efficient and robust

74

three-phase VLS equilibrium calculation algorithm, which performs well at varying

75

pressures, varying gas concentrations and the given reservoir temperature.

76

To build such three-phase VLS equilibrium calculation algorithm, we first need to select a

77

thermodynamic model that can describe the asphaltene-precipitation phenomenon. Over the

78

last 30 years, there are many asphaltene-precipitation models proposed in the literature based

79

on either colloidal theory or solubility theory [5]. In the colloidal theory, asphaltenes are

80

regarded as colloidal particles with resins adsorbed onto their surfaces [6]. If enough resin

81

desorbs, asphaltenes will deposit [7]. In the solubility theory, asphaltenes are assumed to be



82

dissolved in crude oil [8]. If the solubility falls below a certain value, asphaltenes begin to

83

deposit. There are two main approaches in the solubility theory, i.e., regular solution model

84

and equation of state (EOS) [5]. Since a large number of three-phase equilibrium calculations

85

are required to be conducted in compositional simulations, we would like to select an

86

asphaltene-precipitation model with relatively low mathematical complexity that can be

87

readily incorporated into the three-phase VLS equilibrium calculation algorithm. In 1993,

88

Nghiem et al. [1] developed an efficient thermodynamic model for asphaltene precipitation.

89

In their model, vapor and liquid phases are modeled using a cubic EOS, while the

90

precipitated asphaltenes are regarded as a pure dense phase. Moreover, they assumed that the

91

heaviest component in crude oil can be divided into a precipitating component and a

92

non-precipitating component. These two components have identical critical properties and

93

acentric factors, but different binary interaction parameters (BIPs) with respect to light

94

components. A good agreement can be seen between the experimental results and the results

95

yielded by the model with some parameters estimated from experimental data.

96

One challenge in multiphase equilibrium calculations is how to determine the number of

97

present phases in the system. Multiphase equilibrium calculations consist of phase stability

98

test and flash calculations, which are normally conducted in a stage-wise manner [9-11].

99

Phase stability tests are conducted to determine whether the tested mixture is stable or will

100

split into two or more phases [9], while flash calculations are performed at the given number

101

of phases to calculate the phase mole fractions and phase compositions [10]. The most widely

102

applied stage-wise method is alternately conducting phase stability tests and flash


Page 4 of 55

Page 5 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


103

calculations [9]. For example, in a three-phase equilibrium calculation, stability test is first

104

conducted on the feed to check if a new phase needs to be introduced by locating the

105

stationary points on the tangent plane distance (TPD) function [9]. If a negative value of the

106

TPD function is found at any of the stationary points (implying that the feed is unstable), a

107

two-phase flash calculation needs to be conducted; otherwise, the given mixture is considered

108

to be stable as one phase. Then, stability test is again conducted on the phases yielded by the

109

two-phase flash calculation. If those phases are unstable, a three-phase flash calculation is

110

required to split the mixture into three phases; otherwise, the given mixture is considered to

111

split into two phases.

112

However, convergence problems are frequently encountered in phase equilibrium

113

calculations for three or more than three phases, even when one rigorously applies the

114

aforementioned stage-wise procedure. One reason causing the convergence issue is that there

115

may be several stationary points on the TPD surface in three-phase equilibrium calculations.

116

Thus, in three-phase equilibrium calculations, there is a high chance that the stability test fails

117

to detect the negative stationary point(s), since the negative stationary point(s) may be hidden

118

by the non-negative one(s). To increase the possibility of detecting the negative stationary

119

point(s), stability tests are required to be conducted with multiple sets of initial equilibrium

120

ratios [9, 11]. For different types of fluids, different initialization methods might be required.

121

For hydrocarbon fluids mixed with CO2, Li and Firoozabadi [11] suggested that (c+4) initial

122

estimates of equilibrium ratios should be applied for stability tests (where c represents the

123

number of components).



124

Another reason leading to the convergence problems encountered in three-phase equilibrium

125

calculations is that multiple two-phase equilibria may be yielded by two-phase flash

126

calculations that are conducted with different initial equilibrium ratios [12]. This requires one

127

to first determine if the calculated two-phase equilibrium can be used for the subsequent

128

calculations. Gorucu and Johns [13] recommended that the one with the minimum Gibbs free

129

energy should be applied to the subsequent calculations. To yield the two-phase equilibrium

130

with the minimum Gibbs free energy, Gorucu and Johns [13] suggested that two-phase flash

131

calculations should be conducted several times with the use of different initial equilibrium

132

ratios. However, it is time-consuming to repeat two-phase flash calculations for several times.

133

Another way to prevent yielding the trivial two-phase equilibrium and thereby having a

134

chance of yielding a trivial three-phase equilibrium is to perform three-phase equilibrium

135

calculations with an alternative stage-wise method. In this alternative stage-wise method, a

136

three-phase flash calculation is first conducted. If any calculated phase fractions is negative, a

137

two-phase flash calculation needs to be conducted [14, 15]. If any phase fractions calculated

138

by the two-phase flash calculation is negative, the feed is considered to be stable as one

139

phase. Li and Li [16] proposed that the results yielded by the three-phase flash calculation

140

can provide a knowledge of the present phases in the two-phase equilibrium. Negative flash

141

[14, 15, 17] should be allowed in this method. However, it is difficult to provide good initial

142

estimates for equilibrium ratios in three-phase flash calculations [11]. Lapene et al. [18]

143

developed a three-phase vapor-liquid-aqueous (VLA) flash calculation algorithm by

144

assuming that the aqueous phase is pure water phase, i.e., the so-called free-water assumption

145

[19]. With this assumption, the equilibrium ratios of the aqueous phase with respect to the 6 ACS Paragon Plus Environment

Page 6 of 55

Page 7 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


146

reference phase for all the non-aqueous components are zero. Therefore, the number of

147

equilibrium ratio sets (which need to be initialized in the three-phase VLA flash calculation)

148

is reduced, making it possible and easy to provide good initial estimates for equilibrium

149

ratios. They provided efficient initialization method for equilibrium ratios in their three-phase

150

VLA flash calculation algorithm. Li and Li [16] modified this three-phase flash calculation

151

algorithm as well as the alternative stage-wise procedure. By integrating the modified

152

three-phase flash calculation algorithm into the modified alternative stage-wise procedure,

153

they developed an efficient and robust three-phase free-water vapor-liquid-aqueous (VLA)

154

equilibrium calculation algorithm [16]. Since the assumption made on the asphaltene phase in

155

the asphaltene-precipitation model proposed by Nghiem et al. [1] is similar to the free-water

156

assumption, the number of equilibrium ratio sets in the three-phase VLS flash calculation can

157

be also reduced when implementing the asphaltene-precipitation model proposed by Nghiem

158

et al. [1]. Thus, it is possible for us to develop a robust three-phase VLS equilibrium

159

calculation algorithm by taking advantage of the free-solid assumption. To our knowledge,

160

such robust algorithm has not been reported in the literature.

161

In this work, we develop a robust three-phase VLS equilibrium calculation algorithm that is

162

coupled with the asphaltene-precipitation model proposed by Nghiem et al. [1]. The key

163

parameters used in this model are first adjusted based on experimental data. As previously

164

mentioned, such asphaltene-precipitation model assumes that asphaltene phase only contains

165

pure asphaltene. This assumption enables the appearance of the asphaltene phase to be

166

checked simply by comparing the fugacity of asphaltene component in the non-asphaltene



Page 8 of 55

167

phase against that in the asphaltene phase. The appearance of the non-asphaltene phase (i.e.,

168

vapor or liquid phase in this study) is determined by the phase stability test. New

169

initialization methods for equilibrium ratios are provided for both stability test and flash

170

calculations. To develop this three-phase VLS equilibrium calculation algorithm, we first

171

develop a two-phase flash calculation algorithm used to split the mixture into an asphaltene

172

phase and a non-asphaltene phase. Moreover, in order to enhance the robustness and

173

efficiency of our three-phase VLS equilibrium calculation algorithm, two different

174

three-phase VLS flash calculations algorithms are developed and incorporated into our

175

three-phase VLS equilibrium calculation algorithm. The performance of the newly developed

176

algorithm is tested by generating pressure-composition (P-X) diagrams for several real

177

reservoir fluids mixed with pure or impure CO2.

178

2. Mathematical Formulations

179

2.1 Thermodynamic Model

180

Nghiem et al. [1] proposed a thermodynamic model in which vapor and liquid phases are

181

modeled using a cubic EOS, while asphaltene phase is considered as a pure dense phase only

182

containing asphaltene. The fugacity of asphaltene in the asphaltene phase is calculated by the

183

following equation [1],

184

ln f s  ln f  * s

Vs  P  P*  RT

(1)

185

where P and P* represent the actual pressure and the reference pressure, respectively; fs and

186

fs* represent the fugacities of asphaltene in the asphaltene phase at P and P*, respectively; Vs

187

represents the asphaltene molar volume; R represents the universal gas constant; T represents 8 ACS Paragon Plus Environment

Page 9 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


188

temperature; and the subscript s represents asphaltene phase. In this work, the Peng-Robinson

189

(PR) EOS [20] is applied to model vapor and liquid phases.

190

2.2 Phase Stability Test

191

In this study, asphaltene precipitation can be simply checked by the following inequality [21],  f cx  f cs , asphaltene phase does not exist   f cx  f cs , asphaltene phase exists

192

(2)

193

where fcx and fcs represent the fugacity of asphaltene component in the non-asphaltene phase

194

calculated by PR EOS [20] and that in the asphaltene phase calculated by Equation (1),

195

respectively, and the subscripts c, x, and s represent asphaltene component, non-asphaltene

196

phase, and asphaltene phase, respectively. The appearance of vapor or liquid phase is checked

197

by phase stability test [9]. The expression of the TPD function used in the stability analysis is

198

shown as follows [9], c

199

g  y    yi ln i  y   ln yi  ln i  z   ln zi 

(3)

i 1

200

where y and z represent the compositions of the trial phase and the tested phase, respectively;

201

i ( y ) and i ( z ) represent the fugacity coefficients of component i in the trial phase and the

202

test phase, respectively; c is the number of components in the mixture; and the subscript i is

203

component index. If the values of the TPD function are non-negative at any composition of

204

the trial phase, the tested phase is considered to be stable. Michelsen proposed a stability test

205

algorithm [9], suggesting finding all stationary points on the TPD surface. If the values of all

206

stationary points are non-negative, the tested phase can be considered to be stable. The

207

stationary points on the TPD function can be detected by the following equation [9],



ln i  y   ln yi  ln i  z   ln zi  k , i  1,..., c

208

Page 10 of 55

(4)

209

where k is the value of the stationary point (which is a constant). The detailed procedure for

210

conducting stability test is provided by Michelsen [9]. In this work, the quasi-Newton

211

successive-substitution method [22] is used in the stability test to search for the local

212

stationary points. Moreover, multiple estimates of equilibrium ratios are used to initialize the

213

stability test. Section 3 will provide the detailed initialization method used in our stability

214

tests.

215

2.3 Flash Calculation

216

Flash calculations aim to obtain phase fractions and phase compositions for a given mixture

217

equilibrating at a given pressure/temperature condition and a given number of present phases

218

[10]. The convergence of flash calculations is subject to the material balance constraint and

219

the equal-fugacity constraint (i.e., the fugacities of each component in all present phases

220

should be equivalent) [10].

221

2.3.1 Two-Phase Vapor-Liquid (VL) Flash Calculation

222

The equal-fugacity constraint for the two-phase VL flash calculation is shown as follows

223 224

[10],

fix  fiy , i  1,..., c

(5)

225

where fix and fiy represent the fugacities of component i in liquid phase and vapor phase,

226

respectively, and the subscripts x and y represent liquid phase and vapor phase, respectively.

227

Rachford-Rice (RR) equation is developed to solve phase mole fractions based on the

228

material balance [23]. RR equation used for the two-phase VL flash calculation is provided as 10 ACS Paragon Plus Environment

Page 11 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

229


follows [23], c

c

i 1

i 1

RR    yi  xi   

230

zi  K iy  1

1   y  K iy  1

0

(6)

231

where xi, yi, and zi represent the mole fractions of component i in liquid phase, vapor phase,

232

and feed, respectively; βy represents the mole fraction of vapor phase; and Kiy represents the

233

equilibrium ratio of component i in vapor phase with respect to liquid phase. Kiy can be

234

calculated by,

K iy 

235

yi ix  , i  1,..., c x i  iy

(7)

236

where ix and iy represent the fugacity coefficients of component i in liquid phase and

237

vapor phase, respectively. Based on material balance, phase compositions can be calculated

238

by [23],

239

xi 

zi ; yi  zi K iy 1   y  K iy  1

(8)

240

The two-phase VL flash calculation algorithm used in this research is the same as that

241

provided by Michelsen [10]. The two-phase VL flash calculation algorithm is composed of

242

two iteration loops. The outer loop aims to satisfy the equal-fugacity constraint by updating

243

equilibrium ratios, while the inner loop solves RR equation to yield phase mole fractions and

244

phase compositions. In this work, successive substitution iteration (SSI) is used in the outer

245

loop, while Newton method is applied in the inner loop [24, 25].

246

2.3.2 Two-Phase Vapor/Liquid-Asphaltene (V/L-S) Flash Calculation

247

In two-phase V/L-S flash calculation, one of the present phases is the asphaltene phase.

248

Based on the assumption that the asphaltene phase only contains asphaltene, the two-phase 11 ACS Paragon Plus Environment


Page 12 of 55

249

V/L-S flash calculation can be much simpler than the two-phase VL flash calculation. The

250

equal-fugacity constraint for the two-phase V/L-S flash calculation can be simplified to one

251

equation, f cx  f cs

252

(9)

253

where fcx and fcs represent the fugacities of asphaltene in non-asphaltene phase and asphaltene

254

phase, respectively, and the subscripts c, x, and s represent asphaltene component,

255

non-asphaltene phase, and asphaltene phase, respectively. fcs can be calculated based on

256

Equation (1), while fcx can be calculated based on PR EOS with a given composition of

257

non-asphaltene phase. Based on the material balance, composition of non-asphaltene phase

258

can be calculated by the following equation,

259

 z / 1   s  , i  1,..., c  1 xi   i   zi   s  / 1   s  , i  c

(10)

260

where xi and zi represent the mole fractions of component i in non-asphaltene phase and feed,

261

respectively, and βs represents mole fraction of asphaltene phase. Therefore, the

262

equal-fugacity constraint can be satisfied by updating βs. Appendix A provides the detailed

263

procedure of the two-phase V/L-S flash calculation algorithm.

264

2.3.3 Three-Phase VLS Flash Calculation

265

The equal-fugacity constraints for the three-phase VLS flash calculation are given by [1],

266

fix  fiy , i  1,..., c

(11)

267

f cx  f cy  f cs

(12)

268

In this work, two algorithms are developed to conduct the three-phase VLS flash calculations.

269

Both of them are incorporated into the three-phase VLS equilibrium calculation algorithm to 12 ACS Paragon Plus Environment

Page 13 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


270

enhance the robustness and efficiency of the new three-phase VLS equilibrium calculation

271

algorithm. Section 3 will provide the details about how to incorporate these two algorithms.

272

The difference in these two algorithms lies in the method used for solving the phase

273

compositions and phase fractions.

274

Li and Li [16] derived the objective function and constraints (used to solve phase mole

275

fractions) for three-phase free-water VLA flash calculation by simplifying the objective

276

function and constraints developed by Okuno et al. [26] based on the free-water assumption

277

[19]. The objective function and constraints used in Algorithm #1 for three-phase VLS flash

278

calculations is similar to that used in the three-phase free-water VLA flash calculation

279

algorithm developed by Li and Li [16]. The objective function and constraints for three-phase

280

VLS flash calculations can be expressed as follows [16], c 1  min : f   zi ln K iy  K cy  y  K cz     y  i 1  Subject to : K   K   min K   z , K   K z , i  c cy iy y cz i cz iy i 



281

282

283











(13)

where

  1  yc  K cy  1  x  c   K   1  zc  cz 1  xc

(14)

284

where xc, yc, and zc represent the mole fractions of asphaltene component in liquid phase,

285

vapor phase, and feed, respectively. As noted by Li and Li [16], there is only one unknown in

286

the above objective function, i.e., phase mole fraction of vapor phase (βy). Moreover, the

287

constraints in Equation (14) form a reliable feasible region of βy [16, 26]. The mole fractions

288

of asphaltene phase and liquid phase can be solved by the following equations [16], 13 ACS Paragon Plus Environment


s 

289

zc   xc  yc   y  xc 1  xc

;  x  1   y  s

Page 14 of 55

(15)

290

where βx represents the mole fraction of liquid phase. The compositions of liquid and vapor

291

phases can be calculated based on the following equations [16], zi   K  K   K ,i  c  iy cy y cz xi   ; yi  K iy xi 1  ,i  c  K is 

292

293

(16)

where

K is 

294

si ix  , i  1,..., c x i  is

(17)

295

where si represents the mole fraction of component i in asphaltene phase and is represents

296

the fugacity coefficient of component i in asphaltene phase. Based on the assumption that

297

asphaltene phase only contains asphaltene, we can have,

1, i  c si   0, i  c

298

(18)

299

The procedures adopted by Algorithm #1 is similar to that of the three-phase free-water VLA

300

flash calculation algorithm developed by Li and Li [16]. This algorithm can be simply

301

summarized as two iteration loops (which is similar to the two-phase VL flash calculation

302

algorithm). The outer loop is used to update equilibrium ratios for satisfying the

303

equal-fugacity constraint by SSI, while the inner loop is used to calculate phase mole

304

fractions and phase compositions by minimizing Equation (13) using the interior-point

305

method. This algorithm is efficient; but the performance of Algorithm #1 relies on the quality

306

of initial equilibrium ratios. In Section 3, we provide detailed initialization method for

307

equilibrium ratios in Algorithm #1. 14 ACS Paragon Plus Environment

Page 15 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


308

Algorithm #2 for three-phase VLS flash calculations is motivated by two-phase V/L-A flash

309

calculation algorithm. In Algorithm #2, we perform two-phase VL flash calculations in the

310

inner loop to satisfy Equation (11) in the equal-fugacity constraint, while the outer loop is

311

used to update βs to give a new feed used in the inner loop. Equation (12) in the

312

equal-fugacity constraint should be satisfied in the outer loop. The feed used in the inner loop

313

can be updated by,

314

 z / 1   s  , i  1,..., c  1 zi 2   i   zi   s  / 1   s  , i  c

(19)

315

where zi2 represents the mole fraction of component i in the feed used in the inner loop.

316

Algorithm #2 is reliable when three-phase equilibrium actually appears; but it is

317

time-consuming since several times of two-phase VL flash calculations are conducted in the

318

inner loop. The detailed procedure of Algorithm #2 is summarized in Appendix B.

319

3. Three-Phase VLS Equilibrium Calculation Algorithm

320

In this section, we will introduce the robust three-phase VLS equilibrium calculation

321

algorithm. In order to ensure the robustness of the three-phase VLS equilibrium calculation

322

algorithm, it is essential to have good initial values for both stability tests and flash

323

calculations. In general, for different types of fluids, different initialization methods are

324

required. To apply our algorithm to simulate CO2 flooding in light oil reservoirs, we provide

325

new initialization methods for equilibrium ratios in both stability tests and flash calculations.

326

The new initialization method for equilibrium ratios in stability tests is a modification of that

327

proposed by Li and Firoozabadi [11]. This new method is shown as follows,



328



K iST  K iWilson ,1/ K iWilson , 3 K iWilson ,1/ 3 K iWilson , K iCO2  NS , K iCO2  S

Page 16 of 55



(20)

329

where K iST represents the initial equilibrium ratio of component i used for stability tests.

330

The first four terms are the same as those proposed by Li and Firoozabadi [11]. K iWilson

331

represents the equilibrium ratio of component i calculated by Wilson equation [27],

332

K iWilson 

Pci   T exp 5.37 1  i  1  ci P  T 

  

(21)

333

where P and Pci represent pressure and the critical pressure of component i, respectively; T

334

and Tci represent temperature and the critical temperature of component i, respectively; and

335

ωi represents acentric factor of component i. K iCO2  NS and K iCO2  S can be calculated by the

336

following equations, respectively,

337

338

0.8 / zi , i  2  K iCO2  NS    0.2 /  c  1  / zi , i  2

K iCO2  S

0.8 / zi , i  2   108 / zi , i  c    0.2  108 / c  2  / z , i  2 and c    i  

(22)

(23)

339

where 2 and c represent components CO2 and asphaltene, respectively. The new initialization

340

method for equilibrium ratios in flash calculations will be introduced in the detailed

341

procedure of the three-phase VLS equilibrium calculation algorithm. Figure 1 depicts the

342

flow chart of this algorithm. The detailed procedure is listed below,

343

(1) Using the new initialization method for equilibrium ratios as shown in Equation (20),

344

conduct the stability test on the feed (zi) to detect the negative stationary points. Record

345

the compositions of the trial phases at all stationary points with negative TPD values as

346

x_0i. x_01i, x_02i, etc. Note that x_0i corresponds to the smallest TPD value, while the 16 ACS Paragon Plus Environment

Page 17 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

347 348 349 350 351


last trial phase composition corresponds to the largest TPD value. (2) If negative stationary point exists, vapor or liquid phase may appear. If so, go to Step (5). Otherwise, continue to Step (3). (3) Check Inequality (2). If asphaltene phase exists, continue to Step (4); otherwise, output single-phase equilibrium.

352

(4) Conduct two-phase V/L-S flash calculation and output two-phase V/L-S equilibrium.

353

(5) Check Inequality (2). If asphaltene phase exists, go to Step (8); otherwise, continue to

354

Step (6).

355

(6) Conduct two-phase VL flash calculation to yield liquid-phase and vapor-phase

356

compositions (xi2-VL and yi2-VL) and mole fractions of liquid phase and vapor phase (βx2-VL

357

and βy2-VL). The initial equilibrium ratios used in the two-phase VL flash calculation can

358

be generated by the following equation, K iy  zi / x _ 0i , i  1,..., c

359

(24)

360

(7) Check Inequality (2) again. If asphaltene phase does not exist, output two-phase VL

361

equilibrium; otherwise, conduct three-phase VLS flash calculation Algorithm #1 and

362

output three-phase VLS equilibrium. The initial equilibrium ratios used in the three-phase

363

VLS flash calculation are given by,

364 365

366

367

K iy  yi2VL / xi , i  1,..., c

(25)

 xi2VL / 1  0.99  xc2VL  , i  c  xi   2 VL 2 VL 0.01 xi / 1  0.99  xc  , i  c

(26)

where

(8) Conduct three-phase VLS flash calculation Algorithm #1 to yield liquid-phase, 17 ACS Paragon Plus Environment


Page 18 of 55

368

vapor-phase, and asphaltene-phase compositions (xi3-1, yi3-1, and si3-1) and phase mole

369

fractions of liquid phase, vapor phase, and asphaltene phase (βx3-1, βy3-1, and βs3-1). KiWilson

370

are used as the initial equilibrium ratios.

371

(9) If the results calculated by the three-phase VLS flash are physical, output three-phase

372

VLS equilibrium; otherwise, continue to Step (10). This judging criterion is the same as

373

that proposed by Li and Li [16]. As mentioned by Li and Li [16], the results of a

374

three-phase flash calculation can be considered to be unphysical if an open feasible region

375

(i.e., a region where any of the endpoints is a pole) appears in any iteration or any of the

376

ultimately calculated phase mole fractions do not fall into [0, 1].

377 378

(10) If βs3-1 can be calculated and βs3-1tol, update mole fraction of asphaltene phase (βs) and go back to Step (3); otherwise, output the calculated phase mole fractions and phase compositions.



876 877

Figure A-1 Flow chart of the two-phase V/L-S flash calculation algorithm

878 879


Page 52 of 55

Page 53 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


880

Appendix B: Three-Phase VLS Flash Calculation Algorithm #2

881

Figure B-1 depicts the flow chart of the three-phase VLS flash calculation Algorithm #2. The

882

detailed procedure of this algorithms is summarized below,

883

(1) Input the pressure (P), temperature (T), feed (zi), and thermodynamic properties of each

884

component (Tci, Pci, and ωi).

885

(2) Initialize mole fraction of asphaltene phase (βs) and set its tolerance (tol) to be 10-6.

886

(3) Calculate the composition of non-asphaltene phase (zi2) using Equation (19).

887

(4) Conduct two-phase VL flash calculation to calculate phase mole fractions of liquid phase

888

and vapor phase (βx and βy) and compositions of liquid and vapor phases (xi and yi).

889

(5) Calculate the fugacities of each component in liquid and vapor phases (fix and fiy) based on

890

PR EOS.

891

(6) Calculate the absolute relative error (err) between the fugacities of asphaltene component

892

in the non-asphaltene phase (fcx) and asphaltene phase (fcs) using Equation (A-1).

893

(7) If err>tol, update mole fraction of asphaltene phase (βs) and go back to Step (3);

894

otherwise, update phase mole fractions of liquid phase and vapor phase (βx and βy) by the

895

following equation and output the calculated phase fractions and phase compositions,

896

 x   x  1   s  ;  y   y  1   s 


(B-1)


897 898 899 900 901

Figure B-1 Flow chart of the three-phase VLS flash calculation Algorithm #2


Page 54 of 55

Page 55 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

902


Abstract Graphic

903


Robust Three-Phase Vapor-Liquid-Asphaltene Equilibrium

Recommend Documents