Oscillatory Carreau flows in straight channels

The present paper studies the oscillatory flow of Carreau fluid in a channel at different Womersley and Carreau numbers. At high and low Womersley numbers, asymptotic expansions in small parameters, connected with the Womersley number, are developed. For the intermediate Womersley numbers, theoretical bounds for the velocity solution and its gradient, depending on the problem parameters, are proven and explicitly given. It is shown that the Carreau number changes the type of the flow velocity to be closer to the Newtonian velocity corresponding to low or high shear or to have a transitional character between both Newtonian velocities. Some numerical examples for the velocity at different Carreau and Womersley numbers are presented for illustration with respect to the similar Newtonian flow velocity.


Introduction
The non-Newtonian character of blood and other viscoelastic fluids, e.g. polymers, important for some chemical and biochemical engineering applications, are usually described by generalized models of the Newtonian fluids [1]. These models assume the fluids as incompressible and propose a nonlinear dependence of the shear stress on the shear rate, such that the viscosity, which is a constant for the Newtonian fluids, to become a function of the shear rate. For different types of non-Newtonian fluids, this function is empirically determined and represents the rheological model of the fluid. For pseudoplastic or shear-thinning fluids, whose viscosity decreases with the shear rate, the model function is usually a power function (power-law model) or a rational function of the shear rate (Cross model, Carreau model, etc.) [2,3].
The shear-thinning viscosity, for example of blood, is quite well approximated by the Carreau viscosity model, as it has two Newtonian plateaus of constant viscosity at low and high shear rates. These plateaus are connected with a power-law region for the intermediate shear rates. For fluids whose viscosity is described by the Carreau model (Carreau fluids), the dimensionless parameter, Carreau number, is appropriate to be used, which is defined as the ratio between the characteristic shear rate and the transition shear rate [2,4,5]. It is assumed that for low values of the Carreau number the fluid behaviour is localized in the upper Newtonian plateau (low shear rate), while for its higher values the fluid behaviour is essentially defined by the Newtonian flow at the lower viscosity (as shown in figure 1). Usually, in these limiting cases, the flow velocity in straight pipes or channels is similar to that obtained with the Newtonian model based on the higher or lower viscosity. However, this general observation is approximative and depends on the other parameters of the problem. For example, in [6,7], it is shown through simulations in large arteries, that the blood flow, approximated as a Newtonian flow, can lead to errors, particularly in the presence of secondary flows due to curvatures.
Since the blood flow is pulsatile, the problems become additionally complicated if the non-Newtonian character is taken into account [7][8][9]. Apart from the Reynolds number, an additional parameter is introduced for the pulsatile or oscillatory flows, the so-called Womersley number (expressing the ratio of the local pulsating force to the viscous one) [10][11][12]. At low values of the Womersley number, both the Newtonian and Carreau flows in pipes or channels correspond to Poiseuille velocity profiles, while at the high values, they correspond to boundary layers on the walls (usually named Womersley or inertia flow regime). The position of the Poiseuille flow, Womersley flow and flow between these two regime flows is given in figure 1 as a parameter-space diagram of Carreau number and Womersley number. Also, the diagram contains the position of low shear, corresponding to high viscosity, and high shear-low viscosity. Thus, the combination of the Carreau and Womersley numbers can have an interesting influence on the flow characteristics for pulsatile non-Newtonian flows in straight pipe or channel, which is the purpose of the present work.
The pipe flow of Newtonian fluid due to oscillating pressure gradient has been first studied experimentally by Richardson & Tyler [13]. Their observations of the maximum velocity displacement towards the wall is known as the Richardson's annular effect, which is explained also using the analytical solution for the velocity [14,15]. The channel flow has similar behaviour with a slightly different analytical solution as found in [16] and later used by different authors to validate their numerical solutions for non-Newtonian fluid flows, e.g. [17][18][19].
In our previous studies [20][21][22][23][24][25], we have examined the problems of oscillating Carreau blood flow in a straight rigid channel or tube (non-deformable artery) and found numerical solutions for the flow velocity. Also, we have proven that the flow velocity and its gradient are limited from below and above by constants, which depend only on the lower value of the two Newtonian plateau viscosities, amplitude and frequency of the imposed oscillating pressure. Based on these results, we may assume that the blood flow is sufficiently exactly approximated by the Newtonian flow (based on the lower viscosity value) in the larger blood vessels, for example, the abdominal aorta, while in the smaller vessels, like the carotid artery, the non-Newtonian flow character is essential, which is also experimentally accepted [26].
The present paper aims to study more general cases of flows at different Womersley and Carreau numbers and to give solution estimates of the velocity and its gradient for a Carreau flow with respect to the similar Newtonian flow. The high and low Womersley number cases will be studied royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 separately as asymptotic expansions in small parameters, connected with the Womersley number. For the intermediate Womersley numbers, when only numerical solutions can be found, theoretical bounds will be proposed. For their derivation, the comparison principle between sub-and supersolutions of nonlinear uniformly parabolic equations will be used. By means of suitably chosen barrier functions, a priori bounds for the velocity solutions and their gradients, depending on the problem parameters, will be proven and explicitly given. The proven bounds for the Carreau velocity, its gradients and the bounds for the absolute difference between the Newtonian and Carreau velocity solutions will be shown to be valid for every Womersley number, Carreau number and rheological power coefficient n. However, the bound for the absolute difference between the Newtonian and Carreau velocity solutions will be more useful at low values of Carreau number or in the limit n → 1. At high Womersley numbers, it will be shown that the effective Carreau number is responsible for solution type, i.e. if the solution can be approximated with one or the other Newtonian velocity corresponding to low or high shear viscosity or will have a transitional character. The paper is constructed as follows.

