January 4, 2017

Engineering Optimization

root

To appear in Engineering Optimization
Vol. 00, No. 00, Month 20XX, 1–11

1

Modelling and optimization for a wellhead gas flowmeter using
concentric pipes

3

Yana Nec † and Greg Huculak‡

4

(Received 00 Month 20XX; final version received 00 Month 20XX)

5

A novel configuration of a landfill wellhead was analyzed to measure flow rate of gas extracted
from sanitary landfills. The device provides access points for pressure measurement integral
to flow rate computation similarly to orifice and Venturi meters, and has the advantage
of eliminating the problem of water condensation often impairing accuracy thereof. It is
proved that the proposed configuration entails comparable computational complexity and
negligible sensitivity to geometric parameters. Calibration for the new device was attained
using a custom optimization procedure, operating on a quadri-dimensional parameter surface
evincing discontinuity and non-smoothness.

6

Keywords: flow rate, landfill gas, flowmeter, Darcy friction factor, optimization on discontinuous surfaces

1.

2

Background

7
8
9
10
11
12
13
14
15

16

Landfill gas collection often requires installation of flowmeters to comply with environmental regulations. Landfill gas is extracted under vacuum at numerous well points, each
monitored for flow rate and gas composition. Methane, carbon dioxide, nitrogen and oxygen are commonly contained in the gas stream extracted. In the past numerous methods
have been employed to measure wellhead flow rates: orifice plates, Venturi meters as well
as other commercial devices. These devices have been used with some success, however
their accuracy is impaired when wet gases are encountered. At times space requirements
render their use inappropriate. The geometry of the new device addresses these issues.
The operation principle of an orifice flowmeter is briefly reviewed here to facilitate comparison with the proposed configuration infra. The orifice flowmeter comprises a plate
with a centred aperture that is to occlude the fluid conduit, and two sensors to measure
the pressure drop due to the occlusion. The flow rate is computed by means of a theoretical formula based on considerations of momentum and supplemented by an empiric
discharge coefficient accounting for phyiscal phenomena responsible for head loss not
captured by elementary conservation of momentum (Crane 1982). Quondam simplistic
models for this coefficient (Idel’chik 1960) proved unsound, the reason thereto designated
circa 1980s as a marked sensitivity of orifice calibration to the location of pressure sensors.
Since these two entities must be in close conformance to attain adequate accuracy of flow
rate estimation, the sensor locations were standardized to a prescribed distance upstream
and immediately past the plate, wherewith extensive experiements in conjunction with
comprehensive modelling begot the accepted nowadays Reader-Harris / Gallagher dis† Thompson Rivers University 900 McGill rd. Kamloops, British Columbia, Canada. Corresponding author. Email:
cranberryana@gmail.com
‡ GNH Consulting, Delta, British Columbia, Canada. Email: greg@gnhconsulting.ca

1

17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

January 4, 2017

Engineering Optimization

root

65

charge coefficient (ISO 2003). However, somewhat cumbersome form thereof and iterative
calculation process entailed a praxis of commercial orifice flowmeters being calibrated by
the manufacturer in accord with the foregoing principles, providing a formula requiring scarce computational effort, but also delegating the full responsibility for no longer
modifiable installation locations of pressure sensors to the field operator.
The following circumstances render the commercial calibration of the orifice flowmeter
incompatible with the geometry of a landfill well (Nec and Huculak 2015). In contrast to
as universal as it is tacit assumption of horizontal flow used in orifice flowmeters, the well
is vertical. Typical flow rates are low, at times necessitating reduction in the aperture
diameter to obtain a discernible pressure drop for the measurement, whence momentum
loss by gravity is on the same order of magnitude as due to the occlusion by orifice plate.
Thus applying a calibration constant issued for a horizontal flow impairs the accuracy of
flow rate estimation. Furthermore, it is impossible to instal the second pressure sensor
immediately behind the orifice plate due to water vapour, an ever-present component
entrained in the landfill gas, condensing thereon. Therefore the pressure sensor locations
do not conform to standards. In a recent study adequate accuracy of flow rate estimation was achieved by a custom calibration procedure, incorporating an effective relative
roughness parameter to account for turbulence engendered by the presence of the orifice,
as well as discovering a linear dependence on the constriction ratio of plate aperture to
pipe diameter (Nec and Huculak 2015).
The current contribution proposes a novel wellhead geometry, wherein orifice flowmeter
usage is discontinued, eliminating the problem of moisture pooling on the orifice plate
and interference with pressure measurement. The new arrangement is shown to be insensitive to imprecision inevitable in field installation, whilst the computational complexity
involved in flow rate determination is on par with the custom calibrated orifice flowmeter
counterpart (Nec and Huculak 2015). The feasibility of calibration of the new wellhead
is verified through a series of measurements, however experiments comprehending the
geometric and flow parameter space in its entirety are beyond the ambit of this study.

