Non-Newtonian Couette-Poiseuille flow of a dilute gas

The steady state of a dilute gas enclosed between two infinite parallel plates in relative motion and under the action of a uniform body force parallel to the plates is considered. The Bhatnagar-Gross-Krook model kinetic equation is analytically solved for this Couette-Poiseuille flow to first order in the force and for arbitrary values of the Knudsen number associated with the shear rate. This allows us to investigate the influence of the external force on the non-Newtonian properties of the Couette flow. Moreover, the Couette-Poiseuille flow is analyzed when the shear-rate Knudsen number and the scaled force are of the same order and terms up to second order are retained. In this way, the transition from the bimodal temperature profile characteristic of the pure force-driven Poiseuille flow to the parabolic profile characteristic of the pure Couette flow through several intermediate stages in the Couette-Poiseuille flow are described. A critical comparison with the Navier-Stokes solution of the problem is carried out.


1.
Introduction. Two paradigmatic stationary nonequilibrium flows are the plane Couette flow and the Poiseuille flow. In the plane Couette flow the fluid (henceforth assumed to be a dilute gas) is enclosed between two infinite parallel plates in relative motion, as sketched in Fig. 1(a). The walls can be kept at different or equal temperatures but, even if both wall temperatures are the same, viscous heating induces a temperature gradient in the steady state. If the Knudsen number associated with the shear rate is small enough the Navier-Stokes (NS) equations provide a satisfactory description of the Couette flow. On the other hand, as shearing increases, non-Newtonian effects (shear thinning and viscometric properties) and deviations of Fourier's law (generalized thermal conductivity and streamwise heat flux component) become clearly apparent [20]. These nonlinear effects have been derived from the Boltzmann equation for Maxwell molecules [14,28,34,40,57], from the Bhatnagar-Gross-Krook (BGK) kinetic model [6,19,43], and also from generalized hydrodynamic theories [49,51]. A good agreement with computer simulations [16,17,25,30,31,37] has been found. The plane Couette flow has also been analyzed in the context of granular gases [60,65]. In the case of plates at rest but kept at different temperatures, the Couette flow becomes the familiar plane Fourier flow, which also presents interesting properties by itself [4,16,17,24,29,33,40,41,42]. The Poiseuille flow, where a gas is enclosed in a channel or slab and fluid motion is induced by a longitudinal pressure gradient, is a classical problem in kinetic theory [10,35]. Essentially the same type of flow field is generated when the pressure gradient is replaced by the action of a uniform longitudinal body force F = mg x (e.g., gravity), as illustrated in Fig. 1(b). This force-driven Poiseuille flow has received a lot of attention both from theoretical [2,3,15,18,21,32,38,39,44,49,51,55,56,58,59,64,67] and computational [1,12,22,23,27,38,61,62,68] points of view. This interest has been mainly motivated by the fact that the force-driven Poiseuille flow provides a nice example illustrating the limitations of the NS description in the bulk domain (i.e., far away from the boundary layers). In particular, while the NS equations predict a temperature profile with a flat maximum at the center, computer simulations [27] and kinetic theory calculations [55,56] show that it actually has a local minimum at that point.
Obviously, the Couette and Poiseuille flows can be combined to become the Couette-Poiseuille (or Poiseuille-Couette) flow [9,36,50,53]. To the best of our knowledge, all the studies on the Couette-Poiseuille flow assume that the Poiseuille part is driven by a pressure gradient, not by an external force. This paper intends to fill this gap by considering the steady state of a dilute gas enclosed between two infinite parallel plates in relative motion, the particles of the gas being subject to the action of a uniform body force. This Couette-Poiseuille flow is sketched in Fig.  1(c). We will study the problem by the tools of kinetic theory by solving the BGK model for Maxwell molecules.
The aim of this work is two-fold. First, we want to investigate how the fully developed non-Newtonian Couette flow is distorted by the action of the external force. To that end we will assume a finite value of the Knudsen number related to the shear rate and perform a perturbation expansion to first order in the force. As a second objective, we will study how the non-Newtonian force-driven Poiseuille flow is modified by the shearing. This is done by assuming that the shear-rate Knudsen number and the scaled force are of the same order and neglecting terms of third and higher order. In both cases we are interested in the physical properties in the central bulk region of the slab, outside the influence of the boundary layers.

