Turbulent stress measurements with phase-contrast magnetic resonance through tilted slices

Aiming at turbulent measurements in opaque suspensions, a simplistic methodology for measuring the turbulent stresses with phase-contrast magnetic resonance velocimetry is described. The method relies on flow-compensated and flow-encoding protocols with the flow encoding gradient normal to the slice. The experimental data is compared with direct numerical simulations (DNS), both directly but also, more importantly, after spatial averaging of the DNS data that resembles the measurement and data treatment of the experimental data. The results show that the most important MRI data (streamwise velocity, streamwise variance and Reynolds shear stress) is reliable up to at least r¯=0.75\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{r}} = 0.75$$\end{document} without any correction, paving the way for dearly needed turbulence and stress measurements in opaque suspensions.

This is crucial in, for example, the pulp and paper industry where not only do fibre suspensions alter turbulent quantities, but the flowfield in many devices used in the paper making process are complex in nature (MacKenzie et al. 2014).
Methods such as ultrasonic doppler velocimetry (UDV) and electrical impedance tomography (EIT) have conventionally been used for measuring the flowfield of opaque suspensions, whereas particle image velocimetry (PIV) and laser doppler velocimetry (LDV) are more commonly used for measuring the flowfield of transparent fluids. An advantage of MRV over other experimental techniques is that it can measure the full 3D velocity field in a few minutes for both transparent and opaque fluids (Caprihan and Fukushima 1990). The limitation of MRV is that in many cases computational fluid dynamics (CFD) simulations or PIV data are needed to provide details about the turbulent velocity fluctuations and Reynolds stresses (Elkins et al. 2009;Iaccarino and Elkins 2006). If MRV is to become a reliable method for measuring the flowfield of opaque suspensions at high Reynolds numbers, where conventional methods are unsuitable, it is important to extend the capacity towards measuring turbulent quantities. In doing so, it should be shown that MRV can produce rms values and Reynolds shear stress that agree with direct numerical simulations (DNS) in a well-defined geometry. This paper details our efforts in utilising phase-contrast MRV to measure the variance and covariance components of the Reynolds stress tensor for water in pipe flow. It will be shown that MRV provides reliable and accurate data in a large portion of the pipe. Works on turbulent flow in fibre suspensions will follow.
Measurements of diffusion with magnetic resonance has a long history (Stejskal and Tanner 1965) and this ability can be used to quantify turbulence through turbulent

