AN ALGORITHM FOR THE GRADE-TWO RHEOLOGICAL MODEL

. We develop an algorithm for solving the general grade-two model of non-Newtonian fluids which for the first time includes inflow boundary conditions. The algorithm also allows for both of the rheological parameters to be chosen independently. The proposed algorithm couples a Stokes equation for the velocity with a transport equation for an auxiliary vector-valued function. We prove that this model is well posed using the algorithm that we show converges geometrically in suitable Sobolev spaces for sufficiently small data. We demonstrate computationally that this algorithm can be successfully discretized and that it can converge to solutions for the model parameters of order one. We include in the appendix a description of appropriate boundary conditions for the auxiliary variable in standard geometries.


Introduction
Non-Newtonian fluids are found in all aspects of life, from bodily fluids [28] to engine oil [12,27].Rheology (non-Newtonian behavior) plays a significant role in manufacturing, including food [17,19].Thus advances in modeling and simulation of non-Newtonian fluids can have broad impact.
Models of non-Newtonian fluids have been studied extensively for many years, but only recently have there been mathematical advances [10] that allow models for them to be understood more completely.This understanding now allows development of numerical solution methods with a new level of reliability.The grade-two model is the simplest of a family of models proposed by Rivlin and Ericksen [11,15] in which the stress-strain relationship involves derivatives of the fluid velocity.It has been widely studied, but to date no general numerical method has been proposed and analyzed for solving it.
There have been many different approaches to the grade-two model.In two dimensions, certain simplifications can be made if one of the parameters is eliminated, and this allows both rigorous analysis of the system in Lipschitz domains [15] and also extensive numerical analysis of effective discretization schemes [14].However, in [15], it was assumed that the flow velocity was tangential to the boundary.Still in two dimensions, the paper [8] removed that restriction by imposing third-order boundary conditions on the inflow velocity.
However, different approaches were required for general parameters and in three dimensions [1,[3][4][5].Although the method proposed in [1] is quite general, it was developed and analyzed only in the case of tangential flow fields.We propose here a slight variation of the approach [1] to the grade-two model that provides some simplifications both theoretically and computationally, and it applies in full generality.
For certain flow problems, it is critical to allow nontrivial inflow and outflow boundary conditions.For example, simulating many rheometers [21], the basic instruments for measuring fluid properties, requires this.In Figure 1 we depict two different rheometers which involve nontrivial inflow and outflow.The Lodge Stressmeter, depicted in Figure 4.13 of [2], was developed by Arthur Lodge [20].The rheometer measurements are based on the pressure difference between the top of the middle section of the channel and the bottom of "hole" along the bottom of the channel.Contraction rheometers [23] measure the force on the contraction section in the middle between the large inflow channel and the smaller outflow channel.The shape of the contraction section can be chosen differently.Thus a major contribution of this paper is the development of analytical techniques to cover this type of boundary condition.
One issue with the different methods is the requirement for boundary conditions on the inflow boundary.Since the grade-two model is a third-order partial differential equation, we expect there to be another boundary condition in addition to the standard ones for flow problems, such as the Stokes no-slip condition.In [1], this issue was avoided by assuming tangential flow on the boundary.Generalizing [1] to allow an inflow boundary requires a boundary condition on a certain tensor Σ.In the approach proposed here, a condition is posed instead on the vector −∆u + ∇, which is directly related to the divergence of the stress.Thus the additional boundary condition can be viewed as a stress boundary condition.We give examples of what this boundary condition should be for certain geometries.

Rheology models
In all (time-independent) models of fluids, the basic equation can be written as where T is called the extra (also called deviatoric) stress and f represents externally given data.The models only differ according to the dependence of the stress on the velocity u.In the case of a Newtonian fluid where A = ∇u + (∇u)  .Thus, when ∇• u = 0, it follows that ∇• T = ∆u, and we obtain the well known Navier-Stokes equations for Newtonian flow, where  is the kinematic viscosity [18].Typically, the data f is zero, but instead nonhomogeneous boundary conditions are physically relevant.Thus we will assume that (2.1) holds in some domain Ω and that u = g on Ω, where we assume ∫︁ Ω g • n d = 0 to allow divergence-free solutions.Depending on the details of the model, there will also be a need for appropriate boundary conditions for T and other ingredients.

Grade-two fluid model
The grade-two model of Rivlin and Ericksen [11,15] can be expressed as a single equation.The stress tensor for the grade-two fluid model satisfies where A = (∇u) + (∇u)  = 2E and the material derivative and the lower-convected Oldroydian derivative are given by D D f := for any tensor-valued function f .For the steady-state, grade-two fluid model, the stress tensor simplifies to We have used the notation ∘ for tensor multiplication, which here will be just matrix multiplication.Thus the equations of motion (2.1) can be written Here We assume that the boundary data g is defined on all Ω, is divergence free, and sufficiently smooth, to be specified subsequently.

Solving the grade-two model equations
It is helpful to expand the divergence of ̂︀  , defined in (2.4), to get a better sense of what the various terms are in (2.3).Recall (3.2) of [16] that Recall that A is shorthand for A = ∇u + ∇u  .Thus 3) can thus be transformed [1] using (2.5): where and so (2.6) can be further transformed to where Thus Note that  appears at first not to be a symmetric tensor due to the term ∇u  ∘A.However, in two dimensions this is a symmetric matrix if ∇• u = 0.
Proof.Write  as using the fact that the trace of  is zero.Then )︂ .
The latter matrix is evidently symmetric.
The lemma is not true for 3 × 3 matrices as simple examples show.

