Fluid Velocity Profile Reconstruction for Non-Newtonian Shear Dispersive Flow

An inverse problem associated with the mass transport of a material concentration down a pipe where the flowing non-Newtonian medium has a two-dimensional velocity profile is examined. The problem of determining the two-dimensional fluid velocity profile from temporally varying cross-sectional average concentration measurements at upstream and downstream locations is considered. The special case of a known input upstream concentration with a time zero step, and a strictly decreasing velocity profile is shown to be a well-posed problem. This inverse problem is in general ill-posed and mollification is used to obtain a well conditioned problem.


Introduction
We consider an inverse problem associated with the mass transport of a material concentration down a pipe, where the flowing medium has a twodimensional velocity profile.We have previously investigated an inverse problem associated with shear dispersive flow in Newtonian media [1], and in this paper we consider more general inverse problems for certain non-Newtonian fluids.For viscous flow problems, the variation in fluid velocity over the cross section of the pipe leads to shear dispersion of the solute concentration.Molecular diffusion of the solute may also significantly contribute to the dispersive transport of the solute, however in this paper we will assume that shear dispersion is the dominant dispersive effect.
The approach examined in this paper is to use the shear dispersive properties of the concentration of a material tracer down a pipe to deduce the underlying velocity profile of the medium.An early problem of this nature involved estimating the average flow velocity in water mains by injecting salt into the main and measuring the electrical conductivity of the water at some fixed location downstream of the injection point [2].
The problem considered is that of determining the fluid velocity profile over the cross section of the pipe from temporally varying cross-sectional average concentration measurements at upstream and downstream locations.The determination of the two-dimensional fluid velocity profile is in general a difficult nonlinear ill-posed inverse problem.However, the special case of a strictly decreasing velocity profile and a specified continuous upstream concentration with a time zero step is a well-posed problem.It is assumed that no material tracer is input before the time zero step which is required in the input concentration or one of its derivatives.
The problems considered in this paper are applicable to the perifusion apparatus, an in vitro experimental apparatus used in endocrinology investigations [3], [4], and pituitary effluent sampling experiments associated with blood flow in pipes [5].In these experimental systems the observed temporal architecture of the cellular responses to various stimuli is significantly altered by the measurement process.The inverse problem considered in this paper allows the underlying cellular responses to be unmasked.
In section 2 the type of non-Newtonian fluid is specified.In section 3 the equations modelling the flow of a material tracer concentration down a pipe are given, and the mathematical operator mapping the temporally varying cross-sectional average concentration measurements at upstream and downstream locations is specified.In section 4 it is shown that the inverse of this operator is unbounded on the Hilbert space of square integrable functions, L 2 , and that the inverse problem is ill-posed.The special case of known input concentration with a time zero step and a decreasing velocity profile is also examined and is shown to be a well-posed problem.A related problem of finding the relationship between the rate of shear and the shear stress is shown to be ill-posed.In section 5 the mollification method is used to regularize the ill-posed inverse problem and the stability of the regularization is analysed.In section 6 a simple numerical algorithm is constructed, and it is verified by numerical experiment that stable approximations to the velocity profile can be obtained in the presence of moderate amounts of data noise.In section 7 this numerical algorithm is tested on data including small amounts of molecular diffusion.