Introduction
Magnetic resonance velocimetry (MRV) has been gaining popularity in engineering design due to its versatility in measuring laminar and turbulent quantities in complex geometries and opaque suspensions. In many engineering applications, measuring turbulent quantities enables design modifications to enhance production and reduce losses.
1 3 51 Page 2 of 12 diffusion. There have been numerous investigations focusing on deriving a method to measure the mean and rms component of velocity using MRV (Dyverfeldt et al. 2006;Gao and Gore 1991;Gatenby and Gore 1994;Li et al. 1994;Newling et al. 2004;Kuethe 1989). Here, the focus will be on methods where the rms value of the velocity is found by assuming turbulent diffusion to be much greater than the molecular self-diffusion of the moving spins. In doing so, turbulence can be modelled as the loss of the net magnetisation through intravoxel phase dispersion. Spins within a voxel which precess at a different phase will be incoherent to each other, reducing the magnitude of the phase-encoded image.
Depending on method, the rms value of velocity is calculated by modelling the phase-variance for different time intervals, or by assuming a phase correlation to the first gradient magnetic moment (see Eq. 3). Both methods are required to assume a spin velocity distribution of a known form. The works of Seymour and Callaghan (1996), Scheven et al. (2005), Caprihan and Seymour (2000) and Jullien and Lemonnier (2012) shows promise towards measuring turbulent diffusivity where the velocity distribution follows a non-Gaussian behaviour, examples of which include the inner layer in fully turbulent pipe flow. By comparing the evolution of the normalised rms widths, Scheven et al. (2005) found that the second and third moments of the perturbed probability distribution for molecular displacements could be found in terms of the measured moments from the original distribution. Protocols aiming at direct measurements of the velocity distribution has been demonstrated to obtain good agreement with DNS results (Jullien and Lemonnier 2012).
The difficulty inherent in applying this type of analysis in fully turbulent flows surrounds the inability to control the displacement field in regions where the velocity distribution is not Gaussian, largely a result of inadequate spatial resolution near the wall. Similar to Seymour and Callaghan (1996); Seymour et al. (2000), an effective diffusivity can be found by limiting the dispersion coefficient to low q-value data, namely where the attenuation function is logarithmically linear with increasing q-values (i.e. Gaussian approximation). In the work of Gao and Gore (1991), a velocity autocorrelation for the phase-variance and a Gaussian phase distribution was used to derive an expression for signal loss due to phase dispersion.
The more recent work of Dyverfeldt et al. (2006) is similar as that of Gao and Gore (1991), however, they assume that a phase-contrast MRI signal in the presence of a first gradient magnetic moment is governed by the spin velocity distribution, magnetic field inhomogeneities and phase accumulation from the velocity field. With the use of two phase-contrast MRI scans, the effect of field inhomogeneities can be eliminated, and the rms of the velocity can be found. Both works use similar approaches towards measuring an rms value of velocity, the difference coming from the assumption that the correction for gradient amplitude duration, separation time, and rise time cancel when two successive scans, which differ in first gradient moment, are compared.
In this work, we adopt the method of Dyverfeldt et al. (2006), a general form for measuring the rms component of velocity with phase-contrast MRV. In the implementation of 2D phase-contrast MRV available to us, the flow encoding vector is limited to only being defined normal to the slice orientation. When measuring quantities such as the Reynolds stresses, multiple encoding directions would be required to isolate the effects of spanwise perturbations in the flow. Without simultaneously encoding the streamwise and spanwise directions (Elkins et al. 2009) it is difficult to acquire quantitative information about the Reynolds stresses due to coordinate system transforms and slice orientation. We discuss here the methodology required to measure the variance and covariance components of the Reynolds stress tensor in a straight pipe with a conventional 2D phase-contrast MRV sequence. We discuss the reliability of this method by comparing our results to DNS data for a range of Reynolds numbers.

Theory
Magnetic resonance velocimetry (Caprihan and Fukushima 1990;Gladden 1994) works from the basis that hydrogen nuclei are dipoles. In the presence of an external magnetic field, spin packets of hydrogen nuclei are represented by their net magnetisation vector, , which is oriented along the direction of the magnetic field. When the spins are relaxed, hydrogen nuclei will align in the direction of the external magnetic field while precessing at the Larmor frequency, 0 = | |∕2 , in the direction of the external magnetic field. Here, is the magnetic field strength, and is the gyromagnetic ratio. The orientation of the spin magnetisation vector, , can be altered by applying radiofrequency (RF) pulses that match the Larmor frequency of the precessing spins. This effectively tilts the spins away from their original alignment, which upon turning off the RF-pulse, will induce relaxation of the spins towards the external magnetic field direction. This relaxation process is described by the Bloch (1946) equation where T1 and T2 are the spin-lattice and spin-spin relaxation times, respectively. If, for example, the external magnetic field direction is along k , then T1 describes how long it will take for the longitudinal component of to recover to 63% of its original value (M 0 ). T2 is the time it takes for the transverse component of to lose 63% of its original value. As the magnetisation returns to its original alignment, a time varying signal is induced in the receiver coils. This signal is detected in the MR system and is the foundation for magnetic resonance imaging (MRI) and MRV. Note that a short T2 relaxation time, and a long T1 relaxation time would result in a reduction in the signal intensity at shorter times.
To generate 2D images, spatially varying magnetic field gradients are used to encode the phase and frequency of the moving spins. The spins will accumulate phase according to where 0 is the phase accumulation due to magnetic field inhomogeneities, is the magnetic field gradient strength, and is the position. The Larmor frequency of the spins will change according to = (| | + ⋅ )∕2 . If spins are moving at a velocity parallel to a magnetic field gradient, they will acquire a phase shift which is proportional to their velocity . The total phase accumulation of a moving spin can then be written as where, and are the zeroth and first order gradient magnetic moments, respectively. Conventionally, bipolar flow encoding gradients are used to isolate the contribution of velocity on the total phase accumulation of moving spins. Since the accumulation of phase from static spins is independent of gradient polarity, the velocity of the moving spins can be calculated from the net phase shift between the bipolar gradients. The velocity can then be directly calculated by where,

