On the convergence rate of the Kačanov scheme for shear-thinning fluids

We explore the convergence rate of the Kačanov iteration scheme for different models of shear-thinning fluids, including Carreau and power-law type explicit quasi-Newtonian constitutive laws. It is shown that the energy difference contracts along the sequence generated by the iteration. In addition, an a posteriori computable contraction factor is proposed, which improves, on finite-dimensional Galerkin spaces, previously derived bounds on the contraction factor in the context of the power-law model. Significantly, this factor is shown to be independent of the choice of the cut-off parameters whose use was proposed in the literature for the Kačanov iteration applied to the power-law model. Our analytical findings are confirmed by a series of numerical experiments.


Introduction
In this work, we focus on the iterative solution of nonlinear partial differential equations that arise in models of steady flows of incompressible shear-thinning fluids, including models with explicit constitutive relations of Carreau and power-law type. In particular, we consider the following quasi-Newtonian fluid flow problem: find ( , p) such that where Ω ⊂ ℝ d , d ∈ {2, 3} , is a bounded Lipschitz domain, the source term ∈ L 2 (Ω) d is a given external force, is the velocity vector, p denotes the pressure, and e( ) is the d × d rate-of-strain tensor defined by here |e( )| denotes the Frobenius norm of e( ) , and the (real-valued) viscosity coefficient is assumed to satisfy the following structural assumptions: (A1) ∈ C(Ω × ℝ ≥0 ) and it is differentiable in the second variable; (A2) There exist constants m , M > 0 such that (A3) is decreasing in the second variable, i.e., � (x, t) ≤ 0 for all t ≥ 0 and all x ∈ Ω , where ′ denotes the derivative of with respect to the variable t.
The assumption (A3) asserts that the viscosity decreases with increasing strain rate, in line with our assumption that the fluid under consideration is shear-thinning. Moreover, (A2) immediately implies that is bounded from above and below; indeed, by setting s = 0 , we obtain The bounds m and M are, in general, closely related to the infinite and zero shear viscosity plateau, respectively. In the sequel, the dependence of on x ∈ Ω will be suppressed. Upon defining V ∶= { ∈ H 1 0 (Ω) d ∶ ∇ ⋅ = 0} , the weak formulation of (1) is as follows: (1) (2) m (t − s) ≤ (x, t 2 )t − (x, s 2 )s ≤ M (t − s), t ≥ s ≥ 0, x ∈ Ω; for all x ∈ Ω, t ≥ 0.
where e( ) ∶ e( ) denotes the Frobenius inner product of e( ) and e( ) ; we refer to [2] Sect. 2 for more details concerning the weak formulation (4). The space V is endowed with the inner product and the induced norm ||| ||| 2 Ω = ( , ) V , ∈ V . We emphasize that i.e., the norm |||⋅||| Ω is equivalent to the standard norm on H 1 0 (Ω) d ; the first inequality is a special case of Korn's inequality (see, e.g., inequality (1.7) in [17]), while the second can be easily verified by invoking the Cauchy-Schwarz inequality. In particular, V endowed with the inner product of (5) and induced norm |||⋅||| Ω is a Hilbert space.
The weak form (4) of the boundary-value problem under consideration is known to have a unique solution ⋆ ∈ V , which will be shown, nonetheless, in Sect. 2; moreover, this element ⋆ ∈ V is the unique minimiser of the energy functional where Indeed, a straightforward calculation reveals that, for a given ∈ V, where ′ denotes the Gâteaux derivative; we refer to [1,Prop. 2.1] for details. In particular, the weak formulation (4) is the Euler-Lagrange equation for the minimisation of over V.
A prominent iterative solver for the nonlinear problem (4) is Kačanov's scheme, which, in simple terms, fixes the nonlinearity at the previous iterate: for a given n ∈ V find n+1 ∈ V such that where 0 ∈ V is an arbitrary initial guess. Early references concerning this iterative method include [14], where it was used to compute a stationary magnetic field in nonlinear media, and [5], where the convergence of the Kačanov iteration was investigated in the context of Galerkin methods; Fučík, Kratochvíl and Nečas point in their work [5] to pages 369-370 of Michlin's 1966 monograph [15] for a description of the iterative method introduced by Kačanov in [13], in the context of variational methods for plasticity problems. Kačanov's iteration scheme has been, by now, carefully examined; see, e.g., the monographs [16,Sect. 4.5] and [21,Sect. 25.14], or the papers [7,8,10]. More recently, it was shown in the articles [9] and [4] that the energy from (6) contracts along the sequence generated by the Kačanov iteration (8). Indeed, the first of these two papers established the energy contraction for a more general iteration scheme, and the latter focuses on the Kačanov scheme for a 'relaxed p-Poisson problem' involving a truncation of the nonlinearity from below and from above using a pair of positive cut-off parameters − and + . The derived upper bound on the contraction factor depends on the quotient m ∕M involving − and + , and may be extremely close to 1 in certain situations; interestingly, this unfavourable predicted dependence of the contraction factor on the ratio m ∕M has not been observed in numerical experiments. It is this mismatch between the observed behaviour of the method and the rather more pessimistic results of the analysis reported in the literature that motivated the work outlined herein. We will establish an improved upper bound on the contraction factor of the Kačanov iteration for a general class of shear-thinning fluids. The resulting bound will then be further examined for fluids obeying either the Carreau law or a relaxed power-law, to be specified in the lines below. It will be shown that for (finite-dimensional) Galerkin approximations of the relaxed power-law model it is the power-law exponent, rather than the ratio m ∕M , that is responsible for the rate of convergence of the iteration. Specifically, we will show that the contraction factor of the iteration on finite-dimensional spaces is independent of the choice of the lower and upper cut-off parameters featuring in the so-called relaxed Kačanov iteration, where a truncation of the power-law nonlinearity from below and above is carried out by means of these two positive truncation parameters. To the best of our knowledge the proof of such a result was an open question in the literature.
The paper is structured as follows. In Sect. 2 we will show that the weak formulation (4) of the problem under consideration has a unique solution, which, in turn, is the unique minimiser of in V. The proof is based on auxiliary results, which will also be decisive for the derivation of the contraction factor in Sect. 3. In Sect. 4 we will perform a series of numerical experiments, which confirm our theoretical results. The paper closes with concluding remarks recorded in Sect. 5.

