In this paper, we introduce a new technique based on X-ray radiography with high temporal and spatial resolutions to study fast suspension flows regardless of optical access. We benefit from the axial symmetry of our flow configuration, a wide gap Couette setup, to extract a 3D solid volume fraction field from a single X-ray projection image. We propose a mathematical algorithm based on the inversion of Abel transform in conjunction with H1 regularization and data denoising to measure the solid volume fraction field in suspensions in a fraction of a second. We show that the results are in excellent agreement with those obtained from micro Computed Tomography (CT scan) in one hour. This significant reduction in the data acquisition time opens a new avenue in the field of suspensions. As a proof of concept, we study the kinetics of shear-induced migration for suspensions of particles in both Newtonian and yield stress suspending fluids. The latter experiments include two different conditions: With and without a plug region. In both cases, we are able to capture in detail the kinetics of migration. In the presence of a plug region, we manage to accurately describe the particle accumulation at the interface between the sheared and the static regions. Remarkably, even in the absence of sedimentation, the concentration profiles show a complex 2D structure, with no z-invariant region, which illustrates the strong impact of top and bottom boundary effects on migration. We also show the importance of boundary effects on the shear induced migration of particles in a Newtonian suspending fluid. This further shows the necessity of developing techniques that give access to the full spatial concentration field, as the one we present here.
I. INTRODUCTION
An extremely large number of industrial processes involve suspensions such as: Food production, waste disposal, fresh concrete casting, papermaking, drilling muds and cuttings transport, paper coating and ceramic injection molding, as well as natural phenomena like flows of slurries, mudslide, and lava. These materials show all possible non-Newtonian behaviors: Plasticity, thixotropy, shear-thinning, shear-thickening, normal stress differences, etc. At a macroscopic scale, their flows display complex features such as shear-banding, and they tend to develop strong particle concentration inhomogeneities.
In this context, there has been considerable work on the rheology and modeling of noncolloidal Newtonian suspensions, in which the background fluid is Newtonian and the complexity comes from the presence of particles [1–7]. Scholars have observed and studied unexpected aspects of Newtonian suspensions such as irreversible migration of the particles in a reversible—Stokes—regime (shear induced migration [8]), shear thickening viscosity [9], and shear induced normal stress difference [10]. Perhaps, the most famous feature is shear induced migration that has been studied since the remarkable work of Leighton and Acrivos [8], which explains the migration of particles in an inhomogeneous shear field from a high shear rate zone to a low shear rate zone. Although non-Newtonian suspensions such as noncolloidal yield stress suspensions have broad application in nature and industry, they have been much less studied and not much is known about their characteristics [11–13]. In most industrial applications and natural phenomena, the suspensions are either far from equilibrium or strongly non-Newtonian, meaning that the flow kinetics are not only strain-dependent but also strain-rate dependent. Therefore, experimental techniques must be developed to analyze the flows of these complex suspensions in both steady and transient regimes.
Optical techniques such as Particle Imaging/Tracking Velocimetry (PIV/PTV) and Laser Doppler Velocimetry (LDV) have been widely used because of their high spatial and temporal resolutions to study suspensions that have optical access [14]. These techniques need a laser sheet to illuminate the suspension, which makes them incapable to study opaque systems [15]. This limits these techniques to transparent suspensions, where for dense suspensions accurate Refractive Index Matching (RIM) between the fluid and the particles remains a substantial step [16,17]. When optical access is limited, there have been several advances to visualize multiphase flows and different methods such as ultrasonic, Nuclear Magnetic Resonance (NMR), MRI, and X-ray have been widely adopted by scholars. In-depth reviews of available methods indicate that none of them fully satisfies the requirements to study concentration and velocity profiles in multiphase flows [18–20].
Since 1970, ultrasonic velocimetry has been considered as one of the most capable methods to measure the velocity profiles in fluids and especially in studying multiphase flows (e.g., gas and liquid); 2D velocity profiles can be obtained with high spatial and temporal resolutions, but it is not commonly employed in the study of dense systems (e.g., dispersed systems with ) or highly viscous flows mainly because of the presence of multiple scatterings and high attenuation of ultrasonic wave [18,21,22]. Moreover, ultrasonic techniques are not yet adapted to the measurement of volume fraction profiles in concentrated suspensions.
NMR/Imaging (NMR/I) has also been used to study suspensions velocity and particle concentration profiles for a wide range of concentrations [23,24]. This technique has been adopted in different studies such as shear-induced migration of particles [13,25,26] and droplets [27], shear-induced sedimentation [28], measuring the size of droplets in emulsion [29], and measuring the acceleration field independent of the velocity field [18]. Whereas 1D velocity profiles are measured within 1 s, the main challenge for this method remains that the data sampling time scales with the number of scanning points, meaning longer time is required for higher spatial resolution and for 2D maps [18]. Furthermore, volume fraction measurements during flows are tricky and can be difficult or impossible to interpret as the material needs to remain in the same position during the sequence of radiofrequency (RF) and the applied field gradient pulses to measure the liquid content locally. These shortcomings finally make the method capable of measuring volume fractions mostly in static conditions with high spatial resolution [20].
X-ray radiography and Computed Tomography (CT) have been developed mainly in the medical field [30,31]. X-ray radiography is the 2D projection of a 3D object on a detector, which is an attenuation map averaged over the object thickness. X-ray computed tomography is a 3D construction of the object from different 2D projections, which is based on the Radon transformation [32]. CT has high spatial resolution; however, it generally has low temporal resolution. Temporal resolution is typically 1 h for a standard lab X-ray facility, while it can be of the order of 1 s in a synchrotron. Since it needs the projection of the same object at different directions, CT has not been used in studying dynamic flows [19].
Here, we introduce a new technique using an X-ray/CT-scan system to study volume fraction changes in flowing suspensions. In order to increase the temporal resolution, we focus on a system that is axially symmetric, which allows us to get a 2D map of the absorption leading to a 2D map of the volume fraction from a single radiography projection using a simplified version of Radon transformation known as Abel transformation. This is a well-studied problem mathematically [33–37], although it has not been used to study dynamic suspensions [19]. Since our method works for any axially symmetric problem, it is perfectly suitable for viscosimetric flows, such as Couette flows, flows in a parallel plate and a cone-and-plate device, and Poiseuille flows. Consequently, this allows one to study in-depth well-defined problems from the point of view of the stress field (Poiseuille, Couette, cone-and-plate) or the shear rate field (parallel plate), in which the volume fraction evolves during a flow.
In order to implement and test this technique, we study the flow of a dense suspension of particles in a Newtonian and a yield stress suspending fluid yield stress fluid that is loaded in a wide gap Couette cell [13]. As in any suspension of non-Brownian hard spheres, the particles are prone to shear-induced migration, i.e., the volume fraction field is expected to evolve during time, which we expect to be able to capture with our method. The use of the yield stress fluid ensures that the particles do not move when the fluid is left at rest, which enables us to perform a 3D tomography to directly get the final volume fraction field by counting the particles, and allows us to validate the technique by comparing the 2D map extracted from a single X-ray projection image to the actual 2D map extracted from 3D tomography.
The geometry has axial symmetry satisfying the requirement to use Abel transformation. This opaque system cannot be studied by conventional methods such as PIV and PTV due to its limited optical access; the dense suspension of particles () makes the ultrasonic method inappropriate; and the fact that the volume fraction field evolves in time makes it vital to capture the features of the flow in the transient regime that makes NMR/I methods incapable of studying this flow as a result of their low temporal resolution.
We present our method by explaining the principle of the study in Sec. II and our experimental setup, the employed materials, and the governing equations of the X-ray system in Sec. III. These sections set the stage for the analysis of radiography images in Sec. IV, in which we explain the Abel transformation and the common problems in dealing with Abel inversion. We suggest and implement possible ways to address these problems such as regularization of the Abel transformation and denoising the output signal using appropriate filtering techniques. In Sec. V, we show that our methodology toward extracting information from a single 2D radiography projection gives very satisfactory results for a nonlinear dynamic system by comparing its steady state results with the results extracted from 3D tomography. In addition, the acquisition time with the proposed method is reduced from ∼1 h to a fraction of a second, which opens a new avenue in the field of suspensions. As an illustration, we study the kinetics of shear-induced particle migration in a yield stress fluid in two different conditions: With and without a plug region. In both cases, we capture in detail the whole kinetics. In the presence of a plug region, we manage to accurately describe the particle accumulation at the interface between the sheared and the static regions. Remarkably, even in the absence of sedimentation, the concentration profiles show a complex 2D structure, with no z-invariant region, which illustrates the strong impact of top and bottom boundary effects on migration. Moreover, we perform an experiment showing a typical example of shear induced migration of particles in a Newtonian suspending fluid. The results show the significant role of boundary effects on the flow dynamics. This further shows the necessity of developing techniques that give access to the full spatial concentration field, as the one we present here.
II. PRINCIPLE OF THE STUDY
Extracting a concentration map from a 2D image implies that concentration is a function of only two space variables in the problem under study. This information would be straightforward to obtain for a problem invariant in the direction orthogonal to the detector plane. However, if such a problem exists theoretically, this does not correspond to any well-defined experimental situation. On the other hand, in viscosimetric flows such as Couette flows, cone-and-plate flows, parallel-plate flows, and Poiseuille flows, which have axial symmetry, the volume fraction is expected to be invariant along the azimuthal angle θ. In cylindrical coordinates, it is then expected to depend on the radial variable r and on the axial variable z. Possible physical origins of the volume fraction inhomogeneities in a suspension are then shear-induced migration [8], sedimentation [28,38], and—in a Poiseuille flow—self-filtration [39,40].
The challenging problem is now to obtain volume fraction fields in cylindrical coordinates from intensity images I(y, z) in Cartesian coordinates. As we will show in Sec. III, intensity measurements performed at various y positions contain information averaged over different intervals of radial positions (Fig. 1), which makes the problem solvable. We will also show in Sec. IV that accurate information on volume fraction can be obtained provided a series of mathematical operations is applied to the data; this includes appropriate Fourier decomposition, regularization, and filtering.
In order to implement and validate the technique, we have chosen to study flows of a density-matched noncolloidal suspension in a wide-gap Couette flow. In such a situation, particles are expected to migrate toward the region of low shear, i.e., toward the outer cylinder, resulting in an r-dependent volume fraction field. Possible z-dependence of the volume fraction field is then expected, due to the top and bottom boundary conditions.
To validate the technique, we need to perform 3D-tomography on the same sample: This allows us to directly locate all of the particles in space, and to extract the actual volume fraction field from direct particle counting. The method we propose from 2D images thus has to provide quantitatively the same information. The particles have to stay in the same position for a long time at rest (typically 1 h) to allow for a 3D measurement. This is why we have chosen to use a yield stress fluid as a suspending fluid: After a given shear history, the particles are then gelled at rest even in the presence of slight density mismatch or mechanical noise. Using a yield stress, fluid has another interest for the purpose of illustration: Its flow properties depend on the shear rate, and different volume fraction fields are expected to be developed at low and high shear rates.
In our recent work [13], we have shown that similar to the case of Newtonian fluid, we have shear induced migration of particles in the same yield stress suspension. This shows that the migration is not coming from the elastic stresses. Moreover, it is noteworthy to mention that the goal of [13] was mainly to study the effects of the microstructure. We have shown a typical example of shear induced migration, in the absence of unyielded regions, via a different experimental technique involving NMR. In the NMR method, the time required for data acquisition is large, making it incapable of studying fast changes or monitoring continuously the evolution of solid volume fraction. Therefore, it was essential to make use of a new experimental technique to overcome this barrier and to tackle this problem in complex fluids, which is the main goal of our current research.
One of the strengths of the method we propose is that a standard rheometer can be used to perform the study: The only thing is to ensure that the rheometer body is not on the X-ray path. This is not possible for 3D measurements since the sample has to be rotated over to get a series of different 2D projections; the presence of the rheometer body on the X-ray path for some of the angles would not allow for a proper 3D reconstruction. Since we want to compare 2D and 3D measurements for validation of the technique, this led us to design a home-made Couette rheometer that has axial symmetry (Fig. 2).
III. EXPERIMENTS
A. Materials
The experiments are carried out in a wide gap cylindrical Couette cell (Rin = 5 mm and mm), with serrated surfaces to avoid slippage at the walls. Poly Methyl Methacrylate (PMMA) is chosen to build the cell to ensure as low X-ray absorption as possible; cast PMMA is preferred to extrude PMMA as it offers the lowest absorption. The whole cell is then machined from a full block of cast PMMA. The cell is mounted on a home-made velocity controlled rheometer allowing for 3D tomography since it is not situated on the X-ray path: A rotation stage (Newport RGV100BL) that has an angular resolution of , a maximum torque of , a maximum normal force of , and a maximum velocity of 120 rpm is put below the cell and is used to rotate the outer cylinder for a fixed inner cylinder, thus shearing the material. The whole set-up (rotation stage + Couette cell) is rigidly fixed to the microtomography rotating stage, allowing for 3D tomography of the sample (Fig. 2).
The setup is mounted on an X-ray CT scanner Ultratom model by RX-Solutions having 230 kV micro-focus X-ray tube and Paxscan Varian 2520 V flat panel detector (1920× 1560 pixels, each pixel size 127 μm) with spatial resolution down to . During the rotation of the outer cylinder, a real time high resolution digital radioscopy (0.1 s per X-ray image of 25 micron resolution) is performed to get the time-resolved 2D images from which the volume fraction field evolution in time will be computed. It is noteworthy to mention that the frame rate can be larger than 10 Hz. At the end of each experiment, rotation is stopped, and CT-scan of the object at rest is performed with a resolution of 14 μm voxel (volumetric pixel) to get a 3D image of the sample and to identify all particles' location, and thus reconstruct the exact final volume fraction field.
We have designed two model suspensions with Newtonian and yield stress suspending fluids in this study. The model suspensions are density matched, meaning that the fluid phase has the same density as the solid phase. The model suspensions consist of radiopaque suspending fluids and radiolucent suspended particles. The preparation and the rheological properties of the model suspension with a yield stress background fluid are as follows. We use spherical polystyrene particles with an average diameter and a density of (Dynoseeds TS 140, Microbeads SA). The suspending yield stress fluid is a concentrated water-in-oil emulsion. The volume fraction of the dispersed phase is 70%, and the diameters of droplets are approximately 1 μm (measured by the ZEISS-Lab.A1 microscope) with polydispersity of 20% [41]. The emulsion is prepared by dispersing a sodium iodide (NaI) water solution ( NaI in water) in a solution of Span 80 emulsifier in dodecane oil ( Span 80 in dodecane oil) with a Silverson L4RT emulsifier, as in [3,4,13]. The water phase is gradually added to the oil phase while the emulsifier is rotating at 1000 rpm for 20 min, then the rotating speed is increased by 1000 rpm every 20 min until it has reached 6000 rpm. NaI is used here to enhance X-ray absorption and ensure that there is a good contrast between the particles and fluid. In addition, adding NaI allows us to match the density of emulsion to that of the polystyrene particles. The fluid rheological properties are well fitted to a Herschel-Bulkley model () in which (Fig. 3). It has already been shown that this system of particles and emulsion is a yield stress model suspension, and the large particles see the emulsion as a continuum visco-plastic matrix with the bulk rheology that can be explained via classical empirical formulations for non-Brownian suspensions [3,13,42,43]. Therefore, the behavior of a yield stress suspending fluid does not change as we add polystyrene particles and the only impact of adding particles is to increase the values of the yield stress and consistency. In this work, the solid volume fraction for suspensions of particles in yield stress fluids is .
A Newtonian model suspension is also designed to study the shear induced migration of particles suspended in a Newtonian fluid, see Sec. V A. Here, we use rigid spherical particles made of PMMA with a diameter of and a density of . The Newtonian radiopaque suspending liquid is composed of Triton X-100 (), Zinc Chloride (), Water (), and Hydrochloric Acid (HCL) (), with a density of and a viscosity of Pa s at 24 °C [44,45]. This system of Newtonian suspending fluid and PMMA is a model suspension with the viscosity following the Kreiger-Dougherty empirical expression [46]. In this work, the solid volume fraction for suspensions of particles in yield stress fluids is . In Sec. III B, we explain the method that we have employed to calculate the steady-state particle concentration by CT scan.
B. Exact volume fraction field from 3D tomography
3D micro-tomography (CT scan) is performed at rest, after steady state is established, to determine directly the concentration field of the particles from the identification of their location in space. Details on the experiments and on the method developed to identify all particles and their position can be found in [13] and [41]. In [41], an in-depth study of this technique, which is a highly promising technique in studying the microstructure of complex suspensions stable at rest, is presented. This paper is under review for publication and will be available via preprint servers. Briefly, in CT scan tomography, we take X-ray projections of a stationary object at different angles, for example, in this case we need more than 2000 images, which takes about 1 h to build a 3D computed object using Radon transform. Figure 4 shows an image slice at a cross section of the computed object. In each of these slices, we determine the location of the particles following [41]. Since we have the location of each particle, we can select cylindrical volume elements of having , in which d is the diameter of particle, and move this element over the entire cell for , where Ri and Ro are inner and outer radii of the cell, respectively, and H is the height of the cell. This enables us to calculate the particle concentration profile for each volume element by counting the number of particles inside each segment, n, and calculating as follows:
C. Integrated volume fraction field from 2D X-ray radiography
Figure 1 is a schematic that shows the X-ray penetrating the object in a linear path, in which the X-ray is absorbed by the various media in the beam (here, the PMMA Couette cell, the PS beads, the NaI solution droplets, and the oil phase). In the experiments, the camera's response is first calibrated to provide a homogeneous intensity response I0 over all pixels when there is no object between the X-ray source and the detector. When an object is studied by radiography, the detector then provides a field I(y, z) of measured intensity, which follows the Beer-Lambert law [47,48]:
where μi is the attenuation coefficient of the material i, is length of the material i crossed by the beam at position (y, z).
For the problem under study, we can replace the subscripts i = 1, 2, and 3 by “cell,” “beads,” and “fluid.” When the particle volume fraction field ϕ(r, z, t) evolves in time, the lengths and measured at a given position (y, z) do not vary independently: They can be written as and , where is the size of the sheared gap crossed by the beam at (y, z), and Φ(y, z, t) is the averaged particle volume fraction seen by the beam along its path at this position. Finally, the intensity field that we measure during an experiments is
Our goal is to get ϕ(r, z, t). It is easily shown that ϕ(r, z, t) and Φ(y, z, t) are related through (see Sec. IV A)
All the information we need to get the particle volume fraction field is thus encoded in the measured intensity field . Below, we show the procedure we follow to get Φ(y, z, t); Sec. IV will be devoted to all image analysis steps necessary to finally obtain ϕ(r, z, t).
Of course, geometrical functions and can in principle be computed, and all absorption coefficients can either be computed from the knowledge of the material properties or measured independently, which would allow us to extract directly Φ(y, z, t) from . In order to be as accurate as possible we have instead chosen to rely on measurements performed directly within our system in order to remove all parameters from Eq. (3) and finally get Φ(y, z, t). The first step consists in measuring the intensity field when the Couette cell is filled with the pure background fluid only (i.e., ). From Eq. (3)
Then for the suspension, we can write
The next step consists of using the information that the material is initially homogeneous, and was prepared at a known volume fraction . The initial intensity field is
It is worth noting that the initial material homogeneity has been checked from 3D tomography measurements. In the absence of 3D measurement, this initial homogeneity can also be checked from the 2D measurements. Indeed, by measuring the intensity for an empty Couette cell
then, the ratio
has to be constant for all beam paths crossing the suspension if the material is indeed initially homogeneous.
Finally, note that: (i) For long experiments, the source intensity may slowly drift (i.e., I0 might slightly change in time), which is compensated at each time by normalizing the intensity field by its current value in a zone of the image where the X-ray source has crossed no object; (ii) in the experiments and are averaged over 10 s to minimize their contribution to the noise in the analysis. We have averaged over 10 s, which is large enough to minimize the temporal noise and it is small enough to preserve the information. We have performed averaging using different time spans, e.g., 5, 10, 20, and 30 s and 10 s seems to be an optimal choice.
IV. IMAGE ANALYSIS
The work of Johann Radon formed the foundation for one of the important transformations in Mathematics known as the Radon transform. The Radon transform deals with the problem of reconstructing the 3D structure of an object from its projections (i.e., 2D line integrals). Applications of the Radon transform include X-ray radiography, 3D computed tomography, barcode scanners, electron microscopy, optics, etc.
Our objective is to calculate the volume fractions of the particles, ϕ(r, z, t), in the cylindrical Couette cell from an X-ray projection image in the case of a problem with azimuthal invariance. Due to the cylindrical symmetry of the setup, X-ray projections at different angles are identical. Therefore, a single projection is sufficient and contains all the information needed to reconstruct the distribution of ϕ(r, z, t) in the cell. The cylindrical symmetry benefits the calculations in two ways. First, it allows us to determine ϕ at each time step because only one projection is needed while also studying the dynamic evolution of the solid volume fraction. Without this symmetry thousands of projections would be needed at each time step therefore increasing time to determine ϕ. Second, we can use a simplified version of the Radon transform called the Abel transform; from the work of Niels Henrik Abel.
A. Abel transform
Here, we explain the Abel transform, its properties and its application in our problem of interest. The Abel transform is the integral transform which maps a physical quantity of a radially symmetric object, here, for example, in Eq. (8), the solid volume fraction ϕ(r, z, t), to a function Φ(y, z, t) with values at each y equal to the line integral of the function ϕ(r, z, t), see Fig. 1
The Abel transform has a well-defined inverse [49]
To test the validity of Eq. (12), one can calculate from Eq. (11) using integration by parts and the Leibniz integral rule to substitute into Eq. (12) and obtain ϕ(r, z, t). A discrete version of Eq. (11) can be formulated as , where is the matrix of the Abel integral operator, u is the vector of the radially symmetric ϕ(r, z, t), and b is the vector of Φ(y, z, t) extracted from a 2D X-ray image. Several problems exist in finding a solution u to the inverse problem.
First, reconstructing a continuous function from a finite number of measurements is an ill-posed problem. However, the discretization of a model and the determination of an unknown function only in a finite number of points make the problem even-determined.
Second, an output signal b is composed of a true value and a noise, . The noise signal originates from a variety of sources such as shot noise, which is due to the quantum nature of light or setup instability [50]. Due to the noise, an output signal is not in the column space of the matrix of the Abel transform , and therefore, a unique solution u does not exist [51]. To handle this difficulty, one can use the pseudo-inverse of to seek the best solution u. This is equivalent to finding an inexact solution u that minimizes the following objective function:
where the first term and second term correspond to the energy (square of the Euclidean length of a vector) of the vector of the error, and the solution, respectively. Here, λ is a regularization parameter that can be chosen as follows; in a general case, λ = 0 in the overdetermined case (i.e., a full rank column matrix or when left inverse of exists), and λ = 1 in the underdetermined case (i.e., a full rank row matrix or when right inverse exists).
Finally, the matrix of the Abel transform is ill-conditioned, meaning that an infinitesimal change or noise in the right hand side or in the coefficient matrix results in a large change in the solution vector u. We can calculate upper bounds of the normalized change in u to show the following (see [51,52] and Sec. IV A 1 for more details):
Here, is the condition number of . The maximum and minimum singular values of are denoted by ηmax and ηmin, respectively. is the pseudoinverse of [51,52]. An ill-conditioned matrix has a large condition number which amplifies a small perturbation of the right hand side () or the matrix coefficient () leading to a wrong solution. Several approaches exist to modify an ill-conditioned matrix to avoid an erroneous solution.
1. Truncated singular value decomposition
Any rectangular matrix, , can be decomposed into , where is the matrix of eigenvectors for is the matrix of eigenvectors for , and is a rectangular matrix with r diagonal elements (ηi) that are the square roots of the nonzero eigenvalues of or [51]. Using the singular value decomposition (SVD), the solution of can be written as
In case of an ill-conditioned matrix with scaled singular values of , some components of u will be singular. In truncated SVD, we eliminate all the terms in Eq. (15) with singular values (ηi) smaller than a threshold ϵ [53]. This is equivalent to increasing the null space of . The truncated SVD resolves singularity of u, however, it is subject to one drawback, we cannot determine solutions in the null space of which are invisible and associated with no output signal (). To overcome this drawback, Eq. (15) can be regularized as follows:
where λ is a regularization parameter. This method is known as Tikhonov regularization [54]. This is also equivalent to the general minimization of the objective function J(u) in Eq. (13) [53,55]. For a generalized Tikhonov Regularization, we minimize the objective function , where Γ is a Tikhonov matrix that can be a filter (e.g., low pass filter) forcing a condition on the solution u [56].
2. Regularization
There are other regularization techniques such as the H1 regularization and the Total Variation (TV) regularization [34–36,57] that have been used to deal with the inversion of the Abel transform (see works of Asaki and co-workers [33,35,58]). In the H1 regularization, the goal is to obtain a nonsingular and smooth solution u via minimizing
The last term includes the L2 norm of the first derivative of the solution u. This is because the smoothness of u can be measured by the energy of its derivative. Therefore, minimization of J(u) in Eq. (17) gives the solution that minimizes the residual error and the derivative of u. The latter forces the solution u to be smooth. The H1 regularization can be regarded as a generalized Tikhonov regularization. It is noteworthy to mention that minimization of the L2 norm terms generate linear differential equations.
The TV regularization is appropriate for situations in which the solution u is discontinuous. This technique retrieves edges or discontinuities without smoothing the solution. To this end, one should minimize the following objective function [35,59]:
Although the L1 norm of the gradient is the appropriate space in retrieving the discontinuities of u, minimization of an L1 norm term results in solving a nonlinear differential equation. Therefore, computational time of the TV regularization is longer than that of the H1 and the Tikhonov. Since, minimization of the L2 norm terms is equivalent to solving a linear differential equation. It is noteworthy to mention that TV regularization is beneficial for cases in which both yielded and unyielded regions coexist in the flow. We will see later in Sec. V that these cases correspond to a discontinuous jump in the solid volume fraction which requires the minimization of the gradient of u.
3. Basis functions
Using Basis functions is another regularization technique to handle the inversion of the Abel transform [49,57]. The idea is to decompose the solution u in terms of a basis function and calculate the associated coefficients. This method reduces the dimension of the problem via disregarding high-frequency terms in the decomposition. Hypothetically, this technique gives similar results as techniques such as the H1 regularization, Tikhonov and low pass filtering (see above for more details).
4. Output denoising
Noise reduction is the process of removing noise from a signal which is popular in image denoising [60,61]. The main objective is to achieve both noise reduction and feature preservation. This technique has also been applied in the Abel inversion problem. A classic approach is to fit a smooth function to the output signal and then directly use the inverse of the Abel transform (12) [62,63]. More recent ways involve signal denoising via using appropriate filters or regularization matrices [64]. Here, selecting an appropriate filter requires a priori knowledge about the type of noise in the output signal. A satisfactory noise reduction gives a right hand side b close to btrue. The upperbound of Eq. (14) can be written as follows [52]:
We have when . Therefore, the upperbound (19) is a small value, meaning that the linear system is well-conditioned for this very particular choice of right hand side .
B. Analyzing radiography images
We seek a solution u for the linear system , where is the matrix of the Abel transform, is the solid volume fraction and is an X-ray projection or an output signal. Our methodology involves using the methods of basis functions Sec. IV A 3 and output denoising Sec. IV A 4. The validity of our method will be tested via comparing the solution u to an exact solution obtained from 3D scanning. Sections IV B 1–IV B 3 explain all the steps in detail.
1. Fourier expansion
We follow [49] and use harmonic basis functions to calculate the solid volume fraction. We expand ϕ(r, z, t) in terms of a Fourier cosine series
Having Eqs. (20)–(23), the process of finding ϕ(r, z, t) is now straight forward. First, we integrate Eq. (23) to calculate . Then we solve for via minimizing
We substitute an(z, t) into Eq. (20) to obtain ϕ(r, z, t). The most important feature of this method is that we use Eq. (11) instead of Eq. (12) to calculate ϕ(r, z, t). Therefore, the method is robust since it is not required to differentiate our noisy X-ray images. However, it is evident that we cannot compute an infinite number of multipliers an(z, t) and an approximation is required. Here, we cut off all frequencies above an Upper-Frequency (UPF) limit. This resembles using a low pass filter or the H1 regularization.
We take the following steps to determine a cutoff frequency that leads to a satisfactory approximation of ϕ(r, z, t). We obtain the exact distribution in the middle of the cell (z = z0) at the final time step (t = tf) from the 3D scanning §III B. We find the coefficients of a polynomial of degree 17 that fully fits the data of . Equation (11) is used to calculate the Abel transform of the analytic polynomial function fitted to . Finally, we follow the basis function approach to calculate the Abel inversion on this data at different cutoff frequencies. Figure 6 shows that we determine in a good agreement with the exact data (i.e., from 3D scanning) for UPF limits beyond 50. The discrepancy between the final data and the exact data is associated with filtering the high frequency terms. It is noteworthy to mention that due to the discontinuity of at the wall of the inner cylinder, Gibbs phenomenon is present and results in very minor overshoots and undershoots that cannot be resolved by increasing the UPF limits [65,66].
We adopt UPF = 50 and apply the basis function analysis to the single X-ray projection at the final time step t = tf to determine . Figures 7(a) and 7(b) show, respectively, the colormap of determined in (r, z) plane and that of the exact computed from the 3D scanning. The colormap of Fig. 7(b) is noisy, although it crudely holds general features of the exact distribution of , e.g., the “S” shape curve close to the outer cylinder. In the following, we improve this method by denoising the X-ray projections ensuring that the upper bound of Eq. (19) decreases and makes the transformation less problematic.
2. Radio noise reduction
Image denoising requires the knowledge about the nature of noise. There are several sources of noise in our X-ray projections such as shot noise, which is an inherent noise due to the nature of light [50] or the possible instability of the setup [67]. Here, we process temporal and spatial signals to obtain noise distributions in the X-ray projections. Figure 8 shows an X-ray projection of the Couette cell filled with the pure fluid. We select vertical lines between the inner and outer cylinders, lines A, B, and C as an example, and plot its pixel's intensity. The mean value of the intensity is almost constant and any change in the mean value of the intensity is associated with the noise (Fig. 8). This is mainly because the fluid properties are constant in this region. Also, the inhomogeneity of the material used to build the setup is negligible.
We obtain the noise distribution by subtracting the intensity from its mean value () and plotting its histogram in Fig. 9(a). Also, we overlap a continuous line in red that shows the best theoretical Probability Distribution Function (PDF) that fits the data, which is a normal distribution with the variance of and zero mean. A Quantile-Quantile plot (Q-Q plot) in Fig. 9(b) provides a graphical view of the discrepancy between the empirical and theoretical PDF of the noise [68], which shows the noise follows a normal distribution in the selected region. We follow this procedure for at least 100 vertical lines between inner and outer cylinders of the Couette cell to calculate the noise variance for each line and we show that the noise distribution depends on the mean value of the image density. In Fig. 10, we show that the noise variance increases linearly with the image intensity, which is the signature of Poisson noise [50,69–71]. In fact, it has been shown that the noise of X-ray images has Poisson nature since the probability that an emitted photon be counted by the detector has Poisson distribution [71].
We adopt the same procedure as above to analyze the temporal noise distribution of the radially located points. Figure 8(c) shows the image intensity of the selected point in time and Fig. 11(a) shows the noise distribution Q-Q plot. All the results for a number of points located in the gap are summarized in Fig. 11(b), which indicates linear dependency between the mean value of the image intensity and the temporal noise variance. We conclude that our system includes the Poisson noise distribution both temporally and spatially, which is likely due to shot noise [50] known as quantum noise as well [72,73]. Therefore, we need to seek filters capable of denoising the Poisson noise.
3. Filters
There are a large number of studies suggesting different filters to denoise X-ray images such as [64,71,73–78]. We have implemented Savitzky-Golay (S-G), Total Variation (TV), and mean filter on our X-ray data. Here, we discuss these filters and adopt the best case for our problem of interest.
Savitzky-Golay (S-G) filter has been employed to filter noisy signals [79] and has also been used to filter X-ray images [76]. Here, we divide the signal to small segments. The length of the segment depends on the signal, but is usually about 5% of the signal length. The method fits a polynomial Pn of degree n to the data points in each segment using least square method. The method is based on the minimization of
where is the number of pixels in the defined segment, bj is the value of each pixel in the segment, and is the desired polynomial of degree n (we used n = 2) [80].
Total variation (TV) filter has been employed to denoise data (for example, [59] and [81]), and also to denoise X-ray images [78,82]. We used the method presented in [83] to apply TV filter on X-ray images. Consider u to be the original data and b to be the noisy data. TV filter is based on the minimization of the following statement:
where is the regulation parameter and p(u) is the penalty function, which contains signal properties (see Sec. IV A). The selection of the regularization parameter and the penalty function is an active topic of research and has been studied in literature [81]. In here, we have used the gradient of the input signal, explained in Sec. IV A 2, as the penalty function.
Mean filter is one of the most widely used filters to denoise temporal and spatial noise of X-ray images as well as Poisson noise [71,75,77,84]. We adopt centered Moving Box Average (MBA) as the mean filter to denoise X-ray images [85,86]. The moving average principal is to define a window (box) that moves over the image in which the value of the center pixel is replaced by the mean value of all the pixels in the box. For example, in y-direction we will have
where yc is the center of the box and 2k is the number of pixels in the box. The same applies for other spatial and temporal dimensions (t and z).
We take the following step to examine the above filters. For the case that we have homogeneous particle distribution (i.e., at the start of the experiment), we calculate Lgap at certain height, here in the middle of the Couette cell z = z0, based on Eq. (7)
based on the definition of Lgap, which is shown in Eq. (29), its Abel inversion, , should result in constant value, c, in the gap (). We should note that due to the serrated surface of the rotating cylinders, we expect imperfection in the calculated value of c
Figure 12(a) shows this inversion on the X-ray data when denoising filters have not been applied. The value is oscillatory and not constant. Figure 12 shows that the results are improved via implementing different filters on the X-ray images before applying the Abel inversion. We see that the MBA filter gives almost constant value in the gap. Therefore, we adopt MBA filter to denoise our X-ray images hereafter having the size ratio of the MBA filter box to the image equal to 4%. It is noteworthy to mention that the oscillatory behavior of the calculated value as is due to the singularity of the Abel inverse transform close to 0. However, we neglect this because this region is outside the gap, i.e., our domain of interests.
In order to validate MBA filter for our data we consider the case where the particles are in the domain as well. We use the Abel transformation of particle concentration profile extracted from tomography (Sec. IV A) as an input signal. We have added three different noise types in order to simulate X-ray noise and to validate MBA filter. Based on the literature, one of the best ways to simulate X-ray noise is to apply Gaussian noise to the input signal [35–37]. In the first case, a Gaussian noise with the variance () equal to 1.5% of the largest value of the signal is imposed [35–37]. The second case is selected to be a Gaussian noise with the variance of the same order of the radiography noise, which is for the average intensity of (Fig. 10). We have also added Poisson noise with the variance proportional to the intensity () to simulate our radiography noise for the third case, in which based on Fig. 10. The left column of Fig. 13 shows these noise types that have been added to the input signal. The right column of Fig. 13 shows that the Abel inversion of denoised signal using MBA filter (solid red line) is in a good agreement with the exact solution (a solid black line).
Furthermore, it is of great importance to validate the filter and procedure for the case that there is a discontinuity in the solid volume fraction, which is likely to happen in the case of non-Newtonian suspensions, such as noncolloidal particles in complex fluid, due to shear banding or shear localization effects [4,13]. In order to simulate this case, we use the Abel transformation of a step function affected by different noise types as the input signal (left column of Fig. 14). The Abel inversion after denoising the input signal by MBA filter is shown in the right column of Fig. 14, which shows a good agreement with the expected result.
V. RESULTS
A. Shear induced migration of particles suspended in a Newtonian fluid
First, we illustrate the interest of our technique on the shear induced migration of particles suspended in a Newtonian fluid, which is a well studied problem. We also compare our results with the prediction of the Suspension Balance Model (SBM) [87,88]. The Newtonian model suspension (Sec. III A) with is loaded in a wide gap Couette cell. The dimensions of the cell are as follows: . The inner cylinder is rotating with the speed of . Figure 15 shows the colormaps of the solid volume fractions at different times. Over time, the particles migrate from the inner cylinder to the outer cylinder, and a steady state is achieved after about, 3000 s or equivalently, 2400 revolutions. The steady state time is consistent with the experimental result of Corbett et al. [25] explained in Sec. V B, taking into the account the difference in the size of suspended particles. The 1D assumption of the evolution of ϕ is valid for , whereas for and , the concentration field is complex and 2D due to the effects of the top and bottom boundaries. The proposed X-ray technique allows us to visualize the boundary effects which are not addressed in previous studies due to the experimental limitations, i.e., insufficient temporal and spatial resolutions (see Sec. I).
Here, we compare our experimental measurements to predictions obtained from the SBM, solved using the rheological laws of Morris and Boulay [87] and Boyer et al. [89]. The prediction of both rheological laws was almost the same, thus we just included the comparison with the Morris and Boulay law [87]. In the literature, Couette flow is always considered as a 1D problem (as expected for flow between infinite cylinders) and all tests of SBM in this geometry are 1D. Therefore, we need to select a part of our domain in which the variation of ϕ is one dimensional to be able to perform any comparison with the SBM model. Figure 16 illustrates the radially averaged solid volume fraction of 2D concentration fields shown in Fig. 15, , vs z. The boundary effects are negligible when , giving an almost constant along z. Therefore, we averaged over z for and compared it with the prediction of the 1D SBM model (Fig. 17). We have solved the equations of SBM with multipliers similar to that of those given by Morris [87] with the maximum packing fraction of the particles equal to 0.585. Figure 17 shows a general agreement between the experiments and the prediction of SBM model, however, there exists a discrepancy, commented on below.
The SBM model attributes the shear induced migration of particles to the spatial variations of the particle phase stress which includes contact or interparticle contributions as well as hydrodynamic contributions, see [90] and [91]. Therefore, SBM requires constitutive relations for the particle phase stress which, cannot be easily obtained experimentally, see [89]. The available modeling is one dimensional (as expected between infinitely long cylinders). Available experimental measurements in the literature are also 1D (see [8,23,92–97]), mostly because of experimental limitations, assuming that 2D effects are negligible. However, our concentration field is clearly 2D and complex. This is likely due to the top and bottom boundary effects. This requires a full modeling in 2D via a tensorial constitutive modeling of the suspension behavior. The proposed technique in this paper is thus a significant advance for properly testing the modeling of the flows of suspensions.
It is also noteworthy to mention that the exact modeling of the shear induced migration via SBM is still under debate [90,91,98] and further studies are required. For example, Snook recently experimentally studied the dynamics of the migration process in an oscillatory pipe flow and showed that SBM predicts a faster rate of migration than what is observed experimentally, see [98]. The value of solid volume fraction measured at the pipe centerline (i.e., where the shear rate vanishes) was found to decrease with decreasing volume fraction while SBM predicts a centerline concentration at maximum packing.
B. Shear induced migration of particles suspended in a yield stress fluid
We now illustrate the interest of our technique on the shear-induced migration of particles in a visco-plastic fluid. Our experimental setup is a wide gap Couette cell (), in which the shear stress (τ) decreases radially and is proportional to the applied torque (T) via . The suspending fluid in our tests follows the Herschel-Bulkley law: It flows only when the shear stress is more than a value called the yield stress (). In order to enforce flow over the whole gap the applied torque should thus be large enough to exert a shear stress larger than the yield stress at any radial position (Fig. 18, case 1). When the applied torque is not large enough, the shear stress is less than the yield stress for r above a given radius Ry, which defines a yield surface. This results in the solid rotation of the fluid at radii larger than Ry (Fig. 18, case 2). The yield torque () and the critical torque to enforce flow over the whole region () are equal to and . If the applied torque T is between these two values, there is a yield surface in the gap at the position .
In our experiments, the suspension of ϕ = 35% particles in the yield stress fluid (presented in Sec. III) is loaded into the Couette cell, while the outer cylinder rotates at Ω = 2 rps (revolution per second) to produce case 1 and at Ω = 0.5 rps to produce case 2. These values are calculated based on the rheology of the suspension [4,13,99].
We have eliminated irreversible effects such as sedimentation or inertia effects; however, we expect to observe shear-induced migration of particles from high to low shear stress regions [8] that is toward the outer cylinder in our Couette cell. Shear-induced migration in wide gap Couette cells has been well studied for Newtonian suspending fluids [93,94,100]. Corbett et al. [25] have, in particular, studied shear-induced migration of a 40% suspension of particles of same size as ours (140 μm) in a Newtonian fluid, in a geometry of same stress inhomogeneity as our (cup to bob radius ratio equal to 2) and of gap size 2 times larger than our; they report that steady state is reached within 2000 revolutions of the inner cylinder. Based on the diffusive scaling of migration in Newtonian fluids [101], we would thus expect steady state to be achieved after around 500 revolutions in our geometry for our particles in a Newtonian fluid. The case of particles in a yield stress fluid is not yet known, but migration strainscale should be of same order of magnitude. A fast imaging technique is thus a priori required to study the kinetics of migration in such situation; our technique then seems to be perfectly adapted.
Here, we present the migration of particles toward the outer cylinder when the suspending fluid is a yield stress fluid. It is noteworthy that in case 2, where there is a plug region, particles will not penetrate the unyielded region and we expect the particles to accumulate at the yield surface. We show the evolution of the particle concentration (ϕ) fields over time using the Abel inversion of a single shot of radiography at each time. We validate our results by comparing the steady state results of radiography with the 3D tomography results.
The evolution of ϕ over time is shown in Fig. 19 for the case 1 (outer cylinder rotating at ). The tomography result at the end of the experiment is shown in Fig. 19(e), to be compared to the radiography result at the same time, shown in Fig. 19(d). The excellent agreement between these two 2D maps (see also below for ) is a strong validation of our technique. The discrepancy error is calculated by L2 norm () over the gap region. We recall that a CT scan requires the object to be in a stationary position and is able to capture only the steady state results due to its low temporal resolution, while the high resolution of radiography allows us to study the kinetics of migration over the transient regime as well as the steady-state regime.
In Fig. 19, we clearly see the particles migrate toward the outer cylinder over time: The volume fraction gets smaller near the inner cylinder and higher near the outer cylinder. It is worth noting that the volume fraction field has a more complex 2D structure than expected: Indeed, the radial migration is accompanied with a slight vertical migration (from bottom to top), which seems to stop near the top surface, at a distance of the order of the gap size. This leads to a maximum of the volume fraction at the outer cylinder at from the top free surface. The pattern shown by the volume fraction field is also affected near the bottom of the inner cylinder, along a height of the order of the gap size. We checked that this is not a gravity effect that would result from a slight density mismatch by changing the density of the emulsion. The complex 2D structure shown by the concentration profiles, with no z-invariant region, illustrates the strong impact of top and bottom boundary effects on migration. In an infinite cell (or with periodic boundary conditions), a z-invariant profile would be expected; in an experimental set-up, boundary conditions cannot be perfect: Velocity has to decrease down to zero at the bottom of the cup, and the top free surface is deformed due to normal stress differences [102]. This naturally leads to z-variations along the whole gap in order to reach a macroscopic equilibrium. These variations would likely be minimized by increasing the height to gap ratio, which is 10 here. These observations further show the necessity of developing techniques that give access to the full spatial concentration field, as the one we present here.
In order to visualize more quantitatively the variations of the volume fraction and get an idea about the timescale of migration, we show the evolution of the radial volume fraction profile ϕ over time in the middle of the Couette cell in Fig. 20. We also show in Fig. 21 the evolution of the volume fraction in time in three specific radial positions shown in Fig. 20. It is striking that the kinetics is apparently much faster near the inner cylinder than near the outer cylinder, which remains to be understood; boundary conditions and z-migration might play a role in this feature.
In the second case (rotating speed of the outer cylinder ), there is a plug region that likely remains intact over time and we observe the migration of particles toward the yield surface. The particle accumulation at the yield surface results in a maximum of particle concentration (Fig. 22) at a well defined position in the gap of the geometry. Figure 22(e) shows the tomography result, in excellent agreement ( error) with the map extracted from 2D radiography [Fig 22(d)]; it confirms the existence of the plug region, and further validates the technique we have developed.
In order to have a better view of this phenomenon we show the evolution of the particle concentration in the middle of the Couette cell in Fig. 23, which shows that at radii larger than Ry (i.e., in the plug region) the particle concentration remains constant. The peak in the concentration at the yield surface (at ) shows that particles are not able to penetrate into the plug region, which results in a sharp discontinuity in the particle concentration that we were able to detect with our technique. A closer view on the kinetics of migration and on the evolution of volume fraction in the yielded region (point A), at the yield surface (point B) and in the unyielded region (point C) is provided in Fig. 24. The thick black line in Fig. 23 shows the tomography results, which is well matched to the final radiography result.
Finally, note that the 2D concentration map and the yield surface show again a complex spatial structure, which further illustrates the impact of the boundary conditions on the overall response of the suspension.
VI. CONCLUSION
We have introduced a new technique to study fast volume fraction changes in suspensions using an X-ray radiography technique. This method has high spatial and temporal resolutions (down to 25 μm and more than 10 Hz, respectively) and is based on the Abel transformation, which enables us to extract 3D information from a single 2D radiography projection of axial symmetric objects. The Abel inversion technique is a well-studied problem and is considered as an ill-conditioned inverse problem. In order to deal with the Abel inversion for our problem of interest we have employed the H1 regularization in conjunction with denoising techniques to minimize the input signal perturbations. We have implemented a variety of filters and have adopted MBA filter.
We have performed experiments in a wide gap Couette cell filled with a suspension of spherical particles in a yield stress fluid and we have been able to offer a new technique to evidence and characterize the shear-induced migration of particles and its kinetics. In near future, we will present an in-depth study of this open problem. We have shown the evolution of the particle concentration over time for two cases, one case where there is a plug region in the cell and the other one where all region is yielded. For both cases, a complete CT scan was performed at steady state for the sake of comparison with the radiography results, providing a strong validation of our technique.
In the presence of a plug region, we manage to accurately describe the particle accumulation at the interface between the sheared and the static regions. Remarkably, even in the absence of sedimentation, the concentration profiles show a complex 2D structure, with no z-invariant region, which illustrates the strong impact of top and bottom boundary effects on migration. Also, we have shown the importance of boundary effects on the shear induced migration of particles suspended in a Newtonian fluid. This further shows the necessity of developing techniques that give access to the full spatial concentration field, as the one we present here.
As compared to standard CT scans, the proposed method reduces down the processing time from an hour to a fraction of a second. It is adapted to the study of most suspensions (which contrasts with optical techniques) and it is much more convenient and rapid than NMR techniques. Therefore, this technique opens a new avenue to study the dynamics of various kinds of suspensions. In near future, we will present methods capable of measuring both velocity fields and volume fraction fields via X-ray.
ACKNOWLEDGMENTS
This work was supported by a grant from Agence Nationale de la Recherche (ANR 2010 JCJC 0905 01-SUSPASEUIL). The Laboratoire Navier microtomograph has been acquired with the financial support of Region Ile-de-France. S.H. acknowledges financial support of NSF (Grant No. CBET-1554044-CAREER) and ACS PRF (Grant No. 55661-DNI9).