Velocity variance in the direction of the slice
In MRV, intravoxel phase-dispersion within a spin packet reduces the amplitude of the acquired signal from the MR system, which in turn can be modelled as an rms value in velocity. In this work, the method of Dyverfeldt et al. (2006) was adopted, which first starts with considering the signal S( ) within a voxel for given flow and encoding parameters as where s( ) is the spin velocity distribution, defined in this work to be Gaussian, namely as where is the standard deviation. Equation 6 shows that for phase-contrast MRV, the signal acquired from the MR system within a voxel is governed by magnetic field inhomogeneities and intravoxel phase accumulation. Previous works (Dyverfeldt et al. 2006;Gao and Gore 1991) found that if a velocity distribution is assumed, such as the Gaussian distribution shown in Eq. 7, then S( ) = g(u rms ). Substituting Eq. 7 into Eq. 6 results in By combining the results of two scans with different first gradient moments, the effect of magnetic field inhomogeneities can be eliminated. Since the variant term, 2 , is located in the real part of the exponent, the magnitude portion of S( ) gives The flow encoding vector can be defined to measure the velocity in any direction on a 3D cartesian coordinate system. If the encoding direction is defined as: = xî + yĵ + zk, where k is along the streamwise direction z, and ̂j and ̂i are along the spanwise directions y and x, respectively, then Fig. 1 shows the resulting slice orientation for flow encoding vectors = 0î + 0ĵ + 1k and = √ 2∕2î + 0ĵ + √ 2∕2k (remember that our flow encoding gradient can only be normal to the slice).
The tilting of the orientation of the slice to flow encode a straight pipe along a direction in the ̂ik or ̂jkplane results in an elliptic image as shown in Fig. 1 extrinsic matrix rotations the coordinate system defined with (i.e. normal to the streamwise direction) can be derived from , and the Reynolds stress tensor can thus be computed using the measured velocity variances from multiple encoding directions.

Reynolds stresses
The method adopted here for determining the Reynolds stresses was originally used for hot-wire anemometry or single probe LDV measurements (e.g. Inasawa et al. 2003), and has been extended in this work to account for the use of multiple flow encoding directions. To compare the results obtained from MRV to DNS, a relationship between the variance of the measured velocity, flow encoded in a Cartesian system, to the Reynolds stress tensor in cylindrical coordinates can be achieved by decomposing the measured velocity variance to the averaged variance and covariance components of the stress tensor. Table 1 summarises the flow encoding vector directions and the resulting variance and covariance components found from selectively decomposing lines from the 2D images.
To exemplify Table 1, the instantaneous velocity decomposition for a flow encoding vector throughout the crosssection of the pipe is detailed. The tangential position within the pipe is defined as where y and x define the position of a voxel on the Cartesian coordinate system. The instantaneous measured velocity for a flow encoding vector can then be shown to be where u z , u , and u r are the fluctuating components of velocity, and êi and êk are the components of the encoding vector as described in Table 1. A relationship between the variance of the measured velocity, ⟨u 2 M ⟩, and the Reynolds stress components can be found by expanding Eq. 11, squaring the result and averaging to get To calculate the variance and covariance components in Eq. 12 from the measured MRV data ⟨u 2 M ⟩, the set of linear equations, defined by the flow encoding vectors in the ̂ik -plane and k -plane (i.e. ⟨u 2 z ⟩) shown in Table 1, are solved using a least-squares fit. If a line of data along Φ o = 0 is selected, the linear system to be solved can be expressed as where, If a line of data is selected along Φ o = , the terms ⟨u 2 r ⟩ and ⟨u z u r ⟩ in Eq. 13 would be replaced with ⟨u 2 ⟩ and ⟨u z u ⟩, respectively. In this work, the covariant turbulent stresses ⟨u z u ⟩ and ⟨u r u ⟩ were not investigated.