66

2.

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

67
68
69
70
71
72
73
74
75
76

77

Geometry and flow equations

The wellhead studied herein consists of a tube of inner radius r and wall thickness δ
inserted concentrically and secured in a well pipe of radius R > r, whose upper end is
blocked (see figure 1). The upstream pressure sensor is located by the annulus wall and
measures pressure within the nearly stagnant eddy zone. The tube ends with a sharp
elbow. A small horizontal recess holds the downstream pressure sensor, measuring the
static pressure at that point. The elbow outlet connects to pipework collecting gas from
nodes throughout the landfill.
For analysis of mass and momentum balance hereunder consider the control volume
from the upstream plane to the elbow inlet plane. Assume the flow steady, axisymmetric
and incompressible. Integral mass conservation equation reads (Batchelor 1990)

∮

∂V

78
79

80

81

ρu ⋅ ds = 0,

(1)

with ∂V denoting the surface of the control volume V, ρ being fluid density, u –velocity
vector and s – area vector with the normal directed outwards. Completing the integration,

−πR2 uup + πr2 ui⌞ = 0,

(2)

wherein uup and ui⌞ are upstream and elbow inlet velocities respectively. Defining the

2

January 4, 2017

Engineering Optimization

root

elbow outlet plane

elbow recess wall
downstream sensor
elbow inlet plane

ℓout
annulus wall
upstream sensor
ℓ
annulus entry

L
effective flow crosssection adjustment

upstream plane

Figure 1.: Wellhead flow geometry. Eddy region is shaded grey.

constriction ratio β = r/R, equation (2) becomes

82

uup = β 2 ui⌞ .

(3)

Integral momentum conservation equation reads (Batchelor 1990)

∫

∂V

(ρu ⋅ ds)u + ∫

∂V

84

pds = ∫ ρ fbody dV + fsurf ,

(4)

V

with p denoting pressure and f( ⋅ ) – force vectors, remaining quantities defined heretofore.
Thus

⎧

⎫

⎪
1⎪
⎪
⎪
2
πR2 u2up − πr2 u2i⌞ + ⎨πR2 pup − π (R2 − (r + δ ) )paw − π ((r + δ )2 − r2 )pae − πr2 pi⌞ ⎬ = (5)
ρ⎪
⎪
⎩

⎪
⎪
⎭

⎫
⎧
⎪
⎪
⎞
⎛
⎪
⎪
2
g ⎨πR2 L − πℓ( (r + δ ) − r2 ) + πr2 ℓout ⎬ + u2up cf πRL + π (r + δ )ℓ + πr(ℓ + ℓout ) ,
⎪
⎪
⎠
⎝
⎪
⎪
⎭
⎩

wherein as before the subscripts ( ⋅ )up and ( ⋅ )i⌞ refer to upstream and elbow inlet

3

83

85

86
87

88

89

January 4, 2017

Engineering Optimization

90
91
92
93
94
95
96
97
98
99
100

root

planes, and ( ⋅ )aw and ( ⋅ )ae designating respectively annulus wall and entry planes.
Observe that the outlet area of the control volume is threefold: annulus wall, tube wall
ring and tube interior, responsible for the negative terms within the first set of braces.
Gravity is the sole body force, g denoting the gravity constant, the second set of braces
containing the volume of fluid within the pipe and tube compound, with the lengths
L, ℓ, ℓout marked in figure 1. The surface force accounts for wall shear, by dimensional
analysis equalling a product of dynamic pressure 21 ρu2up , friction coefficient cf and area
affected thereby. Without loss of generality uup is a representative velocity for the purpose
of friction modelling, whence cf is in accord with that choice.
The pressure immediately beneath the tube ring is related to the pressure at the
annulus wall through a simple fluid column, since both values are for stagnant fluid:
pae = paw + ρgℓ,