Theoretical assumptions 2.1. Shear thinning
There exist different rheological models for non-Newtonian fluids describing the rheology of such fluids, i.e. describing the relation between the stresses and shear rates. It occurs that the Carreau model [3,23,27] is one of the most appropriate models for the so-called shear-thinning fluids: fluids whose viscosities gradually decrease with the increase of the angular deformation rate (shear rate). Moreover, the viscosities reach two limiting values, in the form of two different plateaus, corresponding to the higher, μ 0 , and lower viscosity, μ ∞ , obtained at the low and high shear rate _ g, respectively. Apart from these two viscosities, the Carreau model also includes the characteristic time λ, which is equal to the inverse of the transition shear rate _ g t , l ¼ 1= _ g t , and the power coefficient n, where 0 < n < 1. Then the Carreau viscosity μ c is given by The values of the physical constants for some shear-thinning fluids, whose viscosity is approximated by equation (2.1), and further named as Carreau fluids, can be found in the literature [1,3,21]. Usually, the Carreau fluids are considered as incompressible at isothermal conditions, i.e. with constant density, ρ, that is assumed in our present work.

Oscillatory
A two-dimensional straight infinite channel (−∞ < x < ∞), with width equal to H (0 ≤ y ≤ H), is considered. The flow is supposed as unsteady laminar, driven by an oscillatory pressure gradient A cos ωt along the channel axis Ox, with pulse amplitude A and angular frequency ω, such that ∂p/∂x = −A cos ωt, where p is the pressure. The equations of continuity and motion in vector form are and r @v @t where v = (v x , v y ) is the velocity vector, T ¼ 2m c ( _ g)E is the viscous stress tensor, E is the shear rate tensor and _ g is the shear rate: _ g 2 ¼ 2tr(E 2 ). From the assumption of an infinite channel and from equations (2.2) and (2.3), it follows that ∂v x /∂x = 0, v y = 0 and ∂p/∂y = 0. The only non-zero terms of T are t xy where μ c = μ c (∂v x /∂y), with no-slip boundary conditions along the channel walls, i.e. v x = 0 at both y = 0 and y = H.

Scaling analysis and dimensionless groups
Using H as a characteristic length (y = HY), 1/ω as a characteristic time (t = T/ω) and B as a characteristic velocity (v x = BU), the dimensionless form of equation (2.4), together with (2.1), becomes where -the Womersley number [11], n ∈ (0, 1), c = 1 − (μ ∞ /μ 0 ) (1 ≥ c > 0), μ 0 -the characteristic viscosity and Cu = λB/H-the Carreau number [2]. The dimensionless no-slip boundary conditions on the channel walls are The introduction of the Carreau number is appropriate to distinguish the different cases of shearthinning: Newtonian flow with the higher viscosity at Cu = 0; low shear thinning at Cu ≪ 1; medium shear thinning at Cu ∼ O(1); high shear thinning at Cu ≫ 1. In the next sections, we shall discuss the influence of Cu on the velocity solution for the different cases with respect to β, namely three different cases of β → 0, β ∼ O(1) and β → ∞.