Existence and uniqueness of the solution
We will show in this section that the weak formulation (4) has a unique solution. To this end we define, for given ∈ V , the linear operator [ ] ∶ V → V ⋆ , where V ⋆ denotes the dual space of V, by and the linear form ∈ V ⋆ by In terms of these, the weak formulation (4) can be restated in the following equivalent form: and the Kačanov iteration (8) takes the form: given 0 ∈ V, By (7) and the definitions of and we further have that We will now show that the operator ↦ [ ]( ) is Lipschitz continuous and strongly monotone, since, in that case, the theory of strongly monotone operators implies that the weak equation (10) has a unique solution ⋆ ∈ V ; see, e.g., [16,Sect. 3.3] or [21,Sect. 25.4]. For the proof of Lipschitz continuity and strong monotonicity of the operator ↦ [ ]( ) we require the following result, which, as well as its proof, is largely borrowed from [2, Lemma 3.1]. However, we place emphasis on sharp bounds, since these will be crucial for our convergence analysis below, leading to the improved factors appearing in (15) and (16).

Proposition 2.2 Let be defined as in
The existence and uniqueness of a solution to the Eq. (4) now follows from the theory of monotone operators, cf. [16,Sect. 3.3] or [21,Sect. 25.4]. ◻ (4) is the Euler-Lagrange equation of the minimisation problem the above proposition yields that ⋆ ∈ V is the unique minimiser of the functional ; we emphasise that is strictly convex thanks to the strong monotonicity (23), cf. [21,Prop. 25.10]. Moreover, the Kačanov scheme (11) is well defined thanks to Proposition 2.2 (a) and the Lax-Milgram theorem.

Energy contraction
In this section, we will show that the energy error, given by ( n ) − ( ⋆ ) , contracts along the sequence { n } generated by the Kačanov iteration (11). To this end, we need an auxiliary result. (9) and (6), respectively, with satisfying (A1)-(A3). Then This result is well-known for the Kačanov iteration in the given setting, and the proof can be found, e.g., in [21,Sect. 25.12] or [16,Sect. 4.5]. However, as it is stated in a slightly different form in those references, and also for the sake of completeness, we will include the proof of this statement nonetheless.