101

102

103

104
105
106
107
108
109

110

111
112

allowing to replace pae in (5) and simplify, yielding

1
u2up − β 2 u2i⌞ + (pup − (1 − β 2 ) paw − β 2 pi⌞ ) =
ρ

115
116
117
118
119
120
121

122

(7)

⎧
⎫
⎪
δ ℓ
ℓ ℓout ⎪
⎪L
⎪
g (L + β 2 ℓout ) + u2up cf ⎨ + (β + ) + β ( +
)⎬.
⎪
R
R R
R
R ⎪
⎪
⎪
⎩
⎭

The brace delimited term in (7) will infra prove essential to support the negligible sensitivity of the studied flow geometry to variations of ℓ and ℓout , two parameters prone to
installation imprecision.
Identical considerations of mass and momentum in conjunction with dimensional analysis allow to introduce the head loss coefficient ζp due to a projection of one concentric
conduit within another
1
pup − pi⌞ = ρu2i⌞ (1 − β 2 ) ζp + ρg (L + ℓout ).
2

(8a)

Empiric studies show ζp to depend on thickness ratio δ /r (Idel’chik 1960). Analogously
for a flow through a sharp elbow with a recess
1
pi⌞ − po⌞ = ρ u2i⌞ ζ⌞ ,
2

113

114

(6)

(8b)

wherein the subscript ( ⋅ )o⌞ refers to elbow outlet plane, with the head loss coefficient ζ⌞
being a function of relative roughness of the tube inner surface ε and Reynolds number
(Idel’chik 1960). Hereunder qualitative use only is made of the dimensionless coefficients
ζp and ζ⌞ , rendering the accuracy and reproducibility of the specific values given in
Idel’chik (1960) of little moment. Since (8b) captures the head loss due to the presence
of the sharp bend alone, for all practical purposes po⌞ = pr⌞ , to wit the pressure measured
in the recess pr⌞ might be without loss of generality be deemed equal to po⌞ , the possible
differences absorbed in ζ⌞ . Combining (8a) and (8b) with (3) gives
pup = pr⌞ +

ρu2up ⎛
2β 4 ⎝

⎞
(1 − β 2 ) ζp + ζ⌞ + ρg (L + ℓout ).
⎠

4

(9a)

January 4, 2017

Engineering Optimization

root

Similarly from (8b) and (3)

123

pi⌞ = pr⌞ +

ρu2up
2β 4

ζ⌞ .

(9b)

Utilizing (9) to replace pup and pi⌞ in (7), upon elementary algebraic manipulation one
obtains
u2up
2β 4

C=

paw − pr⌞
− gℓout ,
ρ

(10)

wherein the coefficient C formally equals

124

125
126

127

128

⎫
⎧
2β 4 cf ⎪
δ ℓ ℓout⎞⎪
⎪
⎪L
2⎛
⎬.
⎨ + β (2 + ) +
C = ζp + ζ⌞ − 2β −
2
⎪
1 − β ⎪R
r r
r ⎠⎪
⎝
⎪
⎭
⎩
2

(11)

129

Expression (11) establishes negligibility of geometric minutiae’s impact on the hydraulic
resistance coefficient C, as is easily seen from the ascending powers of β, the leading
order given by ζp and ζ⌞ , both O(1). Interestingly, only even powers of β appear. Since
0 < β < 1, the power of β 6 , for instance, implies that the variation of ℓ and ℓout must be
O(β −6 ) for that term to bear on the value of C, incontrovertibly exceeding conceivable
adjustments made in the course of installation and operation manyfold. This insensitivity
renders the reliability of the flow rate to be derived from (10) hereinafter preferable to
that of the corresponding orifice estimate evincing marked sensitivity to the location of
the pressure gauges. Here the locations are such as to make incorrect installation virtually
impossible, involving no measurements and none of the deftness and experience required
for an orifice. Therefore only the qualitative dependence C (β, ε, Re) is of import, akin
to the Reader-Harris / Gallagher discharge coefficient for the orifice (ISO 2003) and the
modified coefficient developed in Nec and Huculak (2015).
From (10)