Extracting measured data from DNS
In order for the data obtained from Eq. 13 to be accurate, the artefacts induced by spatial averaging in combination  The encoding vector is defined as: = xî + yĵ + zk, where k is along the streamwise direction. The selected angle, defined as "Φ o ", is the angle about the pipe's origin a line of measured velocity variance data is decomposed. The "components" indicates which components of the Reynolds stress tensor, in cylindrical coordinates, the measured variance decomposes into along the selected line with tilting of the slice must be limited. Data from two published DNS studies (Ahn et al. 2015;Khoury et al. 2013) is used to evaluate these effects. First, a flowfield of the measured quantities in Eq. 11 for multiple flow encoding vectors is created. An example of the velocity variance, ⟨u 2 M ⟩ 1∕2 is shown in Fig. 2. As a comparison, the velocity variance measured directly from MRV with an identical flow encoding vector is shown. It can be seen that qualitatively DNS and MRI agree well, where both exhibit the same asymmetric variance over the cross-section of the pipe. The asymmetry is a result of slice obliqueness and averaging along the slice thickness direction.
Using elemental matrix rotations, the 2D field of ⟨u 2 M ⟩ from DNS was rotated such that the xy-plane is normal to a defined flow encoding vector. An illustration of the rotation is shown in Fig. 3. To investigate the effect of slice thickness and averaging, the variance of the measured quantity ⟨u 2 M ⟩ in the DNS was averaged normal to the plane, similar to the operation of the MRI. This can be expressed as where x ′ c and y ′ c are the coordinates along the centre of the slice, and h(2 ) is the number of positions along z ′ that were averaged. For all numerical simulations, Δz � was set to Δy � .
The DNS data was interpolated to match the MRI grid in the x ′ y ′ -plane by averaging according to where FOV is the MRI's field-of-view, and IPR y is the y-component of the MRI's in-plan matrix resolution (see Sect. 2.4). The spatially matched DNS data can then be rotated back into the xy-plane by where z is zero ∀ (x, y). The variance and covariance components are finally calculated using Eq. 13. Note that a linear interpolation along Φ o was used to match the spatial resolution of the rotated data to that acquired with êk = 1.

MR system
All of the experiments in this work were performed on a 1T MR system (Aspect Imaging, Shoham, Israel) with a 30 G/cm peak gradient strength. The MR system composes of 3 gradient coils and 1 radio frequency coil. Data was collected using a variety of flow encoding directions where the field-of-view (FOV) ranged from 38 to 56 mm, depending on the orientation of the slice. All measurements were made with water with no contrast agent. The T1 and T2 relaxation times are assumed to be approximately 4 and 2 s, respectively. The in-plane matrix resolution (IPR) was chosen to be 128 × 128 pixels. The PC-MRV sequence,  Fig. 3 Illustration of coordinate rotation. x ′ and z ′ denote the rotated coordinate system, where ± is 1/2 of the slice thickness in front of or behind the centre of the slice. Similarly to the MRI, the data is averaged normal to the rotated coordinate system, along z ′ 51 Page 6 of 12 shown in Fig. 4, was set with a slice thickness of 5 mm or D pipe ∕7 and a repetition time (TR) defined as either 10 ms for Re ≤ 321 or 6.1 ms for Re ≥ 597. The echo time (TE) for all Reynolds numbers was 3 ms and the flip angle was 30 • . Note that in this work, Re = D pipe u ∕ , where u is the friction velocity, found experimentally via the pressure gradient relation u 2 = w = D pipe |dP∕dz|∕4. Details on the measurements can be found in Sect. 2.5. The density and viscosity of water, and respectively, where used for all calculations.
To increase the signal-to-noise ratio, five signal acquisitions, or combined images, were prescribed for each experiment where, depending on Reynolds number (see Table 2), between 25 and 120 experiments were repeated. With a TR time of 10 ms the total scan time for 25 experiments at one bipolar gradient amplitude was 6.2 min. The velocity encode (VENC) for measuring the mean component of flow was set such that the maximum velocity in the phaseunwrapped image was always less than the VENC. The VENC defined for the turbulent stress measurements was chosen to be 15% of the velocity encode used to measure the mean component of flow at the same Re . That particular value was chosen because the DNS data showed that no variant or covariant rms components exceeded 14% of the centreline velocity, therefore, the resolution in phase variance was increased by choosing a VENC to be 15% of what was used to measure the mean flow.

