NMR Measurement and Lattice-Boltzman Simulation of the non-local dispersion tensor
Available from
Andrew Jackson's profile on Mendeley.
Page 1
NMR Measurement and Lattice-Boltzman Simulation of the non-local dispersion tensor
NMR Measurement and Lattice-Boltzmann
Simulation of the non-local dispersion tensor
M.W. Hunter a A.N. Jackson a;b P.T. Callaghan a
aMacDiarmid Institute for Advanced Materials and Nanotechnolgy, School of
Chemical and Physical Sciences, Victoria University of Wellington, Wellington,
New Zealand
bThe British Library, Boston Spa, West Yorkshire, United Kingdom
Abstract
The non-local dispersion tensor DNL provides a fundamental description of veloc-
ity correlations and displacement information in a dispersive system. It is shown
that Pulsed Gradient Spin Echo NMR (PGSE-NMR) can be used to measure this
tensor, and we present here the ¯rst measurement of DNL in a complex °ow by
this or any other method. These measurements are complemented by simulations
based on a lattice-Boltzmann calculation of the °uid °ow. For dispersive °ow in
a random beadpack of monosized spheres, six non-zero, independent components
remain. These components have been measured at three times less than ¿v, the
time to °ow one bead diameter. It is shown here that the various elements of DNL
provide insights regarding the dispersive °ow which are extremely sensitive to the
details of local correlations.
Key words: Non-local dispersion tensor, PGSE-NMR, Lattice-Boltzmann, Porous
Media
Preprint submitted to Elsevier 13 July 2009
Simulation of the non-local dispersion tensor
M.W. Hunter a A.N. Jackson a;b P.T. Callaghan a
aMacDiarmid Institute for Advanced Materials and Nanotechnolgy, School of
Chemical and Physical Sciences, Victoria University of Wellington, Wellington,
New Zealand
bThe British Library, Boston Spa, West Yorkshire, United Kingdom
Abstract
The non-local dispersion tensor DNL provides a fundamental description of veloc-
ity correlations and displacement information in a dispersive system. It is shown
that Pulsed Gradient Spin Echo NMR (PGSE-NMR) can be used to measure this
tensor, and we present here the ¯rst measurement of DNL in a complex °ow by
this or any other method. These measurements are complemented by simulations
based on a lattice-Boltzmann calculation of the °uid °ow. For dispersive °ow in
a random beadpack of monosized spheres, six non-zero, independent components
remain. These components have been measured at three times less than ¿v, the
time to °ow one bead diameter. It is shown here that the various elements of DNL
provide insights regarding the dispersive °ow which are extremely sensitive to the
details of local correlations.
Key words: Non-local dispersion tensor, PGSE-NMR, Lattice-Boltzmann, Porous
Media
Preprint submitted to Elsevier 13 July 2009
Page 2
PACS:
1 Introduction
The science of °uid dispersion in porous media is important to a wide range
of applications as diverse as chromatography, ¯ltration, oil recovery, catalysis,
environmental waste management, groundwater °ows, geothermal venting and
processes relevant to animal and plant physiology.
Dispersion is the process whereby molecules that start together in the same
vicinity become separated as a result of translational motions. In thermal equi-
librium, and in the absence of °uid °ow, Brownian motion alone induces molec-
ular separation, although, in a porous medium, the presence of °uid/matrix
interfaces impedes this motion so that apparent di®usion rates depend on
the length and time scales used in making the measurement. Once °ow oc-
curs in a porous medium, other mechanisms for separating initially adjacent
molecules take over and the rate of dispersion rises signi¯cantly above the dif-
fusion \baseline". Like di®usion, dispersion involves stochastic processes that
necessitate the language of statistical physics.
Remarkably, Nuclear Magnetic Resonance turns out to be an e®ective tool in
the measurement of dispersion processes. At the heart of the NMR method is
the use of the spin-echo, a time-reversal process whereby the phases of nuclear
spins are sensitive to random Brownian translations of spin-bearing molecules
when pulses of magnetic ¯eld gradients are applied. Such Pulsed Gradient Spin
Echo (PGSE) methods are not only of equal e®ectiveness in measuring the
enhanced dispersion rate associated with °ow. They also allow, in principle,
2
1 Introduction
The science of °uid dispersion in porous media is important to a wide range
of applications as diverse as chromatography, ¯ltration, oil recovery, catalysis,
environmental waste management, groundwater °ows, geothermal venting and
processes relevant to animal and plant physiology.
Dispersion is the process whereby molecules that start together in the same
vicinity become separated as a result of translational motions. In thermal equi-
librium, and in the absence of °uid °ow, Brownian motion alone induces molec-
ular separation, although, in a porous medium, the presence of °uid/matrix
interfaces impedes this motion so that apparent di®usion rates depend on
the length and time scales used in making the measurement. Once °ow oc-
curs in a porous medium, other mechanisms for separating initially adjacent
molecules take over and the rate of dispersion rises signi¯cantly above the dif-
fusion \baseline". Like di®usion, dispersion involves stochastic processes that
necessitate the language of statistical physics.
Remarkably, Nuclear Magnetic Resonance turns out to be an e®ective tool in
the measurement of dispersion processes. At the heart of the NMR method is
the use of the spin-echo, a time-reversal process whereby the phases of nuclear
spins are sensitive to random Brownian translations of spin-bearing molecules
when pulses of magnetic ¯eld gradients are applied. Such Pulsed Gradient Spin
Echo (PGSE) methods are not only of equal e®ectiveness in measuring the
enhanced dispersion rate associated with °ow. They also allow, in principle,
2
Page 3
a precise and adjustable measurement window over relevant length and time
scales. Note that the PGSE NMR method takes an ensemble average over the
°uid sample, distinguishing it from NMR imaging measurements where the
°ow ¯eld is visualized. The advantage of trading away spatial localization in
such ensemble averaging is a signi¯cant gain in the available length and time
scales over which the motion may be probed.
Pioneering work in the use of ensemble-averaged pulsed gradient spin echo
NMR to study laminar °ow in a pipe was carried out by Hayward, Packer
and Tomlinson in 1972 [8]. In 1996 Lebon et al [13] used PGSE NMR meth-
ods to measure the displacements of molecules in the stochastic °ow as °uid
was forced through a random porous medium and in 1997, Seymour and
Callaghan [17] carried out an extensive PGSE NMR investigation of timescale
dependence of dispersion. These authors made comparisons, with literature
data from other methods, of Peclet number dependence of non-dimensionalized,
asymptotic dispersion coe±cients, measured transverse and longitudinal to the
°ow direction, while subsequent studies by other groups [23,16,14] extended
this work. More recently, the PGSE method has been enhanced to allow inde-
pendently encoding for dispersion at two separated time intervals [11,3]. These
so-called Double PGSE methods allow for the direct measurement of the ve-
locity autocorrelation function. What is missing, in all these approaches, is in-
formation regarding spatial correlations in °ow. What is called for is a PGSE
NMR method that allows the full spatio-temporal structure of dispersive °ow
to be characterized, at the same time preserving information regarding dif-
ferent directional components. Such an NMR method, involving measurement
of the non-local dispersion tensor [12], was ¯rst demonstrated in 2007 [9].
The present paper reports on the ¯rst full measurement of DNL in a model
3
scales. Note that the PGSE NMR method takes an ensemble average over the
°uid sample, distinguishing it from NMR imaging measurements where the
°ow ¯eld is visualized. The advantage of trading away spatial localization in
such ensemble averaging is a signi¯cant gain in the available length and time
scales over which the motion may be probed.
Pioneering work in the use of ensemble-averaged pulsed gradient spin echo
NMR to study laminar °ow in a pipe was carried out by Hayward, Packer
and Tomlinson in 1972 [8]. In 1996 Lebon et al [13] used PGSE NMR meth-
ods to measure the displacements of molecules in the stochastic °ow as °uid
was forced through a random porous medium and in 1997, Seymour and
Callaghan [17] carried out an extensive PGSE NMR investigation of timescale
dependence of dispersion. These authors made comparisons, with literature
data from other methods, of Peclet number dependence of non-dimensionalized,
asymptotic dispersion coe±cients, measured transverse and longitudinal to the
°ow direction, while subsequent studies by other groups [23,16,14] extended
this work. More recently, the PGSE method has been enhanced to allow inde-
pendently encoding for dispersion at two separated time intervals [11,3]. These
so-called Double PGSE methods allow for the direct measurement of the ve-
locity autocorrelation function. What is missing, in all these approaches, is in-
formation regarding spatial correlations in °ow. What is called for is a PGSE
NMR method that allows the full spatio-temporal structure of dispersive °ow
to be characterized, at the same time preserving information regarding dif-
ferent directional components. Such an NMR method, involving measurement
of the non-local dispersion tensor [12], was ¯rst demonstrated in 2007 [9].
The present paper reports on the ¯rst full measurement of DNL in a model
3
Page 4
porous media using a non-invasive method in which all elements of the tensor
have been measured and compared with an independent lattice-Boltzmann
simulation.
Our experimental method is based on nuclear magnetic resonance (NMR),
in which special magnetic ¯eld gradient pulses provide the necessary signal
encoding to track the velocities and spatial displacements of water molecules
labelled by their proton spin states. We demonstrate this method experimen-
tally for dispersive low Reynolds number (Re) °ow in a random bead pack of
mono-sized spheres. A lattice-Boltzmann simulation of the °ow ¯eld through
an independently generated bead-pack has also be performed. This °ow ¯eld
is used to simulate dispersive °ow by allowing virtual tracer particles to °ow
and di®use through the pore space. These tracers can then be used to estimate
the non-local dispersion tensor in the simulated °ow, thus allowing a detailed
comparison with the experimental results.
2 The non-local dispersion tensor
We begin by de¯ning a steady state Eulerian °ow ¯eld vE(r; t) = vE(r) and
stationary Lagrangian °ow ensemble vL(t) with mean °ow hvi. The °uctuating
(zero mean) parts of the velocities, uE(r) and uL(t) are thus de¯ned by
vE(r) = uE(r) + hvi (1)
and
vL(t) = uL(t) + hvi (2)
The asymptotic dispersion tensor, D¤, is described in terms of the velocity
4
have been measured and compared with an independent lattice-Boltzmann
simulation.
Our experimental method is based on nuclear magnetic resonance (NMR),
in which special magnetic ¯eld gradient pulses provide the necessary signal
encoding to track the velocities and spatial displacements of water molecules
labelled by their proton spin states. We demonstrate this method experimen-
tally for dispersive low Reynolds number (Re) °ow in a random bead pack of
mono-sized spheres. A lattice-Boltzmann simulation of the °ow ¯eld through
an independently generated bead-pack has also be performed. This °ow ¯eld
is used to simulate dispersive °ow by allowing virtual tracer particles to °ow
and di®use through the pore space. These tracers can then be used to estimate
the non-local dispersion tensor in the simulated °ow, thus allowing a detailed
comparison with the experimental results.
2 The non-local dispersion tensor
We begin by de¯ning a steady state Eulerian °ow ¯eld vE(r; t) = vE(r) and
stationary Lagrangian °ow ensemble vL(t) with mean °ow hvi. The °uctuating
(zero mean) parts of the velocities, uE(r) and uL(t) are thus de¯ned by
vE(r) = uE(r) + hvi (1)
and
vL(t) = uL(t) + hvi (2)
The asymptotic dispersion tensor, D¤, is described in terms of the velocity
4
Page 5
autocorrelation function (VACF) of the Lagrangian velocities by [12,1]
D¤ = limt!1 sym
Z t
0
d¿huL(0)uL(¿)i (3)
where h: : :i represents the Lagrangian ensemble average. Note that D¤ may
also be de¯ned in Einsteinian terms involving the dyadic of mean square dis-
placements, ¾2(t) by [1]
D¤ = limt!1
1
2
d¾2
dt (4)
The VACF link with the Eulerian ¯eld may be made via a propagator P (rjr0; ¿)
which describes the conditional probability that a °uid element initially at
r will migrate to r0 at a later time ¿ . P (rjr0; ¿) is governed by the mi-
croscale advection-di®usion equation for the system. Given a starting proba-
bility P (r; 0), the velocity autocorrelation function becomes
huL(0)uL(¿)i =
Z
dr0
Z
druE(r; 0)P (r; 0)P (rjr0; ¿)uE(r0; ¿) (5)
Clearly, just as the VACF contains details of the temporal correlations other-
wise buried in the asymptotic dispersion tensor, so the expression
Z
druE(r; 0)P (r)P (rjr0; ¿)uE(r0; ¿) (6)
contains spatial correlation information buried in the VACF. This integral
has been termed by Koch and Brady [12] the non-local dispersion tensor. It
is a primary quantity of interest in any detailed description of dispersive °ow.
Writing this quantity in terms of relative displacements in time and space we
5
D¤ = limt!1 sym
Z t
0
d¿huL(0)uL(¿)i (3)
where h: : :i represents the Lagrangian ensemble average. Note that D¤ may
also be de¯ned in Einsteinian terms involving the dyadic of mean square dis-
placements, ¾2(t) by [1]
D¤ = limt!1
1
2
d¾2
dt (4)
The VACF link with the Eulerian ¯eld may be made via a propagator P (rjr0; ¿)
which describes the conditional probability that a °uid element initially at
r will migrate to r0 at a later time ¿ . P (rjr0; ¿) is governed by the mi-
croscale advection-di®usion equation for the system. Given a starting proba-
bility P (r; 0), the velocity autocorrelation function becomes
huL(0)uL(¿)i =
Z
dr0
Z
druE(r; 0)P (r; 0)P (rjr0; ¿)uE(r0; ¿) (5)
Clearly, just as the VACF contains details of the temporal correlations other-
wise buried in the asymptotic dispersion tensor, so the expression
Z
druE(r; 0)P (r)P (rjr0; ¿)uE(r0; ¿) (6)
contains spatial correlation information buried in the VACF. This integral
has been termed by Koch and Brady [12] the non-local dispersion tensor. It
is a primary quantity of interest in any detailed description of dispersive °ow.
Writing this quantity in terms of relative displacements in time and space we
5
Page 6
have
DNL(R; ¿) =
Z
druE(r; 0)P (r)P (rjr+R; ¿)uE(r+R; ¿) (7)
and
huL(0)uL(¿)i =
Z
dRDNL(R; ¿) (8)
and
D¤ = limt!1 sym
Z t
0
d¿
Z
dRDNL(R; ¿) (9)
Nuclear magnetic resonance (NMR) provides a nearly ideal tracer method
for the measurement of dispersion in that every single molecule is labelled
non-invasively by its local precession frequency in a non-uniform magnetic
¯eld [2,17], the tracer being distributed uniformly within the pore space. Of
course, in principle, magnetic resonance imaging (MRI) allows us to measure
the complete Eulerian velocity ¯eld, to some speci¯ed resolution, limited by
the intrinsic sensitivity of the NMR method. If however an ensemble-averaged
signal from the entire sample is acquired, then resolution in time or with re-
spect to molecular displacements may be optimized [2]. Foremost among such
ensemble-averaged methods for the measurement of molecular displacements
is pulsed gradient spin echo nuclear magnetic resonance (PGSE NMR) [19,2].
This class of experiment not only allows simple analysis of ensemble-averaged
mean-squared displacements over a well-de¯ned time interval [17], but has also
been shown to be e®ective in measuring the VACF for porous media °ow [11].
The tensor DNL(R; ¿) should also be amenable to direct measurement using
PGSE NMR. One earlier such conjecture, by Ding and Candela [6], indicated
the need to encode the NMR signal with information concerning the displace-
ment propagator. However, the measurement of DNL requires that the ex-
periment also be sensitive to velocities separated in space and time. We here
6
DNL(R; ¿) =
Z
druE(r; 0)P (r)P (rjr+R; ¿)uE(r+R; ¿) (7)
and
huL(0)uL(¿)i =
Z
dRDNL(R; ¿) (8)
and
D¤ = limt!1 sym
Z t
0
d¿
Z
dRDNL(R; ¿) (9)
Nuclear magnetic resonance (NMR) provides a nearly ideal tracer method
for the measurement of dispersion in that every single molecule is labelled
non-invasively by its local precession frequency in a non-uniform magnetic
¯eld [2,17], the tracer being distributed uniformly within the pore space. Of
course, in principle, magnetic resonance imaging (MRI) allows us to measure
the complete Eulerian velocity ¯eld, to some speci¯ed resolution, limited by
the intrinsic sensitivity of the NMR method. If however an ensemble-averaged
signal from the entire sample is acquired, then resolution in time or with re-
spect to molecular displacements may be optimized [2]. Foremost among such
ensemble-averaged methods for the measurement of molecular displacements
is pulsed gradient spin echo nuclear magnetic resonance (PGSE NMR) [19,2].
This class of experiment not only allows simple analysis of ensemble-averaged
mean-squared displacements over a well-de¯ned time interval [17], but has also
been shown to be e®ective in measuring the VACF for porous media °ow [11].
The tensor DNL(R; ¿) should also be amenable to direct measurement using
PGSE NMR. One earlier such conjecture, by Ding and Candela [6], indicated
the need to encode the NMR signal with information concerning the displace-
ment propagator. However, the measurement of DNL requires that the ex-
periment also be sensitive to velocities separated in space and time. We here
6
Page 7
demonstrate such a method.
3 NMR implementation
3.1 The pulse sequence
Shown in ¯gure 1 is the pair of two-dimensional pulse sequences needed to
extract the components of the non-local dispersion tensor.
Fig. 1. The `compensated' (a) and `uncompensated' (b) versions of the pulse se-
quence. A superposition is required when post-processing the signal.
7
3 NMR implementation
3.1 The pulse sequence
Shown in ¯gure 1 is the pair of two-dimensional pulse sequences needed to
extract the components of the non-local dispersion tensor.
Fig. 1. The `compensated' (a) and `uncompensated' (b) versions of the pulse se-
quence. A superposition is required when post-processing the signal.
7
Page 8
The analysis and implementation of each dimension is di®erent. As will be
shown below, the analysis in the displacement dimension requires the dis-
placement encoding to be in the propagator sense, so that full attenuation is
achieved and the propagator information can be recovered via fourier trans-
form. The encoding is implemented using a Bi-Polar PGSE [5] with each pulse
having a strength g and a duration of ±. By contrast, in order to get the neces-
sary velocity encoding in the second dimension, the experiments must be per-
formed and analyzed in the low-q limit. This encoding is commonly referred as
a double-PGSE experiment [11]. Each pair of gradients, with a strength of gu
and duration of ± are separated by a time ¢u. Our de¯nition of q will be dif-
ferent for the two dimensions: in the displacement dimension qD = °2±g such
that q is the Fourier conjugate of displacement, and in the velocity encoding
dimension qu is the phase conjugate of velocity such that qu = 2¼°±¢ugu.
In this de¯nition the magnetogyric ratio ° has units of HzT¡1. The two di-
mensions of gradient pulses are generally implemented concurrently with the
same pulse duration, ±, such that when the displacement and velocity encoding
directions are the same, the gradient pulses superpose.
The time ¿ , describes the separation of the double PGSE pulses that provide
the velocity correlation information, in this case ¿ can be thought of to as
a `mixing-time'. While in this experiment ¿ also refers to the single PGSE
encoding for displacement, we will refer to it as an `encoding-time' in that
context.
8
shown below, the analysis in the displacement dimension requires the dis-
placement encoding to be in the propagator sense, so that full attenuation is
achieved and the propagator information can be recovered via fourier trans-
form. The encoding is implemented using a Bi-Polar PGSE [5] with each pulse
having a strength g and a duration of ±. By contrast, in order to get the neces-
sary velocity encoding in the second dimension, the experiments must be per-
formed and analyzed in the low-q limit. This encoding is commonly referred as
a double-PGSE experiment [11]. Each pair of gradients, with a strength of gu
and duration of ± are separated by a time ¢u. Our de¯nition of q will be dif-
ferent for the two dimensions: in the displacement dimension qD = °2±g such
that q is the Fourier conjugate of displacement, and in the velocity encoding
dimension qu is the phase conjugate of velocity such that qu = 2¼°±¢ugu.
In this de¯nition the magnetogyric ratio ° has units of HzT¡1. The two di-
mensions of gradient pulses are generally implemented concurrently with the
same pulse duration, ±, such that when the displacement and velocity encoding
directions are the same, the gradient pulses superpose.
The time ¿ , describes the separation of the double PGSE pulses that provide
the velocity correlation information, in this case ¿ can be thought of to as
a `mixing-time'. While in this experiment ¿ also refers to the single PGSE
encoding for displacement, we will refer to it as an `encoding-time' in that
context.
8
Page 9
3.2 The NMR signal
The normalized signal acquired at the ¯nal spin echo from the sequence shown
in ¯gure 1a may be written
E(qD;qu) =
Z Z
exp(iqD¢R) exp(¡iqu ¢ (uE(r; 0) + hvi))
£P (r)P (rjr+R; ¿) exp(iqu¢(uE(r+R; ¿) + hvi))dRdr
(10)
taking the inverse Fourier transform in the qD dimension, F¡1qDf: : :g, gives us
an expression we call S as a function of displacement R and qu,
S(R;qu) = F¡1qDfE(qD;qu)g =
Z
exp(¡iqu¢uE(r; 0))P (r)P (rjr+R; ¿)
£ exp(iqu¢uE(r+R; ¿))dr
(11)
This integral over starting position R : : : dr can be also be represented by an
ensemble average notation h:::i. Since this quantity involves an average over
all starting positions, the r dependance is dropped. The encoding due to the
second motion encoding pulse occurs when the particle is at a displacement R
from the original position, at a later time ¿ . The displacement information is
implied through the average propagator P (R; ¿). Remembering that the strict
de¯nition is given in equation 11, we can write S as
9
The normalized signal acquired at the ¯nal spin echo from the sequence shown
in ¯gure 1a may be written
E(qD;qu) =
Z Z
exp(iqD¢R) exp(¡iqu ¢ (uE(r; 0) + hvi))
£P (r)P (rjr+R; ¿) exp(iqu¢(uE(r+R; ¿) + hvi))dRdr
(10)
taking the inverse Fourier transform in the qD dimension, F¡1qDf: : :g, gives us
an expression we call S as a function of displacement R and qu,
S(R;qu) = F¡1qDfE(qD;qu)g =
Z
exp(¡iqu¢uE(r; 0))P (r)P (rjr+R; ¿)
£ exp(iqu¢uE(r+R; ¿))dr
(11)
This integral over starting position R : : : dr can be also be represented by an
ensemble average notation h:::i. Since this quantity involves an average over
all starting positions, the r dependance is dropped. The encoding due to the
second motion encoding pulse occurs when the particle is at a displacement R
from the original position, at a later time ¿ . The displacement information is
implied through the average propagator P (R; ¿). Remembering that the strict
de¯nition is given in equation 11, we can write S as
9
Page 10
S(R;qu) =
D
exp(¡iqu¢uE(0))P (R; ¿) exp(iqu¢uE(¿))
E
(12)
As we have changed the notation from explicit de¯nition of the non-local
dispersion given in equation 7, it is worth repeating the de¯nition of DNL
here using the more concise h: : :i notation.
DNL(R; ¿) =
D
uE(0)P (R; ¿)uE(¿)
E
(13)
The compact notation of equation in equation 13 will be used henceforth
although the explicit de¯nition is given in equation 7. The subscripts denoting
Eulerian velocities will be omitted.
3.3 The low-q limit
We now return to the expression for the Fourier transformed normalized echo
attenuation in equation 12. To obtain the desired expression for DNL we need
to expand out the echo attenuation, in the qu dimension, in the low q limit.
In the example of the displacement encoding and velocity encoding in the Z
direction, which we choose to be the direction of the tube velocity, we get
Scomp(Z; qu) =
*Ã
1¡ iquuz(0)¡ q2u
uz(0)2
2 + iO(q
3
u) +O(q4u)
!
P (Z; ¿)
£
Ã
1 + iquuz(¿)¡ q2u
uz(¿)2
2 ¡ iO(q
3
u) +O(q4u)
!+
10
D
exp(¡iqu¢uE(0))P (R; ¿) exp(iqu¢uE(¿))
E
(12)
As we have changed the notation from explicit de¯nition of the non-local
dispersion given in equation 7, it is worth repeating the de¯nition of DNL
here using the more concise h: : :i notation.
DNL(R; ¿) =
D
uE(0)P (R; ¿)uE(¿)
E
(13)
The compact notation of equation in equation 13 will be used henceforth
although the explicit de¯nition is given in equation 7. The subscripts denoting
Eulerian velocities will be omitted.
3.3 The low-q limit
We now return to the expression for the Fourier transformed normalized echo
attenuation in equation 12. To obtain the desired expression for DNL we need
to expand out the echo attenuation, in the qu dimension, in the low q limit.
In the example of the displacement encoding and velocity encoding in the Z
direction, which we choose to be the direction of the tube velocity, we get
Scomp(Z; qu) =
*Ã
1¡ iquuz(0)¡ q2u
uz(0)2
2 + iO(q
3
u) +O(q4u)
!
P (Z; ¿)
£
Ã
1 + iquuz(¿)¡ q2u
uz(¿)2
2 ¡ iO(q
3
u) +O(q4u)
!+
10
Page 11
where qu represents the magnitude of a qu vector applied along the z-axis.
In this version of the pulse sequence the phase shift due to the bulk °ow is
compensated, hence the subscript Scomp. Expanding gives
Scomp(Z; qu) = P (Z; ¿) + iqu
D
uz(0)P (Z; ¿)
E
¡q
2
u
2
D
uz(0)2P (Z; ¿)
E
¡ iqu
D
P (Z; ¿)uz(¿)
E
¡q2u
D
uz(0)P (Z; ¿)u(¿)
E
¡ q
2
u
2
D
P (Z; ¿)uz(¿)2
E
+iO(q3u) +O(q4u)
(14)
The truncation errors in q3u and q4u are notated separately since the real and
imaginary parts of the data are treated separately. The desired encoding is
present in the above as the partial coe±cient of q2u. The velocities in the corre-
lation term in equation 12 are denoted as being separated in time by virtue of
the encoding from the propagator. The other terms, while still containing the
propagator encoding, no longer contain any velocity correlation so for these
terms the apparent time dependence of uz may be dropped.
If a second experiment is performed with the sign of the second velocity en-
coding pulse reversed, in accordance with ¯gure 1b, we get
11
In this version of the pulse sequence the phase shift due to the bulk °ow is
compensated, hence the subscript Scomp. Expanding gives
Scomp(Z; qu) = P (Z; ¿) + iqu
D
uz(0)P (Z; ¿)
E
¡q
2
u
2
D
uz(0)2P (Z; ¿)
E
¡ iqu
D
P (Z; ¿)uz(¿)
E
¡q2u
D
uz(0)P (Z; ¿)u(¿)
E
¡ q
2
u
2
D
P (Z; ¿)uz(¿)2
E
+iO(q3u) +O(q4u)
(14)
The truncation errors in q3u and q4u are notated separately since the real and
imaginary parts of the data are treated separately. The desired encoding is
present in the above as the partial coe±cient of q2u. The velocities in the corre-
lation term in equation 12 are denoted as being separated in time by virtue of
the encoding from the propagator. The other terms, while still containing the
propagator encoding, no longer contain any velocity correlation so for these
terms the apparent time dependence of uz may be dropped.
If a second experiment is performed with the sign of the second velocity en-
coding pulse reversed, in accordance with ¯gure 1b, we get
11
Page 12
Suncomp(Z; qu) = exp (i2qu hvzi)
Ã
P (Z; ¿) + iqu
D
uzP (Z; ¿)
E
¡q
2
u
2
D
u2zP (Z; ¿)
E
+ iqu
D
P (Z; ¿)uz
E
+q2u
D
uz(0)P (Z; ¿)uz(¿)
E
¡ q
2
u
2
D
P (Z; ¿)u2z
E
+iO(q3u) +O(q4u)
!
(15)
The extra phase factor, exp (i2qu hvzi) at the beginning of the expression arises
because this version of the double PGSE experiment does not compensate for
the bulk °ow. This experiment is usually referred to as `uncompensated'. The
phase factor is easily determined and corrected by examining the data E(0; qu).
Performing this phase correction and implementing a di®erence superposition
results in
Scomp(Z; qu)¡ exp(¡i2qu hvzi)Suncomp(Z; quz) = i2qu
D
P (Z; ¿)uz
E
+2q2u
D
uz(0)P (Z; ¿)uz(¿)
E
+iO(q3u) +O(q4u)
(16)
A ¯t to q2u to the real part of the above superposition will give the desired
encoding for the non-local dispersion tensor. Further details will be published
elsewhere.
12
Ã
P (Z; ¿) + iqu
D
uzP (Z; ¿)
E
¡q
2
u
2
D
u2zP (Z; ¿)
E
+ iqu
D
P (Z; ¿)uz
E
+q2u
D
uz(0)P (Z; ¿)uz(¿)
E
¡ q
2
u
2
D
P (Z; ¿)u2z
E
+iO(q3u) +O(q4u)
!
(15)
The extra phase factor, exp (i2qu hvzi) at the beginning of the expression arises
because this version of the double PGSE experiment does not compensate for
the bulk °ow. This experiment is usually referred to as `uncompensated'. The
phase factor is easily determined and corrected by examining the data E(0; qu).
Performing this phase correction and implementing a di®erence superposition
results in
Scomp(Z; qu)¡ exp(¡i2qu hvzi)Suncomp(Z; quz) = i2qu
D
P (Z; ¿)uz
E
+2q2u
D
uz(0)P (Z; ¿)uz(¿)
E
+iO(q3u) +O(q4u)
(16)
A ¯t to q2u to the real part of the above superposition will give the desired
encoding for the non-local dispersion tensor. Further details will be published
elsewhere.
12
Page 13
3.4 Extracting DNL
To inspect the complete non-local tensor we need to consider all possible
encoding directions, and the implications on the superposition. The choices of
directions for the initial and ¯nal velocity encoding need not be parallel. The
expression above can be generalized using the dimension subscripts ® and ¯
for the initial and ¯nal velocity and qu respectively, and °, for the direction
of the displacement encoding. The displacement direction will be denoted by
X with the appropriate subscript. The terms describing the truncation errors
will be omitted. This gives the general expression for the superposition as
exp(iqu® hv®i ¡ iqu¯ hv¯i)Scomp(X°; qu®;¯)
¡ exp(iqu® hv®i+ iqu¯ hv¯i)Suncomp(X°; qu®;¯)g = iqu®
D
u®P (X°; ¿)
E
+iqu¯
D
P (X°; ¿)u°
E
+2qu®qu¯
D
u®(0)P (X°; ¿)u¯(¿)
E
(17)
where qu®;¯ means the ¯rst pair of velocity encoding gradients are along ® and
the second along ¯.
In some cases the velocity encoding is in a direction that has no bulk °ow and
as such, no phase correction is necessary. In contrast if the initial and ¯nal
motion encodings are orthogonal and only one is component is parallel to the
direction of the main °ow, z, the phase correction factor is exp(¡iqu hvi) and
appears in both Scomp(X°; qu®;z) and Sun(X°; qu®;z). Despite the fact the nei-
13
To inspect the complete non-local tensor we need to consider all possible
encoding directions, and the implications on the superposition. The choices of
directions for the initial and ¯nal velocity encoding need not be parallel. The
expression above can be generalized using the dimension subscripts ® and ¯
for the initial and ¯nal velocity and qu respectively, and °, for the direction
of the displacement encoding. The displacement direction will be denoted by
X with the appropriate subscript. The terms describing the truncation errors
will be omitted. This gives the general expression for the superposition as
exp(iqu® hv®i ¡ iqu¯ hv¯i)Scomp(X°; qu®;¯)
¡ exp(iqu® hv®i+ iqu¯ hv¯i)Suncomp(X°; qu®;¯)g = iqu®
D
u®P (X°; ¿)
E
+iqu¯
D
P (X°; ¿)u°
E
+2qu®qu¯
D
u®(0)P (X°; ¿)u¯(¿)
E
(17)
where qu®;¯ means the ¯rst pair of velocity encoding gradients are along ® and
the second along ¯.
In some cases the velocity encoding is in a direction that has no bulk °ow and
as such, no phase correction is necessary. In contrast if the initial and ¯nal
motion encodings are orthogonal and only one is component is parallel to the
direction of the main °ow, z, the phase correction factor is exp(¡iqu hvi) and
appears in both Scomp(X°; qu®;z) and Sun(X°; qu®;z). Despite the fact the nei-
13
Page 14
ther pulse sequence is uncompensated or compensated, the opposite encodings
still provide the appropriate superposition.
A further superposition using Scomp(X°;¡qu®;¯) and Sun(X°;¡qu®;¯) can be
performed to eliminate the q term or alternatively the real and the imaginary
parts of the data can be treated separately. In the latter case, this superposi-
tion also allows the two ensemble averages hu®P (X°)i and hu2®P (X°)i to be
independently extracted. See section 6.6 for more details.
The non-local dispersion tensor, written in this subscript form is
DNL®¯ (X°; ¿) =
D
u®(0)P (X°; ¿)u¯(¿)
E
(18)
It is worth noting that there are some non-local components that can be
measured without the superposition described in equation 17. If the correlation
function hu®(0)u¯(¿)i is known, in advance, to be equal to zero, then
Z
DNL®¯ (X°; ¿)dX° = 0 (19)
independent of the choice of X°. A non-zero DNL®¯ (X°; ¿) must then be an odd
function. For a transverse displacement encoding, the other terms in q2 will
be even functions, thus the data from one pulse sequence can be superposed
with itself, °ipped along the X° direction thus
14
still provide the appropriate superposition.
A further superposition using Scomp(X°;¡qu®;¯) and Sun(X°;¡qu®;¯) can be
performed to eliminate the q term or alternatively the real and the imaginary
parts of the data can be treated separately. In the latter case, this superposi-
tion also allows the two ensemble averages hu®P (X°)i and hu2®P (X°)i to be
independently extracted. See section 6.6 for more details.
The non-local dispersion tensor, written in this subscript form is
DNL®¯ (X°; ¿) =
D
u®(0)P (X°; ¿)u¯(¿)
E
(18)
It is worth noting that there are some non-local components that can be
measured without the superposition described in equation 17. If the correlation
function hu®(0)u¯(¿)i is known, in advance, to be equal to zero, then
Z
DNL®¯ (X°; ¿)dX° = 0 (19)
independent of the choice of X°. A non-zero DNL®¯ (X°; ¿) must then be an odd
function. For a transverse displacement encoding, the other terms in q2 will
be even functions, thus the data from one pulse sequence can be superposed
with itself, °ipped along the X° direction thus
14
Page 15
exp(iqu® hv®i § iqu¯ hv¯i)S(X°; qu®;¯)¡
exp(iqu® hv®i§iqu¯ hv¯i)S(¡X°; qu®;¯)g = iqu®
D
u®P (X°; ¿)
E
+iqu¯
D
P (X°; ¿)u°
E
+ 2qu®qu¯
D
u®(0)P (X°; ¿)u¯(¿)
E
where S results from either the `compensated' or `uncompensated' version of
the pulse sequence. In practice the full superposition is usually used.
The three choices of directions available for the three di®erent encodings give
a total of 27 terms. Hereafter discussions relating to the non-zero and indepen-
dent terms will have ®, ¯ and ° denoted by k, > and ? for encoding parallel to
°ow direction, and the two transverse directions respectively. Inspection of the
data, both experimental and simulated will replace ® and ¯ with the chosen
encoding directions. The displacement encoding direction X° will be denoted
simply by X, Y or Z. In the case of displacement encoding in the direction
parallel to the main °ow, i.e. along the Z axis, it is often more informative
to express the propagators and non-local dispersion terms as a function of
velocity, if this is the case the displacement encoding will be denoted as Vz
with appropriate units.
4 Experimental method
NMR experiments were performed on an Bruker AVANCE 400 MHz spectrom-
eter, with a microimaging gradient set capable of providing 1.45Tm¡1. Dis-
tilled, degassed water was pumped (BVP-Z Ismatec) through a random bead
pack of 500¹m latex monodisperse spheres. A column of randomly packed
15
exp(iqu® hv®i§iqu¯ hv¯i)S(¡X°; qu®;¯)g = iqu®
D
u®P (X°; ¿)
E
+iqu¯
D
P (X°; ¿)u°
E
+ 2qu®qu¯
D
u®(0)P (X°; ¿)u¯(¿)
E
where S results from either the `compensated' or `uncompensated' version of
the pulse sequence. In practice the full superposition is usually used.
The three choices of directions available for the three di®erent encodings give
a total of 27 terms. Hereafter discussions relating to the non-zero and indepen-
dent terms will have ®, ¯ and ° denoted by k, > and ? for encoding parallel to
°ow direction, and the two transverse directions respectively. Inspection of the
data, both experimental and simulated will replace ® and ¯ with the chosen
encoding directions. The displacement encoding direction X° will be denoted
simply by X, Y or Z. In the case of displacement encoding in the direction
parallel to the main °ow, i.e. along the Z axis, it is often more informative
to express the propagators and non-local dispersion terms as a function of
velocity, if this is the case the displacement encoding will be denoted as Vz
with appropriate units.
4 Experimental method
NMR experiments were performed on an Bruker AVANCE 400 MHz spectrom-
eter, with a microimaging gradient set capable of providing 1.45Tm¡1. Dis-
tilled, degassed water was pumped (BVP-Z Ismatec) through a random bead
pack of 500¹m latex monodisperse spheres. A column of randomly packed
15
Page 16
beads was contained in a poly(ether ether ketone) (PEEK) cylinder of in-
ner diameter 10mm. The beadpack was contained by plugs of porous plastic,
outside the sensitive region of the rf coil, to a length of 50mm. A `birdcage'
resonator with a diameter of 15mm was used. The porosity of the beadpack
was around 40%.
In order to dampen out high-frequency pulsations from the pump, a length
of rubber tubing was used in series with the otherwise te°on tubing. Any
remaining variation in the °ow was not detrimental to the signal quality. The
liquid was pumped at a tube velocity of around 10.2mms¡1, giving a ¿v of
49.0ms and a Peclet number of approximately 1600 and a Reynolds number of
approximately 3.3. Experiments reported here are performed at displacement
encoding times of 10ms, 21.5ms and 46.3ms. To reduce the delay between
scans, caused by the need to allow for near-complete T1 relaxation the water
was doped by 0.025wt% with GdCl3 to give a T1 of approximately 110ms.
The time for one measurement is around 8 hours for transverse encodings and
12 hours for longitudinal.
5 Lattice Boltzmann simulation
The simulation code consists of three components, the bead pack generator,
the lattice Boltzmann (LB) °ow generator, and the virtual tracer particle
simulation component. These components are combined in order to build a
suitable beadpack, model the °ow through it, and then model the motions of
virtual tracer particles through that °ow ¯eld.
16
ner diameter 10mm. The beadpack was contained by plugs of porous plastic,
outside the sensitive region of the rf coil, to a length of 50mm. A `birdcage'
resonator with a diameter of 15mm was used. The porosity of the beadpack
was around 40%.
In order to dampen out high-frequency pulsations from the pump, a length
of rubber tubing was used in series with the otherwise te°on tubing. Any
remaining variation in the °ow was not detrimental to the signal quality. The
liquid was pumped at a tube velocity of around 10.2mms¡1, giving a ¿v of
49.0ms and a Peclet number of approximately 1600 and a Reynolds number of
approximately 3.3. Experiments reported here are performed at displacement
encoding times of 10ms, 21.5ms and 46.3ms. To reduce the delay between
scans, caused by the need to allow for near-complete T1 relaxation the water
was doped by 0.025wt% with GdCl3 to give a T1 of approximately 110ms.
The time for one measurement is around 8 hours for transverse encodings and
12 hours for longitudinal.
5 Lattice Boltzmann simulation
The simulation code consists of three components, the bead pack generator,
the lattice Boltzmann (LB) °ow generator, and the virtual tracer particle
simulation component. These components are combined in order to build a
suitable beadpack, model the °ow through it, and then model the motions of
virtual tracer particles through that °ow ¯eld.
16
Page 17
5.1 Beadpack generation
The beadpack is established, bead-by-bead, using a discrete random-walk to
move each sphere's center around within the con¯ning cylinder. To simplify
the construction process, this random walk is constrained so that the sphere
centers always lie upon a node on the underlying grid that will be used during
the LB °uid simulation. Each sphere is ¯rst random-walked along the bottom
of the tube, to generate an initial position in the X-Y plane. The code then
random-walks the sphere upward through the beadpack, until the current
sphere is above all of the preceding spheres. The sphere then random-walks
downward, avoiding any overlaps with neighboring spheres, until the sphere
can no longer be moved downwards or horizontally (i.e. until it reaches a local
minimum).
Once complete, the program outputs a ¯le containing a list of the locations
of the sphere-centers. As this packing algorithm tends to stack the beads
unevenly, the list of sphere-centers is post-processed to remove those for which
of which any part lies above the given height. This allows one to `over¯ll' the
tube, and remove the unevenly-¯lled upper layers to produce a more uniformly
¯lled cylinder of the desired height. Finally, the list of sphere centers is used
to generate a set of lattice boundary conditions suitable for simulating the
°uid °ow via the Lattice Boltzmann method. These boundary conditions are
expressed as a three-dimensional array where each node indicates whether that
point should be considered as °uid or solid, and in the latter case whether the
node is on the surface of the solid or in the bulk. As each sphere has been
constrained to lie on a grid-node, the discretization process can use exactly
the same representation for each bead.
17
The beadpack is established, bead-by-bead, using a discrete random-walk to
move each sphere's center around within the con¯ning cylinder. To simplify
the construction process, this random walk is constrained so that the sphere
centers always lie upon a node on the underlying grid that will be used during
the LB °uid simulation. Each sphere is ¯rst random-walked along the bottom
of the tube, to generate an initial position in the X-Y plane. The code then
random-walks the sphere upward through the beadpack, until the current
sphere is above all of the preceding spheres. The sphere then random-walks
downward, avoiding any overlaps with neighboring spheres, until the sphere
can no longer be moved downwards or horizontally (i.e. until it reaches a local
minimum).
Once complete, the program outputs a ¯le containing a list of the locations
of the sphere-centers. As this packing algorithm tends to stack the beads
unevenly, the list of sphere-centers is post-processed to remove those for which
of which any part lies above the given height. This allows one to `over¯ll' the
tube, and remove the unevenly-¯lled upper layers to produce a more uniformly
¯lled cylinder of the desired height. Finally, the list of sphere centers is used
to generate a set of lattice boundary conditions suitable for simulating the
°uid °ow via the Lattice Boltzmann method. These boundary conditions are
expressed as a three-dimensional array where each node indicates whether that
point should be considered as °uid or solid, and in the latter case whether the
node is on the surface of the solid or in the bulk. As each sphere has been
constrained to lie on a grid-node, the discretization process can use exactly
the same representation for each bead.
17
Page 18
Note that the radius of the spheres is reduced by 1 lattice spacing, a, before
mapping them into the cubic space. This is done because the LBGK algorithm
does not resolve the hydrodynamics correctly in channels that are less than 4
lattice units wide [20], and therefore reducing the bead radius helps to ensure
that the °ow simulation is stable. This has a signi¯cant a®ect on the overall
porosity of the discretized beadpack, which for coarse grids can be as high as
0:7, much larger than the » 0:4 expected for random bead packs. However,
the porosity of the discrete system is calculated by counting numbers of solid
and °uid lattice sites. This is likely to be something of an under-estimate on
such a coarse mesh, as we would ideally use the e®ective hydrodynamic radius
of the particles which is di±cult to estimate directly. having explored a range
of discretization algorithms, we have found this radius-reduction method to
be a simple way of constructing a computable system while still producing the
overall °ow and dispersion behavior characteristic of this type of system [10].
5.2 The °ow generation
The Lattice-Boltzmann code simulates the °uid °ow around the beadpack.
It takes the boundary conditions described above and simulates the °ow of
a Lattice-Boltzmann °uid through that system. The code uses the Lattice
Bhatnagar-Gross-Krook (LBGK) algorithm on a D3Q19 lattice, with 2nd-
order mid-grid bounce-back boundary conditions. A Poiseuille bulk force is
applied along the °ow axis, along with periodic boundary conditions, in order
to generate the desired °ow. The magnitude of the bulk force was chosen so
as to induce a low Reynolds number °ow of about Re = 0:05.
The code is written in C++ can be run in serial (one machine), or on any
18
mapping them into the cubic space. This is done because the LBGK algorithm
does not resolve the hydrodynamics correctly in channels that are less than 4
lattice units wide [20], and therefore reducing the bead radius helps to ensure
that the °ow simulation is stable. This has a signi¯cant a®ect on the overall
porosity of the discretized beadpack, which for coarse grids can be as high as
0:7, much larger than the » 0:4 expected for random bead packs. However,
the porosity of the discrete system is calculated by counting numbers of solid
and °uid lattice sites. This is likely to be something of an under-estimate on
such a coarse mesh, as we would ideally use the e®ective hydrodynamic radius
of the particles which is di±cult to estimate directly. having explored a range
of discretization algorithms, we have found this radius-reduction method to
be a simple way of constructing a computable system while still producing the
overall °ow and dispersion behavior characteristic of this type of system [10].
5.2 The °ow generation
The Lattice-Boltzmann code simulates the °uid °ow around the beadpack.
It takes the boundary conditions described above and simulates the °ow of
a Lattice-Boltzmann °uid through that system. The code uses the Lattice
Bhatnagar-Gross-Krook (LBGK) algorithm on a D3Q19 lattice, with 2nd-
order mid-grid bounce-back boundary conditions. A Poiseuille bulk force is
applied along the °ow axis, along with periodic boundary conditions, in order
to generate the desired °ow. The magnitude of the bulk force was chosen so
as to induce a low Reynolds number °ow of about Re = 0:05.
The code is written in C++ can be run in serial (one machine), or on any
18
Page 19
parallel machine that supports MPI (Message Passing Interface). As a rough
outline of performance, we note that simulations on grids of around 200£200£
400 lattice sites can be performed on a single machine in a few hours, whereas
simulations at around 400£400£1000 will take a few tens of hours on a parallel
machine. The code outputs the velocity ¯eld as three binary ¯les corresponding
to the velocity at each point in the system in each of the [x, y, z] directions.
Note also that the code can be compiled to use either single-precision (32-bit)
or double-precision (64-bit) °oating point numbers for storing the discrete-
velocity populations and °ow velocities. We use single-precision numbers to
save on memory, but must occasionally use double-precision to ensure that
using low-precision storage is not damaging our results 1 . A number of test
systems have been run at both levels of precision, and no signi¯cant variation
is results has been discovered [10].
5.3 Tracer particles and calculation of DNL
The virtual tracer code takes the velocity ¯eld determined using LBGK, and
simulates the movements of tracer particles through that velocity-¯eld. Tracer
motion is both mechanical (due to the °ow-¯eld, following streamlines) and
di®usive (crossing streamlines). The code takes the Peclet number Pe as a
parameter, and determines the di®usion constant D0 for the simulation from
that via Pe = dÁhvzi(1¡Á)D0 where d is the bead diameter, Á is the porosity and
hvzi is the mean rate of °ow in the z-direction. The measurement system is
structured like the experimental set up, following exactly each time interval
1 Note, however, that all integration over the whole system, e.g. the calculation of
overall °ow rate, are accumulated using double-precision variables.
19
outline of performance, we note that simulations on grids of around 200£200£
400 lattice sites can be performed on a single machine in a few hours, whereas
simulations at around 400£400£1000 will take a few tens of hours on a parallel
machine. The code outputs the velocity ¯eld as three binary ¯les corresponding
to the velocity at each point in the system in each of the [x, y, z] directions.
Note also that the code can be compiled to use either single-precision (32-bit)
or double-precision (64-bit) °oating point numbers for storing the discrete-
velocity populations and °ow velocities. We use single-precision numbers to
save on memory, but must occasionally use double-precision to ensure that
using low-precision storage is not damaging our results 1 . A number of test
systems have been run at both levels of precision, and no signi¯cant variation
is results has been discovered [10].
5.3 Tracer particles and calculation of DNL
The virtual tracer code takes the velocity ¯eld determined using LBGK, and
simulates the movements of tracer particles through that velocity-¯eld. Tracer
motion is both mechanical (due to the °ow-¯eld, following streamlines) and
di®usive (crossing streamlines). The code takes the Peclet number Pe as a
parameter, and determines the di®usion constant D0 for the simulation from
that via Pe = dÁhvzi(1¡Á)D0 where d is the bead diameter, Á is the porosity and
hvzi is the mean rate of °ow in the z-direction. The measurement system is
structured like the experimental set up, following exactly each time interval
1 Note, however, that all integration over the whole system, e.g. the calculation of
overall °ow rate, are accumulated using double-precision variables.
19
Page 20
of the NMR experiments.
The simulation moves a set of Nt independent tracers through the system
simultaneously. At every step, the simulation advances by a ¯xed timestep
¢t, during which the position of each particle is updated. This time step is
chosen so as to ensure that the °ow ¯eld is always fully explored, by ensuring
that the corresponding spatial step size is no more than 0.1a. We use the
Heun (mid-step) predictor corrector algorithm, which is known to be weakly
second-order convergent for our system [21]. Thus, for each spatial coordinate
(X) of every tracer particle, the update is determined using:
X(t+¢t) = X(t) + 1=2fv(X(t)) + v(Xp(t+¢t))g¢t+
q
2D0¢W (t) (20)
Where v(X) is the point velocity determined from the LB °ow ¯eld, ¢W (t)
is a Gaussian random variable with variance ¢t, and the predicted position
Xp determined using the Euler method:
X(t+¢t) = X(t) + v(X(t))¢t+
q
2D0¢W (t) (21)
The point velocity is calculated from the LB °ow ¯eld using simple linear
interpolation over the cubic cell containing the tracer position. As for the
LB simulation, periodic boundary conditions are imposed along the °ow axis.
Any tracer update that would cause the particle to cross the discretized bead
surfaces and leave the pore space is rejected. This approach has been used to
calculate the non-dimensionalized dispersion coe±cient as a function of Peclet
number and agrees well with the published data over the range Pe = 10¡1 to
20
The simulation moves a set of Nt independent tracers through the system
simultaneously. At every step, the simulation advances by a ¯xed timestep
¢t, during which the position of each particle is updated. This time step is
chosen so as to ensure that the °ow ¯eld is always fully explored, by ensuring
that the corresponding spatial step size is no more than 0.1a. We use the
Heun (mid-step) predictor corrector algorithm, which is known to be weakly
second-order convergent for our system [21]. Thus, for each spatial coordinate
(X) of every tracer particle, the update is determined using:
X(t+¢t) = X(t) + 1=2fv(X(t)) + v(Xp(t+¢t))g¢t+
q
2D0¢W (t) (20)
Where v(X) is the point velocity determined from the LB °ow ¯eld, ¢W (t)
is a Gaussian random variable with variance ¢t, and the predicted position
Xp determined using the Euler method:
X(t+¢t) = X(t) + v(X(t))¢t+
q
2D0¢W (t) (21)
The point velocity is calculated from the LB °ow ¯eld using simple linear
interpolation over the cubic cell containing the tracer position. As for the
LB simulation, periodic boundary conditions are imposed along the °ow axis.
Any tracer update that would cause the particle to cross the discretized bead
surfaces and leave the pore space is rejected. This approach has been used to
calculate the non-dimensionalized dispersion coe±cient as a function of Peclet
number and agrees well with the published data over the range Pe = 10¡1 to
20
Page 21
Pe = 105 [10].
At prede¯ned times in the tracer simulation, the ¯nal position and velocity
of each tracer particle is recorded. This enables each particles contribution to
the VACF to be calculated, for example the jth particle u®(0)ju¯(¿)j. This
contribution is then added to a histogram bin dictated by the displacement X i°.
The ensemble average of these particles gives the desired quantity DNL®¯ (X°; ¿).
Further tensors, discussed in section 6.6, are similarly calculated.
5.4 Dimensionality of the non-local dispersion tensor
The statistically calculated non-local dispersion tensor allows for convenient
examination of the non-zero and independent terms. Each encoding can be
either parallel or transverse to the main °ow, however, for a full survey we
will treat the two orthogonal encodings independently. We use the subscripts
>, ? and k to describe the two transverse directions and the parallel direc-
tion, respectively. Terms are equal when all of the transverse subscripts are
swapped, that is
D
u?(0)P (X?; ¿)u?(¿)
E
6=
D
u?(0)P (X>; ¿)u?(¿)
E
(22)
but
D
u?(0)P (X?; ¿)u?(¿)
E
=
D
u>(0)P (X>; ¿)u>(¿)
E
(23)
21
At prede¯ned times in the tracer simulation, the ¯nal position and velocity
of each tracer particle is recorded. This enables each particles contribution to
the VACF to be calculated, for example the jth particle u®(0)ju¯(¿)j. This
contribution is then added to a histogram bin dictated by the displacement X i°.
The ensemble average of these particles gives the desired quantity DNL®¯ (X°; ¿).
Further tensors, discussed in section 6.6, are similarly calculated.
5.4 Dimensionality of the non-local dispersion tensor
The statistically calculated non-local dispersion tensor allows for convenient
examination of the non-zero and independent terms. Each encoding can be
either parallel or transverse to the main °ow, however, for a full survey we
will treat the two orthogonal encodings independently. We use the subscripts
>, ? and k to describe the two transverse directions and the parallel direc-
tion, respectively. Terms are equal when all of the transverse subscripts are
swapped, that is
D
u?(0)P (X?; ¿)u?(¿)
E
6=
D
u?(0)P (X>; ¿)u?(¿)
E
(22)
but
D
u?(0)P (X?; ¿)u?(¿)
E
=
D
u>(0)P (X>; ¿)u>(¿)
E
(23)
21
Page 22
5.4.1 Longitudinal displacement encoding
When the displacement encoding is parallel to the °ow, all the o®-axis terms
are zero, the non zero terms being DNL??(Xk; ¿) and DNLkk (Xk; ¿). The latter
component perhaps shows the most structure at long times, which is indicative
of the persistence of the parallel velocity auto-correlation function.
5.4.2 Transverse displacement encoding
The diagonal transverse displacement encodings show di®erent behavior, and
one o®-axis element is present, giving four independent componentsDNL??(X?; ¿),
DNL>>(X?; ¿), DNLkk (X?; ¿) and the o®-axis component DNL?k (X?; ¿). The o®-
axis component describes the correlations between transverse and longitu-
dinal velocities. The correlation term itself is zero but when resolved by a
transverse displacement aligned with the transverse velocity direction, as in
the term DNL?k (X?; ¿), some structure is seen. This is not seen when resolved
in the other transverse direction DNL?k (X>; ¿), or in the longitudinal direction
DNL?k (X?; ¿).
The three terms describing transverse correlations show di®erent structure
but when integrated over displacement will give the same results, namely the
transverse velocity auto-correlation function hu?(0)u?(¿)i.
22
When the displacement encoding is parallel to the °ow, all the o®-axis terms
are zero, the non zero terms being DNL??(Xk; ¿) and DNLkk (Xk; ¿). The latter
component perhaps shows the most structure at long times, which is indicative
of the persistence of the parallel velocity auto-correlation function.
5.4.2 Transverse displacement encoding
The diagonal transverse displacement encodings show di®erent behavior, and
one o®-axis element is present, giving four independent componentsDNL??(X?; ¿),
DNL>>(X?; ¿), DNLkk (X?; ¿) and the o®-axis component DNL?k (X?; ¿). The o®-
axis component describes the correlations between transverse and longitu-
dinal velocities. The correlation term itself is zero but when resolved by a
transverse displacement aligned with the transverse velocity direction, as in
the term DNL?k (X?; ¿), some structure is seen. This is not seen when resolved
in the other transverse direction DNL?k (X>; ¿), or in the longitudinal direction
DNL?k (X?; ¿).
The three terms describing transverse correlations show di®erent structure
but when integrated over displacement will give the same results, namely the
transverse velocity auto-correlation function hu?(0)u?(¿)i.
22
Page 23
6 Experiment and simulation results
6.1 De¯nition of Peclet number
The random beadpacks of the experiment and simulations di®er. First the
ratio of bead diameter to cylinder is 20 to 1 experimentally and 10 to 1 in
simulations. Second, the porosity of the experimental beadpack is around 0.40
whereas the simulated beadpack has a porosity of around 0.71. This has im-
plications for the matching of Peclet number and the choice of encoding time.
The common experimental de¯nition of Peclet number is given as [17]
Pe =
à Á
1¡ Á
! hvi dbead
D0 (24)
Where Á is the porosity, D0 is the self di®usion coe±cient and hvi is the
interstitial velocity. This includes a de¯nition of the `characteristic length'
as l = (Á=(1 ¡ Á))dbead. Because of the di®erence in porosity between the
experimental and simulated beadpacks there is a marked di®erence in the
characteristic length, when expressed as a fraction of the bead diameter. For
the experimental beadpack we have lexp = 0:66dbead and in simulation, lsim =
2:4dbead. Note that a characteristic time ¿v may be de¯ned, corresponding to
the time to °ow a bead diameter
¿v = dbeadhvi : (25)
In this case, to enable meaningful comparisons between experiment and simu-
lation it may be more natural to look at encoding times expressed in terms of
23
6.1 De¯nition of Peclet number
The random beadpacks of the experiment and simulations di®er. First the
ratio of bead diameter to cylinder is 20 to 1 experimentally and 10 to 1 in
simulations. Second, the porosity of the experimental beadpack is around 0.40
whereas the simulated beadpack has a porosity of around 0.71. This has im-
plications for the matching of Peclet number and the choice of encoding time.
The common experimental de¯nition of Peclet number is given as [17]
Pe =
à Á
1¡ Á
! hvi dbead
D0 (24)
Where Á is the porosity, D0 is the self di®usion coe±cient and hvi is the
interstitial velocity. This includes a de¯nition of the `characteristic length'
as l = (Á=(1 ¡ Á))dbead. Because of the di®erence in porosity between the
experimental and simulated beadpacks there is a marked di®erence in the
characteristic length, when expressed as a fraction of the bead diameter. For
the experimental beadpack we have lexp = 0:66dbead and in simulation, lsim =
2:4dbead. Note that a characteristic time ¿v may be de¯ned, corresponding to
the time to °ow a bead diameter
¿v = dbeadhvi : (25)
In this case, to enable meaningful comparisons between experiment and simu-
lation it may be more natural to look at encoding times expressed in terms of
23
Page 24
the time ¿l to °ow a characteristic length. For completeness, both comparisons
will be made.
6.2 Propagator measurement
The 6 independent components and propagators were measured at three dis-
placement encoding times 10ms, 21.5ms and 46.3ms, corresponding to the
approximate logarithmic distribution of 0.20¿v, 0.44¿v and 0:94¿v. In terms of
the time to °ow a characteristic length the encoding times are 0.31¿l, 0.66¿l
and 1.43¿l. It is clear that the mixing is well within the pre-asymptotic regime
where it is expected that most structure will be shown in the non-local dis-
persion tensor [18]. Once we have these 6 measurements we can begin to infer
characteristics of the °ow.
It is clear from the measured and simulated longitudinal propagators in ¯g-
ure 2 that there are signi¯cant di®erences between the two systems. Note the
use of a displacement axis for the transverse direction and a velocity axis for
the longitudinal. The use of the velocity axis makes for easier comparison over
widely di®erent encoding times where longitudinal displacements may sub-
stantially di®er. The simulations shows a large number of particles at slow,
or zero °ow and shows little mixing over a time of ¿ = ¿v. In contrast the
measured data show a Gaussian tendency at ¿ = ¿v. The simulated propaga-
tors when encoding times are expressed as a function on ¿l are more favorable,
¯gure 3, and show the distribution moving towards gaussian for times » 1:5¿l.
Nevertheless, the simulations show a large subset of particles with zero °ow.
Other groups [15] have had success in using a lattice Boltzmann simulation
on a lattice generated directly from a MRI image of a real bead-pack, using a
24
will be made.
6.2 Propagator measurement
The 6 independent components and propagators were measured at three dis-
placement encoding times 10ms, 21.5ms and 46.3ms, corresponding to the
approximate logarithmic distribution of 0.20¿v, 0.44¿v and 0:94¿v. In terms of
the time to °ow a characteristic length the encoding times are 0.31¿l, 0.66¿l
and 1.43¿l. It is clear that the mixing is well within the pre-asymptotic regime
where it is expected that most structure will be shown in the non-local dis-
persion tensor [18]. Once we have these 6 measurements we can begin to infer
characteristics of the °ow.
It is clear from the measured and simulated longitudinal propagators in ¯g-
ure 2 that there are signi¯cant di®erences between the two systems. Note the
use of a displacement axis for the transverse direction and a velocity axis for
the longitudinal. The use of the velocity axis makes for easier comparison over
widely di®erent encoding times where longitudinal displacements may sub-
stantially di®er. The simulations shows a large number of particles at slow,
or zero °ow and shows little mixing over a time of ¿ = ¿v. In contrast the
measured data show a Gaussian tendency at ¿ = ¿v. The simulated propaga-
tors when encoding times are expressed as a function on ¿l are more favorable,
¯gure 3, and show the distribution moving towards gaussian for times » 1:5¿l.
Nevertheless, the simulations show a large subset of particles with zero °ow.
Other groups [15] have had success in using a lattice Boltzmann simulation
on a lattice generated directly from a MRI image of a real bead-pack, using a
24
Page 25
Fig. 2. (Color online) The simulated (a) and measured (b) propagators. The sim-
ulated longitudinal propagators clearly show a larger proportion of slower moving
particles than the measured. The greater mixing shown in the experimental mea-
surements over simulation highlights the di®erence between the two systems.
monte-carlo simulation to produce agreement with measured propagators. It
is not the goal of this work to be able to fully predict NMR results with a sim-
ulation technique, merely to develop two independent methods for inspecting
the non-local dispersion tensor. Hence the preferred method of independently
generating the simulated beadpack.
25
ulated longitudinal propagators clearly show a larger proportion of slower moving
particles than the measured. The greater mixing shown in the experimental mea-
surements over simulation highlights the di®erence between the two systems.
monte-carlo simulation to produce agreement with measured propagators. It
is not the goal of this work to be able to fully predict NMR results with a sim-
ulation technique, merely to develop two independent methods for inspecting
the non-local dispersion tensor. Hence the preferred method of independently
generating the simulated beadpack.
25
Page 26
−0.4 −0.2 0 0.2 0.4
0
5
10
15
X (mm)
P
(X
,τ
)(
m
m
−
1
)
0 20 40 60
0
0.02
0.04
0.06
0.08
0.1
0.12
Vz (mms−1)
P
(V
z,
τ)
(m
m
−
1
s) 0.37
0.80
1.72
τ/τl
Fig. 3. (Color online) Simulated propagators with encoding times in terms of ¿l.
When expressed in fractions of the time to °ow a characteristic length, the experi-
mental propagators have encoding times of a similar fraction.
−0.4 −0.2 0 0.2 0.4
0
200
400
600
800
1000
X (mm)
D
N
L
zz
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−40
−20
0
20
40
60
X (mm)
D
N
L
xx
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
0
20
40
60
80
100
X (mm)
D
N
L
yy
(X
,τ
)(
m
m
s−
2
)
0 20 40 60
0
2
4
6
8
10
Vz (mms−1)
D
N
L
zz
(V
z,
τ)
(m
m
s−
1
)
0 20 40 60
−0.5
0
0.5
1
Vz (mms−1)
D
N
L
xx
(V
z,
τ)
(m
m
s−
1
)
−0.4 −0.2 0 0.2 0.4
−50
0
50
X (mm)
D
N
L
xz
(X
,τ
)(
m
m
s−
2
)
0.19
0.42
0.90
τ/τv
Fig. 4. (Color online) The six independent non-local dispersion tensor components
of for °ow in a porous medium, simulated. See text for discussion. The compo-
nents in the direction of the main °ow are plotted as average velocity, rather than
displacement.
26
0
5
10
15
X (mm)
P
(X
,τ
)(
m
m
−
1
)
0 20 40 60
0
0.02
0.04
0.06
0.08
0.1
0.12
Vz (mms−1)
P
(V
z,
τ)
(m
m
−
1
s) 0.37
0.80
1.72
τ/τl
Fig. 3. (Color online) Simulated propagators with encoding times in terms of ¿l.
When expressed in fractions of the time to °ow a characteristic length, the experi-
mental propagators have encoding times of a similar fraction.
−0.4 −0.2 0 0.2 0.4
0
200
400
600
800
1000
X (mm)
D
N
L
zz
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−40
−20
0
20
40
60
X (mm)
D
N
L
xx
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
0
20
40
60
80
100
X (mm)
D
N
L
yy
(X
,τ
)(
m
m
s−
2
)
0 20 40 60
0
2
4
6
8
10
Vz (mms−1)
D
N
L
zz
(V
z,
τ)
(m
m
s−
1
)
0 20 40 60
−0.5
0
0.5
1
Vz (mms−1)
D
N
L
xx
(V
z,
τ)
(m
m
s−
1
)
−0.4 −0.2 0 0.2 0.4
−50
0
50
X (mm)
D
N
L
xz
(X
,τ
)(
m
m
s−
2
)
0.19
0.42
0.90
τ/τv
Fig. 4. (Color online) The six independent non-local dispersion tensor components
of for °ow in a porous medium, simulated. See text for discussion. The compo-
nents in the direction of the main °ow are plotted as average velocity, rather than
displacement.
26
Page 27
−0.4 −0.2 0 0.2 0.4
0
200
400
600
800
1000
X (mm)
D
N
L
zz
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−40
−20
0
20
40
60
X (mm)
D
N
L
xx
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
0
20
40
60
80
100
X (mm)
D
N
L
yy
(X
,τ
)(
m
m
s−
2
)
0 20 40 60
0
2
4
6
8
10
Vz (mms−1)
D
N
L
zz
(V
z,
τ)
(m
m
s−
1
)
0 20 40 60
−0.5
0
0.5
1
Vz (mms−1)
D
N
L
xx
(V
z,
τ)
(m
m
s−
1
)
−0.4 −0.2 0 0.2 0.4
−50
0
50
X (mm)
D
N
L
xz
(X
,τ
)(
m
m
s−
2
)
0.20
0.44
0.94
τ/τv
Fig. 5. (Color online) The six independent non-local dispersion tensor components
for °ow in a porous medium, measured by NMR. See text for discussion. The com-
ponents in the direction of the main °ow are plotted as average velocity, rather than
displacement. The error bars show the estimate of error in the ¯t to q2.
6.3 DNL components from simulation and experiments
The non-local dispersion components clearly show a good agreement for en-
coding times matched via the time to °ow a bead (¯gures 4 and 5) in the cases
when transverse correlations are expressed. The agreement between the terms
describing longitudinal correlations is acceptable. Interestingly the greatest
di®erence in structure is highlighted by the cross-correlation term, DNLxz (X; ¿).
Here the strong measured correlation describing particles moving a small dis-
tance in ¡X, is not present in simulation.
While the agreement between the simulated and experimental propagators was
27
0
200
400
600
800
1000
X (mm)
D
N
L
zz
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−40
−20
0
20
40
60
X (mm)
D
N
L
xx
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
0
20
40
60
80
100
X (mm)
D
N
L
yy
(X
,τ
)(
m
m
s−
2
)
0 20 40 60
0
2
4
6
8
10
Vz (mms−1)
D
N
L
zz
(V
z,
τ)
(m
m
s−
1
)
0 20 40 60
−0.5
0
0.5
1
Vz (mms−1)
D
N
L
xx
(V
z,
τ)
(m
m
s−
1
)
−0.4 −0.2 0 0.2 0.4
−50
0
50
X (mm)
D
N
L
xz
(X
,τ
)(
m
m
s−
2
)
0.20
0.44
0.94
τ/τv
Fig. 5. (Color online) The six independent non-local dispersion tensor components
for °ow in a porous medium, measured by NMR. See text for discussion. The com-
ponents in the direction of the main °ow are plotted as average velocity, rather than
displacement. The error bars show the estimate of error in the ¯t to q2.
6.3 DNL components from simulation and experiments
The non-local dispersion components clearly show a good agreement for en-
coding times matched via the time to °ow a bead (¯gures 4 and 5) in the cases
when transverse correlations are expressed. The agreement between the terms
describing longitudinal correlations is acceptable. Interestingly the greatest
di®erence in structure is highlighted by the cross-correlation term, DNLxz (X; ¿).
Here the strong measured correlation describing particles moving a small dis-
tance in ¡X, is not present in simulation.
While the agreement between the simulated and experimental propagators was
27
Page 28
−0.4 −0.2 0 0.2 0.4
0
200
400
600
800
1000
X (mm)
D
N
L
zz
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−40
−20
0
20
40
60
X (mm)
D
N
L
xx
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
0
20
40
60
80
100
X (mm)
D
N
L
yy
(X
,τ
)(
m
m
s−
2
)
0 20 40 60
0
2
4
6
8
10
Vz (mms−1)
D
N
L
zz
(V
z,
τ)
(m
m
s−
1
)
0 20 40 60
−0.5
0
0.5
1
Vz (mms−1)
D
N
L
xx
(V
z,
τ)
(m
m
s−
1
)
−0.4 −0.2 0 0.2 0.4
−5
0
5
X (mm)
D
N
L
xz
(X
,τ
)(
m
m
s−
2
)
0.37
0.80
1.72
τ/τl
Fig. 6. (Color online) The six independent non-local dispersion tensor components
of for °ow in a porous medium, simulated at times to match the experimental
measurement times when expressed as ¿l. With the exception of the component
DNLxz (X; ¿) the components are plotted on the same scale as ¯gure 5. See text for
discussion. The components in the direction of the main °ow are plotted as average
velocity, rather than displacement.
best with encoding times expressed as ¿l, the agreement of the non-local disper-
sion components are mixed (¯gures 5 and 6). The agreement between compo-
nents describing parallel velocity correlations is very good, but the transverse
terms are poor. In this case the measured structure of the cross-correlation
term DNLxz (X; ¿) is more clearly seen in the simulations, albeit with a di®erent
magnitude.
The clear di®erences are not seen as a failing of either method but perhaps
shows the capacity of non-local dispersion measurements to highlight subtle
28
0
200
400
600
800
1000
X (mm)
D
N
L
zz
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−40
−20
0
20
40
60
X (mm)
D
N
L
xx
(X
,τ
)(
m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
0
20
40
60
80
100
X (mm)
D
N
L
yy
(X
,τ
)(
m
m
s−
2
)
0 20 40 60
0
2
4
6
8
10
Vz (mms−1)
D
N
L
zz
(V
z,
τ)
(m
m
s−
1
)
0 20 40 60
−0.5
0
0.5
1
Vz (mms−1)
D
N
L
xx
(V
z,
τ)
(m
m
s−
1
)
−0.4 −0.2 0 0.2 0.4
−5
0
5
X (mm)
D
N
L
xz
(X
,τ
)(
m
m
s−
2
)
0.37
0.80
1.72
τ/τl
Fig. 6. (Color online) The six independent non-local dispersion tensor components
of for °ow in a porous medium, simulated at times to match the experimental
measurement times when expressed as ¿l. With the exception of the component
DNLxz (X; ¿) the components are plotted on the same scale as ¯gure 5. See text for
discussion. The components in the direction of the main °ow are plotted as average
velocity, rather than displacement.
best with encoding times expressed as ¿l, the agreement of the non-local disper-
sion components are mixed (¯gures 5 and 6). The agreement between compo-
nents describing parallel velocity correlations is very good, but the transverse
terms are poor. In this case the measured structure of the cross-correlation
term DNLxz (X; ¿) is more clearly seen in the simulations, albeit with a di®erent
magnitude.
The clear di®erences are not seen as a failing of either method but perhaps
shows the capacity of non-local dispersion measurements to highlight subtle
28
Page 29
di®erences in the behavior of systems where particles have moved over short
times or short distances. Further discussions will involve encoding times ex-
pressed in the time to °ow on bead diameter ¿v.
6.4 Analysis of longitudinal displacement encodings
The non-local dispersion tensor contains correlation information that is sup-
pressed in the VACF. Considering the longitudinal VACF huz(0)uz(¿)i, we see
that in ¯gure 5a and d we have two ways of describing the correlations, ¯rst
displaced by X and second by Z, or as shown, Vz. The component DNLzz (Vz; ¿)
shows three distinct regions, a correlation, decorrelation and a further corre-
lation. That at lowest Vz highlights the slow and stationary part as having a
very highly correlated uz. Since these particles are moving slowly with respect
to the mean °ow, they will have a large, negative, uz which will be relatively
unchanged at short ¿ . The unusual feature of this peak is that while this
peak will contribute quite strongly to the VACF, indicating a high degree of
correlation, the °ow in the laboratory frame is very small. The second cor-
relation peak arises from particles that have moved a considerable distance
in the encoding time, and yet retained a high degree of correlation. While
the longitudinal propagator shown in ¯gure 2(a) does not show a signi¯cant
fraction of particles moving that distance, the non-local dispersion component
DNLzz (Vz; ¿) clearly reveals this and shows its contribution to the VACF.
It is much more e±cient to investigate the full position and time dependance
in simulation. As an example the non-local component DNLzz (Vz; ¿) is shown
for times up to 20 ¿v in ¯gure 7. The time evolution of the two correlation
regions and the decorrelation region can be seen. Furthermore, these three
29
times or short distances. Further discussions will involve encoding times ex-
pressed in the time to °ow on bead diameter ¿v.
6.4 Analysis of longitudinal displacement encodings
The non-local dispersion tensor contains correlation information that is sup-
pressed in the VACF. Considering the longitudinal VACF huz(0)uz(¿)i, we see
that in ¯gure 5a and d we have two ways of describing the correlations, ¯rst
displaced by X and second by Z, or as shown, Vz. The component DNLzz (Vz; ¿)
shows three distinct regions, a correlation, decorrelation and a further corre-
lation. That at lowest Vz highlights the slow and stationary part as having a
very highly correlated uz. Since these particles are moving slowly with respect
to the mean °ow, they will have a large, negative, uz which will be relatively
unchanged at short ¿ . The unusual feature of this peak is that while this
peak will contribute quite strongly to the VACF, indicating a high degree of
correlation, the °ow in the laboratory frame is very small. The second cor-
relation peak arises from particles that have moved a considerable distance
in the encoding time, and yet retained a high degree of correlation. While
the longitudinal propagator shown in ¯gure 2(a) does not show a signi¯cant
fraction of particles moving that distance, the non-local dispersion component
DNLzz (Vz; ¿) clearly reveals this and shows its contribution to the VACF.
It is much more e±cient to investigate the full position and time dependance
in simulation. As an example the non-local component DNLzz (Vz; ¿) is shown
for times up to 20 ¿v in ¯gure 7. The time evolution of the two correlation
regions and the decorrelation region can be seen. Furthermore, these three
29
Page 31
0 5 10 15 20
0
0.2
0.4
0.6
0.8
1
τ/τv
〈u
(0
)u
(τ
)〉
(d
2 be
ad
τ−
2
v
)
1st correlation
decorrelation
2nd correlation
VACF
Fig. 9. (Color online) The contributions from the three regions as de¯ned in ¯gure 8.
At very short ¿ , i.e. ¿ < 0:1¿v, there is no decorrelation and hence no distinct regions
exist.
regions, shown for clarity in ¯gure 8 can be independently integrated and
monitored as a function of ¿ . Figure 9 shows the contribution from each region
to the VACF . It is interesting to note that the two correlation regions show
di®erent characteristic times, the faster correlation decay of the second region
is perhaps due to the mechanical mixing of the °owing particles. In contrast
the particles in the 1st region, which are stationary or moving slowly, show
a longer timescale, perhaps dominated by Taylor or hold-up dispersion. The
mean position in dimensionless units of the three regions are also shown in
¯gure 10 to a longer encoding time. The maximum contribution from the
decorrelation occurs at around 2:5¿v.
6.5 Analysis of transverse displacement encodings
The non-local dispersion tensor components DNLxx (X; ¿) and DNLzz (Z; ¿) de-
scribe the same correlation function when integrated over displacement. We
can therefore inspect the two components and infer the behavior of the parti-
31
0
0.2
0.4
0.6
0.8
1
τ/τv
〈u
(0
)u
(τ
)〉
(d
2 be
ad
τ−
2
v
)
1st correlation
decorrelation
2nd correlation
VACF
Fig. 9. (Color online) The contributions from the three regions as de¯ned in ¯gure 8.
At very short ¿ , i.e. ¿ < 0:1¿v, there is no decorrelation and hence no distinct regions
exist.
regions, shown for clarity in ¯gure 8 can be independently integrated and
monitored as a function of ¿ . Figure 9 shows the contribution from each region
to the VACF . It is interesting to note that the two correlation regions show
di®erent characteristic times, the faster correlation decay of the second region
is perhaps due to the mechanical mixing of the °owing particles. In contrast
the particles in the 1st region, which are stationary or moving slowly, show
a longer timescale, perhaps dominated by Taylor or hold-up dispersion. The
mean position in dimensionless units of the three regions are also shown in
¯gure 10 to a longer encoding time. The maximum contribution from the
decorrelation occurs at around 2:5¿v.
6.5 Analysis of transverse displacement encodings
The non-local dispersion tensor components DNLxx (X; ¿) and DNLzz (Z; ¿) de-
scribe the same correlation function when integrated over displacement. We
can therefore inspect the two components and infer the behavior of the parti-
31
Page 32
0 5 10 15 20
0
0.5
1
1.5
2
2.5
3
3.5
4
τ/τv
Co
rr
ela
tio
n
po
sit
io
n
(d
be
ad
τ−
1
v
)
1st correlation
decorrelation
2nd correlation
Fig. 10. (Color online) The movement of the mean of each of the three regions as
de¯ned in ¯gure 8. The decorrelation lobe is relatively unchanged but the 1st corre-
lation region is slowly populated with slightly faster moving particles as ¿ increases,
while the 2nd correlation is rapidly populated with slower moving particles.
cles in X and Z that contribute to that correlation. With reference to ¯gure 11,
the term DNLxx (X; ¿) at an encoding time of 0:44¿v shows a strong decorrelation
for particles that have not moved very far, if at all, in the X direction (re-
gion A). The region of DNLxx (Z; ¿) labelled D show a decorrelation for particles
that have moved a larger distance than the average. We might surmise that
the subset of particles that are decorrelated correspond to particles similar to
that labelled i that move round a bead as shown in ¯gure 11. The particles
with a large correlation have moved a short distance in Z. In the case of region
C we might guess that these particles are on a slow moving stream-line j, and
as a result have moved in the same direction as the ux component of velocity,
region B.
To con¯rm the above, a simultaneous encoding in X and Z would be neces-
sary, this would give a two-dimensional non-local dispersion tensor. The extra
displacement dimension can easily be added to the NMR experiment with
time as the only cost. In simulation the `binning' as discussed in section 5.3
32
0
0.5
1
1.5
2
2.5
3
3.5
4
τ/τv
Co
rr
ela
tio
n
po
sit
io
n
(d
be
ad
τ−
1
v
)
1st correlation
decorrelation
2nd correlation
Fig. 10. (Color online) The movement of the mean of each of the three regions as
de¯ned in ¯gure 8. The decorrelation lobe is relatively unchanged but the 1st corre-
lation region is slowly populated with slightly faster moving particles as ¿ increases,
while the 2nd correlation is rapidly populated with slower moving particles.
cles in X and Z that contribute to that correlation. With reference to ¯gure 11,
the term DNLxx (X; ¿) at an encoding time of 0:44¿v shows a strong decorrelation
for particles that have not moved very far, if at all, in the X direction (re-
gion A). The region of DNLxx (Z; ¿) labelled D show a decorrelation for particles
that have moved a larger distance than the average. We might surmise that
the subset of particles that are decorrelated correspond to particles similar to
that labelled i that move round a bead as shown in ¯gure 11. The particles
with a large correlation have moved a short distance in Z. In the case of region
C we might guess that these particles are on a slow moving stream-line j, and
as a result have moved in the same direction as the ux component of velocity,
region B.
To con¯rm the above, a simultaneous encoding in X and Z would be neces-
sary, this would give a two-dimensional non-local dispersion tensor. The extra
displacement dimension can easily be added to the NMR experiment with
time as the only cost. In simulation the `binning' as discussed in section 5.3
32
Page 33
Fig. 11. (Color online) Inferring the behavior of particles moving round a bead
with experimental measurements. By inspecting the DNLxx (X; ¿) and DNLxx (X; ¿) at
an encoding time such that the mean displacement is approximately half a bead
diameter. A guess is made that particles similar to i make up that parts of the
graphs A and D, whereas the correlations shown as B and C are from particles
similar to j. See text for complete discussion.
can include a second dimension. A large number of tracer particles, » 107, is
required in order to produce statistically reliable estimates. Indeed the exper-
imental data (¯gure 12) shows three distinct regions, the particles that move
simultaneously in X and Z showing a strong correlation in ux. The particles
that show a strong decortication occupy a broad region and are clearly not lim-
ited to those which have moved a small X. This two-dimensional measurement
is consistent with our interpretation but perhaps highlights a little naivety. In
33
with experimental measurements. By inspecting the DNLxx (X; ¿) and DNLxx (X; ¿) at
an encoding time such that the mean displacement is approximately half a bead
diameter. A guess is made that particles similar to i make up that parts of the
graphs A and D, whereas the correlations shown as B and C are from particles
similar to j. See text for complete discussion.
can include a second dimension. A large number of tracer particles, » 107, is
required in order to produce statistically reliable estimates. Indeed the exper-
imental data (¯gure 12) shows three distinct regions, the particles that move
simultaneously in X and Z showing a strong correlation in ux. The particles
that show a strong decortication occupy a broad region and are clearly not lim-
ited to those which have moved a small X. This two-dimensional measurement
is consistent with our interpretation but perhaps highlights a little naivety. In
33
Page 35
6.6 Other tensors
The di®erence superposition and ¯t to q2 in equation 16 is designed to extract
the correlations that yield the non-local dispersion tensor. An additive super-
position and ¯t to q2 will give a term that can be interpreted as the second
moment of u resolved by displacement. This term is potentially interesting
as it gives some structure to the second moment, and because this structure
shows temporal information.
Additionally, a ¯t in q will give the ¯rst moment of u resolved by displacement.
This term is interesting since the ¯rst moment is, by de¯nition, zero but with
displacement resolution we have some displacement and temporal structure.
Figure 13 shows the compete set of 8 extra tensors available from the same
experiments that produced the 6 non-local elements in ¯gure 5. The extra
tensors do not contain any information correlating velocities in space and time
and the superposition from the two pulse sequences used here is not required to
measure them. We show these tensors simply because the encoding necessary
for non-local dispersion allows their extraction.
In experiments where the motion encoding is parallel to the main °ow, the
necessary phase correction due to hvi can be omitted, thus giving tensors
not in u but v. It is assumed that tensors in u give more structure, since
the distinction between velocities above and below the mean are clear, but
nevertheless, these quantities are available for examination. In total, from an
experiment at a single encoding time ¿ , we can obtain the quantities shown
in table 1, where ®, ° and ¯ refer to the chosen direction of initial motion
encoding, displacement encoding and ¯nal motion encoding respectively, while
35
The di®erence superposition and ¯t to q2 in equation 16 is designed to extract
the correlations that yield the non-local dispersion tensor. An additive super-
position and ¯t to q2 will give a term that can be interpreted as the second
moment of u resolved by displacement. This term is potentially interesting
as it gives some structure to the second moment, and because this structure
shows temporal information.
Additionally, a ¯t in q will give the ¯rst moment of u resolved by displacement.
This term is interesting since the ¯rst moment is, by de¯nition, zero but with
displacement resolution we have some displacement and temporal structure.
Figure 13 shows the compete set of 8 extra tensors available from the same
experiments that produced the 6 non-local elements in ¯gure 5. The extra
tensors do not contain any information correlating velocities in space and time
and the superposition from the two pulse sequences used here is not required to
measure them. We show these tensors simply because the encoding necessary
for non-local dispersion allows their extraction.
In experiments where the motion encoding is parallel to the main °ow, the
necessary phase correction due to hvi can be omitted, thus giving tensors
not in u but v. It is assumed that tensors in u give more structure, since
the distinction between velocities above and below the mean are clear, but
nevertheless, these quantities are available for examination. In total, from an
experiment at a single encoding time ¿ , we can obtain the quantities shown
in table 1, where ®, ° and ¯ refer to the chosen direction of initial motion
encoding, displacement encoding and ¯nal motion encoding respectively, while
35
Page 36
−0.4 −0.2 0 0.2 0.4
−20
−10
0
10
20
X (mm)
〈
u x
P
(X
,τ
)〉
(s
−
1
)
−0.4 −0.2 0 0.2 0.4
−40
−30
−20
−10
0
10
20
X (mm)
〈
u z
P
(X
,τ
)〉
(s
−
1
)
0 20 40 60
−0.6
−0.4
−0.2
0
0.2
Vz (mms−1)
〈
u z
P
(V
z
,τ
)〉
−0.4 −0.2 0 0.2 0.4
−100
0
100
200
300
X (mm)
〈
u2 y
P
(X
,τ
)〉
(m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−50
0
50
100
150
X (mm)
〈
u2 x
P
(X
,τ
)〉
(m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−500
0
500
1000
X (mm)
〈
u2 z
P
(X
,τ
)〉
(m
m
s−
2
)
0 20 40 60
−2
0
2
4
6
Vz (mms−1)
〈
u2 z
P
(V
z
,τ
)〉
(m
m
s−
1
)
0 20 40 60
−1
0
1
2
3
Vz (mms−1)
〈
u2 x
P
(V
z
,τ
)〉
(m
m
s−
1
)
0.20
0.44
0.94
τ/τv
Fig. 13. (Color online) The 8 tensors also available from the experimental encoding
necessary for DNL as measured by NMR. The error bars show the estimate of error
in the ¯t to q or q2.
by integrating over the displacement we obtain the quantities shown in table 2.
Of course some of these quantities will be zero or dependent on others. Inspec-
tion shows that, for °ow in a bead pack we have 21 independent, non-zero,
displacement resolved quantities, including 6 non-local elements. In addition,
two propagators, two VACF measurements and the moments hvzi hu2zi, hu2xi,
hv2zi and hv2xi. All of this information can be obtained from 6 experiments.
36
−20
−10
0
10
20
X (mm)
〈
u x
P
(X
,τ
)〉
(s
−
1
)
−0.4 −0.2 0 0.2 0.4
−40
−30
−20
−10
0
10
20
X (mm)
〈
u z
P
(X
,τ
)〉
(s
−
1
)
0 20 40 60
−0.6
−0.4
−0.2
0
0.2
Vz (mms−1)
〈
u z
P
(V
z
,τ
)〉
−0.4 −0.2 0 0.2 0.4
−100
0
100
200
300
X (mm)
〈
u2 y
P
(X
,τ
)〉
(m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−50
0
50
100
150
X (mm)
〈
u2 x
P
(X
,τ
)〉
(m
m
s−
2
)
−0.4 −0.2 0 0.2 0.4
−500
0
500
1000
X (mm)
〈
u2 z
P
(X
,τ
)〉
(m
m
s−
2
)
0 20 40 60
−2
0
2
4
6
Vz (mms−1)
〈
u2 z
P
(V
z
,τ
)〉
(m
m
s−
1
)
0 20 40 60
−1
0
1
2
3
Vz (mms−1)
〈
u2 x
P
(V
z
,τ
)〉
(m
m
s−
1
)
0.20
0.44
0.94
τ/τv
Fig. 13. (Color online) The 8 tensors also available from the experimental encoding
necessary for DNL as measured by NMR. The error bars show the estimate of error
in the ¯t to q or q2.
by integrating over the displacement we obtain the quantities shown in table 2.
Of course some of these quantities will be zero or dependent on others. Inspec-
tion shows that, for °ow in a bead pack we have 21 independent, non-zero,
displacement resolved quantities, including 6 non-local elements. In addition,
two propagators, two VACF measurements and the moments hvzi hu2zi, hu2xi,
hv2zi and hv2xi. All of this information can be obtained from 6 experiments.
36
Page 37
Description Quantity
Propagator P (X° ; ¿)
®
Non-local dispersion tensor DNL®¯ (X° ; ¿)
First moment of velocity u®P (X° ; ¿)
® ; P (X° ; ¿)u¯
®
Second moment of velocity 1=2
³u2®P (X° ; ¿)
®+
D
P (X° ; ¿)u2¯
E´
Non-local dispersion in v v®(0)P (X° ; ¿)v¯(¿)
®
First moment in v v®P (X° ; ¿)
® ; P (X° ; ¿)v¯
®
Second moment in v 1=2
³v2®P (X° ; ¿)
®+
D
P (X° ; ¿)v2¯
E´
Table 1
A table showing all of the displacement resolved quantities available from a single
experiment with the initial motion encoding in the direction ®, ¯nal motion encoding
in the direction ¯ and displacement encoding in the direction °. Not all of these
quantities will be non-zero.
6.7 Higher dimensionality
The component DNL?k (X?; ¿) clearly shows structure that is otherwise hid-
den since the correlation function,
D
u?(0)uk(¿)
E
is equal to zero. This same
correlation is described in the DNL component DNL?k (X?; ¿) but this is ob-
served to be zero. Since some structure is present in the former quantity, per-
haps a 2-dimensional displacement encoding could provide more information,
indeed the component
D
u?(0)P (X?; Xk; ¿)uk(¿)
E
shows this structure. Fur-
ther investigation reveals 8 independent 2-dimensional components, namely
DNL??(X>; X?; ¿), DNLkk (X>; X?; ¿), DNL>?(X>; X?; ¿) and DNLk? (X>; X?; ¿) for
37
Propagator P (X° ; ¿)
®
Non-local dispersion tensor DNL®¯ (X° ; ¿)
First moment of velocity u®P (X° ; ¿)
® ; P (X° ; ¿)u¯
®
Second moment of velocity 1=2
³u2®P (X° ; ¿)
®+
D
P (X° ; ¿)u2¯
E´
Non-local dispersion in v v®(0)P (X° ; ¿)v¯(¿)
®
First moment in v v®P (X° ; ¿)
® ; P (X° ; ¿)v¯
®
Second moment in v 1=2
³v2®P (X° ; ¿)
®+
D
P (X° ; ¿)v2¯
E´
Table 1
A table showing all of the displacement resolved quantities available from a single
experiment with the initial motion encoding in the direction ®, ¯nal motion encoding
in the direction ¯ and displacement encoding in the direction °. Not all of these
quantities will be non-zero.
6.7 Higher dimensionality
The component DNL?k (X?; ¿) clearly shows structure that is otherwise hid-
den since the correlation function,
D
u?(0)uk(¿)
E
is equal to zero. This same
correlation is described in the DNL component DNL?k (X?; ¿) but this is ob-
served to be zero. Since some structure is present in the former quantity, per-
haps a 2-dimensional displacement encoding could provide more information,
indeed the component
D
u?(0)P (X?; Xk; ¿)uk(¿)
E
shows this structure. Fur-
ther investigation reveals 8 independent 2-dimensional components, namely
DNL??(X>; X?; ¿), DNLkk (X>; X?; ¿), DNL>?(X>; X?; ¿) and DNLk? (X>; X?; ¿) for
37
Page 38
Description Quantity
VACF hu®(0)u¯(¿)i
Second moment of velocity 1=2
³u2®
®+
D
u2¯
E´
First moment of velocity hu®i ; hu¯i
VACF in v hv®(0)v¯(¿)i
First moment in v hv®i ; hv¯i
Second moment in v 1=2
³v2®
®+
D
v2¯
E´
Table 2
A table showing the remaining quantities available from a single experiment with
the initial motion encoding in the direction ®, ¯nal motion encoding in the direction
¯. Not all of these quantities will be non-zero.
fully transverse displacements andDNL??(X?; Xk; ¿),DNL>>(X?; Xk; ¿),DNLkk (X?; Xk; ¿)
and DNLk? (X?; Xk; ¿) which include a parallel displacement encoding. The eight
components and the 2 two-dimensional propagators obtained by simulation are
shown in ¯gures 14 and 15. The enhanced structure in the non-local dispersion
is clear.
The quantity DNL>?(X>; X?; ¿) is interesting since none of the DNL showed any
structure for orthogonal transverse correlations. The symmetry ofDNL>?(X>; X?; ¿)
shows that the integrals over each dimension will indeed give zero. Also of note
is that a one-dimensional representation of the transverse orthogonal correla-
tions can be seen if the displacement encoding directed at an angle between the
two transverse directions, this quantity can be called
D
u>(0)P (X45; ¿)u?(¿)
E
and is non-zero. It is somewhat di±cult to interpret correlations where trans-
38
VACF hu®(0)u¯(¿)i
Second moment of velocity 1=2
³u2®
®+
D
u2¯
E´
First moment of velocity hu®i ; hu¯i
VACF in v hv®(0)v¯(¿)i
First moment in v hv®i ; hv¯i
Second moment in v 1=2
³v2®
®+
D
v2¯
E´
Table 2
A table showing the remaining quantities available from a single experiment with
the initial motion encoding in the direction ®, ¯nal motion encoding in the direction
¯. Not all of these quantities will be non-zero.
fully transverse displacements andDNL??(X?; Xk; ¿),DNL>>(X?; Xk; ¿),DNLkk (X?; Xk; ¿)
and DNLk? (X?; Xk; ¿) which include a parallel displacement encoding. The eight
components and the 2 two-dimensional propagators obtained by simulation are
shown in ¯gures 14 and 15. The enhanced structure in the non-local dispersion
is clear.
The quantity DNL>?(X>; X?; ¿) is interesting since none of the DNL showed any
structure for orthogonal transverse correlations. The symmetry ofDNL>?(X>; X?; ¿)
shows that the integrals over each dimension will indeed give zero. Also of note
is that a one-dimensional representation of the transverse orthogonal correla-
tions can be seen if the displacement encoding directed at an angle between the
two transverse directions, this quantity can be called
D
u>(0)P (X45; ¿)u?(¿)
E
and is non-zero. It is somewhat di±cult to interpret correlations where trans-
38
Page 39
Non-local dispersion Tensor DNLkk (Xk; ¿),DNL??(Xk; ¿),DNLkk (X?; ¿)
DNL??(X?; ¿), DNL>>(X?; ¿), DNL?k (X?; ¿)
First moment of velocity u?P (X?; ¿)
®, ukP (X?; ¿)
® , ukP (Xk; ¿)
®
Second moment of velocity u2>P (X?; ¿)
®, u2?P (X?; ¿)
®,
D
u2kP (X?; ¿)
E
D
u2kP (Xk; ¿)
E
, u2?P (Xk; ¿)
®
Non-local dispersion in v vk(0)P (X?; ¿)vk(¿)
®, vk(0)P (Xk; ¿)vk(¿)
®
v?(0)P (X?; ¿)vk(¿)
®
First moment in v vkP (X?; ¿)
®, vkP (Xk; ¿)
®
Second moment in v
D
v2kP (X?; ¿)
E
,
D
v2kP (Xk; ¿)
E
Table 3
A table showing the independent, non-zero, displacement-resolved tensors available
for the non-local dispersion measurements.
verse velocity encodings are used since it is unknown whether the positive
correlation results from particles moving from a region where the velocity is
in one direction to a region where the velocity is in the second direction or
the particle has continued to moved in a way where the velocity is non-zero in
both directions. Nevertheless, these quantities are available to be measured.
A natural extension to the two-dimensional encoding would be to add the
third dimension. Indeed the three terms containing transverse correlations all
show di®erent structure and a three-dimensional term describing this correla-
tions with three dimensions of displacement encoding would be illuminating.
However, the considerable time needed for either, the experimental data to
be acquired, or to perform a simulation long enough to get reliable statistics,
39
DNL??(X?; ¿), DNL>>(X?; ¿), DNL?k (X?; ¿)
First moment of velocity u?P (X?; ¿)
®, ukP (X?; ¿)
® , ukP (Xk; ¿)
®
Second moment of velocity u2>P (X?; ¿)
®, u2?P (X?; ¿)
®,
D
u2kP (X?; ¿)
E
D
u2kP (Xk; ¿)
E
, u2?P (Xk; ¿)
®
Non-local dispersion in v vk(0)P (X?; ¿)vk(¿)
®, vk(0)P (Xk; ¿)vk(¿)
®
v?(0)P (X?; ¿)vk(¿)
®
First moment in v vkP (X?; ¿)
®, vkP (Xk; ¿)
®
Second moment in v
D
v2kP (X?; ¿)
E
,
D
v2kP (Xk; ¿)
E
Table 3
A table showing the independent, non-zero, displacement-resolved tensors available
for the non-local dispersion measurements.
verse velocity encodings are used since it is unknown whether the positive
correlation results from particles moving from a region where the velocity is
in one direction to a region where the velocity is in the second direction or
the particle has continued to moved in a way where the velocity is non-zero in
both directions. Nevertheless, these quantities are available to be measured.
A natural extension to the two-dimensional encoding would be to add the
third dimension. Indeed the three terms containing transverse correlations all
show di®erent structure and a three-dimensional term describing this correla-
tions with three dimensions of displacement encoding would be illuminating.
However, the considerable time needed for either, the experimental data to
be acquired, or to perform a simulation long enough to get reliable statistics,
39
Page 40
X (dbead)
Y
(d
be
ad
)
P (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLxx (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLxy (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLzz (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLzx (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
Fig. 14. (Color online) The transverse two dimensional propagator and non-local
dispersion components obtained by simulation for an encoding time of 0.09 ¿v. The
propagator and the term DNLzz (X;Y; ¿) are plotted with a saturated color map to
show some of the structure. Lighter colors correspond to a positive correlation,
darker colors correspond to a decorrelation.
means such a measurement has not been attempted.
Without further exploration we note that when a two-dimensional experiment
is performed, the qu and q2u ¯ts and superposition would give rise to other
tensors describing the ¯rst and second moments, as in the one-dimensional
case.
40
Y
(d
be
ad
)
P (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLxx (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLxy (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLzz (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
X (dbead)
Y
(d
be
ad
)
DNLzx (X, Y, τ )
−0.5 0 0.5
−0.5
0
0.5
Fig. 14. (Color online) The transverse two dimensional propagator and non-local
dispersion components obtained by simulation for an encoding time of 0.09 ¿v. The
propagator and the term DNLzz (X;Y; ¿) are plotted with a saturated color map to
show some of the structure. Lighter colors correspond to a positive correlation,
darker colors correspond to a decorrelation.
means such a measurement has not been attempted.
Without further exploration we note that when a two-dimensional experiment
is performed, the qu and q2u ¯ts and superposition would give rise to other
tensors describing the ¯rst and second moments, as in the one-dimensional
case.
40
Page 41
X (dbead)
V z
(d
be
ad
τ−
1
v
)
P (X, Vz, τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLxx (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLyy (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLzz (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLzx (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
Fig. 15. (Color online) The transverse-longitudinal two dimensional propagator and
non-local dispersion components for an encoding time of 0.09 ¿v. The propagator
and the term DNLzz (X;Vz; ¿) are plotted with a saturated color map to show some
of the structure. Lighter colors correspond to a positive correlation, darker colors
correspond to a decorrelation.
7 Conclusion
We have demonstrated the ability to measure the full non-local dispersion ten-
sor for dispersive °ow in a model porous medium. The six non-zero indepen-
dent terms show remarkable structure that can be used to infer characteristics
about the °ow. Through our NMR measurements we have also been able to
41
V z
(d
be
ad
τ−
1
v
)
P (X, Vz, τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLxx (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLyy (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLzz (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
X (dbead)
V z
(d
be
ad
τ−
1
v
)
DNLzx (X, Vz , τ )
−0.5 0 0.5
0
1
2
3
4
5
Fig. 15. (Color online) The transverse-longitudinal two dimensional propagator and
non-local dispersion components for an encoding time of 0.09 ¿v. The propagator
and the term DNLzz (X;Vz; ¿) are plotted with a saturated color map to show some
of the structure. Lighter colors correspond to a positive correlation, darker colors
correspond to a decorrelation.
7 Conclusion
We have demonstrated the ability to measure the full non-local dispersion ten-
sor for dispersive °ow in a model porous medium. The six non-zero indepen-
dent terms show remarkable structure that can be used to infer characteristics
about the °ow. Through our NMR measurements we have also been able to
41
Page 42
extract a range of other tensors add to another displacement dimension to the
measurements
Simulations of the non-local dispersion tensor show broad agreement with
experiment and will provide a useful tool for further investigations. Simula-
tions naturally enable a more e±cient means of investigating, over a range of
encoding times, the large-dimensional space involving multiple cartesian di-
rections. The direct comparison between the simulation and experiment can
be improved in three ways: by matching the Reynolds number of the °ows,
by more accurately resolving the ¯ne structure of the porous media, and by
producing more precise estimates of the various observables. However, such
simulations are computationally expensive due to the very large lattices and
a large number of tracers required to achieve this.
A remarkable feature of both experiment and simulation is the highly detailed
structure of DNL, as well as the signi¯cant di®erences between the various
tensor components. A simple picture involving °ow around a single bead gives
¯rst order insight. But the structural details available suggest that by mea-
suring DNL we gain a highly sensitive probe of the mechanism of dispersion.
Already, performing experiments on porous rocks show signi¯cantly di®erent
behavior. We anticipate the DNL measurements will provide a wealth of in-
formation.
References
[1] Brenner, H. Dispersion resulting from °ow through spatially periodic porous-
media. PHILOSOPHICAL TRANSACTIONS OF THE ROYAL SOCIETY OF
LONDON SERIES A-MATHEMATICAL PHYSICAL AND ENGINEERING
42
measurements
Simulations of the non-local dispersion tensor show broad agreement with
experiment and will provide a useful tool for further investigations. Simula-
tions naturally enable a more e±cient means of investigating, over a range of
encoding times, the large-dimensional space involving multiple cartesian di-
rections. The direct comparison between the simulation and experiment can
be improved in three ways: by matching the Reynolds number of the °ows,
by more accurately resolving the ¯ne structure of the porous media, and by
producing more precise estimates of the various observables. However, such
simulations are computationally expensive due to the very large lattices and
a large number of tracers required to achieve this.
A remarkable feature of both experiment and simulation is the highly detailed
structure of DNL, as well as the signi¯cant di®erences between the various
tensor components. A simple picture involving °ow around a single bead gives
¯rst order insight. But the structural details available suggest that by mea-
suring DNL we gain a highly sensitive probe of the mechanism of dispersion.
Already, performing experiments on porous rocks show signi¯cantly di®erent
behavior. We anticipate the DNL measurements will provide a wealth of in-
formation.
References
[1] Brenner, H. Dispersion resulting from °ow through spatially periodic porous-
media. PHILOSOPHICAL TRANSACTIONS OF THE ROYAL SOCIETY OF
LONDON SERIES A-MATHEMATICAL PHYSICAL AND ENGINEERING
42
Page 43
SCIENCES 297, 1430 (1980), 81{133.
[2] Callaghan, P. Principles of Nuclear Magnetic Resonance Microscopy. Oxford
University Press, 1991.
[3] Callaghan, P., and Codd, S. Flow coherence in a bead pack observed using
frequency domain modulated gradient nuclear magnetic resonance. PHYSICS
OF FLUIDS 13, 2 (FEB 2001), 421{427.
[4] Codd, S., Manz, B., Seymour, J., and Callaghan, P. Taylor dispersion
and molecular displacements in Poiseuille °ow. PHYSICAL REVIEW E 60, 4,
Part A (OCT 1999), R3491{R3494.
[5] Cotts, R., Hoch, M., Sun, T., and Markert, J. Pulsed ¯eld gradient
stimulated echo methods for improved NMR di®usion measurements in
heterogeneous systems. JOURNAL OF MAGNETIC RESONANCE 83, 2 (JUN
15 1989), 252{266.
[6] Ding, A., and Candela, D. Probing nonlocal tracer dispersion in °ows
through random porous media. PHYSICAL REVIEW E 54, 1 (JUL 1996),
656{660.
[7] Goossens, M., Mittelbach, F., and Samarin, A. The LATEX Companion.
Addison-Wesley, Reading, Massachusetts, USA, 1994.
[8] Hayward, R., Tomlinso, D., and Packer, K. Pulsed ¯eld-gradient spin-
echo NMR studies of °ow in °uids. MOLECULAR PHYSICS 23, 6 (1972),
1083{&.
[9] Hunter, M. W., and Callaghan, P. T. NMR measurement of nonlocal
dispersion in complex °ows. PHYSICAL REVIEW LETTERS 99, 21 (NOV 23
2007).
[10] Jackson, A. N. and Hunter, M. W., and Callaghan, P. T. Simulating
dispersion in porous media. To be published.
43
[2] Callaghan, P. Principles of Nuclear Magnetic Resonance Microscopy. Oxford
University Press, 1991.
[3] Callaghan, P., and Codd, S. Flow coherence in a bead pack observed using
frequency domain modulated gradient nuclear magnetic resonance. PHYSICS
OF FLUIDS 13, 2 (FEB 2001), 421{427.
[4] Codd, S., Manz, B., Seymour, J., and Callaghan, P. Taylor dispersion
and molecular displacements in Poiseuille °ow. PHYSICAL REVIEW E 60, 4,
Part A (OCT 1999), R3491{R3494.
[5] Cotts, R., Hoch, M., Sun, T., and Markert, J. Pulsed ¯eld gradient
stimulated echo methods for improved NMR di®usion measurements in
heterogeneous systems. JOURNAL OF MAGNETIC RESONANCE 83, 2 (JUN
15 1989), 252{266.
[6] Ding, A., and Candela, D. Probing nonlocal tracer dispersion in °ows
through random porous media. PHYSICAL REVIEW E 54, 1 (JUL 1996),
656{660.
[7] Goossens, M., Mittelbach, F., and Samarin, A. The LATEX Companion.
Addison-Wesley, Reading, Massachusetts, USA, 1994.
[8] Hayward, R., Tomlinso, D., and Packer, K. Pulsed ¯eld-gradient spin-
echo NMR studies of °ow in °uids. MOLECULAR PHYSICS 23, 6 (1972),
1083{&.
[9] Hunter, M. W., and Callaghan, P. T. NMR measurement of nonlocal
dispersion in complex °ows. PHYSICAL REVIEW LETTERS 99, 21 (NOV 23
2007).
[10] Jackson, A. N. and Hunter, M. W., and Callaghan, P. T. Simulating
dispersion in porous media. To be published.
43
Page 44
[11] Khrapitchev, A., and Callaghan, P. Reversible and irreversible dispersion
in a porous medium. PHYSICS OF FLUIDS 15, 9 (SEP 2003), 2649{2660.
[12] Koch, D., and Brady, J. A nonlocal description of advection di®usion
with application to dispersion in porous-media. JOURNAL OF FLUID
MECHANICS 180 (JUL 1987), 387{403.
[13] Lebon, L., Leblond, J., Hulin, J., Martys, N., and Schwartz, L. Pulsed
¯eld gradient NMR measurements of probability distribution of displacement
under °ow in sphere packings. MAGNETIC RESONANCE IMAGING 14,
7-8 (1996), 989{991. 3rd International Meeting on Recent Advances in MR
Applications to Porous Media, LOUVAIN, BELGIUM, SEP 03-06, 1995.
[14] Manz, B., Alexander, P., and Gladden, L. Correlations between
dispersion and structure in porous media probed by nuclear magnetic resonance.
PHYSICS OF FLUIDS 11, 2 (FEB 1999), 259{267.
[15] Manz, B., Gladden, L., and Warren, P. Flow and dispersion in porous
media: Lattice-Boltzmann and NMR studies. AICHE JOURNAL 45, 9 (SEP
1999), 1845{1854.
[16] Sederman, A., Johns, M., Alexander, P., and Gladden, L. Structure-
°ow correlations in packed beds. CHEMICAL ENGINEERING SCIENCE 53,
12 (JUN 1998), 2117{2128.
[17] Seymour, J., and Callaghan, P. Generalized approach to NMR analysis
of °ow and dispersion in porous media. AICHE JOURNAL 43, 8 (AUG 1997),
2096{2111.
[18] Stapf, S. NMR investigations of correlations between longitudinal
and transverse displacements in °ow through random structured media.
CHEMICAL PHYSICS 284, 1-2, Sp. Iss. SI (NOV 1 2002), 369{388.
44
in a porous medium. PHYSICS OF FLUIDS 15, 9 (SEP 2003), 2649{2660.
[12] Koch, D., and Brady, J. A nonlocal description of advection di®usion
with application to dispersion in porous-media. JOURNAL OF FLUID
MECHANICS 180 (JUL 1987), 387{403.
[13] Lebon, L., Leblond, J., Hulin, J., Martys, N., and Schwartz, L. Pulsed
¯eld gradient NMR measurements of probability distribution of displacement
under °ow in sphere packings. MAGNETIC RESONANCE IMAGING 14,
7-8 (1996), 989{991. 3rd International Meeting on Recent Advances in MR
Applications to Porous Media, LOUVAIN, BELGIUM, SEP 03-06, 1995.
[14] Manz, B., Alexander, P., and Gladden, L. Correlations between
dispersion and structure in porous media probed by nuclear magnetic resonance.
PHYSICS OF FLUIDS 11, 2 (FEB 1999), 259{267.
[15] Manz, B., Gladden, L., and Warren, P. Flow and dispersion in porous
media: Lattice-Boltzmann and NMR studies. AICHE JOURNAL 45, 9 (SEP
1999), 1845{1854.
[16] Sederman, A., Johns, M., Alexander, P., and Gladden, L. Structure-
°ow correlations in packed beds. CHEMICAL ENGINEERING SCIENCE 53,
12 (JUN 1998), 2117{2128.
[17] Seymour, J., and Callaghan, P. Generalized approach to NMR analysis
of °ow and dispersion in porous media. AICHE JOURNAL 43, 8 (AUG 1997),
2096{2111.
[18] Stapf, S. NMR investigations of correlations between longitudinal
and transverse displacements in °ow through random structured media.
CHEMICAL PHYSICS 284, 1-2, Sp. Iss. SI (NOV 1 2002), 369{388.
44
Page 45
[19] Stejskal, E., and Tanner, J. Spin di®usion measurements - spin echoes
in presence of a time-dependent ¯eld gradient. JOURNAL OF CHEMICAL
PHYSICS 42, 1 (1965), 288{&.
[20] Succi, S. The lattice Boltzmann equation for °uid dynamics and beyond.
Oxford University Press, 2001.
[21] Szymczak, P., and Ladd, A. Boundary conditions for stochastic solutions of
the convection-di®usion equation. PHYSICAL REVIEW E 68, 3, Part 2 (SEP
2003).
[22] Tallarek, U., van Dusschoten, D., Van As, H., Bayer, E., and
Guiochon, G. Study of transport phenomena in chromatographic columns
by pulsed ¯eld gradient NMR. JOURNAL OF PHYSICAL CHEMISTRY B
102, 18 (APR 30 1998), 3486{3497.
[23] Tessier, J., and Packer, K. The characterization of multiphase °uid
transport in a porous solid by pulsed gradient stimulated echo nuclear magnetic
resonance. PHYSICS OF FLUIDS 10, 1 (JAN 1998), 75{85.
45
in presence of a time-dependent ¯eld gradient. JOURNAL OF CHEMICAL
PHYSICS 42, 1 (1965), 288{&.
[20] Succi, S. The lattice Boltzmann equation for °uid dynamics and beyond.
Oxford University Press, 2001.
[21] Szymczak, P., and Ladd, A. Boundary conditions for stochastic solutions of
the convection-di®usion equation. PHYSICAL REVIEW E 68, 3, Part 2 (SEP
2003).
[22] Tallarek, U., van Dusschoten, D., Van As, H., Bayer, E., and
Guiochon, G. Study of transport phenomena in chromatographic columns
by pulsed ¯eld gradient NMR. JOURNAL OF PHYSICAL CHEMISTRY B
102, 18 (APR 30 1998), 3486{3497.
[23] Tessier, J., and Packer, K. The characterization of multiphase °uid
transport in a porous solid by pulsed gradient stimulated echo nuclear magnetic
resonance. PHYSICS OF FLUIDS 10, 1 (JAN 1998), 75{85.
45
Sign up today - FREE
Mendeley saves you time finding and organizing research. Learn more
- All your research in one place
- Add and import papers easily
- Access it anywhere, anytime
Start using Mendeley in seconds!
Readership Statistics
1 Reader on Mendeley
by Discipline
by Academic Status
100% Researcher (at a non-Academic Institution)
by Country
100% United Kingdom
Groups
 |