COUETTE-POISEUILLE FLOW 3
The organization of the paper is as follows. The Boltzmann equation for the Couette-Poiseuille flow is presented in Sec. 2. Section 3 deals with the NS description of the problem. The main part of the paper is contained in Sec. 4, where the kinetic theory approach is worked out. Some technical calculations are relegated to Appendix A. The results are graphically presented and discussed in Sec. 5. The paper ends with some concluding remarks in Sec. 6.
2. The Couette-Poiseuille flow. Symmetry properties. Let us consider a dilute monatomic gas enclosed between two infinite parallel plates located at y = ±L/2. The plates are in relative motion with velocities U ± along the x axis and are kept at a common temperature T w . The imposed shear rate is therefore ω = (U + − U − )/L. Besides, an external body force F = mg x, where m is the mass of a particle and g is a constant acceleration, is applied. The geometry of the problem is sketched in Fig. 1(c). In the absence of the external force (g = 0) this problem reduces to the plane Couette flow [see Fig. 1 In the steady state only gradients along the y axis are present and thus the Boltzmann equation becomes where f is the one-particle velocity distribution function and J[f, f ] is the Boltzmann collision operator [7,8,11,13,26,48], whose explicit expression will not be written down here. The notation f (y, v|ω, g) emphasizes the fact that, apart from its spatial and velocity dependencies, the distribution function depends on the independent external parameters ω and g. As said above, g = 0 and ω = 0 correspond to the Couette and Poiseuille flows, respectively. In general, Eq. (1) must be solved subjected to specific boundary conditions, which can be expressed in terms of the kernels K ± (v, v ′ ) defined as follows. When a particle with velocity v ′ hits the wall at y = L/2, the probability of being reemitted with a velocity v within the range dv is K + (v, v ′ )dv; the kernel K − (v, v ′ ) represents the same but at y = −L/2. The boundary conditions are then [13] where Θ(x) is Heaviside's step function. In the case of boundary conditions of complete accommodation with the walls, so that K ± (v, v ′ ) = K ± (v) does not depend on the incoming velocity v ′ , the kernels can be written as where ϕ ∓ (v) represents the probability distribution of a fictitious gas in contact with the system at y = ∓L/2. Equation (3) can then be interpreted as meaning that when a particle hits a wall, it is instantaneously absorbed and replaced by a particle leaving the fictitious bath. Of course, any choice of ϕ ∓ (v) must be consistent with the imposed wall velocities and temperatures, i.e.
In these equations n is the number density, u is the flow velocity, is the peculiar velocity, T is the temperature, k B is the Boltzmann constant, p is the hydrostatic pressure, P ij is the pressure tensor, and q is the heat flux. Taking velocity moments in both sides of Eq. (1) one gets the following exact balance equations ∂ y P yy = 0, (15) ∂ y P xy = mng, (16) ∂ y q y + P xy ∂ y u x = 0.
(17) Henceforth, without loss of generality, we will assume u x (0) = 0. In other words, we will adopt a reference frame solidary with the flow at the midpoint y = 0.
The symmetry properties of the Couette-Poiseuille flow imply the following invariance properties of the velocity distribution function: As a consequence, if χ(y|ω, g) denotes a hydrodynamic variable or a flux, one has χ(y|ω, g) = S g χ(−y|ω, −g) Table 1. Parity factors S g and S ω for the hydrodynamic fields and the fluxes [see Eq. (19)].
where S g = ±1 and S ω = ±1. The parity factors S g and S ω for the non-zero hydrodynamic fields and fluxes are displayed in Table 1. In general, if χ is a moment of the form then S g = (−1) kx+ky and S ω = (−1) ky . The general solution to the stationary Boltzmann equation (1) with the boundary conditions (6) can be split into two parts [45,46,47]: Here, f H represents the hydrodynamic, Hilbert-class, or normal contribution to the distribution function. This means that f H depends on space only through a functional dependence on the hydrodynamic fields, i.e.
The contribution f B represents the boundary-layer correction to f H , so that f = f H +f B verifies the specified boundary conditions. The correction f B is appreciably different from zero only in a thin layer (the so-called boundary layer or Knudsen layer), adjacent to the plates, of thickness of the order of the mean free path. Consequently, if the separation L between the plates is much larger than the characteristic mean free path, there exists a well defined bulk region where the boundary correction vanishes and the distribution function is fully given by its hydrodynamic part.
In the boundary layers the hydrodynamic profiles are much less smooth than in the bulk domain. The values of the flow velocity near the walls are different from the velocity of the plates (velocity slip phenomenon), i.e. u x (y = ±L/2) = U ± . Besides, the extrapolation of the velocity profile in the bulk to the boundaries, u x,H (y = ±L/2), is also different from both the actual values u x (y = ±L/2) and the wall velocities U ± . Of course, an analogous temperature jump effect takes place with the temperature profile. The boundary contribution f B for small Knudsen numbers has been analyzed elsewhere [8,47].
In the remaining of this paper we will focus on the hydrodynamic part f H (and will drop the subscript H), with special emphasis on the corresponding hydrodynamic contributions to the momentum and heat fluxes. In order to nondimensionalize the problem, we choose quantities evaluated at the central plane y = 0 as units: In the above equations we have found it convenient to introduce the dimensionless scaled spatial variable where ν(y) is an effective collision frequency. For the sake of concreteness, we choose it as where η is the NS shear viscosity. The change from the boundary-imposed shear rate ω to the reduced local shear rate a is motivated by our goal of focusing on the central bulk region of the system, outside the boundary layers. Note that a represents the Knudsen number associated with the velocity gradient at y = 0. Likewise, g * measures the strength of the external field on a particle moving with the thermal velocity along a distance on the order of the mean free path. The relationship (27) can be inverted to yield The invariance properties (18) translate into Given the symmetry properties (30), we can restrict ourselves to a > 0 and g * > 0 without loss of generality.
3. Navier-Stokes description. To gain some insight into the type of fields one can expect in the Couette-Poiseuille flow, it is instructive to analyze the solution provided by the NS level of description. In the geometry of the problem, the NS constitutive equations are P xx = P yy = P zz = p, where η is the shear viscosity, as said above, and κ is the thermal conductivity. Inserting the NS approximate relations (31)-(34) into the exact conservation equations (15)-(17) one gets p = const, where Pr = (5k B /2m)η/κ ≃ 2 3 is the Prandtl number. In dimensionless form, Eqs. (36) and (37) can be rewritten as For simplicity, let us assume that the particles are Maxwell molecules [7,11,63], so ν(y) ∝ n(y) and ν * (s) = n * (s). In that case, Eqs. (38) and (39) allow for an explicit solution: Here we have applied the Galilean choice u x (0) = 0 and the symmetry property ∂ y T | y=0 = 0. Equation (40) shows that, according to the NS approximation, the velocity field in the Couette-Poiseuille flow is simply the superposition of the (quasi) linear Couette profile and the (quasi) parabolic Poiseuille profile. In the case of the temperature field, however, apart from the (quasi) parabolic Couette profile and the (quasi) quartic Poiseuille profile, a (quasi) cubic coupling term is present. Here we use the term "quasi" because the simple polynomial forms in Eqs. (40) and (41) refer to the scaled variable s. To go back to the real spatial coordinate y one needs to make use of the relationship (27), taking into account that for Maxwell molecules ν ∝ n. Instead of expressing s as a function of y it is more convenient to proceed in the opposite sense by using Eq. (29). Since 1/ν * = T * one simply has For further use, note that, according to Eq. (41), Thus, the NS temperature profile presents a maximum at the midpoint y * = 0. Before closing this section, let us write the pressure tensor and the heat flux profiles provided by the NS description: q * y (s|a, g * ) = s a 2 − ag * s + 4. Kinetic theory description. Perturbation solution. Now we want to get the hydrodynamic and flux profiles in the bulk domain of the system from a purely kinetic approach, i.e., without assuming a priori the applicability of the NS constitutive equations. To that end, instead of considering the detailed Boltzmann operator J[f, f ] we will make use of the celebrated BGK kinetic model [5,7,66].
where ν is the effective collision frequency defined by Eq. (28) and is the local equilibrium distribution function. In terms of the dimensionless variables introduced in Eqs. (23)-(27), Eq. (48) can be rewritten as Its formal solution is The formal character of the solution (51) is due to the fact that f * appears on the right-hand side explicitly and also implicitly through M * and ν * . The solvability (or consistency) conditions are Let us assume now that g * is a small parameter so the solution to Eq. (50) can be expanded as Likewise, where χ * denotes a generic velocity moment of f * . The expansions of n * , u * , and T * induce the corresponding expansion of M * . The expansion in powers of g * allows the iterative solution of Eq. (51) by a scheme similar to that followed in Ref. [54] in the case of an external force normal to the plates. where These are just the equations corresponding to the pure Couette flow. The complete solution has been obtained elsewhere [6,20,25] and so here we only quote the final results. The hydrodynamic profiles are where the dimensionless parameter γ(a) is a nonlinear function of the reduced shear rate a given implicitly through the equation [6,20] where the mathematical functions F r (x) are defined by K 0 (x) being the zeroth-order modified Bessel function. Equation (62) clearly shows that F r (x) has an essential singularity at x = 0 and thus its expansion in powers of x, is asymptotic and not convergent. However, the series representation (63) is Borel summable [6,25], the corresponding integral representation being given by Eq. (62). The functions F r (x) with r ≥ 3 can be easily expressed in terms of F 0 (x), F 1 (x), and F 2 (x) as It is interesting to compare the hydrodynamic profiles with the results obtained from the Boltzmann equation at NS order (see Sec. 3). We observe that Eq. (58) agrees with Eq. (35) and Eq. (59) agrees with Eq. (40) for g * = 0. On the other hand, Eq. (41) with g * = 0 differs from Eq. (60), except in the limit of small shear rates, in which case γ(a) ≈ 1 5 a 2 (Note that Pr = 1 in the BGK model). The relevant transport coefficients of the steady Couette flow are obtained from the pressure tensor and the heat flux. They are highly nonlinear functions of the reduced shear rate a given by [6,19,20,30] q * y,0 (s|a) = a 2 F 0 (γ)s.
(71) Notice that, although the temperature gradient is only directed along the y axis (so that there is a response in this direction through q * y ), the shear flow induces a nonzero x component of the heat flux [19,20,30,37]. Furthermore, normal stress differences (absent at NS order) are present. Equations (69) and (71) can be used to identify generalized nonlinear shear viscosity and thermal conductivity coefficients. In general, the velocity moments of degree k of f * 0 are polynomial functions of the spatial variable s of degree k − 2. An explicit expression for the velocity distribution function f * 0 has also been derived [20,25].

4.1.2.
Limit of small shear rates. The coefficient γ(a) characterizing the profile of the zeroth-order temperature T * 0 is a complicated nonlinear function of the reduced shear rate a, as clearly apparent from Eq. (61). Obviously, the zeroth-order pressure tensor and heat flux given by Eqs. (66)-(71) inherit this nonlinear character.
It is illustrative to assume that the reduced shear rate a is small so one can express the quantities of interest as the first few terms in a (Chapman-Enskog) series expansion. From Eqs. (61)-(71) one obtains q * y,0 (s|a) = a 2 1 − The terms of order a 2 , a, and a 2 in Eqs. (72), (76), and (78), respectively, agree with the corresponding NS expressions, Eqs. (41), (45), and (47). On the other hand, as noted above, the normal stress differences (P * xx − P * yy and P * zz − P * yy ) and the streamwise heat flux component q * x reveal non-Newtonian effects of orders a 2 and a 3 , respectively. where and we have already specialized to Maxwell molecules, so that ν * = p * /T * . In order to apply the consistency conditions (52) in the derivation of the hydrodynamic fields p * 1 , u * x,1 , and T * 1 , it is convenient to define the moments Therefore, the consistency conditions are The evaluation of Φ n1n2n3 is carried out in Appendix A. The firstorder profiles are In the above equations the functions F r are understood to be evaluated at x = γ. As shown in Appendix A, the moment Φ n1n2n3 is a polynomial function of s of degree n 1 + n 2 + n 3 − 1. In particular, the non-zero elements of the first-order pressure tensor are P * ij,1 (s|a) = P xy,1 (a) = 1.

5.1.
Finite shear rates. First order in g * . The nonlinear dependence on the reduced shear rate a of the zeroth-order quantities has been analyzed elsewhere [6,19,20], so that here we focus on the first-order corrections. Figure 2(a) shows the coefficients associated with the hydrodynamic profiles, i.e., p  x,1 (a), and T x,1 , T The normal-stress coefficients P (1) xx,1 (a) and P (1) zz,1 (a) are plotted in Fig. 2(b). Both coefficients vanish in the NS description. Again, the truncated series (112) and (113) are reliable only for a 0.1. We observe that the xx element has a much larger magnitude than the zz element. The other relevant coefficients of the pressure tensor are not plotted because they are identically P (1) yy,1 = 0 and P (1) yy,1 = 1, as a consequence of the exact momentum balance equations (15) and (16). The coefficients associated with the heat flux vector are plotted in Fig. 3. According to the NS equations, q (0) x,1 = q (2) x,1 = q (0) y,1 = 0 and q (2) y,1 = −a, what strongly contrasts with the nonlinear behavior observed in Fig. 3. This is especially dramatic in the case of the streamwise coefficient q (2) x,1 , which deviates from zero even in the limit a → 0. It is interesting to note that, while q (2) x,1 , q (0) y,1 , and q (2) y,1 have definite signs (at least in the interval 0 ≤ a ≤ 1), the coefficient q (0) x,1 changes from negative to positive around a ≃ 0.42.
In order to illustrate how the Couette-flow profiles are distorted by the action of the external force, we will take a = 1 with g * = 0 and g * = 0.1. While the latter value is possibly not small enough as to make the first-order calculations sufficient, our aim here is to highlight the trends to be expected when the fully nonlinear Couette flow coexists with the force-driven Poiseuille flow. The profiles are shown in Fig. 4. In the pure Couette flow (g * = 0) the temperature and pressure profiles are symmetric, while the velocity and heat flux profiles are antisymmetric. The application of the external force breaks these symmetry features since the firstorder terms have symmetry properties opposite to those of the zeroth-order terms, in agreement with the signs of S g in Table 1. As a consequence, the temperature gradient increases across the channel with respect to that of the Couette flow, as shown by Fig. 4(a). The elements of the pressure tensor are no longer uniform but exhibit negative gradients, especially in the case of the normal stress P * xx [see Fig. 4(b)]. An exception is P * yy , which is exactly uniform as a consequence of the momentum balance equation (15). To first order in g * the value of P * yy is the same as in the Couette flow [see Eq. (98)], but this situation changes when terms of order g * 2 are added [40,55,56]. We observe from Fig. 4(c) that the flow velocity (in the reference frame moving with the midplane y = 0) is decreased by the action of the external force. A similar behavior is presented by q * y , what qualitatively correlates with the increase in the temperature gradient observed in Fig. 4(a). Regarding the component q * x , it takes larger values in the Couette-Poiseuille flow than in the pure Couette flow. In both cases (g * = 0 and g * = 0.1) the shearing is so large (a = 1) that the two components of the heat flux have a similar magnitude, i.e., |q * x | ∼ |q * y |. An interesting effect induced by the external field is the existence of a non-zero heat flux at y * = 0, even though the temperature gradient vanishes at that point. More specifically, q * x (0) = q (0) x,1 > 0 and q * y (0) = q (0) y,1 < 0. The sign of the former quantity changes, as noted before, at a ≃ 0.42.

5.2.
Small shear rates. Second order in g * . Thus far, all the results are valid for arbitrary values of the reduced shear rate a but are restricted to first order in the reduced external field g * . One could continue the perturbation scheme devised in Sec. 4 to further orders in g * but the analysis becomes extremely cumbersome if one still wants to keep a arbitrary. Furthermore, the perturbation expansion in powers of g * is expected to be only asymptotic [56]. On the other hand, we can combine the results obtained here [see Eqs. (72)-(78) and (109)-(117)] with those of Ref. [56] to obtain the hydrodynamic and flux profiles to second order in g * , by assuming that both the reduced shear rate and the reduced force are of the same order, i.e., a ∼ g * . The results are p * (s|a, g * ) = 1 − 12 5 ag * s + 6 5 g * 2 s 2 + · · · , u * x (s|a, g * ) = as − q * y (s|a, g * ) = a 2 s − ag * 19 5 + s 2 + 1 3 g * 2 s 3 + · · · .
In these equations the ellipses denote terms that are at least of orders a 3 , g * 3 , a 2 g * or ag * 2 . As for the relationship between y * and s, from Eqs. (27) or (29) we get s = y * + 1 15 a 2 y * 3 − 1 30 ag * y * 2 36 + y * 2 + 1 150 g * 2 y * 3 22 + y * 2 + · · · . (127) Therefore, we can safely replace s by y * in Eqs. (118)-(126). Again, it is instructive to compare Eqs. (118)-(126) against the NS predictions worked out in Sec. 3 (with Pr = 1, in consistency with the BGK value of the Prandtl number). As can be seen from Eqs. (35), (40), (41), and (44)- (47), the NS expressions do not contain terms of order higher than a 2 , g * 2 , or ag * and therefore their comparison with the kinetic-theory results (118)-(126) is not biased. We see that only the NS results for u * x and P * xy are supported by kinetic theory. This does not mean that Newton's law for the shear stress, Eq. (32) is satisfied, since η = p/ν and the hydrostatic pressure p is not actually uniform. As said before, the NS constitutive equations also fail to account for the existence of normal stress differences (typical non-Newtonian effects) as well as of a streamwise component of the heat flux (failure of Fourier's law). Perhaps the most interesting and subtle differences refer to the presence in Eqs. (120) and (126) of the extra terms 19 25 g * 2 s 2 and − 19 5 ag * , respectively, which are absent in their NS counterparts, Eqs. (41) and (47). The extra term in q * y implies that q * y (0) = 0, what represents a clear violation of Fourier's law (34). The term 19 25 g * 2 s 2 present in the temperature field (120) has dramatic consequences on the curvature of the temperature profile at the midpoint y * = 0. From Eq. (120) we get Therefore, while the NS temperature profile presents a maximum at y * = 0 [see Eq. (43)], Eq. (128) shows that the profile has actually a local minimum if a 2 < 19 5 g * 2 . To analyze this feature in more detail, let us rewrite Eq. (120) as Figure 5(a) shows the scaled temperature difference (T * − 1)/g * 2 , as given by Eq.
(129), for a/g * = 0, 1, 19/5, 3, 2 19/5, and 5. In the case a/g * = 0 (pure Poiseuille flow) the temperature profile has a minimum at y * = 0 surrounded by two symmetric maxima at y * = ±2 19/5. When a/g * = 0 (Couette-Poiseuille flow) several possibilities arise. If 0 < a/g * < 19/5 one still has a local minimum at y * = 0 but the two maxima are no longer symmetric: the one with y * < 0 moves to the center, while the one with y * > 0 departs from it. This is represented by the case a/g * = 1 in Fig. 5(a). At a/g * = 19/5, the left maximum and the central minimum merge to become an inflection point of zero slope. Next, in the range 19/5 < a/g * < 2 19/5 the temperature presents a local maximum at y * = 0 followed by a minimum and an absolute maximum, both with y * > 0. This situation is illustrated by the case a/g * = 3 in Fig. 5(a). At a/g * = 2 19/5 the minimum and maximum with y * > 0 merge to create an inflection point of zero slope. Finally, if a/g * > 19/5 [see case a/g * = 5 in Fig. 5(a)] only the central maximum remains and the profile becomes more and more symmetric as a/g * increases. In the limit a/g * → ∞ (or, equivalently, g * → 0) one recovers the pure Couette flow. This rich phenomenology is absent in the case of the NS temperature profile, as shown in Fig.  5(b).
Given the physical interest of Eq. (120) or, equivalently, Eq. (129), it is convenient to rewrite it in real units. This yields where C T = 19 25 in the BGK model, while C T ≃ 1.0153 in the Boltzmann equation for Maxwell molecules [40,44,55]. Equation (130) still holds in the NS description, except that C T = 0.
6. Concluding remarks. In this paper we have studied the stationary Couette-Poiseuille flow of a dilute gas. As illustrated by Fig. 1(c), the gas is enclosed between two infinite parallel plates in relative motion (Couette flow) and at the same time the particles feel the action of a uniform longitudinal force (force-driven Poiseuille flow) along the same direction as the moving plates. Our main goal has been to assess the limitations of the NS description of the problem and highlight the importance of non-Newtonian properties.
In order to get explicit results, the complicated Boltzmann collision operator has been replaced by the mathematically much simpler BGK model with a temperatureindependent collision frequency (Maxwell molecules). The kinetic model has been solved to first order in the reduced force parameter g * for arbitrary values of the reduced shear rate a. Moreover, complementing these results with those obtained in previous works for the pure Poiseuille flow to second order in g * , we have been able to get the solution to second order in both a and g * .
Starting from the pure nonlinear Couette flow, we have studied the influence of a weak external force, as measured by the nonlinear shear-rate dependence of the nine coefficients p x,1 (a), T x,1 (a), q (2) x,1 (a), q (0) y,1 (a), and q (2) x,1 (a). These functions are plotted in Figs. 2 and 3. A more intuitive picture on the distortion produced by the force on the Couette profiles for the hydrodynamic fields (p, u x , and T ), the pressure tensor (P xx , P yy , P zz , and P xy ), and the heat flux (q x and q y ) is provided by Fig. 4.
Complementarily, we have obtained the quantities of interest [cf. Eqs. (118)-(126)] when the shear rate and the force are treated on the same footing, both to second order. This has allowed us to analyze [see Fig. 5(a)] how, by starting from the pure Poiseuille flow, the symmetric bimodal temperature profile is strongly distorted by the shearing until arriving at the symmetric parabola characteristic of the pure Couette flow.
To put the paper in a proper context, it must be stressed that in our treatment we have focused on the nonlinear bulk properties of the Couette-Poiseuille flow, while ignoring linear effects like Knudsen boundary layers. The latter effects, however, may be relevant in applications like slow MEMS flows where the nonlinear effects considered here might be less important.
Considering the great current interest in the force-driven Poiseuille flow as a playground to test hydrodynamic theories and theoretical approaches, we expect that the work presented here may contribute to motivate further studies, both theoretical and computational, on the Couette-Poiseuille flow.