Flow loop
A centrifugal pump (Flygt model no. 040-6716260) provided flow, ranging from 0.11 to 2.5 l/s. The volumetric flow rate was measured using a Krohne aquaflux 070 flowmeter and the pressure gradient was measured using a Fuji electric FCX-AII differential pressure transducer. Both the flow meter and the differential pressure transducer were connected to a NI-9205 DAQ device for data logging. Flow exits the pump through approximately 10 m of flexible 5 cm I.D. tubing, before splitting into two 5 cm I.D. PVC pipes. The split flow enters a single 5 cm I.D. PVC with a 90 • bend before constricting to the 3.5 cm I.D. acrylic test section. The MRV is located more than 140 pipe diameters from the constriction. Figure 5 shows the friction factor data for a range of Reynolds numbers measured in the test rig. It can be seen that a transition from laminar flow occurs near Re ≈ 2000, where the data agrees well to the Colebrook friction factor line (Colebrook and White 1936), with an assumed roughness height of 1.6 μm, from Re = 4000 to 10,000. Note that the Reynolds number was experimentally calculated as Re = 4 Q flow ∕ D pipe , where Q flow is the volumetric flow rate measured from the Krohne aquaflux flowmeter. The data suggests that the experimental rig agrees with the typical range for transition from laminar to turbulent flow. The   position of the differential pressure transducer was found to have no effect on the measured friction factor, indicating the flow is fully developed for all Reynolds numbers within the MRI test section.

Results
The results presented here are broken into two sections: evaluation of the errors inherent in the MRI measurements by utilisation of the DNS data, and results from the MRI. The evaluation of the error was made for Re = 180, 360, 550, 1000, and 3000, whereas eight Reynolds numbers were measured using MRV. A summary of the measurements can be found in Table 2.

