Modeling ingress losses into a horizontal landfill gas well:
reconciling 2D turbulent flow and quasi-1D approximation
Istafa Qureshi
Mitacs Globalink Internship
Summer 2023

Thompson Rivers University
Department of Mathematics
Supervisor: Yana Nec

1

1

Introduction

This document details on the steps that the author has followed to try to incorporate turbulent ingress head
losses with the Darcian flow head loss through porous media in a horizontal landfill gas well to produce a
model for calculating the head loss characteristics across the gas well for efficient extraction of landfill gas[2].
The horizontal well is a pipe with equidistant apertures penetrating the length of a landfill, with a vacuum
inducing pump at the outlet and the far end blocked. The well is supported by a thin lamina of gravel and
drains a cylindrical volume of gas generating medium around it. In [2], a quasi-1D model for this system was
developed and implemented in MATLAB/Octave, but a computational study in COMSOL [1] has shown
that the quasi-1D construction overlooks the highly localized turbulent structures at the ingress apertures.
The purpose of the current project was to attempt to incorporate these local head losses into the quasi-1D
model.

2

COMSOL Model Review

To generate the simplified MATLAB/Octave model, a COMSOL model was created to compare against.
The gas well can be simplified to a 2D axisymmetric analogue. The circular apertures are replaced by
axisymmetric slits of equivalent area that wrap around the well, essentially splitting it into individual pipe
segments. The specifications of the nominal model are given in Table 1.
Symbol
µ
Rgas
rp
rg
rw
t
l
E
L
rh
nh
n
T
Patm
Pout
Pout
C
kg
kw
ϕ
k
ϵ

Variable
Dynamic viscosity
Gas constant
Pipe radius
Outer radius gravel domain
Outer radius waste domain
Pipe wall thickness
Pipe segment length
Pipe wall roughness
Wall length
Hole radius
Holes per section
Pipe segments
Temperature
Atmospheric pressure
Outlet pressure (Open BC)
Outlet pressure (Closed BC)
Generation rate
Gravel permeability
Waste permeability
Porosity
Turbulent kinetic energy
Turbulent dissipation rate

Value
0.000014291
287
0.0762(3”)
1
10
0.00711(0.28”)
15
1.5 × 10−6
420
0.00635(0.25”)
6
28
288.15
101325
98825
95075
1 × 10−6
1 × 10−5
1 × 10−8
0.55
0.00016
0.005

Units
Pa · s
J/kg/K
m
m
m
m
m
m
m
m
K
Pa
Pa
Pa
kg/s/m3
m2
m2
m2 /s2
m2 /s3

Table 1: Model Specifications

3

COMSOL Java Interface

The Java API for COMSOL was used to automate the data generation, extraction, and sorting. The
procedure for this is detailed in [4]. The rationale for using this interface over the COMSOL Parametric
Sweep was to export data after each parameter shift, as opposed to exporting after the complete scan.
Additionally, it allows for the entire parameter space to be run at once. However, the COMSOL .mph files
are not saved using the Java API (requires libraries that are not available, refer to model.save ), so all
necessary data must be explicitly extracted.
COMSOL uses a parametric solver, which can utilize previous solutions for faster convergence. Since
this solver is not designed for changing geometries, it runs into convergence issues if the involved parameters change the geometry. COMSOL parametric sweep is an alternative to this which works for changing

2

geometries and converges as fast as a Java scan. For initial tests and setup, it is highly recommended to use
the Parametric Sweep with the parametric solver1 that is available within COMSOL.
Some errors with the COMSOL GUI have been known to occur, one notable error is the GUI being
inconsistent to the Java code generated from the same file. This phenomenon was observed for one file and
the source of the disparity is yet to be identified. To ensure the credibility of the .mph file it is advisable to
export its associated Java file and check for any discrepancies.

4

Parametric Scan

