Thin-Film Flow of an Inhomogeneous Fluid with Density-Dependent Viscosity

In this paper, we study the pressure-driven thin film flow of an inhomogeneous incompressible fluid in which its viscosity depends on the density. The constitutive response of this class of fluids can be derived using a thermodynamical framework put into place to describe the dissipative response of materials where the materials’ stored energy depends on the gradient of the density (Mechanics of Materials, 2006, 38, pp. 233–242). Assuming a small aspect ratio for the channel, we use the lubrication approximation and focus on the leading order problem. We show the mathematical problem reduce to a nonlinear first order partial differential equation (PDE) for the density in which the coefficients are integral operators. The problem is solved numerically and plots that describe the evolution of the density in the fluid domain are displayed. We also show that it was possible to determine an analytical solution of the problem when the boundary data are small perturbations of the homogeneous case. Finally, we use such an analytical solution to validate the numerical scheme.


Introduction
Most bodies are inhomogeneous, but in many of them the inhomogeneity is mild, and they can be approximated as homogeneous bodies. In this study, we consider bodies wherein the inhomogeneity cannot be neglected. Some of these inhomogeneous bodies could be approximated as being incompressible. Before providing a rigorous definition of what is meant by an inhomogeneous fluid, we will define the concept in loose terms. The homogeneity of a body is not decided on the basis of the properties of a body in its present configuration but, rather, is determined by whether there is some configuration wherein its material points are indistinguishable in terms of their response. For example, a homogeneous fluid which is shear thinning, on being sheared could, in its current flowing state, have its property, say viscosity, vary at different particles; but this notwithstanding the fluid is homogeneous if in its reference state, and this might be close to its static state, all particles of the body would respond in an identical fashion. A rigorous definition of the notion of homogeneity can be found in Reference [1], but here we provide an elementary treatment of the issues. Two material points, P 1 and P 2 , belonging to an abstract body, are said to be materially uniform within the context of a purely mechanical purview, if there exist neighborhoods of the particles which, when placed by different placers into a three dimensional Euclidean space, respond in an indistinguishable manner with regard to their mechanical response. If all the particles belonging to the body are pairwise materially uniform, the body is said to be materially uniform. A body is said to be homogeneous if all the material points belonging to the body are materially uniform with respect to a single placer. A body that is not homogeneous is said to be inhomogeneous. Thus, based on the fact that properties of a body vary in a specific configuration, say the current deformed configuration, does not deem the body inhomogeneous. Whether a body is homogeneous or not cannot be decided within the context of an Eulerian specification. However, the statement that one ought to find if there exists some configuration in which all neighborhoods of material points respond in an identical fashion is an impractical definition, as we can never exhaust all possible configurations in which the body might be placed and, in this sense, never be in a position to decide whether the body is homogeneous. We were interested in studying bodies which, in their rest state, are inhomogeneous.
Anand and Rajagopal [2] have shown how inhomogeneous fluids with properties that vary mildly in the mean value may lead to differences in the response between the inhomogeneous body and its homogenized approximation, which can be larger than one order of magnitude. This means that one has to be very careful when approximating a fluid as being homogeneous because small variations of the properties can lead to significant errors in the computation of global and local quantities. Inhomogeneous fluids can be described within the theory, developed by Rajagopal and co-workers [3], that appeals to the notion of a body having multiple natural configurations. This theory allows one to derive the constitutive equations of a large disparate class of dissipative materials [4][5][6][7][8].
The central idea of the theory is the fact that the body can exist in multiple stress free configurations, and its evolution is regulated by thermodynamical considerations. In particular, the evolution of the natural configuration is determined by the maximization of the rate of entropy production which, in the case of a pure mechanical setting (isothermal case), is equivalent to the maximization of the rate of dissipation.
Malek et al. [9] have used this thermodynamical framework to derive the response of inhomogeneous incompressible fluid-like bodies where the stored energy (henceforth quantities with superscript * denote a dimensional quantity) (Helmoltz potential ψ * ) depends on the density ρ * and on the gradient of the density (∇ρ) * and where the rate of dissipation ξ * is a function of the mean normal stress (which we shall denote by the pressure p * ), of the density, and of the symmetric part of the velocity gradient, namely the tensor D * , that is: Using the maximization criterion under the constraint of incompressibility and imposing isotropy for the Helmoltz potential, one can derive the generic form of the Cauchy stress tensor that now depends on the pressure, the density and its gradient, and the symmetric part of the velocity gradient, i.e., (see Reference [9] for all the details): where ψ * (∇ρ) * means differentiation with respect to (∇ρ) * and where µ * has the dimension of a viscosity. The specification of the stored energy function ψ * and of the viscosity µ * , i.e., the way energy is stored in the continuum and the way the system dissipates energy, allows one to determine the response of the material. Examples of inhomogeneous fluids described in Reference [9] are studied in Reference [9,10].
In this paper, we are interested in materials described by a constitutive equation of the form given by Equation (1). In particular, we consider a fluid-like material in which ψ * did not depend on ∇ρ * and µ * = µ * (ρ * ). In this specific case, Equation (1) reduced to: Hence, we consider an inhomogeneous incompressible fluid where the viscosity changes with the density ρ * and where volumes were preserved, i.e., tr (D * ) = 0. In practice, we are interested in materials with mass distribution that is not uniform in the reference configuration but with motion that could be reasonably approximated as isochoric (i.e., volume-preserving). In other words, the density did not vary as long as the particle is fixed.
The paper investigates the pressure-driven flow between parallel plates of a material whose constitutive equation is described by Equation (2) (see Figure 1). The height of the channel is assumed to be far smaller than its length so that the lubrication approximation applies because of the smallness of the aspect ratio. The leading order momentum equation can be solved analytically, producing a velocity field with components that depends on the density in a functional way. Substitution into the mass balance equation finally provides a very particular integro-differential equation for the density that is solved numerically. A numerical code is implemented for such an equation, and the evolution of the density is plotted as a function of time and space. To investigate the validity of the numerical scheme, we finally consider the case in which the density (and, hence, viscosity) is a small perturbation of the homogeneous case (Navier-Stokes system). In this particular situation, we are able to find an explicit analytical solution that could be compared with the one obtained numerically. The comparison shows good agreement, verifying the validity of the numerical scheme.