Lemma 3.1 Let and be defined as in
Proof It can be shown that see, e.g., [10, Sect. 5.1], and therefore for any , ∈ V . Hence, by the definition of , cf. (6), we find that which completes the proof of the claim. ◻ Now we are in a position to prove the contraction of the energy along the sequence generated by the Kačanov scheme (11). We note that similar results can be found, e.g., in [9, Thm. 2.1] or [4,Cor. 19]. (6). Then, the energy error contracts along the sequence { n } generated by the Kačanov iteration (11) in the sense that where and (t) = (t 2 )t for t ≥ 0.

Theorem 3.2 Assume that (A1)-(A3) hold and let be defined as in
Proof We largely proceed along the lines of [4]. However, as we want to improve the contraction factor from that reference and, in addition, remove any unknown constants, some non-trivial modifications are necessary in the second part of the proof.
Let us define the real-valued function (t) ∶= ( ⋆ + t( n − ⋆ )) , t ∈ [0, 1] . Then, by invoking the fundamental theorem of calculus, we obtain We will first show that � (t), t ∈ [0, 1] , is increasing. A straightforward calculation reveals that Recall that [ n ]( n+1 )( ) = ( ) for all ∈ V , cf. (11), which leads to and hence, by (24), Next, we will take care of the second summand in (26). As was done in [4], we want to bound this summand in terms of the energy difference ( n ) − ( ⋆ ) . However, in order to improve the contraction factor whilst removing all unknown constants, some modifications to the argument presented in [4] are necessary. For The triangle inequality yields that

and thus
This further implies that, for any s ∈ [0, 1] , we have and consequently, in regard of (30), which immediately implies (29). Combining the equalities (28) and (29)  It is straightforward to verify that the contraction factor is minimal for = 1 ∕2Q(n) , and, in that case, one has that which proves the claim. ◻

Remark 3.3
Since m ≤ (t) ≤ M as well as m ≤ � (t) ≤ M for all t ≥ 0 , we get the following crude uniform bound on the contraction factor: We note that, in the context of the relaxed power-law model, cf. Sect. 3.2, this bound, in principle, coincides with the contraction factor from [4].
We note that the contraction factor (25) is not computable as it involves ⋆ , and the uniform upper bound from Remark 3.3 is rather pessimistic. In the following, we will establish an improved computable bound, up to higher order error terms, for the contraction factor on finite-dimensional subspaces. (6)

Proposition 3.6 Let be defined as in (6), with satisfying (A1)-(A3), and let ⋆ be the unique minimiser of in V; then,
An analogous result holds on any finite-dimensional subspace W ⊂ V, with V replaced by W in the assertion above.

Application to the Carreau model
A widely used model for the flow of incompressible non-Newtonian fluids is the Carreau law, cf. [3]. In that case the viscosity coefficient in (1) is of the form where for shear-thinning fluids, r ∈ (1, 2) , > 0 is the relaxation time, and 0 < ∞ < 0 < ∞ denote the infinite and zero shear rate, respectively. The function from (35) is smooth, decreasing since r ∈ (1, 2) , and satisfies the structural assumption (A2), cf. (2), thanks to the following lemma.