Evaluation of errors in the measurements
The error evaluation using DNS data is used to determine the optimal combination of encoding vectors to minimize the effect of spatial averaging along the rotated streamwise direction, z ′ , and secondly, the effect of slice thickness on the deviation from the un-averaged ⟨u z u r ⟩ and ⟨u r u r ⟩ components. Nine pairs of encoding vectors were considered, ranging from to in increments of ∕36 (see Table 1). The covariance and variance components, ⟨u z u r ⟩ and ⟨u r u r ⟩, were calculated using Eq. 13 for multiple combinations of encoding pairs. Figure 6 shows the L2-norm of ⟨u z u r ⟩ and ⟨u r u r ⟩, normalised to the streamwise centreline velocity U ∞ , for every combination of flow encoding vector pairs. This can be expressed as, for example, In Eq. 15 the difference in the covariance component ⟨u z u r ⟩ is defined as, where ⟨u z u r ⟩ measured is obtained by averaging DNS data and is a function of N. Note that r in Eq. 15 is the nondimensional radius and N is the number of encoding pairs ranging from 1 to 9. The L2-norm here is used as a quantitative measure of the deviation from the original DNS = ± sin 36̂i + 0ĵ + cos 36k = ± sin 4̂i + 0ĵ + cos 4k , Δ⟨u z u r ⟩ = ⟨u z u r ⟩ measured − ⟨u z u r ⟩ DNS , data due to spatial averaging of the MRV for each combination of encoding vectors. It can be seen from Fig. 6 that the L2-norm for the value of ⟨u z u r ⟩ is always less than 3% of U ∞ , regardless of the number of encoding pairs considered in Eq. 13. The L2-norm in the rms value of ⟨u r u r ⟩, however, was found to deviate up to 43% of U ∞ if only ± 9 = ± sin 4î + 0ĵ + cos 4k was used. This is a result of the gradient in the rms values along z ′ , which increases with decreasing size of the outer layer (i.e. (1 − r) + > 50) or Re . Note that (1 − r) + is the distance from the wall measured in viscous lengths, calculated only with DNS as (1 − r) + = u (R − r)∕ .
In Fig. 7 the effect of slice thickness is shown in the diagnostic plot for Re = 180 and 1000. The result of only spatially condensing the data to match the resolution of the MRI (refer to Sect. 2.3) is illustrated with 2 = 0, where it can be seen that �Δ⟨u z u r ⟩ 1∕2 � and �Δ⟨u r u r ⟩ 1∕2 � for Re = 1000 is greater than with Re = 180 near U∕U ∞ = 0.65. This correlates well to (1 − r) + = 50 or r = 0.95 for Re = 1000, which effectively means that the MRI will likely be unable to measure beyond the peak in ⟨u z u z ⟩ for high Reynolds number flows purely due to spatial resolution. Increasing the slice thickness to 2 = D∕7 and 2D/7, it can be seen that the effect of averaging along z ′ resulted in Δ⟨u z u r ⟩ 1∕2 values always below 1% of U ∞ for either Re equal to 180 or 1000 in the outer layer.
The most profound effect of the slice thickness is on the radial variance ⟨u r u r ⟩ for Re = 180, where �Δ⟨u r u r ⟩ 1∕2 � Fig. 6 L2-norm of the variance and covariance turbulent stress components from DNS data at Re = 550 as a function of number of velocity encode pairs, N, used in Eq. 13. The slice thickness for this case was defined as 2 = D∕7. The circles show the resulting L2-norm when using a random combination of flow encoding vector pairs, whereas the square markers show the optimal combination of encoding pairs. The light green square is the combination which minimizes the L2-norm in ⟨u r u r ⟩, whereas the blue square shows the combination which minimizes the L2-norm in ⟨u z u r ⟩ 51 Page 8 of 12 varied between 2 and 8% of U ∞ within U∕U ∞ = 0.8 and 0.6, respectively. As mentioned, this is a result of the outer layer starting at r = 0.72 or U∕U ∞ = 0.83. The radial variance as well diverged from zero near (1 − r) + = 50 or U∕U ∞ ≈ 0.65 for Re = 1000. The error in both ⟨u z u r ⟩ and ⟨u r u r ⟩ near U∕U ∞ = 1 was found to exist for all Reynolds numbers, and is a result of the steep gradient in all rms components within U∕U ∞ ∈ (0.92, 0.96).