The data stored from the parametric scans have been limited to pressure values at our points of interest at
the ingresses. Two probing points for each ingress have been selected at seven pipe thicknesses (7 × t) from
the center of the ingress on either side of the ingress; one within the turbulent pipe domain and the other
within the porous lamina. The parameters being changed are shown in Table 2.
Symbol
Pout
kw
l
nh
n
C

Parameter
Outlet pressure
Waste permeability
Pipe Segment length
Number of holes
Number of pipe segments
Generation Rate

Test range
88825 : 1250 : 100075
10−6:−1:−12
7.5 : 7.5 : 30
1, 2, 4, 6
6 : 27
10−5:−1:−7

Unit
Pa
m2
m
kg/s/m3

Table 2: Testing Ranges (The colon operator represents the MATLAB/Octave linspace function.)
For the closed boundary case, the outlet pressure was fixed at the median value of our testing range,
while varying other parameters. However, for the open boundary, lower pressure tends to cause convergence
problems, so a pressure value for stable convergence was used.

Mass collected for changing pressures

Mass collected for changing l

0.97

quasi-1D
COMSOL
predicted

0.52

Mass output(kg/s)

Mass output(kg/s)

0.81

0.34

0.3
0.22
0.19
0.13
8.8825

10.0075

Pressure (Pa)

10

7.5

4

30

l (m)

Figure 1: Output Mass for Open Boundary

1 Parametric solver is available in the Advanced Settings of Parametric Sweep in the COMSOL GUI

3

5

COMSOL Model Setup

The geometry of the COMSOL file can be understood from the .mph files, ergo, requires no detailed instructions.
The COMSOL model was created in two parts; the turbulent flow regime of the pipe and the Darcian
flow in the porous media. In the turbulent regime, the model was set to be compressible (M a < 0.3) with
an outlet having Pout set as fixed pressure value at one end of the well. The pressures at the ingresses of the
turbulent regime were equated to the corresponding pressure values at those ingresses obtained from Darcy’s
law. The radial and axial velocities were first initialized to zero for the closed boundary condition which
revealed that the radial velocities are minimal throughout the entire pipe domain excluding the vicinity of
the ingress. Moreover, the axial velocities showed a linear variation within the pipe domain which supported
previous evidence[2]. Hence, the radial velocities for subsequent setups were initialized to zero for the entire
domain and the axial velocities were initialized from the predicted linear Reynolds number variation given
in [2]. The initial value for pressure for the entire pipe domain was set to the atmospheric pressure. Within
the porous media, the pressure values have been initialized by a linear function. The pressure increases from
Pout to Patm moving radially outwards from the pipe exterior (rp + t) to the exterior of the waste layer (rw ).
For the closed boundary case, the outer radius of the waste layer at rw is set to have no flow across it
and the pressure at this point varies with the applied outlet pressure. An important characteristic of this
setup is that the absolute pressure within the domain never approaches atmospheric levels and the entire
domain is essentially under vacuum throughout most of our test range. The mass collected from the domain
is always equal to the mass generated in the waste layer.
The open boundary model has been created identically with the exception that the outer radius of the
waste layer is set to have a fixed inlet pressure equal to atmospheric pressure (Patm ). The mass collected
in this case varies directly with the applied pressure as shown in Figure 1. The available data suggest that
there are fixed differences in the mass collected from the domain due to the ingresses with varying pressures,
which proposes that it relates to the geometry of the ingress rather than the applied pressures,s further
supported by observing the mass output for changing l in Figure 1.

6

Juxtaposition with the quasi-1D model in MATLAB/Octave

6.1

Closed Boundary

When compared to the equivalent MATLAB/Octave setup the results, they have a similar pressure profile as
shown in Figure 2, which indicates a fair degree of accuracy. This can be attributed to the mass collected in
both systems being the same, as both are equal to the mass generated in the systems. For both cases, it was
noteworthy that for the entire test range of the outlet pressure, the pressure at the pipe’s closed end never
exceeded the outlet pressure by more than a fixed pressure difference which we will call Pm . Pm changes
with varying lengths of the pipe segments and is much lower than the applied pressure difference. For our
nominal case, Pm is about 500 Pa.