Lemma 3.7
Let r ∈ (1, 2) , > 0 , and 0 < ∞ < 0 < ∞ . Then, the following inequalities hold: Proof Define (t) ∶= (t 2 )t , t ≥ 0 . The mean value theorem yields and thus we need to show that ∞ = inf ≥0 � ( ) and 0 = sup ≥0 � ( ) . A straightforward calculation reveals that �� ( ) ≠ 0 for all ≥ 0 , i.e., ′ has no local extrema in the interval (0, ∞) . Since, in addition, lim →0 � ( ) = 0 and lim →∞ � ( ) = ∞ , the lemma is established. ◻ In particular, we may apply Theorems 3.2 and 3.4 to the Carreau model. In this case, the computable contraction factor from (33) reads as follows, with n ∈ W: here X ∶= { ∈ W 1,r 0 (Ω) d ∶ ∇ ⋅ = 0} and ∈ X ⋆ , where, for shear-thinning fluids, r ∈ (1, 2) . In particular, the viscosity coefficient is given by Clearly, ∶ ℝ >0 → ℝ >0 is neither bounded away from zero nor bounded from above, i.e., (A2) is not satisfied. Therefore, as was proposed in the work [4], we consider a relaxed version of : for 0 < − < + < ∞ we define the viscosity coefficient The function is decreasing, strictly positive, bounded, globally Lipschitz continuous, and satisfies (A2) with it is, furthermore, differentiable at all t ∈ [0, ∞) ⧵ { 2 − , 2 + } and has finite leftand right-derivatives at t = 2 − and t = 2 + , respectively. Hence, even though is not continuously differentiable on [0, ∞) , Theorem 3.2 can, nevertheless, be applied in the given setting. Moreover, in the generic case when the set Ω n S ∶= {x ∈ Ω ∶ |e( n (x))| ∈ { − , + }} , for every n ≥ 0 , has Lebesgue measure zero, the operator from the proof of Theorem 3.4 is Fréchet differentiable at n ∈ W . Thus, in turn, Theorem 3.4 can then also be applied to the relaxed powerlaw model 1 . A simple calculation reveals that the computable contraction factor from (33) can again be bounded; indeed, Moreover, one even has that q A (n) = 1 − 4 −1 (r − 1) if the set {x ∈ Ω ∶ − ≤ |e( n (x))| ≤ + } is of positive Lebesgue measure. We further remark that since r ∈ (1, 2) . This shows that the bound (39) is, for every value r ∈ (1, 2) , sharper than the bound from Remark 3.3. Furthermore, this bound predicts that it is the physical parameter r that affects the convergence rate of the iteration, in the finite-dimensional setting at least, rather than the quotient r−2 + ∕ r−2 − implied by existing bounds on the contraction factor, cf. [4,Cor. 19]. Significantly, the upper bound (r − 1) on the contraction factor appearing of the right-hand side of (39) is independent of the relaxation parameters ± . This is of importance as we are interested in the power-law model (37) and we thus need to let − → 0 and + → ∞ . We note that the existence of a bound independent of ± on the contraction factor of the relaxed Kačanov iteration applied to the power-law model with r ∈ (1, 2) was stated in the infinite-dimensional case as an open problem in [4,Ex. 20].
We further note that the energy functional corresponding to the viscosity from (38) coincides with the energy functional J from [4] up to a constant shift depending on − . To be precise, one has that and thus the results established in [4] may be directly applied in our setting. In particular, this implies that the sequence of unique minimisers ⋆ ∈ V of converges in W 1,r 0 (Ω) d to the unique minimiser ⋆ ∈ X of cf. [4,Cor. 10].

Remark 3.8
The relaxed power-law model could also be solved by using a (damped) Newton method, cf. [10,Prop. 5.3]. However, it is unclear whether and how the convergence rate will deteriorate as − → 0 and + → ∞ . For an application of Newton's method to the power-law model with a different regularisation approach we refer to [12]; however, the convergence rate in relation to the choice of the regularisation parameter is not examined in that work.

Remark 3.9
We emphasise that our analysis does also apply to a variable (measurable) exponent r ∶ Ω → (1, 2) for both the relaxed power-law model as well as the Carreau model. Then, in (39) and (36), respectively, we need to replace 1 − 1 ∕4(r − 1)

Experiments
In this section, we will perform some numerical tests to assess our findings. To this end, we consider the simplified problem where Ω ∶= (−1, 1) 2 ⧵ [0, 1] × [−1, 0] ⊂ ℝ 2 is an L-shaped domain, f ∈ L 2 (Ω) , and the coefficient either obeys the Carreau law (35) or the relaxed power-law (38). We remark that the theory derived before equally applies to this simpler case. In all our experiments below, we use a conforming P1-finite element discretisation, where the mesh consists of O(10 6 ) triangles, except where explicitly stated otherwise.