MRV results
The experimental rig used for this work is primarily used for analysing multiphase cellulose fibre suspension flows, and as a result does not contain a settling chamber to remove swirl from pipe bends and pumping, or oscillations from connections in the piping system. This is because fibre suspensions tend to flocculate and clog settling chambers. Upon measuring the mean tangential velocity profile it was found that the flow of water exhibited relatively high magnitudes of swirl. Figure 8 shows the magnitude of the swirl measured from the MRI in a format similar to the diagnostic plot. A swirl number was defined as, S = U ,w ∕U ∞ , where U ,w is the maximum swirl component measured nearest to the wall. Conventionally, a swirl number is used for rotating pipe flows, however, it has been adopted in this work to account for the effects of the swirl on the mean streamwise velocity component. Of the Reynolds numbers investigated in this work, it was found that the swirl number ranged from 0.1 to 0.7. This indicates that near the wall there is shear layer from the swirling component, leading to an increase in the friction velocity at lower bulk velocities. To account for the swirl when comparing the rms values measured from MRV to DNS, a scaling factor was introduced on the centreline velocity measured from the MRI. This scaling is justified with the work of Örlü and Alfredsson (2008), who studied the mean and rms of the streamwise component in a swirling jet emanating from a fully developed axially rotating pipe flow. It was found that a rotating pipe of swirl number 0.5 reduced the centreline velocity by almost 12% at the pipe exit, whereas the rms component was unaffected.
We adopt the outer layer scaling law proposed by Oberlack (1999) for a rotating pipe flow, which can be expressed as  Table 2 for the marker definitions where is a function of the velocity ratio and is a constant. Since the form of is unknown, and the swirl cannot be controlled for a given Re , a first order empirical swirlcorrection for U ∞ was used. Assuming that = 1 and is a non-uniform expansion, chosen to be with limits, Substituting and into Eq. 16 and rearranging yields For simplicity, the corrected centreline velocity will be defined as Figure 9 compares the MRV data in the diagnostic plot, with and without scaling in U ∞ , to DNS. It can be seen that without the scaling, the normalised rms values deviate further from the DNS data for Reynolds numbers which display strong swirl (e.g. Re = 160). Applying the swirl correction to the rms value improves the agreement to DNS for Re ≥ 321, however, the correlation is still limited to either the resolution of the MRI, or the onset of the viscous wall region, (1 − r) + < 50. It must be noted that we make no attempt to model the swirl in either the viscous wall region, or the outer layer. It is assumed that the effect of the swirling component on the distribution of the mean flow is negligible in the outer layer, and that the swirl correction fails for Re = 160 and 2283 as a result of dU ∕dr being highly non-linear in the outer layer.
To establish the effect of encoding angles on the calculated variance and covariance components from Eq. 13, the L2-norm was calculated for all possible combinations of encoding pairs ranging from 4 to 9 . Figure 10 presents the minimum L2-norm for each Reynolds number considered in this work, compared against DNS data at an Re closest to that measured with MRI. It can be seen that, depending on encoding angles, the L2-norm of ⟨u z u r ⟩ was below 0.06 for all Reynolds numbers. The encoding angles which minimised the L2-norm for three Reynolds numbers, Re = 321, 1115, and 1827 are highlighted. It was found that ±45 • or ± ∕4 was a common encoding pair between the Reynolds numbers shown in Fig. 10. It appears that the encoding  Table 2 for the summary of experiments vector 9 gives leading order accuracy for the covariance component, where higher order accuracy is achieved by acquiring data at encoding angles which approach the direct measurement of the streamwise variance (i.e. 3 ). The L2-norm of the radial variance was found to range between 0.18 and 0.46 for Re = 160-2283. It was found that out of all possible combinations, the minimum L2-norm was acquired with only 9 . This is not surprising considering that an encoding angle of 90 • or ∕2 would directly measure the radial variance, however, only at the centre of the pipe. Based upon the corrected DNS, measuring the radial variance will be a challenge as steep encoding angles result in increased error but approach a direct measurement of the radial variance as êi → 1. A reduction in slice thickness would undoubtedly reduce the averaging effects, but decreases the signal-to-noise ratio. It also must be noted that the effects of the swirl on the radial variance are not quantified in this work, which may well be contributing to the divergence from DNS. Limiting the swirl and optimising the MR measurement is required before a quantitative analysis on ⟨u r u r ⟩ is possible. Figure 11 presents the DNS and MRV (see Table 2 for marker summary) covariance rms component in the diagnostic format for all the Reynolds numbers considered in this work. The data shown was calculated using the encoding angles which minimised the L2-norm, the values of which can be seen in Fig. 10. The MRV data measured at Re = 321 and 597 (top figure) was found to agree to DNS results at Re = 360 (green) and 550 (red), respectively. The middle plot in Fig. 11 shows good agreement for Re = 922 and 1115 when compared to DNS acquired at 1000. The normalised covariance measured from MRI at Re = 1600 is slightly discontinuous when compared to DNS near U∕U ∞ = 0.95, a possible result of the correction factor failing to account for the swirl. Similarly, Fig. 11 (bottom) indicates that the swirl correction factor failed for Re = 180, 1827, and 2283. The covariance rms measured at Re = 2283 was significantly lower than the DNS data both with and without the swirl correction factor, an analogous result to the streamwise variance shown in Fig. 9. Figure 12 presents the radial distribution of the normalised turbulent stress for the MRV cases with minimal swirl. The turbulent stress and position were defined as where u is the friction velocity determined experimentally from the axial pressure gradient measured along the test section. The MRV results are in excelled agreement to the corrected DNS where the discrepancies are largely a result of the uncertainty in wall position. To overcome this uncertainty, it is suggested to diagnose turbulent stress data when using titled slices in the format shown in Figs. 9 and 11.