Consideration of Non-Newtonian Flow
In our previous investigation [1], the fluid was assumed to be Newtonian, and therefore the fluid velocity profile in the pipe was governed by Poiseuille's equation [6], (p 57 et.seq.).It is natural to consider associated problems for non-Newtonian fluid flow.Our motivation for this generalisation arose from investigating blood flow in endocrinology experiments [5], which is known to behave in a non-Newtonian manner [7].
Suppose that the velocity profile of a viscous, incompressible fluid flowing through a rigid pipe of circular cross section, with a no-slip boundary condition is given by v(r), where r is the circular polar radial coordinate associated with the cylindrical coordinate system so that the x-axis is aligned with the axis of the pipe.The radius of the pipe is R, and the maximum flow velocity is v m .The fluid is assumed to have a suitably low Reynolds number 1 , so that the flow is laminar.The dispersion of matter with high Reynolds number can be analysed by considering a virtual coefficient of diffusion [8].We will assume that the pipe cross-sectional geometry is circular.The theory is applicable to more general pipe geometries [9], however we do not pursue this aspect here.
Many fluids display a marked shear-dependent viscosity.Solutions such as blood have a reduced viscosity when the shear rate is large and are termed pseudoplastic.Conversely some fluids such as concentrated solutions of sugar in water exhibit an increase in viscosity with the rate of shear and are termed dilatant.Materials requiring a finite yield stress before flow can commence are called plastic and are not considered here.For steady fully developed flow in a circular pipe we assume that the rate of shear is a function of the shear stress (τ ) only, that is where r is the distance from the centre of the pipe.Equation ( 1) is termed the velocity gradient or the rate of shear.This assumption excludes rheopectic, thixotropic and viscoelastic materials.For developed flow, the relation between the shear stress in the fluid (τ ) and the radial position (r) is where τ 0 is the mean stress on the fluid at the pipe boundary [10], (p 150 et.seq.).The functional form of f (τ ) and v(r) for the different types of non-Newtonian fluids are indicated schematically in figure 1.It follows that there is less shear dispersion for pseudoplastic fluids than for dilatant fluids.This model of fluid flow is well known, and is discussed further in [12], (chapter 1).

The Mass Transport Equations
The mass transport of the concentration of a material tracer, c(x, r, t), down the pipe is described by the simple transport equation where v(r) is the fluid velocity profile.Because the injection of the tracer into the pipe at x = 0 is assumed to be independent of r, the boundary condition takes the simple form c(0, r, t) = c 0 (t).If the injection of the tracer is radially dependent then unless some extra system information is available, solutions of the inverse problem are not unique and therefore the problem is highly ill-posed.Assuming the initial concentration in the pipe is zero, i.e., c(x, r, 0) = 0, then the solution of (3) is where H denotes the Heaviside function.If radial dependence in the downstream concentration is available, then reconstruction of the velocity profile is a simple well-posed problem [11].However in applications radial concentration measurements are unavailable, and the concentration measurement at the station x = is assayed over the whole pipe.Simple analysis then shows that the cross-sectional average concentration at a fixed station x = is The nature of the smoothing operator defined in ( 4) is dependent on the functional form of v(r).A reasonable physical assumption is to assume that the velocity profile is a strictly decreasing function of the radial distance from the centre of the pipe, i.e., v (r) < 0. With this assumption the maximum fluid velocity is at the centre of the pipe and is denoted by v m .The nature of the smoothing operator can then be seen after the transformation s = t − v(r) to the integral in equation ( 4), so yielding where the kernel has the functional form The parameter a represents the time for fluid travelling at speed v m to reach the station at x = , and use has been made of the identity r = v −1 /(t + a) .It is important to note that the strict monotonicity of the flow is used in rewriting (4) in the form (5) The lower limit of integration in equation ( 5) would be −∞ if we relax the assumption that the tracer is injected into the pipe at t = 0.This change makes equation ( 5) a singular integral equation, altering the character of the problem.If the upstream and downstream concentrations are measured over finite intervals, then since the concentration in the pipe may be nonzero, the downstream measurement may be contaminated with an unknown concentration.Thus the determination of the velocity profile becomes more difficult.