6.2

Open Boundary

When compared with the MATLAB/Octave solution, this case gives results as shown in Figure 3. This data
is obtained from a generic case and the pressure discrepancy persists for all cases. As we would expect, the
data from COMSOL approaches the maximum possible pressure within the well, while there is a substantial
pressure difference in the MATLAB/Octave solution. The disparity in both curves can be attributed to the
output masses being different in both cases as shown in Figure 1. This disparity may be resolved by either
forcing the ingress pressure losses or the masses collected through the ingresses. Both options are explored
hereunder.

7

COMSOL Files and Datasets

To create a model for the ingress head loss, the datasets were generated from COMSOL Java API of the test
variables given in Table 2, excluding the parameters that were deemed unrealistic or showed convergence
complications. For future use, the simulations have also been run on COMSOL using parametric scan via
the parametric solver. The incomplete or missing scans are given in Table 3.

4

Pressure profile along pipe length

P - Pout

Pm

quasi-1D
COMSOL

0
0

L

Distance from outlet
Figure 2: Closed BC Pressure Profile
Parameter
C
n
n
nh

Case
Closed Boundary
Open Boundary
Closed Boundary
Open Boundary

Deficient Files
.mph & .csv
.mph & .csv
.mph
.mph

Error Type
Unrealistic Conditions
Convergence Issues
Convergence Issues
Convergence Issues

Table 3: COMSOL Deficient Files & Datasets

8

Ingress Loss Modeling

The incongruity between the quasi-1D and 2D models stems from the head loss at the ingresses. For the
quasi-1D solution, the losses over the ingresses are not accounted for, so the system opts to disperse the head
loss over the porous domain as shown in Figure 4. Consequently, we need to incorporate the ingress losses
into the quasi-1D solution, which is done by forcing the model to use the correct pressure predicted in the
gravel layer by COMSOL to solve for the integration constant aa defined in appendix B in [2].
This pressure is obtained from the COMSOL data by fitting a quartic polynomial
ping = c1 i4 + c2 i3 + c3 i2 + c4 i1 + c5 i0 ,

(1)

where cn are constants and i is the ingress position. The values of cn were obtained from the parametric scan
data and are calculated based on the input conditions. Depending on the input conditions, the polynomial
for pressure loss was created and added to the corresponding quasi-1D section pressure to simulate the
pressure outside the ingress, namely in the gravel layer; resulting in the following calculation:




1 2
rA
rA
ka rA
µ
2
2
2
2
(r − rA ) + rA ln
ln
−
ln
,
(2)
aa = pB − (pi + ping ) + RT Cb
kb
2 B
rB
rP
kb
rB
where the term (pi + ping ) is supposed to simulate the pressure inside the gravel layer. This condition is
only applicable for the open boundary case in MATLAB/Octave, as the closed boundary case is modeled
independently of the ingress losses. The characteristics observed for the variables are detailed hereinafter.

5

Pressure profile along pipe length

Patm-P out

P - Pout

quasi-1D with losses
quasi-1D without losses
COMSOL

0
0

L

Distance from outlet
Figure 3: Open BC Pressure Profile

8.1

Closed Boundary

As expected, changing the geometry alters the resistance to the flow within the pipe. Reducing nh increases
the resistance to flow at the ingress, resulting in higher losses. Increasing l results in greater mass generation
in porous laminae and higher pipe wall friction fostering higher losses.
As the ingress loss is a very localized phenomenon, there is no significant change that occurs with changing
kw as expected[1].
Changing the outlet pressure does not produce an appreciable difference in the ingress pressure drop. This
can be attributed to the mass generated in the system being constant.
Increasing the generation rate of the system tends to increase the pressure throughout, resulting in some
regions in the domain exceeding atmospheric pressure, yielding impractical conditions.