Error decay in dependence on r
First, we will examine how the convergence rate of the error depends on the exponent r−2 2 ; recall that the norm error is equivalent to the energy error, cf. Proposition 3.6. This will be done for both the Carreau and the relaxed power-law model, for smooth and irregular solutions.

Error decay for the Carreau model
Let the function obey the Carreau law (35), with ∞ = 1 , 0 = 100 , = 2 , and varying values of r ∈ (1, 2) . The source term f is chosen so that the unique solution is given by (a) The smooth function u ⋆ (x, y) = sin( x) sin( y) , where (x, y) ∈ ℝ 2 denote the Euclidean coordinates; (b) The function where R and are polar coordinates, which exhibits a singularity at the origin (0, 0). In the smooth case (a) the mesh is uniform, and in the singular case (b) the mesh is increasingly refined in the vicinity of the singularity point (0,0). In Fig. 1 we plot the error ‖ ‖ ∇u n − ∇u ⋆‖ ‖L 2 (Ω) against the number of iterative steps n. We can clearly see that the convergence rate deteriorates with decreasing r, as was predicted in Sect. 3. We further note that the irregularity of the solution in (b) does not affect the convergence rate, as was conjectured in Sect. 3.1.

Error decay for the relaxed power-law model
Now consider the relaxed power-law model, cf. (38), with − = 10 −6 and + = 10 6 . As before, the source term f is chosen so that (a) u ⋆ is smooth, and (b) u ⋆ exhibits a singularity at the origin (0,0). In Fig. 2, the error ‖ ‖ ∇u n − ∇u ⋆‖ ‖L 2 (Ω) is plotted against the number of iterative steps n. We observe that for the power-law model the dependence of the convergence rate on the exponent is even stronger than for the Carreau model.

Error decay for close to constant viscosity
In the experiments before we had that the ratio of the infinite and zero shear rates was much smaller than (r − 1) , cf. (36). Now we choose the parameters so that the 'shear stress' depends almost linearly on the 'shear rate', and we further let the source term f be such that the unique solution of (4) is given by the smooth function u ⋆ (x, y) = sin( x) sin( y) . For the Carreau law we set ∞ = 1 , 0 = 2 , = 2 , and take again varying values of r ∈ (1, 2) ; we emphasize that, in Relaxed power-law model this test, we consider even smaller values of r than in the experiments before. In view of the a posteriori computable contraction factor (36) we expect that the convergence rate will not deteriorate drastically for r close to one, which is confirmed by our numerical experiment, cf. Fig. 3 (left). In the case of the relaxed power-law model, we set − = 1 , + = 2 , and test the same values r ∈ (1, 2) as for the Carreau model. We note that these choices of the relaxation parameters ± are in practice of no interest, as one is, rather, interested in − → 0 and + → ∞ . Nonetheless, we still presume that the convergence deteriorates for r close to one by our analysis in Sect. 3.2. This is indeed the case, as illustrated in Fig. 3 (right).

Error decay in dependence on the zero and infinite shear rates
Next, we will show that, for fixed r ∈ (1, 2) , the convergence rate does not essentially deteriorate when the ratio of the infinite and zero shear rates decreases. As in the experiment before, we choose the source term f so that the unique solution is given by the smooth function u ⋆ (x, y) = sin( x) sin( y) . For the Carreau model we set = 2 , r = 1.5 , 0 = 10 a , and ∞ = 10 −a for a ∈ {1, 2, 3, 4, 5} . As we can see from Fig. 4 (left), the convergence rate is (almost) independent of a, i.e., the convergence does not deteriorate for a decreasing quotient ∞∕ 0 . For the relaxed power-law model we set r = 1.5 , − = 10 −a , and + = 10 a for a ∈ {1, 2, 3, 4, 5} . Even though the plots differ for the various values of a, the convergence rate is almost the same for all of them; indeed, no significant deterioration of the convergence rate can be observed in Fig. 4 (right) for increasing a.

Energy decay and the contraction factor
We now focus on the energy decay, and compare the exact contraction factor, cf. (40), the a posteriori computable factor (41), and the worst case factor from Remark 3.3, cf. (42). Again, this will be done for the Carreau and the relaxed powerlaw models. In our figures below, we plot the energy decay (u n ) − (u ⋆ ) , as well as the aforementioned factors against the number of iteration steps n.