Solutions of the Newtonian flow velocity for different Womersley number cases
From equation (2.5), it is clear that the Womersley number β, significantly changes this equation. For low β, the velocity profiles are Poiseuille like, while for high β, in the limit β → ∞, they are with boundary layer character, but still symmetric with respect to Y = 1/2. As mentioned above, at Cu = 0, the fluid is Newtonian. From equations (2.5) and (2.6), its velocity, further denoted by V(T, Y ), is found explicitly for the first time in the classical book of Landau & Lifshitz [16] and also used in our previous papers [20][21][22][23][24][25] where If the lower viscosity μ ∞ is used as a characteristic viscosity and B 1 ¼ AH 2 m 1 ¼ B=(1 À c) as characteristic velocity, the corresponding Newtonian flow velocity f (T, Y ) satisfies the equation royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 . The solution f(T, Y ) has the same form as (3.1)-(3.3) at β replaced by β ∞ [21]. However, for many polymer fluids, μ ∞ is assumed as zero [1] and then its usage is limited.
In order to give some further estimates of the Carreau fluid velocity in connection with the Newtonian solution (3.1), the latter is developed into series of β → 0 and denoted by V 0 (T, Y ) The Newtonian solution f (T, Y ) (corresponding to μ ∞ ) at β → 0 and if b 1 ¼ b= ffiffiffiffiffiffiffiffiffiffi ffi 1 À c p ! 0, can be developed also as an asymptotic expansion in b 2 1 , similarly to (3.5): where f 00 = V 00 and f 01 = V 01 .
In the limit β → ∞, the solution (3.1) is developed into series of the small parameter 1/β 2 and denoted by V ∞ (T, Y ) The solution above is found to be a uniformly valid solution, after applying the perturbation theory [28]. The first term corresponds to the solution in the interior region (between the two boundary layers), the second and third ones-to the solutions in the boundary layers near the walls Y = 0 and Y = 1, respectively. Moreover, in the limit β → ∞, also β ∞ → ∞ and then the solution f(T, Y ) can be developed into series of 1=b 2 1 . It will have the same form as (3.8) and (3.9), but with β replaced by β ∞ In the general case of a Carreau fluid, at Cu ≠ 0, the velocity satisfying equations (2.5) and (2.6), can be found only numerically. In the limiting case of low Cu → 0, there exists an asymptotic solution, given in [22]. For Cu/β → ∞, the solution of equations (2.5) and (2.6) is found numerically to be close to the solution f (T, Y )/(1 − c). In the following section, we shall present some bounds of the solution U(T, Y ), its derivatives and its reference with respect to the Newtonian flow solution V(T, Y ).

Bounds for the Newtonian velocity and its derivatives
Below we shall list some bounds of the Newtonian solution (3.1), which will be further used to estimate its divergence from the Carreau solution. First, we turn back to equation (2.5) with boundary conditions (2.6) satisfied by V(T, Y ) (at Cu = 0) and where the derivatives with respect to T and Y are denoted as subscripts. Some a priori bounds for the derivatives of in the next lemma (its proof is given in [21]). These bounds are valid for all values of β. and (4:6)