Velocity Profile Determination
The problem of determining fluid velocity profiles has generated a large number of different experimental procedures [12], (chapter 8).Methods used to measure the flow properties can be characterised as invasive or non-invasive, depending on whether the flow properties are disturbed by the measurement process.The techniques employed utilize the electrical, optical, chemical, thermal and kinetic properties of the medium.The approach considered in this paper is to use the shear dispersive properties of the concentration of a material tracer down a pipe to deduce the underlying velocity profile of the medium.
Consider the problem of determining the fluid velocity profile over the cross section of the pipe from temporally varying cross-sectional average concentration measurements at upstream and downstream locations.In terms of equation ( 4) if we measure c 0 and Q, can we determine v.This inverse problem is a nonlinear Fredholm integral equation of the first kind, and is in general highly ill-posed.The ill-posedness of the problem is strongly dependent on the functional form of c 0 .However, if we make the assumption that the velocity profile is a strictly decreasing function of the radial distance from the centre of the pipe, then we can use equation (5) to determine v(r).Observe that equation ( 5) is a linear Volterra integral equation of the first kind, which are well known ill-posed problems [13].However, we are really more interested in finding the velocity profile, v(r), than the function k(t).Therefore on integrating equation ( 5) by parts we obtain the Volterra integral equation of the second kind where with K(0) = 0.Because the variable of interest is K(t), the derivative of the downstream cross-sectional average concentration does not appear in equation (7).It follows from equation (8) and the assumption v (r) < 0, that K(t) is a strictly increasing function such that lim t→∞ K(t) = 1.Volterra equations of the second kind are known to be well-posed problems [13], (p 45 et.seq.).To avoid the degenerate case of equation (7) reducing to a Volterra equation of the first kind we shall assume that c 0 (0) = 0.With c 0 (0) vanishing the conditioning of the Volterra equation of the first kind will be directly related to the smoothness properties of the input function c 0 .However the experimentalist can adjust the input concentration to improve the conditioning of the velocity profile reconstruction.It is obvious from equation ( 7) that the best conditioning is achieved with a c 0 having discontinuities, generating wavefronts of concentration travelling down the pipe.This enables a "layer stripping" type solution to the velocity reconstruction problem.
A unique continuous solution K(t) to equation ( 7) exists if Q and c 0 lie in the space of continuous functions, C[0, ∞] [13], (p 30 et. seq.).The velocity profile can then be determined from by a simple inverse functional relationship (see [14] for another inverse problem in nonlinear optics leading to a similar relation).Thus the following theorem follows.
Theorem 1 The inverse problem of determination of the strictly decreasing velocity profile has a well-posed solution provided that c 0 (0) = 0 and Thus for example if c 0 (t) = c 0 2 is constant, then equation ( 7) reduces to Another useful special case is when the upstream input concentration is the pulse function, c 0 (t) = c 0 (H(t) − H(t − L)), L > 0. Equation ( 7) then reduces to the simple recursive relationship and v can be determined from equation (9).However in a practical application, the computation of the derivative of the upstream concentration profile, c 0 , renders equation ( 7) ill-posed.This is because realistic measurement data can generally only be placed in the function spaces L 2 , or C, and in these function spaces differentiation operators are unbounded.It follows that the operator mapping the measured data to the velocity profile is unbounded.The regularization of this operator mapping is examined in section 5.
A related inverse problem is determining the functional relationship between the rate of shear and the shear stress specified in equation ( 1), which provides a fuller understanding of the constitutive laws underlying the flow mechanism.It follows from equations ( 1) and ( 2) that the functional form of f can be determined by differentiating the constructed velocity profile.Thus the determination of f is an ill-posed problem.The ill-posed nature of the problem is immediately apparent upon differentiating equation (5).The operator mapping the measured data to the velocity profile is now associated with the derivative of both data inputs, and again is unbounded.
The determination of the velocity profile assumes that v(r) is a strictly decreasing function.A similar theory applies to plastic materials which display a finite yield stress.However because plastic materials can be suitably approximated by pseudoplastic materials the theory offers little extra insight.