Conclusions
A method for measuring the turbulent stresses using phasecontrast magnetic resonance velocimetry with the flow encoding gradient normal to the slice has been detailed. Data from DNS was used to quantify the effects of spatial averaging from the MRI, where it was found that slice thickness and encoding angle had a much more profound effect on the radial variance rather than the covariance component. Results from MRV were presented in a diagnostic format for the streamwise variance and a covariance component of the turbulent stress tensor. It was found that the streamwise variance agreed well to DNS once a correction factor was applied to account for the swirl within the test section. The swirl is a deficiency in the flow loop since it is optimised to handle suspensions, not to create a perfect single phase flow.
The covariance rms ⟨u z u r ⟩ 1∕2 calculated from the MRV data was in good agreement to DNS for Reynolds numbers with minimal swirl without any corrections. The discrepancies between DNS and MRV measurements are likely a result of the additional shear caused by the swirl. Quantifying the angular and axial momentum along the test section would likely rectify swirl effects, however, is outside the scope of this work. + T = ⟨u z u r ⟩ u 2 , Fig. 10 Minimum L2-norm of the variance and covariance turbulent stress components from MRI data. The L2-norm was calculated using Eq. 15 where Δ⟨u i u j ⟩ = ⟨u i u j ⟩ MRI − ⟨u i u j ⟩ DNS . The DNS data which was used to compare against the MRI data was chosen such that the difference in Re was at a minimum. Recall that the MRV data is compared against DNS acquired at Re = 180, 360, 550, 1000, and 3000. The encoding vector arrows are described in Table 1 Page 11 of 12 51 The use of tilted slices was found to be effective in measuring the covariance component ⟨u z u r ⟩ and the streamwise variance of the turbulent stress tensor for r∕R < 0.75 with phase-contrast MRV. Excellent agreement with DNS was obtained in this region. Whether A particularly interesting application for this method will be measurements in opaque suspensions, where turbulent stress data is crucial for modelling but hard to obtain with optical or intrusive methods. Fig. 11 Diagnostic plot of the covariance rms, ⟨u z u r ⟩ 1∕2 , from MRV and DNS. The top and middle figure shows the DNS and MRI data where the effect of swirl was found to be minimal, whereas the bottom figure shows the data for Re = 160, 1827, and 2283. The solid lines are DNS data for Re = 180, 360, 550, 1000, and 3000 whereas the dashed lines containing markers are the corrected DNS data. As before, the filled stars markers indicate where (1 − r) + = 50 and the filled diamonds markers indicate where r = 1-10 Δy � ∕D pipe Fig. 12 Turbulent stress data from MRV and DNS. The Reynolds numbers shown from MRV measurements were found to have the least amount of swirl, namely Re = 597, 922, and 1115. The solid lines are DNS data from calculations at Re = 550 and 1000 whereas the dashed lines containing markers are the corrected DNS data. As before, the filled stars markers indicate where (1 − r) + = 50 and the filled diamonds markers indicate where r = 1-10 Δy � ∕D pipe