A priori bounds for the Carreau flow
In this subsection, we shall prove some bounds for the Carreau fluid velocity and its gradients. For this purpose, equation (2.5) is rewritten in a non-divergence form in the following manner: with boundary and initial conditions The case c = 1 will not be discussed in the present work. It seems more complicated, because equation (4.7) is not uniformly parabolic and the problem (4.7) and (4.8) has no more a global classical solution for T ∈ [0, ∞) (cf. [29], where the discussed problem corresponds to ours at n = 0), since the gradient of the solution blows up on the boundary and the solution detaches from the boundary data after a finite time. However, the cases of c → 1 are included in the present study.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 Theorem 4.3 (Global gradient bounds). Suppose U(T, Y ) is the solution of (4.7)-(4.9). Then U Y (T, Y ) attains its maximum and minimum on the parabolic boundary Proof. Differentiating (4.7) and (4.8) with respect to Y, we get that U Y (T, Y ) satisfies the boundary value problem From (4.10) and the regularity of U(T, Y ), it follows that the operator P is uniformly parabolic. According to the strong maximum principle for uniformly parabolic equations (see Theorem 2, Section 3 in [30]) U Y (T, Y ) attains its maximum and minimum on Γ and the bound (4.11) holds. □ and : is a supersolution to the boundary value problem (4.7)-(4.9). Indeed, for the operator we have from the choice of K 0 1 and from (4.3) From the comparison principle, we get Analogously, the function ÀH 0 1 (T, Y) is a subsolution to (4.7)-(4.9) and the opposite inequality holds. The bounds (4.14) are consequence of (4.13). □ From (4.11), (4.14) and (4.4), we get the following: is the solution of the Carreau flow problem (4.7)-(4.9). Then the bound However, for c & 1 (c close to 1, but still c < 1), the solution U(T, Y ) of (4.7) is bounded. This means that the constant K 0 1 given by (4.16) is not appropriate to estimate the gradient |U Y (T, Y )|. In order to improve (4.16) for c close to 1, we repeat iteratively the proofs of theorem 4.4 and corollary 4.5, starting with the initial iteration K 0 1 .
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 7 Theorem 4.6 (Improved boundary gradient bounds). Suppose U(T, Y ) is the solution of (4.7)-(4.9). Then the bounds and : Proof. Using (4.16), the function H 1 is a supersolution of (4.7)-(4.9), similarly to the proof of theorem 4.4. Hence from the comparison principle we get Analogously, the function ÀH 1 1 (T, Y) is a subsolution of (4.7)-(4.9) and the opposite inequality holds, which proves (4.17). The bounds (4.18) are consequences of (4.17). □ From the expressions for K 0 1 and K 1 1 given by equations (4.15) and (4.19), it is obvious that K 0 1 ! K 1 1 for all values of c, Cu and n. Moreover, K 1 1 tends to K 0 1 at n → 0 and/or Cu → ∞. Here, we include a plot to show this tendency given in figure 2.
As a consequence of theorems 4.3 and 4.6, we get the following corollary: For c ≥ 0.5, simple computations give the inequality the proof of which is presented in appendix A. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 Repeating the proofs of theorem 4.6 and corollary 4.7, we get inductively the following bounds: : For c ≥ 0.5 and m ≥ 2, similarly to (4.21), the following inequality is inductively obtained The proof of (4.26) is given in appendix B.

Relations between the Carreau flow and Newtonian flow velocity
In this subsection, we shall give bounds of the difference between the Newtonian flow solution and the Carreau flow solution as well as of the difference between their gradients.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 Proof. The function Z(T, Y ) = U(T, Y ) − V(T, Y ) satisfies the boundary value problem where the solution V(T, Y ) of the Newtonian flow is given by (3.1). From lemma 4.2, it follows that Φ(η) is a decreasing function of η. Then from (4.10) and (4.29), we obtain the inequality (1 À n)Cu 2 (K 1 ) 2 : (4:37) The auxiliary function H 2 (T, Y ) = K 2 Y(1 − Y ) is a supersolution for (4.36). Indeed, simple computations give us from (4.5) and (4.37) the inequalities and using the comparison principle, we get Analogously, by means of −H 2 (T, Y ) the opposite bound holds. Bounds (4.34) are a consequence of (4.33). □ 5. Carreau flow velocity estimates at low and high β 5.1. β → 0 In the limit β → 0, similarly to the Newtonian solution (3.5), the Carreau fluid velocity U 0 (T, Y ) is developed into series of β 2 Then U 00 (T, Y ) satisfies the uniformly elliptic equation for every fixed T > 0  it can easily be shown that U 00 (T, Y ) is a concave function of Y for T ∈ ( − (π/2), π/2) and a convex one for T ∈ (π/2, 3π/2). From equations (5.3) and (3.6), it follows that À @ 2 U 00 @Y 2 ! cos T ¼ À Thus from the strong interior maximum principle for elliptic equations, we get and Since in the a priori bounds for U(T, Y ) the barrier function is independent of T (similar to the one in the proof of theorem 4.8), the same bounds hold for U 00 (T, Y ). Thus (4.33) holds also for U 00 (T, Y ) It is seen from equation (5.2) that the solution U 00 strongly depends on Cu. If Cu → 0, since (1 − n) < 1, K 2 → 0 and the solution U 00 is the Newtonian solution V 00 , given with equation (3.6).
Here, we have to note that the coefficient 1/(1 − c) appears because of the relation between the characteristic velocities and because the Carreau dimensional velocity is the same at both characteristic velocities, i.e. BU ≡ B ∞ u, where u is the dimensionless Carreau velocity calculated with μ ∞ as characteristic viscosity [21].