Problem Regularization
As we have seen in section 4, the ill-posed nature of reconstructing the strictly decreasing velocity profile from noisy input concentration measurements is equivalent to differentiation.Because the measured data lie in function spaces where differentiation operators are unbounded it is therefore necessary to restore continuity with respect to the data to solve the problem.The regularization procedure for solution of equation ( 7) is simply to compute c 0 in a stable manner.It is well known that numerical differentiation can be made a well-posed problem, and there are a number of regularization methods available.Any of these regularization methods can be used to yield stable solutions to equation (7).In this paper we use the method of mollification [15] based on the treatment of Murio [16], in a similar manner to our use in [1].Suppose that due to measurement difficulties an ideal data function g has been corrupted by noise n and is measured as g m .That is where the functions are defined on some interval I = [0, T ], for some T > 0. Let the extension of the data function g m to the interval I δ = [−3δ, T + 3δ] be defined by T < t ≤ T + 3δ (11) and define the Gaussian convolution, or mollification of g by where δ is the radius of mollification.Then if we define the function norm g(t) ∞ = sup t∈I |g(t)|, the following results can be obtained [16], (p 6 et.seq.).
This consistency result shows that J δ g → g as δ → 0.

Theorem 3 With the extended noisy measurement function g
The mollification method provides the differentiation operator with a Lipschitz continuity result provided that the data g m ∈ C, and δ > 0 is fixed.Furthermore as g m − g → 0, δ can be reduced, and the consistency error then decreases provided that g ∈ C 2 .The regularized solution to the inverse problem consists of first mollifying the measured upstream concentration data, which we denote c 0 , by forming J δ c 0 (t), and then solving the well-posed inverse problem (7) with this regularized data.The measured downstream cross-sectional average concentration is denoted by Q m .