Fluid
Fluid Figure 1. Sketch of the system.

The Mathematical Model
It is important to recognize that we are not dealing with a fixed set of particles in the flow domain at any instant of time, as particles are constantly entering and leaving the domain. As we are dealing with an inhomogeneous body, we could then have particles in the flow domain that have different response characteristics associated with them and, more importantly, the governing equations would apply to a different set of particles. That is, in order to study inhomogeneous bodies it is necessary to use the governing equations in the Lagrangian form and follow the partlcle and not use an Eulerian approach as the same point in space will be occupied by different particles. However, using a Lagrangian approach is very cumbersome and hence we make an approximation that will allow us to use an Eulerian approach. However, it is important to recognize that this approximation is only valid for a very small class of flows and also for materials whose inhomogeneity is mild. We make the following assumption concerning the particles that enter the flow domain that ensures we can indeed apply the Eulearian form of the governing equations for the different particles that are occupying the flow domain, as though they are the same set of particles. The assumption that we make is rather strong, but it is applicable for a reasonable class of problems. We assume that, at any instant of time, the properties associated with the particles that enter the flow domain at a fixed point at the boundary are identical, that assures that they have the same response characteristics. We further assume that they also satisfy the same boundary condition at x * = 0 at all times. This is a critical assumption since this allows our domain, in which we did not have a fixed set of particles (as different particles entered the domain at x * = 0), to be equivalent to a fixed set of particles, as those that enter the fluid domain satisfy the same conditions at x * = 0. This allows us to approximate the problem in which we, in reality, had a different set of particles, as though they were the same set of particles, using an Eulerian framework.
Let us consider a fluid in which its Cauchy stress is given by Equation (2), where µ * (ρ * ) is a prescribed positive smooth function. The mass balance is expressed as: where v * (x * , t * ) is the velocity field. Denoting the material derivative as: Equation (3) can be rewritten as: We say that the motion is volume-preserving (isochoric) if: for all t ∈ R and x * ∈ Ω * ⊂ R 3 . For an incompressible continuum, Equation (4) yields: i.e., the density of a fixed particle does not change, even though the density can change from particle to particle in the body.
Using the standard procedure in continuum mechanics, we denote the Lagrangian density byρ * , the coordinate of the particles in the reference configuration by X * , and the motion by χ * . If the body is inhomogeneous, then the density varies in the reference configuration, i.e.,ρ * =ρ * (X * ) and also in the current configuration, since: Following an Eulerian approach ρ * (x * , t * ) is an unknown in the problem, even thoughρ * (X * ) is given. Indeed, the field ρ * (x * , t) is determined only when the mapping of χ * (X * , t * ) has been explicitly found, solving the equations of motion. Hence, according to the Eulerian approach, we have to solve Equation (3) coupled with Equation (5) and with the balance of linear momentum which, in the absence of body forces, is given by: Equations (5)-(7) constitute the mathematical formulation of the model to which one must add appropriate initial and boundary conditions. In particular, we consider a velocity field of the form v * (x * , y * ) = (u * (x * , y * ), v * (x * , y * )) and assume that the fluid is confined in a channel of width 2H * and length L * (see Figure 1).
We notice that the problem as specified by Equations (5)- (7) is given in terms of Eulerian quantities, while the inhomogenerity, which is a consequence of the density varying in the reference configuration, is a Lagrangian specification.
Concerning the boundary conditions for the velocity, because of symmetry w.r.t. the y * = 0 plane, we limit ourselves to study the upper part of the channel. We shall write: where u * y * is the derivative w.r.t. y * of the longitudinal component of the velocity. The pressure drop that drives the flow is taken as: where, for simplicity, we have rescaled the outlet pressure to zero and where ∆p * may depend on time.
We also assign the initial condition for the velocity as v * (x * , y * , 0) = v * s (x * , y * ). With regard to the density, we impose the boundary condition: and the initial condition: where ρ * (density at the lateral inlet) and ρ * s (initial density) are given functions. It is convenient to introduce: Looking at the symmetry condition expressed by Equation (8) and Equation (3), we see that ρ * m must satisfy the following first order partial differential equation (PDE): where the boundary and initial conditions are obtained from the compatibility conditions: Therefore ρ m is an auxiliary unknown of the problem that can be determined solving Equations (9) and (10).
We suppose that the aspect ratio of the domain of the fluid is sufficiently small, so that we write: We rescale the variables as follows: where we choose: and where ρ * o , µ * o > 0 are characteristic values for the density and for the viscosity, respectively, and where U * is the mean velocity in the channel. The time t * o represents the characteristic transit time (i.e., the time taken by a fluid particle to cross the entire channel), while p * o comes from the classical Poiseuille law. After little algebra, the system Equations (5)-(7) can be rewritten in the following non-dimensional form: where: is the Reynolds number.