Transformed grade-two model equations
Let  be defined by solving with suitable inflow boundary conditions [5,8].Then Thus (2.9) transforms to Note that  is not a symmetric tensor due to the term ∇u  .The incompressibility condition ∇ Now consider a coupled system that looks initially like a problem potentially different from (2.3), which is a slight variant of the one proposed in [1]: where Much of the paper will be devoted to proving this system is well posed and provides an equivalent formulation for solution of (2.3). Theorem which is the same as (2.9).Reversing the derivation of (2.9) proves (u, ) satisfies (2.3) with ̂︀  defined by (2.4).The statement about w just involves replacing −∆u in (2.3) by the indicated expressions.
The difference between (3.8) and equation (2.6) of [1] is that w replaces ∇•  for a certain tensor , and a transport equation is posed for the full tensor  as opposed to the vector w.Using (3.8) gives a smaller system to solve.The issue of inflow boundary conditions [6,7] did not arise in [1] which was restricted to tangential flows.Thus in the general case, some suitable expression for  on the inflow boundary would be required.

An algorithm for the transformed equations
The system (3.8) is analogous to the reduced system in [15], and the algorithm in that paper suggests an algorithm for solving (3.8): start with some w 0 , then solve for  ≥ 1 For definiteness, we will take w 0 = w  .We prove convergence of this iteration for small data (g and w  ) in Section 4.3.To begin with, let us establish a basic bound.We collect details on the Lebesgue and Sobolev spaces and norms used in Appendix A. Consider the Sobolev inequalities where ,  is a constant that depends only on the dimension , and   is the Sobolev constant in (3.10).
Proof.Applying (3.10) to (3.6), we get for a constant  that depends only on the dimension .Note that Combining these estimates yields (3.11). where ,  is a constant that depends only on the dimension , and   is the Sobolev constant in (3.10).
Proof.One consequence of the assumption  > /2 is the Sobolev inequality for any tensors T, U.
In view of (3.6), we get We have from the Sobolev inequalities (3.13) and (3.10) that Similarly, for a constant  that depends only on the dimension , we have The remaining terms are similar.
We can recover the physical pressure  from (3.1), that is One computational challenge is that (2.3) is a third-order PDE due to the presence of the term u • ∇(∆u).Thus we need to be careful about the number of boundary conditions required to get a unique solution.