Theorem 4
The strictly decreasing velocity profile reconstruction problem, as stated in section 4, with mollified measurement data J δ c 0 , has a wellposed solution provided that c 0 (0) = 0 and c 0 , Q m ∈ C[0, T ].
Proof: Before we use the stability result for the differentiation operator we must first bound the Volterra operator in equation (7).Consider the equation (λI where I is the identity operator, K is a compact linear Volterra integral operator in a space of continuous functions defined on an appropriate interval, and λ > 0 is constant.The associated perturbed mapping equation is ((λ + ∆λ)I − (K + ∆K))(h + ∆h) = g + ∆g where g and ∆g are continuous functions, the operators K and K + ∆K have continuous kernels k and k + ∆ k, and (λ + ∆λ) > 0 is constant.The continuity of g and k will guarantee the existence of a unique continuous solution h to equation ( 13 The |∆λ| term in equation ( 13) is obtained by considering a perturbed mapping equation with ∆K = ∆g = 0 [13], (p 41 et.seq.).Hence the solution h has continuous dependence on the data inputs k, g and λ, and it follows that equation ( 13) is well-posed.The bound in equation (13) indicates that the modulus of continuity of the solution on the measurement depends strongly on the magnitude of λ, which is c 0 (0) in our problem.This implies that unless c 0 (0) is much larger than any measurement noise then well-posedness will be lost.It follows that if we define the family of operators T −1 by T −1 ( k, g) = (λI − K) −1 g, then the mollification of the upstream data function, J δ c 0 , and Murio's stability lemma imply the continuity result It then follows that the solution of equation ( 7), and the computation of the velocity profile v(r) through equation ( 10) can be performed in a stable manner.
Consistency of the mollified problem solution to the exact K can be shown from Murio's consistency lemma provided that c 0 ∈ C 2 [0, T ].

Numerical Method and Results
In practice the upstream and downstream cross-sectional average concentration measurements, c 0 (t) and Q(t + a), can only be measured over some finite interval 0 ≤ t ≤ T .This means we cannot reconstruct v(r) completely, the minimum resolvable velocity is /T .The parameters in this section are taken from the perifusion apparatus [17], [18], they are v m = 6 × 10 −3 m/s, = 2 m, R = 5 × 10 −4 m, T = 3000 s, so that a = 1/3 × 10 3 s.
We now describe a simple numerical algorithm to solve the inverse problem of velocity reconstruction.Consider a uniform time mesh {t i } N i=0 with t 0 = 0, t i = t i−1 + h, 1 ≤ i ≤ N , where h = T /N .Then the kernel K(t) can be approximated by a B-spline of degree n where M + 1 is the cardinality of the B-spline basis {b i } M i=0 .If n = 1 so that the B-spline basis functions are roof functions, and the spline knot points are chosen to coincide with the time mesh points, then collocation of equation ( 7) provides the finite dimensional equation The properties of B-splines and their support imply that this matrix system is lower triangular.If the integral is approximated by the trapezoidal rule then the solution for the coefficients where the prime on the summation signifies that the last value in the summation is to be halved.To ensure the algorithm is well-posed the coefficients {c 0 (t i )} N i=0 must be replaced by their mollified derivatives.This is achieved by constructing {J δ c 0 (t i )} N i=0 and then using central differencing to estimate the numerical derivatives.The extension and mollification of the data functions utilizes equations ( 11) and (12).
One can reason that due to noise, J δ c 0 (0), and J δ Q(t + a) are better approximations to the true data functions c 0 (0) and Q(t + a) respectively.Due to the nature of the smoothing operator in equation ( 5), Q(t + a) is expected to contain less high frequency information than c 0 .Therefore the mollified functions J δ c(0), and J δ Q(t + a) are used in equation ( 14) for improved accuracy.The numerical approximation to K(t) then allows the reconstruction of the velocity profile v(r) through use of equation (9).
The measurement data is simulated through computation with equation (5), the results of which are corrupted with noise n having a normal distribution with zero mean and a standard deviation of 0.05.Because concentrations must be positive, simulated negative concentrations are set to zero.For the numerical results presented in this paper N = 100, so that the sampling time interval is h = 30 s.
The amount of mollification to be used in the numerical reconstruction is dependent on the perceived amount of noise.There are methods for determining the optimum mollification radius [16], however in this paper we employ moderate smoothing with δ = 75.In general the noise levels in the measured data functions c 0 and Q will be different, and hence their appropriate optimal mollification radii will also differ.
The simulated non-Newtonian fluid is pseudoplastic and the velocity profile obeys the power law relation [19] The upstream concentration input is chosen to be the pulse-like Gaussian function The downstream cross-sectional average concentration profile, Q(t + a), is calculated by equation ( 5).These true data functions, the noise corrupted measurements c 0 and Q m , and the mollified signals are shown in figure 2. The directly estimated, the mollified, and the true derivatives of the upstream input concentration c 0 are shown in figure 3. Note the noise amplification of the directly estimated derivative.
The true velocity profile (---) and the reconstructed velocity profile (-) are shown in figure 4. Given the data noise levels, the reconstruction is reasonably accurate.At lower velocities the noise is more significant and the reconstruction is more difficult.The reconstructed velocity is observed to be not uniformly spaced and the resolution of the higher velocity components is much lower.This velocity spacing is given by equation (8) and is /(t + a).If the input data functions are measured according to this spacing, then the reconstructed velocity spacing is uniform (not shown).
Also shown in figure 4 is the reconstructed velocity profile (-) for a dilatant power law fluid (• • •).This power law fluid is simulated through equation ( 15 concentration input, simulated noise, and mollification radius remain unchanged.Thus, even with moderate data noise levels, the reconstructed velocity profiles allow the different non-Newtonian fluids to be distinguished.The relationship between the rate of shear and the shear stress as specified by equations ( 1) and ( 2), can now be found by differentiating the reconstructed velocity profile.This can be found in a stable manner by computing the mollified derivative of the reconstructed velocity profile.Provided the mollification radius is suitably large, this method yields satisfactory results, but reconstruction algorithms that require a strictly decreasing velocity profile yield better results.However this added constraint will result in a more difficult nonlinear inverse problem.

Consideration of Molecular Diffusion
In this section we test the numerical method presented in section 6 on simulated data including small amounts of molecular diffusion.The molecular diffusion of the solute may significantly contribute to the dispersive transport of the solute.A simple finite difference scheme that incorporated shear dispersion and molecular diffusion was used in the data simulation.This simulated data is not corrupted with noise, and the sampling time interval is h = 30 s.The upstream input concentration is the pulse-like Gaussian function in equation ( 16), and the pseudoplastic fluid is specified by the power law relationship in figure 15.
The simulated downstream cross-sectional average concentration is shown in figure 5 for different molecular diffusion coefficients (D).As diffusion increases the concentration peak increases and shifts to the right.There is also a more rapid concentration decrease after the concentration peak.These effects are predominantly due to the diffusive transport of the material tracer transverse to the fluid flow.
The reconstructed velocity profiles are shown in figure 6 with a linear spline fit through the reconstructed points.The jump in the velocity profiles at v m reflect the fact that molecular diffusion allows some of the material tracer to travel at velocities higher than in pure shear dispersive flow.The velocity reconstruction for the diffusion adjusted data underestimates the high velocities down the centre of the pipe.This is due to the high initial concentration gradient across the pipe and the resulting diffusive transport of the tracer towards the pipe boundary.Similarly, after the bulk of the tracer has travelled down the centre of the pipe the tracer is transported towards the centre of the pipe.This results in the overestimation of the small velocities near the pipe boundary.Hence as the amount of diffusion increases, the downstream cross-sectional average concentration measurement contains less information about the fluid velocity profile, and the velocity profile determination becomes more ill-posed.

Discussion
We have considered the problem of determining the two-dimensional fluid velocity profile from temporally varying cross-sectional average concentration measurements at upstream and downstream locations.Assuming a strictly decreasing velocity profile and a specified continuous upstream concentration with a time zero step, then the problem becomes well-posed.It is assumed that no material tracer is input before the time zero step which is required in the input concentration or one of its derivatives.Failure to use this time zero step protocol produces an ill-posed problem.
In contrast if the upstream concentration also has to be measured, then the problem of determining the velocity profile is an ill-posed deconvolution problem with ill-conditioning equivalent to the problem of numerical differentiation.A simple regularization of the ill-posed deconvolution problem based on the mollification method is developed, and it is verified by numerical experiment that stable approximations to the velocity profile can be obtained in the presence of moderate amounts of data noise and small amounts of molecular diffusion.
The contribution of molecular diffusion is an important factor in perifusion experiments.We will include this aspect into more general inverse problems in a subsequent paper.Interestingly for problems in which diffusion transverse to fluid flow is significant the determination of the upstream concentration profile becomes better posed and the determination of the velocity profile becomes more ill-posed.
The velocity profile in the pipe was determined for a tracer in solution system.A related problem is to deduce the velocity profile for the solution with no material tracer injected.One might also wish to consider problems in which the velocity profile to be determined had some time dependency.Problems of this nature are highly ill-posed.
Some interesting inverse problems arise if one considers developing fluid flow.In this case a component of the fluid velocity is into the centre of the pipe, and there is also an acceleration of the fluid down the pipe.Problems include the reconstruction of an upstream tracer concentration and the determination of the flow development characteristics.These are interesting problems for future work.

Figure 1 .
Figure 1.The functional form of the rate of shear, f (τ ), and the fluid velocity profile, v(r), for the different types of non-Newtonian fluids.

Figure 2 .
Figure 2. The true data functions, the noise corrupted measurements, and the mollified signals.(a) The true upstream input concentration profile c 0 (---), the noise corrupted measurement c 0 (+), and the mollified signal J δ c 0 (-).(b) The true downstream crosssectional average concentration profile Q (---), the noise corrupted measurement Q m (+), and the mollified signal J δ Q m .

Figure 4 .
Figure 4.The true velocity profile (---) and the reconstructed velocity profile with a linear spline fit through the reconstructed points (-).The minimum resolvable velocity is 1/1500 m/s.This reconstructed velocity profile is distinguishable from the reconstructed velocity profile (-) for a dilatant power law fluid (• • •).