Leading Order
Assuming Re O(1) and neglecting all the terms containing ε, we look for the zero order approximation of the general problem (Equation (11)): It immediately follows that p = p(x, t) and, integrating (12) 3 with the symmetry condition u y (x, 0, t) = 0 and no-slip conditions u(x, 1, t) = 0, we obtain: ds. Now, exploiting the balance of mass u x = −v y , we find: Integrating again between y and 1 and exploiting v(x, 1, t) = 0, we obtain: Next, imposing the symmetry condition v(x, 0, t) = 0, we get: namely: where C may depend only on time. Exploiting the boundary conditions (we have set ∆p = ∆p * /p * o , p(0, t) = ∆p, and p(1, t) = 0), we integrate (14) between x = 0 and x = 1 and get: where µ = µ(ρ(x, s, t)). Relations (15), (16) represent the formal expressions for the velocity components. One immediately sees that the velocity components are operators (and not functions) of the density ρ. With an abuse of notation, we shall write u = u(ρ) and v = v(ρ), recalling that the dependence on ρ is through an integral. In conclusion, we have ended up with a non-standard equation for ρ, i.e., with u = u(ρ) and v = v(ρ) given by (15), (16). This equation can be obviously solved only numerically.
Recalling the discussion of Section 2, the problem to be solved is the following: with the compatibility condition: ρ s (0, y) = ρ (y, 0).
In conclusion, we have to solve a system of two first order PDEs for the unknowns ρ m (x, t) and ρ(x, y, t).
We then use ρ (1) (x, y, t) as the new guess and determine ρ (2) (x, y, t) = Ψ(ρ (1) ). Proceeding in this way, we build a sequence of functions: If the map Ψ has a unique fixed point in some function space, then the unique fixed point is the solution to Equation (17). Here, we do not focus on the proof of the existence and uniqueness of a fixed point, but we use such an iterative procedure to evaluate the solution numerically. In particular, the iterative procedure will be stopped when the quantity: is less than a fixed tolerance. The numerical scheme is implemented through a first-order upwind scheme for both Equations (17) 1 and (17) 4 . A schematic representation of the stencils for the finite difference schemes for ρ m and ρ is shown in Figure 2. Let us focus first on Equation (17) 1−3 , Figure 2a. We build the grid: x i = i∆x t k = k∆t i = 0, . . . , n x , k = 0, . . . ., n t . and we approximate the derivatives in the following way: The first-order upwind explicit scheme is given by: where: We build the 3D grid: x i = i∆x y j = j∆y t k = k∆t i = 0, . . . , n x , j = 0, . . . , n y , k = 0, . . . , n t , and we approximate the derivatives: ∆y .
The first-order upwind explicit scheme is given by: where and where ρ (n) is again the guess at the nth iteration. The solution of Equation (21) provides the solution at the (n + 1)th iteration, that is ρ (n+1) . We have set the tolerance with Equation (19) to be 10 −5 , so that, when tol < 10 −5 , the iterations are stopped and the resulting solution is deemed as the solution to the problem.