β → ∞
In the limit β → ∞, boundary layer problem is again obtained as in the Newtonian fluid case. Equation (2.5) can be rewritten as The solution U ∞ is sought in the two boundary layers, adjacent to the walls Y = 0 and Y = 1, with thickness ∼O(β −1 ) and in the interior region between them as a perturbation expansion in 1/β 2 , similarly to (3.8) , instead of the parameter Cu in (2.5), it is better to analyse the solution with respect to the ratio Cu/β. Then in the boundary layers and in the interior region the solution strongly depends on the value of the ratio Carreau number Cu to Womersley number β. It occurs that in the limit Cu/β → 0, the solution U ∞0 → V ∞0 . For Cu ≫ 1 and Cu/β ≫ 1, the solution in the boundary layer can be found only numerically or by some approximate methods. At Cu/β ≫ 1, the solution U ∞0 tends numerically to f ∞0 (T, Y ), which is the first term in (3.10), i.e.

Results
In order to illustrate our results, the problem (2.5) and (2.6) has been solved numerically by the Crank-Nickolson method in finite differences. The time interval has been taken long enough (up to 20π for lower Cu and up to 30π for higher Cu), in order to eliminate the influence of the initial condition (4.8). In the following analysis of the solution, the interval T ∈ [18π, 20π] has been considered for Cu < 1000 and T ∈ [28π, 30π] for Cu * 1000.
For the special case of β → 0, the Carreau velocity solution together with the Newtonian solutions V 00 and f 00 /(1 − c) are shown in figure 5 for c = 0.9, n = 0.5 and T = 20π at Cu ≪ 1 and Cu ≫ 1. In these examples β = 0.01, while Cu = 0.01 and Cu = 10 5 , respectively. Then, β ∞ is also small, i.e. b 1 ¼ b= ffiffiffiffiffiffiffiffiffiffi ffi 1 À c p ¼ 0:0316 ( 1, such that the asymptotic expansion (3.7) holds. It is well seen in figure 5 that the Carreau velocity is very close to V 00 at Cu = 0.01 and to f 00 /(1 − c) at Cu = 10 5 .
For high values of β, the solutions V(T, Y ) and f (T, Y )/(1 − c) are presented in figure 6 for c = 0.9, T = 2π and β = 100 close to the left boundary Y = 0 (the plot close to Y = 1 is mirror image of this one). As  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 pointed out above, the boundary layer of f 10 =(8b 2 1 (1 À c)) ¼ f 10 =8b 2 is thinner than that of V ∞0 /8β 2 . These solutions are used for comparison with the Carreau solution. It occurs that at Cu ≤ 100 the Carreau solution is very close to V ∞ , except in some tiny regions of the boundary layers, which are almost invisible and the plot is not presented here. Analysing the solution U ∞ (T, Y ) of (5.5), we arrive to the fact that the function Φ (∂U ∞ /∂Y ) governs the solution form (if equation (5.5) is rewritten similarly to equation (4.7)), since it takes values between 0 and 1 according to lemma 4.2. If the solution U ∞ is substituted by U ∞0 in Φ, then Φ(∂U ∞0 /∂Y ) differs slightly from 1 for a large range of Cu/β ≪ 1. This means that the solution U ∞ (T, Y ) can be approximated by V ∞0 (T, Y ) up to O (1/β 4 ). Thus, it occurs that, if Cu/β ≪ 1, the solution U ∞0 ≈ V ∞0 .  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 In figure 6, the cases of Cu = 10 k (k = 3, 5, 7) at β = 100, i.e. Cu/β > O(1), are also considered. From the plots, it can be concluded that with the increase of Cu the solution U ∞0 goes closer to the solution f ∞0 instead of to V ∞0 . The calculations have been performed for a longer time interval up to 30π in order to obtain periodicity.
Let us analyse the Newtonian velocity solution f (T, Y ), given by (17) and (18) of [21] (there denoted by v), corresponding to the lower viscosity μ ∞ as characteristic viscosity, and the solution U(T, Y ) of equations (4.7) and (4.8). Since the dimensional velocities are the same, then from equation (27) of [21], the function U(T, Y ) satisfies the bound while the bound given by equation (4.33) with (4.35) concerns the difference between U(T, Y ) and V(T, Y ). From the two bounds, it follows that U(T, Y ) is close to both Newtonian solutions, in the limit c → 0, which is evident since μ ∞ ≈ μ 0 , i.e. there exists only one Newtonian velocity solution. However, in the limit c → 1, but still c ≠ 1, the bound (7.1) is not appropriate, i.e. does not give any valuable information for the relation between the Carreau and Newtonian velocity solution. In this respect, the bound (4.33) is more suitable as it depends on the other parameters. In figure 7, the bound K 2 is plotted for different values of Cu and n, at c = 0.999. The line c/2(1 − c) is also plotted for comparison with K 2 , which is an indication that it is not possible to regard the solution difference using only the parameter c. It is evident that K 2 strongly depends on n and Cu at fixed c: decreasing with n and increasing with Cu. The behaviour of K 2 for other values of c is similar, as increasing with c. Here, it must be noted that the bound K 2 is only a qualitative measure of the solution difference.  the channel wall (Y = 0 or Y = 1) are very important when analysing the wall shear stresses (WSS), and deciding with which Newtonian solution the flow can be eventually approximated. For example, at c = 0.9 the solution U(T, Y ) can be approximated with f (T, Y )/(1 − c), when their gradient difference is small enough, e.g. equal to 0.0149 in the case of Cu = 10 5 . Finally, we could make the following statement, that the effective Carreau number is Cu in the Poiseuille and transition regime and Cu/β in the Womersley flow regime. This leads to the conjecture: basically, the effective Carreau number is responsible for solution type changes, converging to one or the other Newtonian solution. The other parameters: Womersley number β, n and c can only accelerate or delay this convergence process when increasing the effective Carreau number.
In order to support this statement, we reconsider the example cases of our previous work [21] for flows in a channel with width 5 mm. The case of blood, shown in fig. 2a of [21], corresponds to Cu = 1775, c = 0.938, n = 0.357 and β = 0.649. The Carreau velocity profile is closer to the Newtonian velocity at viscosity μ ∞ , i.e. to f (T, Y )/(1 − c). However, the presented case of the polymer solution HEC 0:5% in fig. 2b of [21], corresponding to Cu = 9, c = 0.995, n = 0.5088 and β = 0.327, is closer to the Newtonian velocity with the higher viscosity μ 0 , i.e. to V (T, Y ). These examples show that the big difference in the Carreau number in both cases leads to a big difference of their velocity profiles. The corresponding WSS = μ c B/H|∂U/∂Y| Y=0,Y=1 , which are very important for practical applications, also show the same tendency ( fig. 3a, b in [21]). In the blood case, the obtained peak WSS corresponding to that in a human brachial artery (with diameter 5 mm) is 4.25 Pa, which is close to the experimental limits [31,32]: 3.3 ± 0.7 Pa for an artery of diameter 4.4 ± 0.6 mm. The obtained peak WSS of the Carreau model is slightly higher than the WSS of the Newtonian model calculated with the lower viscosity μ ∞ , which is 4.04 Pa. It is worthwhile to mention that the peak WSS of the Newtonian model, calculated with the higher viscosity μ 0 , is 14.24 Pa, which is far away from the experimental data. This result is very important to support our conjuncture that the fluid can be approximated as a Newtonian fluid with the lower viscosity, if the Carreau number is high enough. In the cases when this viscosity is unknown (hardly to be measured), the flow solution remains non-Newtonian, described by the Carreau model, as done in the present work, or by another appropriate model.
From the obtained results, we can conclude that the flow remains laminar, i.e. the peak Reynolds numbers defined as Re max ¼ r VH=min (m c ) ( Re cr , where V is the maximal mean cross-sectional velocity, corresponding to maximal volume flow rate in time. The experimental observations of Patel & Head [33] show that the approximate value of 1300 may be accepted as the lower critical Reynolds number for steady channel flow. We suspect that for oscillatory channel flow the critical Reynolds number will depend on the Womersley number, as it has been reported for Newtonian flows in tubes [34], and to be higher or around the critical Reynolds number for steady flows. In the two examples of blood and HEC solution flows, as cited above, the obtained peak Re max are 850 and 14.2, respectively.
We expect that the critical Reynolds number, Re cr , of Carreau fluid flows in straight channels will not be very different from 1300. Although that we have not found any experimental confirmation for Re cr of shear-thinning flows in channels, there are many results for circular pipes, that support this idea. For example, the experiments show that the transition to turbulence of shear-thinning flows in pipes may be delayed in comparison to Newtonian fluids [35].

Conclusion
The oscillatory flow of a Carreau fluid in a straight infinite channel has been studied in comparison to the two limiting cases of Newtonian fluids (with higher or lower viscosity). The longitudinal velocity is a solution of a parabolic nonlinear equation, which depends on the Carreau and Womersley numbers. An analysis of the non-linearity of the Carreau problem has been performed with respect to the Womersley number: low, high and intermediate. For the first two cases, asymptotic expansions are proposed for the Carreau flow velocity. Since for the intermediate Womersley numbers the Carreau flow velocity cannot be found in an analytic form, theoretical bounds have been proven with theorems depending on the other parameters of the Carreau viscosity model: Carreau number, Womersley number, power coefficient and the Newtonian viscosity ratio. The theoretical bounds also concern the differences between the Newtonian and Carreau flow velocity and between their gradients on the channel wall but have only a qualitative character. To give some quantitative results for these differences, the Carreau flow velocity problem has been solved numerically. Its solution shows that these differences increase with the Carreau number, for any value of the Womersley number (Reynolds number being in the limits of laminar flow). Therefore, we can state the following royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191305 conjecture: that the effective Carreau number only is responsible for the different behaviour of the velocity solution. For practical purposes, it is very important which one of the Newtonian velocities to be used in the case to simplify the problem, i.e. if it is possible for the Carreau velocity to be approximated with one of the limiting Newtonian velocities, corresponding to the lower and higher viscosity. For example, in the case of a blood flow, which has been studied in our previous works, the Carreau number is high enough for the solution to be approximated by the Newtonian solution corresponding to the lower viscosity μ ∞ , as their peak velocities are very close, and the WSS are in the limits of the experimentally measured ones for the brachial human artery [31,32].
The obtained estimates can serve as an indicator to what extent the considered problem may have one or another asymptotic solution corresponding to a developed flow in a channel. As the asymptotic solutions are given by simple expressions, they can be easily implemented as initial or boundary velocity profiles when solving more complicated problems in complex geometries by professional or home-made software.
The flows in perfectly straight two-dimensional channels are considered in the present work. In fact, small disturbances of the channel width can be added to the model. Then the dimensionless wall position will be given as Y = −ɛg(T ) and Y = 1 + ɛg(T ), where |g(T )| ≤ 1 and ɛ ≪ 1. In this way, the present problem will occur in the zero-th order approximation in ɛ of the more general problem of wall perturbations in time (from elastic or other sources). For example, the proper knowledge of the flow velocity in rigid channels/tubes is a starting point before the introduction of elasticity in the model, as a fluid-structure interaction.
Another further continuation of the present work is to use a more general function G(T ) of the pressure gradient instead of the pure oscillation in equation (2.5). G(T ) must be bounded and smooth enough. The obtained bounds for the Carreau velocity will be similar, but it is necessary to know explicitly the function G(T ). Moreover, for general function G(T ), the solution of the Newtonian velocity cannot be given in a closed analytic form like in the present work for G(T ) = cos T.
Finally, we point out the open problem connected with the special case of c = 1. Then the equation (2.5) is not a uniformly parabolic one and the existence of a classical solution is questionable. In this case it is possible for the gradient ∂U/∂Y on the boundaries Y = 0 and Y = 1 to become infinite for some times T. However, for c ∈ [0, 1), the gradient @U=@Y cannot reach infinite values anywhere inside the region Y ∈ [0, 1] and T ≥ 0 according to corollary 4.5. Our conjecture is that at c = 1 and Cu → 0, the gradient is bounded and the treated problem still has a classical solution.
Data accessibility. This work does not have any additional data. Authors' contributions. All authors contributed equally to all aspects of the research and gave their final approval for publication.
Competing interests. The authors declare no competing interests. Funding. The author N.K. has been partially supported by the grant no. BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program (2014-2020).