8.2

Open Boundary

The geometry change for the open boundary shows the same characteristics as the closed boundary condition;
increasing resistance results in higher ingress head loss. Changing kw tends to have no effect as previously
discussed. Additionally, changing the generation rate has inappreciable impact in the system, which can
be attributed to a significant fraction of the mass collected coming from the open boundary and not inner
generation. Changing outlet pressure in this case, does in fact, change the ingress pressure drop. The lower
the outlet pressure, the greater the pressure drop across the ingresses due to the increased mass being drawn
through the ingress.
Incorporating the ingress losses in this case tends to cause a drop in pressure along the pipe due to less
mass being collected from the porous domain as seen in Figure 3, which proves that forcing ingress losses
is not a viable solution for this case. Forcing ingress losses tends to decrease the mass collected from the
system, resulting in reduced pressure drop across the shared domain boundary and causing the pressure to
drop in the porous domain as shown in Figure 4.

6

Radial Pressure Profile
Patm-P out

COMSOL
quasi-1D

0
0 rp

rw
Figure 4: Open BC Radial Pressure Profile

9

Forced Darcy Pressure Loss

In the open boundary case, the pressure loss across the porous domain remained fairly constant for all
parameters. At the ingress closest to the outlet, the pressure drop across the porous domain was minute
(< 10P a), which reduced to null when moving to the closed end of the well. So, to replicate this in
MATLAB/Octave, the pressure was changed to these values to force the model; resulting in the following
solution for aa




µ
rA
ka rA
1 2
rA
2
2
2
2
aa = pB − pdi + RT Cb
(r − rA ) + rA ln
ln
−
ln
,
(3)
kb
2 B
rB
rP
kb
rB
where


pdi = patm − pc ×

n−i
n


(4)

and pc is the maximum pressure drop across the porous media. However, this results in the reduction in
the overall mass collected from the porous domain resulting in the same outcome as before. It reduces the
pressure drop in the pipe even further resulting in almost no pressure loss over the turbulent domain as
shown in Figure 5.
To check whether this method can yield any results whatsoever, a value for pc was selected which forced
the masses in the quasi-1D solution to match the values known from COMSOL. This enforces the solution
to use a set value of integration constant aa yielding the solution given in Figure 6. The value of pc was
assigned empirically but it showed that the pressure profiles match if the masses in both solution are the
same, which reveals that the solution is possible with the quasi-1D but it is not clear where the inadequacies
lie. This same process can be repeated for all test parameters with the understanding that the variation of
assigned pressures is not always linear axially, yielding that there is a possible solution equating the quasi-1D
model to the COMSOL model.

7

Patm-P out

Radial Pressure Profile

Patm-P out

Pressure profile along pipe length

COMSOL
quasi-1D

0

0

P-P out

P - Pout

quasi-1D
COMSOL

0

0 rp

L

Distance from outlet

rw

Radial Distance

Figure 5: Forced Darcy Pressure Loss Profiles

Pressure profile along pipe length

P - Pout

Patm-P out

quasi-1D
COMSOL

0
0

L

Distance from outlet
Figure 6: Arbitrary pc Pressure Profile

8

References
[1]

Abhishek Baral, Mehrdad Sepehri, and Yana Nec. “Effects of turbulence at the ingress into landfill gas
wells”. In: Journal of Engineering Mathematics 135, 7 (2022).

[2]

Yana Nec and Greg Huculak. “Landfill gas flow: collection by horizontal wells”. In: Transport in Porous
Media 130 (2019), pp. 769–797.

[3]

Yana Nec and Greg Huculak. “Solution of weakly compressible isothermal flow in landfill gas collection
networks”. In: Fluid Dynamics Research 49:065505 (2017).

[4]

Istafa Qureshi. Tutorial for COMSOL API using Eclipse. 2023.

9