Examples
To illustrate the behavior of the solutions of Equation (17), we have selected a pair of boundary conditions: where the compatibility condition ρ (y, 0) = ρ s (0, y) is always satisfied. These functions are taken arbitrarily and do not come from any particular application. They have been chosen just to to illustrate the variation of the density with time and space. In both cases, we have taken the viscosity function as: that is an exponential. The pressure drop is ∆p = 1. We look for a solution in the time interval t ∈ [0, 1] and we solve the problem using the numerical scheme of Equations (20) and (21). In Figure 3 we have plotted the function ρ s and ρ(x, y, 1) for the case (B1), while in Figure 4, the same functions for the case (B2). In both cases, the evolution of the density in the fluid domain from time t = 0 to time t = 1 can be observed. Notice that the boundary conditions are satisfied. Moreover, we observe that, because of the symmetry of the problem w.r.t. the plane y = 0, we have ρ y (x, 0, t) = 0. In the central part of the channel (near y = 0), the density tends to assume the value it has at the inlet. This is due to the fact ρ satisfies a homogeneous transport equation in which the boundary data are "transported" along the characteristics.

Analytical Solution
Recalling Remark 2, we know that ρ = 1 automatically satisfies the system (Equation (17)) and leads to the classical Navier-Stokes system. Suppose we take a "small perturbation" of this very special case, namely: Then, expanding the function µ around ρ = 1, we get (the prime denotes differentiation w.r.t. ρ): Since we are dealing with non-dimensional variables, we can take µ(1) = 1, without loss of generality. Therefore, at the first order, we write: where β = µ (1) and where the o(δ) terms have been neglected. Let us now substitute Equations (22) and (23) into the system of Equation (15): Using Taylor expansion (because of the smallness of δ), we approximate: so that: Using Equation (24) again, we finally find: Neglecting the o(δ) terms in the above, we find: Finally, recalling that: we make use of Equation (16) to find: v = δ∆p (y − 1) 2 (y + 2) 6 · In conclusion, we have obtained the expansion of the velocity field at the first order in δ: with: and: Next, substitute Equation (22) and Equation (25) into Equation (17). The zero order problem, i.e., the one in which all the terms containing δ are neglected, is the classical Navier-Stokes system with constant density. The first order problem becomes: where ζ (y, t) and ζ s (x, y) are obtained, expanding around δ the boundary data ρ (y, t) and ρ s (x, y), respectively, namely: ρ (y, t) = 1 + δζ (y, t) ρ s (x, y) = 1 + δζ s (x, y).
The compatibility of ρ (y, 0) = ρ s (0, y) implies the condition ζ (y, 0) = ζ s (0, y). From Equation (26) 4 , we notice that y plays the role of a parameter for the equation for ζ. Equation (26) 1−3 can be solved by means of the method of characteristics. Indeed, extending ζ s (x, y) for x 0 in the following way: we find that the solution ζ m is given by: Therefore we are left with the problem: Still, exploiting the method of characteristics, we find that the solution of Equation (29) is given by: To verify the validity of the solution, we need to merely substitute Equations (28) and (30) into Equation (26). (27), we notice that the extensionζ s (x, y) of the function ζ s (x, y) to the strip:

Remark 3. From Equation
is "physically meaningful" only if ζ (y, t) is a bounded function of time t. Indeed, the functionζ s (x, y) is defined up to y = 1, implying that, for x < 0, we have to evaluate the limit lim t→∞ ζ (1, t), which, of course, has to be bounded. On the other hand, if ζ does not depend on time, i.e., the density at the inlet does not change with time, the functionζ s (x, y) becomes: and the solution of the problem is still given by (28), (30).
To prove that the numerical scheme is consistent, we compare the solution obtained using the perturbation with the one obtained with the numerical procedure. Towards this aim, we consider the following pair of boundary conditions:        ρ s (x, y) = 1 + δ sin(6x 2 ) + 4 5 cos(4y 2 ) , ρ l (y, t) = 1 + δ cos(4y 2 ) .
We solve the problem using the perturbation of Equation (30) and with the numerical scheme of Equation (21), getting two solutions, ρ p (x, y, t) and ρ u (x, y, t), where ρ p is the analytical solution using the perturbation and ρ u is the solution obtained with the numerical scheme. We plot the difference of θ = ρ p − ρ u for different values of δ. In other words, we show that the perturbation is advantageous when the boundary data are small perturbations of the homogeneous case. Notice that, when δ = 4/5, the numerical solution is exactly the one plotted in Figure 3. In Figures 5 and 6, the function θ (absolute error) is plotted for decreasing values of δ. As one can notice, the error decreases for decreasing δ, as expected.
The good agreement in the comparison suggests that the perturbation is not only useful as an approximation for small density variation around the constant value (homogeneous system), but it also provides some faith that the numerical method is reliable and can be used to infer the density distribution with general boundary data.

Concluding Remarks
We investigated the thin film flow of an inhomogeneous fluid in which its viscosity was a function of the density. The constitutive equation for the class of fluids considered can be derived by assuming that the stored energy that depends on the gradient of the density [9]. We considered a channel flow driven by a given pressure gradient, and we assumed that the aspect ratio between the width and the length of the channel was such that we can use the lubrication approximation. We then focused on the leading order problem. The mathematical problem reduced to a nonlinear first order PDE for the density, where the coefficients are integral operators. The problem was solved numerically and plots that describe the evolution of the density in the fluid domain were displayed.