136

(12)

144

√ 2√
paw − pr⌞
2β
uup = √
− gℓout ,
ρ
C

yielding the volumetric flow rate

130
131
132
133
134
135

137
138
139
140
141
142
143

145

√
√
paw − pr⌞
2 πr2
q= √
− gℓout .
ρ
C

(13)

Invoking the state equation for ideal gas p = ρRT , where R and T refer to gas constant
and temperature respectively, juxtaposition of (13) above and result (13) of Nec and
Huculak (2015) forthwith reveals the conceptual equivalence of the orifice flowmeter and
wellhead geometry suggested herein. Identically for both devices the gravity term is
negligible by comparison to the pressure drop term for high flow rates, becoming coequal
only for nearly stagnant wells.
Albeit the suitability of the geometry analyzed here might not appear surprising given
that the wellhead has been ideated with the distinctive features of the landfill flow regime
in mind, whilst the orifice is a generic device, to attain its full potential, formula (13)
must be furnished with an adequate calculation procedure for the hydraulic resistance coefficient C. Hence it is the authors’ intention to glean the qualitative form of C (β, ε, Re),

5

146

147
148
149
150
151
152
153
154
155
156
157

January 4, 2017

Engineering Optimization

158
159

root

ascertain the feasibility of calibration and verify that overall computational complexity
does not exceed that existing for related models of similar accuracy.

kRe

A
Re×
Re
Figure 2.: Typical functional shape of kRe (Re).

160

161
162
163
164
165
166
167
168
169

3.

Functional form of C

By (11) the resistance coefficient is a function of head loss coefficients ζp and ζ⌞ (ε, Re),
as well as the constriction ratio β and parameters pertaining to longitudinal geometry
together with the friction factor cf . Bearing in mind that formally 0 < β < 1 and in praxis
0 < β < 21 , identically to the orifice, retainment of terms with powers higher than square
is incongruent with the measurement precision expected for the pressure difference in
(13). The generic friction factor cf is of the same order of magnitude as, for instance,
Darcy friction factor f (Moody 1944), and for the case of a simple pipe flow can be shown
to equal cf = 41 f (Nec and Huculak 2016). Therefore cf ∼ o(1) and expression (11) is
henceforth curtailed to read
C = ζp + ζ⌞ − 2β 2 + o (β 4 ) .

170

171
172

173

174
175

(14)

The qualitative dependence of ζp and ζ⌞ (ε, Re) is adopted from empiric studies (Idel’chik
1960) and consequently modified below to suit the current application.
3.1

Coefficient ζp

The head loss coefficient ζp accounts for the disturbance to flow in the main pipe due
to projection of the inner tube thereinto. In general ζp depends on the projection length

6

January 4, 2017

Engineering Optimization

root

ℓ/r as well as the thickness ratio δ /r, however for sufficient δ /r the dependence becomes
trivial (Idel’chik 1960, pg. 98), as is indeed the case here with 0.1 < δ /r < 0.2. Therefore
henceforward ζp is regarded a constant.

3.2

Coefficient ζ⌞ (ε, Re)

176
177
178

179

The head loss coefficient ζ⌞ pertains to the abrupt change in flow direction at the wellhead outlet, its two arguments being relative roughness of the conduit material and
Reynolds number, here ranging 5×10−5 < ε < 20×10−5 and 104 < Re < 105 respectively.
The functional dependence of ζ⌞ comprises three parts (Idel’chik 1960, pg. 215, table
6-11): generic constant, factor due to Re and factor due to ε:
ζ⌞ = c kRe kε .

(15)