Variational formulation
A variational formulation of (3.9) is as follows.The first two equations can be approximated by the iterated penalty method: find u ,ℓ ∈  ℎ + g such that Once this is converged, we set u  = u ,ℓ and define the pressure via [22] ∫︁ Note that   has mean zero if constant functions are in Π ℎ , in view of the divergence theorem: We can pose the transport equation (3.9) via: where w  is posed only on the inflow boundary, that is, Note that it is tempting to integrate by parts to get but there would be boundary terms that would need to be added to the formulation.We can take  ℎ to be continuous, vector-valued, piecewise polynomials of degree  and Π ℎ to be continuous, scalar-valued, piecewise polynomials of degree  − 1.The use of continuous elements in (3.16) is called the unified Stokes algorithm (USA) [22].
Recall from (3.6) that Recalling (2.8), we have We can compute the physical pressure   from (3.14) via but this does not need to be done at each iteration.

Required inflow boundary conditions
One drawback to the proposed method (3.9) is that it requires specification of boundary conditions for w = −∆u + ∇.Although we cannot provide general guidance for this, we can compute boundary conditions for w for typical flow geometries.We present this in Appendix B.

Theoretical details
Here we collect the theoretical details required to prove the validity of our algorithm.We begin with an assumption about the smoothness of the data and domain.First we assume that for some  > , g ∈  2  (Ω), with Further, we assume that there is a constant   such that for any g as above and any w ∈   (Ω) the solution (u, ) of satisfies, for  = 0, 1, where   depends on Ω as well as .For  = 0, we require   > , but for  = 1 we only require   > /2.This follows from Theorem 5.4 in [13], page 88 when Ω is sufficiently smooth.In our computational tests, we will have less smoothness with the polyhedral domains used, but these could be approximated by smooth domains.We will prove the following theorem which establishes the existence of solutions for the grade-two model (2.3) as well as for the equivalent model (3.8).
Theorem 4.1.Assume that (4.2) holds for the Stokes problem (4.1).Suppose that /( + 1) <  ≤   , for  = 0, 1, and that  satisfies Then there exist positive, finite constants  and   such that if the boundary data satisfy and and the initial iterates are sufficiently small, then the iterates (3.9) are bounded for all  > 0: where  is a finite positive constant and  = 0, 1.Moreover, (u  ,   , w  ) converge geometrically in The constraint (4.3) implies  > 2 for  = 2 and  > 6 for  = 3, and thus the constraint  >  is satisfied implicitly.In our computational experiments, we will see that the assumptions on the data size may not be very restrictive in practice.
Let   =     /.Then we have proved that Define  to be any constant such that Then (4.9) implies Then we conclude that   ≤ 1 4 as well.Note that by definition,  ≥ 1, so  ≤ 1/4.Recall that we have taken w 0 = w  .
Therefore, if the boundary data is sufficiently small, e.g., we conclude that in particular that ⃦ ⃦ w 0 ⃦ ⃦   (Ω) ≤ 1 4 , and thus for all  > 0. Thus also for all  > 0. Note that we can take the constant  as large as we like.

𝑊 1 𝑞 bounds on the iterates
Applying (4.2) with  = 1 to the algorithm (3.9), we have In [24], it is proved that the unique solution of (4.7) satisfies Note that Hölder's inequality and (4.13) imply Combining this with (4.16) yields where Then (4.17) implies . Then we have Then By taking  sufficiently large, we have More precisely, this holds when where  0 is defined in (4.10).Assume further that Then (4.18) and (4.20) imply that Note that (4.19) implies that Under the inductive hypothesis that then (4.23) implies that provided that  ≥ 16.Therefore (4.24) implies that for all , provided that where  1 is defined in (4.21).Using (4.14), we find for all  > 0, under the assumptions (4.19) and (4.22),where