Energy contraction for the Carreau model
We consider the Carreau model, cf. (35), for ∞ = 1 , 0 = 100 , = 2 , and r = 1.3 , respectively r = 1.1 . In both cases, we approximate the discrete solution for the source term f from case (a) before by performing seventy steps of the Kačanov iteration (11), and subsequently use this approximation for the determination of the reference energy (u ⋆ ) ; here, u ⋆ denotes the unique minimiser in the finite element space. We can clearly observe in Fig. 5 that, on the one hand, the a posteriori computable factor q A (n) , cf. (41), is much larger than the actual factor q E (n) , cf. (40). On the other hand, however, the factor q A (n) clearly still improves the worst case factor q W (n) from Remark 3.3, cf. (42).

Energy contraction for the relaxed power-law model
Let the coefficient obey the relaxed power-law model with − = 10 −6 , + = 10 6 , and r = 1.3 , respectively r = 1.1 . In this experiment, we approximate the discrete solution u ⋆ by performing fifty and one hundred iteration steps for r = 1.3 and  r = 1.1 , respectively. As before, the a posteriori computable contraction factor q A (n) is noticeably larger than the exact factor q E (n) , however, this is less marked than before; see Fig. 6. Moreover, it considerably improves the worst case contraction factor q W (n) ≈ 1 − 10 −12 .
We now repeat this experiment on a coarser mesh consisting of O(10 5 ) uniform triangles. In Fig. 7 we plot the factors q A (n) , q E (n) , as well as q(n) from (25) against the number of iteration steps. We observe that the (non-computable) factor q(n) from (25) has a similar trend as the exact factor q E (n) , cf. (40), and approximates the computable factor q A (n) , cf. (41), as the number of iteration steps increases.
Finally, we remark that, in the context of fixed point iterations, the contraction factor can be (heuristically) approximated by as n → ∞ , see, e.g., [18]; we emphasize that q H (n) ≥ 0 , for n ≥ 2 , thanks to (27). As can be observed in Fig. 8, the factor q H (n) , cf. (43), does indeed approximate the exact factor q E (n) from (40) well for sufficiently large n. However, in contrast with the bound q A (n) from Theorem 3.4, the computable factor q H (n) does not provide any guaranteed a priori information.

Energy decay for different mesh sizes
We conclude this section with a comparison of the energy decay for different mesh sizes. For the Carreau model, cf. (35), we set ∞ = 1 , 0 = 100 , = 2 , and r = 1.3 . In the case of the relaxed power-law model, let − = 10 −6 , + = 10 6 , and r = 1.3 . In each case we approximated the discrete solution, and, in turn, the corresponding energy by performing one hundred iteration steps. As we can see from Fig. 9, the asymptotic convergence rates (almost) coincide for the different mesh sizes.

Conclusion
In this article, we established an a posteriori computable (energy) contraction factor for the Kačanov scheme (on finite-dimensional Galerkin spaces), motivated by applications to quasi-Newtonian fluid flow problems. For the relaxed power-law model, this factor is independent of the relaxation parameters ± ; we also demonstrated that it is, instead, the power-law exponent that affects the convergence rate of the iteration. In contrast, existing bounds on the contraction factor of the relaxed Kačanov iteration depend on the relaxation parameters ± in an unfavourable manner, in the sense that they tend to 1 as − → 0 and/or + → ∞ . A series of numerical tests have confirmed that our a posteriori computable contraction factor improves, on finite-dimensional Galerkin spaces, existing bounds, and that, as predicted by our analysis, for the power-law model it is in fact the closeness of the power-law exponent r ∈ (1, 2) to 1 that influences the convergence rate of the iteration. However, our experiments revealed that the theoretically derived bound on the contraction factor of the Kačanov scheme is still too pessimistic.

Page 26 of 27
The idea is to identify quadratic functions g ± , on [ 2 ± − , 2 ± + ] , respectively, which smoothly connect the constant parts of with the map t ↦ t r−2 2 = (t) on [ 2 − + , 2 + − ] , i.e., A straightforward calculation reveals that the properties (a)-(d) are satisfied for and Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.