The quantitative dependence thereof as given in Idel’chik (1960) is as follows. The parameter kRe is defined by
kRe = {

Re < 40000
,
Re ⩾ 40000

45f
1

(16a)

where f refers to Darcy-Weisbach friction coefficient, for a fully turbulent regime obtained
by solution of Colebrook equation (Colebrook 1939)
2.51
ε
1
√
√ = −2 log10 (
+
3.7 Re f
f

).

(16b)

In (16a) the numeric factor 45 is interlocked with the cross-over Reynolds number 40000,
so that the curve kRe (Re) is a continuous function, but non-differentiable at the cross-over
point. A typical dependence is shown in figure 2. It is desired to preserve the aforesaid
functional dependence whilst converting all fixed empiric constants into variables open
to optimization. Therefore definition (16a) is generalized as

⎧
f
⎪
⎪
⎪A
kRe = ⎨ f×
⎪
⎪
⎪
⎩A

Re < Re×

Re ⩾ Re×

,

(17a)

wherein A supplants unity in (16a) and f× maintains continuity of kRe through the
solution to Colebrook equation with Re× :
2.51
1
ε
√ ) = 0.
√ + 2 log10 (
+
3.7 Re× f×
f×

(17b)

For the geometry at hand Re× is expected to be lower than in (16a) due to turbulence
engendered by the projecting tube before the bend in flow within the elbow, begetting
an earlier cross-over. Furthermore, with the introduction of A in (17a) the multiplicative
constant c in (15) is to be set to equal unity without loss of generality.
The most general functional form of kε is given by
kε = {

1
1 + Aε ε

7

Re < Re×
.
Re ⩾ Re×

(18)

180
181
182
183
184

185

186
187

188

189
190

191

192
193
194
195
196

197

198
199

200

201
202
203
204
205

206

January 4, 2017

Engineering Optimization

root

212

In Idel’chik (1960) the nominal value of Aε , corresponding to Re× = 40000, is Aε = 500,
however here Aε is a degree of freedom.
In light of the generalizations above four parameters are to be determined before
the computation of resistance coefficient C in (14) can be effected: ζp , A, Re× , Aε . This
quadruple is obtained infra by minimization of a cost function based on a set of independent measurements.

213

4.

207
208
209
210
211

214
215
216
217
218
219
220

221
222

223

224

225

226
227
228
229
230
231

232

233
234
235
236
237
238
239
240
241

Optimization

The optimization centres on fitting the four variables spanning the parameter space of
C with the purpose to show that given a set of trustworthy flow rate measurements,
calibration of C can be effected with the stated functional forms for ζp and ζ⌞ . Mathematically, taking a set of measured flow rates {qm } and corresponding estimates {q } as
computed by (13), find a set of parameters {ζp , A, Re× , Aε } such that the norm ∣∣q − qm ∣∣n
is minimal, n designating the norm’s exponent (solved here for n = 1 and n = 2 with
qualitatively similar results):
min ∣∣q − qm ∣∣n

(19a)

√
√
paw − pr⌞
2 πr2
− gℓout ,
s.t.: q = √
ρ
C

⎧
f
⎪
⎪
⎪
A
ζ⌞ = ⎨
f×
⎪
⎪ A(1 + Aε ε)
⎪
⎩

C (ζp , A, Re× , Aε ; ε, β ) = ζp + ζ⌞ − 2β ,
2

ζp , A, Re× , Aε > 0.

(19b)

Re < Re×

Re ⩾ Re×

(19c)

This problem is unusual from several aspects, the combination thereof compounding the
complexity of what at first glance might appear an easily amenable task. First, one of
the optimization variables, the cross-over Reynolds number Re× , is interlocked with the
estimate q by a doubly implicit relation: the computation of each estimate point in the
set q requires the hydraulic resistance coefficient C that depends on the flow Reynolds
number Re, which depends on q by
Re =

4qρ
,
πµd

(20)

and furthermore the discontinuity point of C is exactly Re× , meaning that for any tested
value of Re× it is unknown beforehand whether Re < Re× or Re ⩾ Re× . Second, C is a
function of four variables, the quadruple of optimization parameters, but also depending
on two parameters ε and β. The relative roughness ε can be deemed fixed, as the pipe
material is not expected to change during the optimization procedure. By contrast, the
diameter ratio β, whilst fixed for a given geometric configuration, must perforce change
for wellheads of different typical flow rates. For a sound accuracy the pressure term
(paw − pr⌞ )/ρ in (19b) must significantly exceed the gravity term gℓout , thence for an
individual well the geometry might be adjusted depending on the gas production: for

8

January 4, 2017

Engineering Optimization

root

high flow rates β can be as high as moiety, whilst low flow rates might necessitate
β < 61 . This means that the four arguments of C in (19c) all depend on β. Thence
the minimum in (19a) attained upon solution is in fact β-dependent. In this light, an

additional, difficult to quantify constraint regards the possible solutions {ζp , A, Re× , Aε }:
it is undesirable to have the order of magnitude of these quantities vary significantly
with β, albeit a smooth or even monotonic dependence is not expected due to possibly
qualitatively different turbulent flow regimes. Unreasonable variation, for instance over
several orders of magnitude, is an indicator the physical modelling is wanting. Third,
although (19) appears conceptually to be a classic curve fitting problem, it is not: the
flow rate q depends explicitly on the pressure difference, but also on temperature T
through the density ρ. The viscosity µ that affects Re also depends on temperature. In
reality the measurements {qm } are taken at different temperatures as well. Therefore
there is no curve in the classic sense: the available data correspond to disjoint points on
a multivariable surface q with an explicit dependence on T through ρ, and an implicit
one through Re. One might argue that a carefully controlled experiment will permit to
maintain a fixed temperature, however since this model is to be used in reality, artificially
controlling the temperature will nullify the applicability of the results of this experiment:
then the optimization will only be valid for the particular sub-space conforming to the
chosen temperature out of the entire physical parameter space.
Problem (19) is at an obvious variance with studies dedicated to provision of consistent
experimental data and subsequent comprehensive modelling (e.g. Colebrook equation for
Darcy friction factor and Reader-Harris / Gallagher discharge coefficient for orifice flow
rate). Therefore the authors forbear to explore the physical parameter space and focus in
lieu on flow regimes characteristic to a landfill, whose field data are summarily accessible
(courtesy of GNH Consulting, British Columbia).
Problem (19) was solved by an adaptation of a golden ratio search algorithm to handle
the implicit nature, discontinuity and non-smoothness of the involved functions. The
search procedure was run on the following part of the parameter space

{ζp , A, Re× , Aε } ∶ 0 < ζp < 1,

0 < A < 1.5,

103 < Re× < 105 ,

100 < Aε < 104 .

(21)

A series of local minima of (19a) were stored as the search progressed. None of these
were close to the bounds (21). As the bounds were set based on physically meaningful
quantities, it is unlikely the global minimum is situated outside (21). The constraint that
optimal values be of similar magnitude was exercised only upon completion of the entire
search. The existence of numerous local minima is discussed hereinafter.
Although some problems incorporating non-smooth functions are nonetheless amenable
to reformulations that permit application of gradient based algorithms, it is conjectured
that here this will be impossible due to the implicit nature of the constraints that further
compounds the inherent discontinuity.
4.1

Numerical considerations

242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

270

271
272
273
274
275
276
277
278
279

280

The computation of flow rate q in (13) is iterative: the hydraulic resistance coefficient
C depends on Re, whilst Re depends on q by equation (20). Therefore q is estimated
by (13) with an initial guess for C, followed by determination of Re from (20), then
repeated recomputation of C, q and Re until proper convergence1 . The optimization was
implemented in GNU Octave (open source software), utilizing the function fminbnd. For
1 This is identical to the procedure for an orifice, where Reader-Harris / Gallagher discharge coefficient is employed.

9

281
282
283
284
285

January 4, 2017

Engineering Optimization

286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306

root

each suggested quadruple ζp , A, Re× , Aε the hydraulic resistance coefficient C must be
converged.
The product of non-differentiable kRe , equation (17a), and discontinuous kε , equation
(18), renders ζ⌞ a function discontinuous in Re, the selfsame variable that creates the
implicit dependence of C on q. The discontinuity point is Re× , one of the degrees of
freedom in the optimization. This intricate inter-dependence poses certain numerical
complications, the most notable being dearth of convergence of C if the Reynolds number
is close to the discontinuity point. Even when the likelihood to obtain in the course of the
iterative computation of C the singular equality Re = Re× , minuscule as it might be, is
eliminated entirely by a negligible shift wheresoever necessary, Reynolds numbers falling
sufficiently close to Re× occasion toggling between the two branches of the discontinuous
function. When the computation of C does not require evaluation of ζ⌞ on both sides of
Re× , the convergence to precision of 10−6 is obtained within 3-6 iterations, on par with
Reader-Harris / Gallagher discharge coefficient that converges to said precision in 2-5
iterations (Nec and Huculak 2015).
A further concomitant of the discontinuity of ζ⌞ is related to the optimization, but
will not affect the estimation of flow rate upon completion thereof. It was found that
the parameter space spanned by the quadruple {ζp , A, Re× , Aε } offered multiple local
minima to the cost function (19a), when the fit against a set of independent measurements {qm } was performed. Interestingly, designating as a global optimum the quadruple {ζp , A, Re× , Aε } that gives minimal error, there exist several disparate quadruples
opt

312

{ζp , A, Re× , Aε } that correspond to local minima with negligibly higher error. This
loc
means the field operator will be at liberty to choose any one set of parameters without introducing perceptible inaccuracy into the estimation of flow rate. This situation
is unlikely in the extreme with a continuously differentiable multi-dimensional surface
underlying the optimization and is the direct outcome of the discontinuity in (18) and
non-smoothness in (17a).

313

4.2

307
308
309
310
311

314
315
316
317
318
319
320

Results

Three flow regimes were considered: high flow rate engendered by an active well, medium
flow rate corresponding to moderate production and low flow rate conforming to a slowly
producing well. In all cases a set of reliable measurements {qm } was chosen so as to span
as uniformly as possible the concomitant interval of differential pressure. The optimization was performed separately, verifying the combined results presented a physically
acceptable solution. Figures 3-5 depict the comparison between the measured flow rate
and estimate (13) with {ζp , A, Re× , Aε } .
opt

321
322
323
324
325
326
327
328
329
330
331
332

When a well is active, a typical constriction ratio is β ≈ 13 with ensuing Reynolds
numbers ranging 104 < Re < 6×104 . The corresponding cross-over Reynolds number Re× is
relvatively low, whereby for most of the working range the hydraulic resistance coefficient
C is constant, entailing a rapid convergence. When the Reynolds number is close to Re× ,
it will behove the operator to adjust the wellhead to a smaller constriction ratio so as to
maintain a viable pressure difference. An example of this situation is shown in figure 3.
If gas production within the landfill cavity in the proximity of a well diminishes, the
constriction ratio will be reduced to β ≈ 14 in order to sustain the same working range
of Reynolds numbers, this time Re× found in the upper part thereof, again enabling
effortless convergence in praxis. Figure 4 details this intermediate flow regime.
For slow production wells the suitable constriction ratio might be as low as β ≈ 15 , with
a similar range of Reynolds numbers and Re× falling midmost (figure 5). This is the only

10

January 4, 2017

Engineering Optimization

root

regime, where toggling around the discontinuity point is likely to present difficulty. As a
rule, higher constriction ratios are preferable as long as the pressure difference is aptly
measurable. Therefore a configuration with such low β will be installed for Re < Re× and
changed in favour of a higher β if the cross-over point is approached.
The optimal quadruples ζp , A, Re× , Aε reported in figures 3-5 corresponded to the global
error minimum for each set of reference measurements. Whilst it is possible to choose
one of the quadruples conforming to local minima for any particular constriction ratio,
only the global optima entail a logical adjustment of flow regimes with the fluctuations
in gas production, when the entire spectrum of operation is considered.

5.

Conclusion

333
334
335
336
337
338
339
340
341

342

The novel wellhead configuration presented targeted a long-standing problem of impaired
flow rate metering in landfill wells due to condensation of water vapour entrained in the
fluid upon the orifice plate and interference with the pressure sensor that must be placed
in an immediate proximity thereof. In the proposed setting the fluid proceeds from the
main well pipe into a concentric tube of a smaller diameter, permitting to withdraw the
orifice plate. The resulting flow rate was proved to be conceptually equivalent to the
orifice device by considerations of mass and momentum conservation. The sensitivity of
concomitant pressure measurement to installation geometry was shown to be minor, and
the associated error was evaluated asymptotically in constriction ratio β, an inherent
small parameter of the system.
The hydraulic resistance coefficient was modelled by adopting the functional form of
empiric coefficients in related geometries and performing a custom optimization to attain
a fit against a set of independent measurements. The optimization involved a quadridimensional discontinuous surface with an implcit dependence that required iterative
numerical solution. In field use the computational complexity was shown to be on par
with devices endowed with similar accuracy of flow rate estimation.

References

343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358

359

Batchelor, G.K. 1990. Introduction to fluid dynamics. Cambridge.
Colebrook, C.F. 1939. “Turbulent flow in pipes, with particular reference to the transition region
between the smooth and rough pipe laws.” Journal of the Institution of Civil Engineers 11:
133–156.
Crane, Co. 1982. Flow of fluids through valves, fittings and pipe. Technical paper No. 410M. 300
Park Avenue, New York, NY 10022, appendix A-5.
Idel’chik, I.E. 1960. Handbook of hydraulic resistance. Coefficients of local resistance and of friction (translated from Russian). Published for US Atomic Energy Commission and NSF, Washington DC, AEC-tr-6630, 1966.
ISO. 2003. “Measurement of fluid flow by means of pressure differential devices inserted in circular
cross-section conduits running full - part 2: Orifice plates.” .
Moody, L.F. 1944. “Friction factors for pipe flow.” Transactions of ASME 66: 671–684.
Nec, Y., and G. Huculak. 2015. “On the problem of a vertical gas flow through an orifice with
non-standard pressure tappings locations.” Can.J.Civ.Eng. 42: 563–569.
Nec, Y., and G. Huculak. 2016. “Solution of weakly compressible isothermal flow in landfill gas
collection networks.” submitted to Fluid Dynamics Research .

11

360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375

Engineering Optimization

root

flow rate q [cfm]

35

20

5
1

2
pressure difference paw − pr⌞ [Wc”]

3

1.4
hydraulic resistance coefficient C

January 4, 2017

1.25

1.1
10

20

30

40

50

60

Reynolds number / 1000
Figure 3.: A typical example of flow regime for β = 31 . Top: optimization output (diamond) and independently measured
(cross) flow rate. Reference points included multifarious values of density ρ, viscosity µ and relative roughness ε, rendering
interpolation unfeasible. Bottom: hydraulic resistance coefficient C versus Re. Optimal parameters: ζp = 0.66, A = 0.71,
Re× = 13000, Aε = 7300. Diamond symbols on both panels correspond.

12

Engineering Optimization

root

flow rate q [cfm]

22

16

10

4
0

1

2
3
pressure difference paw − pr⌞ [Wc”]

1.95

hydraulic resistance coefficient C

January 4, 2017

4

5

1.75

1.55
10

20

30

40

Reynolds number / 1000
Figure 4.: A typical example of flow regime for β = 41 . Top: optimization output (diamond) and independently measured
(cross) flow rate. Reference points included multifarious values of density ρ, viscosity µ and relative roughness ε, rendering
interpolation unfeasible. Bottom: hydraulic resistance coefficient C versus Re. Optimal parameters: ζp = 0.72, A = 1.02,
Re× = 31800, Aε = 1000. Diamond symbols on both panels correspond.

13

Engineering Optimization

root

flow rate q [cfm]

15

10

5

5
pressure difference paw − pr⌞ [Wc”]

0

hydraulic resistance coefficient C

January 4, 2017

10

2.5

2.4

2.3
10

20

30

40

Reynolds number / 1000
Figure 5.: A typical example of flow regime for β = 51 . Top: optimization output (diamond) and independently measured
(cross) flow rate. Reference points included multifarious values of density ρ, viscosity µ and relative roughness ε, rendering
interpolation unfeasible. Bottom: hydraulic resistance coefficient C versus Re. Optimal parameters: ζp = 1.06, A = 1.33,
Re× = 23000, Aε = 800. Diamond symbols on both panels correspond.

14