Convergence estimates
Recall the tensor  introduced in (2.10): Thus (3.4) and (3.6) imply To estimate terms involving  , note that for any two sequences   and   , Thus (4.29) implies We examine these four terms separately.First, ( For the second term, there is a constant  depending only on the dimension such that and a similar estimate holds for the third and fourth terms.Thus (4.30) shows that At the last step, we utilized the fact that u  − u −1 is zero on the boundary, so we could apply (4.2) directly.Define e = w  − w −1 .Then from (3.9) (4.34) Applying (4.2) to (4.34), we get Taking  sufficiently large, that is, Therefore the sequence w  converges geometrically in   (Ω)  .Subtracting iterates in (3.9), we find Thus (4.2) implies that and thus the sequence u  converges geometrically in  2  (Ω)  and the sequence   converges geometrically in  1  (Ω).This completes the proof of Theorem 4.1.

Computational experiments
Here we explore computational techniques for implementing the grade-two algorithm.Table 1 presents results for the algorithm (3.9) for  =  1 =  2 = 1, implemented using quadratic and cubic elements using the iteratedpenalty method (IPM) [9], on the unit square domain (B.1) with  = 1 and boundary data (B.11) with  = 1.The exact  is quadratic as indicated in (B.8).Using piecewise degree  elements for  ℎ results in piecewise degree  − 1 elements for the pressure approximation.Thus for  = 2, the pressure approximation is only piecewise linear, and the approximation of  dominates the overall errors.Table 1a indicates that the error   is close to second order.But for  = 3, the exact  is in the pressure space, and we get essentially round-off error.Due to some sort of instability, the errors grow as the mesh size is reduced, but they are significantly smaller than for the case  = 2. Figure 2 shows that there is a localized error that occurs in  1 (which should be identically zero) near the corners of the inflow boundary.This pollutes the component  2 (which is isolated in Tab. 1) and causes errors in u and .Computations for  = 4 yielded similar results as for the  = 3 case.
The error for approximating w in  1 are much worse than for other errors.But we know from Section 4.2 that the transport problem does not have uniform bounds in  1 , so the larger errors are not surprising.Table 1.Grade-two simulations of Poiseuille flow in the domain (B.1) with  = 1 and  = 1 and boundary data (B.11) with  = 1 using (a) piecewise quadratics for  ℎ and (b) piecewise cubics for  ℎ .The mesh consisted of an  × array of Malkus splits (squares divided into four triangles by the bisectors) [26].The algorithm (3.9) was implemented using the iterated penalty method (3.15).The column "iters" indicates the number of iterations of (3.9).The pressure was computed via USA [22] as described in (3.16) with Π ℎ begin continuous piecewise polynomials of degree  − 1. Errors:   1b with  = 16.

Conclusions
We developed an algorithm for solving the general grade-two model of non-Newtonian fluids which for the first time allows nontrivial inflow boundary conditions.The new algorithm couples a Stokes equation for the fluid velocity with a transport equation for an auxiliary vector-valued function.As a third-order partial differential equation, the grade-two model requires an additional boundary condition, and our new formulation leads to a

Thus
We can use the formula (2.2) to compute the stress: The tensor  is given by (2.10):For shear flow,  is linear, so  ′ is constant, and thus T G is constant.Therefore ∇• T G ≡ 0. Similarly, ∆u ≡ 0, and  is constant.Suppose that  0 is this constant.If we specify that | Γ− =  0 , then we conclude that  is also constant ( =  0 ).Thus w = 0 as well.But there could be other solutions for other choices of | Γ− , leading to nonconstant .In that case, w = ∇ ̸ = 0.

B.2. Poiseuille flow
For Poiseuille flow,  is quadratic, and ∇• T G is not even constant.Since u • ∇u = 0, the top equation in (2.3) takes the form Then we get two equations for the pressure: )︂

Figure 1 .
Figure 1.Two different rheometers: (a) the Lodge Stressmeter and (b) a contraction rheometer.The arrows on the left of each rheometer indicate the fluid inflow, and the arrows on the right of each rheometer indicate the fluid outflow.
the solution (u, , w) of (3.8),In view of Theorem 3.1, (u, ) is the solution of the grade-two model (2.3), where  is related to  by (3.1).