{'fraction_non_alphanumeric': 0.05136418866674031, 'fraction_numerical': 0.014525571633712582, 'mean_word_length': 3.925734494015234, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Convective transport in nanofluids: regularity of solutions and error estimates for finite element approximations

We study the stationary version of a thermodynamically consistent variant of the Buongiorno model describing convective transport in nanofluids. Under some smallness assumptions it is proved that there exist regular solutions. Based on this regularity result, error estimates, both in the natural norm as well as in weaker norms for finite element approximations can be shown. The proofs are based on the theory developed by Caloz and Rappaz for general nonlinear, smooth problems. Computational results confirm the theoretical findings.


Introduction
Nanofluids, i.e. a dilute mixture of a conventional base fluid and particles of submicron size, have received much attention for instance as a cooling liquid. This is due to their superior heat transfer properties. The enhanced heat transfer cannot solely be understood by the altered heat conducting coefficients of the mixture, but rather by effects of a heterogeneous distribution of the particles. Among the mathematical models to explain such behavior, the Buongiorno model [8] has become rather popular. By now, many simulations are based on this model, see [1,4,11,14,16,17] for a by far not complete list of applications. In [5] the mechanism of the enhanced heat transfer for laminar flow conditions was revealed: strong temperature gradients at a hot wall lead to reduction of concentration of particles by thermophoresis there and this in turn reduces the concentration dependent viscosity of the dispersion. This then alters the flow profile leading to a stronger convective heat transfer.
To the best of our knowledge, despite its relevance in applications, there is hardly any rigorous mathematical analysis of the Buongiorno model. In [5] existence of weak solutions was shown using energy techniques. It was shown that solutions of a decoupled semi-implicit time-discretization converge to a solution of the continuous system, thereby also suggesting an effective numerical method.
In [6] existence of solutions to the stationary system was shown. Interestingly, the proof is somewhat technically more demanding than for the time-dependent problem.
The objective of the present work is to first show (under some smallness assumptions) regularity for the stationary problem and then use these regularity results to prove quasi-optimal error estimates for finite element approximations of the system. To this end it is shown that the system can be cast into the general framework of Caloz and Rappaz [9] for nonlinear (smooth) problems.
It turns out that the right space for the scalar quantities concentration and temperature is W 1 p with p > d, d the space dimension, whereas for the fluid part we can stay in a Hilbert space setting.
The rest of the paper is organized as follows. In Section 2 we present the mathematical model and set some notation. In Section 3 we present the regularity results for the solutions of the system of PDEs. In Section 4 we present a linearization of the problem which will allow us to prove the error estimates for a finite element discretization in Section 5. We close this article with some numerical experiments in Section 6 where we illustrate the orders of convergence, as well as the interesting effect of thermophoresis as a means to enhance the heat transfer properties.

The mathematical model
We consider the stationary system of a variant of the four equations, two-phase Buongiorno model [8] describing the motion of a nanofluid including concentration transport by thermophoresis. The model has been slightly modified to make it thermodynamically consistent, see [5,6]. In non-dimensional form it reads as follows: Let Ω ⊆ R d , d ∈ {2, 3} be an open, bounded domain with C 2 boundary. We look for a concentration field φ, a temperature T as well as a velocity u and a pressure p fulfilling the following system of equations in Ω (in the distributional sense) ∇· u = 0, with D(u) = ∇u + ∇u T and the particles' flux given by with N BT the ratio of Brownian diffusivity/thermophoretic diffusivity, T 0 a non-dimensional ambient temperature, Re the Reynolds number, Pr the Prandtl number, Sc, Sc f the Schmidt and fluid Schmidt number, respectively and Le, the Lewis number. Buoyancy effects through a Boussinesq approximation are considered with β > 0 and e g denoting a unit vector in the direction of gravity. The above system must of course be supplemented by appropriate boundary conditions.
The flux j is the non convective slip flux consisting on a Brownian part −∇φ and the so called ther- T0 that drives particles from hot to cold. Phenomenologically this can be explained by the fact that collisions of the particles with molecules from the base fluid are stronger on the hot side of the particle than on the cold part, resulting in a net flux from hot to cold. For more on thermophoresis we refer for instance to [13].
The mathematical challenge with the above system lies in the rather strong nonlinearity. The right space for φ, T is therefore W 1 p (Ω) with p > d. However, an energy estimate is only available in H 1 (Ω), see [5,6]. To overcome this problem, in Section 3 we prove regularity estimates in W 2 p (Ω) based on some bootstrap arguments and a smallness assumption.
This paves the way to cast the problem in the general framework developed in [9] for nonlinear problems.
We define a vector-valued cut-off function, for R > 0 as follows Note that σ R is Lipschitz and For R > 0 we now consider the regularized problem in Ω, can be proved as in [6] obtaining also the a priori estimate with C independent of R and 0 ≤ φ R ≤ 1, a.e.
The following lemma will be instrumental in proving our regularity results.
Proof: The regularity for u follows from the proof of [2, Lemma 4]. There, regularity in H 2 (Ω) for j = 0, 1 was shown. However, a closer inspection of the proof shows that the assertion is also valid for M as above and arbitrary j ∈ N 0 . The regularity for T can be shown in the same way with even some simplifications.
The first regularity result for the solution of (3.1) holds under a smallness assumption on the data of the problem. The precise result is the following.
is also a solution to (3.1).
Proof: The proof is based on a couple of bootstrap arguments.
We fix p > d and denote D := φ D 2,p + f 0,p + g 0,p < ∞. For each R > 0, let (φ R , T R , u R ) be a solution to (3.2) satisfying (3.3), and emphasize that the constants involved in the symbols below do not depend on R, D, or the problem data φ D , f , g, but only on Sobolev embeddings and the constant C from (3.3).
By embedding, u R ∈ L 6 (Ω). The right-hand side for the φ R -equation (3.2a) is the sum of a function in L 3/2 (Ω) and the divergence of a function in L ∞ (Ω) (bounded by R). Then by regularity [3,Theorem 3.29], ∇φ R ∈ L 3 (Ω). Now, the right-hand side of (3.2a) is the sum of a function in L 2 (Ω) and the divergence of a function in L ∞ (Ω). This sum is now the divergence of a function in L 6 (Ω), so that by the same regularity result [3, Theorem 3.29], ∇φ R ∈ L 6 (Ω) and repeating this argument, ∇φ R ∈ L M (Ω) for all 1 ≤ M < ∞ with the estimate which thereby implies Let us now turn to the equations for T R and u R . Thanks to (3.2a) and the definition of η, ρ the right-hand sides can be written asf First, we need an intermediate regularity result for u R , namely u R ∈ L M (Ω) for all 1 ≤ M < ∞. To this end, we observe Therefore, by Lemma 3.1 applied to where we have used (3.8).
Finally, using again Lemma 3.1 and the embedding so that under this assumption, ∇T R 0,∞ ≤ R and the first assertion follows with The second assertion is an immediate consequence of (3.9).
Combining the previous theorem with regularity results we obtain the following corollary.

Proof:
We let F 0 be as in Theorem 3.2 and the problem data satisfy D : Also, T is a solution of (3 It remains to show H 2 -regularity for u, which is an immediate consequence of the fact that u is a solution to (3.5) withg = g − (j + (1 + φ)u) · ∇u − T e g , which satisfies Note that in the proof of Theorem 3.2 above we have already shown that u R ∈ W 2 3/2 (Ω) → W 1 3 (Ω).
Once regularity in W 2 p (Ω), H 2 (Ω) is established, it is not difficult to get higher regularity, provided data is more regular. This is stated in the next corollary.
Proof: The proof is based on induction over j. The case j = 0 is the assumption of this corollary, which holds under the smallness assumption of Corollary 3.4. So let us assume that the statement is correct for j. We shall then show that it also holds for j + 1. First, note that by the induction assumption and Let us now check the regularity of ∂ αf for |α| = 1 + j: Clearly, the first term on the right-hand side is in L p (Ω). By Leibniz' formula, the second term can be written as Our first observation yields ∂ α (∇φ · ∇T ) ∈ L p (Ω). The last term can be written by Leibniz' formula as In view of our first observation, the worst summand in the sum above (if 6 < p) is attained for |β 2 | = 1+j with ∂ β2 u ∈ H 1 (Ω) → L 6 (Ω), showing that altogether ∂ α (1+φ)u·∇T ) ∈ L min{6,p} (Ω). Since |α| = 1+j was arbitrary, this shows ∂ αf ∈ L min{6,p} (Ω). As an intermediate result one gets T ∈ W 3+j min{6,p} (Ω).
The velocity u is a solution to (3.5) withg = g + ∇φ + h(φ)∇T − (1 + φ)u · ∇u − T e g . The derivative of the right hand sideg for this momentum equation reads Expanding again the derivative by Leibniz' formula and observing W 1 min{6,p} (Ω) → L ∞ (Ω) the worst term is identified to be for all |α| ≤ 1 + j and then u ∈ H 3+j (Ω), which is already the desired regularity for u.
One more sweep of the above arguments, but now using the intermediate regularity results, concludes the proof. Proof: From standard theory it is clear that there exists a unique pressure p ∈ L 2,0 (Ω) such that the above equation is fulfilled in the distributional sense. Now, shifting the term ∇· (µ(φ)D(u)) to the right hand side, differentiating the right hand side successively up to the desired order and using the regularity of φ and u it follows that ∇p ∈ W j M (Ω) and the estimate is then also immediate.
Due to the properties of the coefficients and since p > 3 ≥ d (which in particular implies W 1 p (Ω) → L ∞ (Ω)) F is Frechet differentiable with derivative given by We state the differentiability and also the Lipschitz continuity of DF in the next lemma.
Proof: The differentiability was already discussed above. To show the Lipschitz continuity we have to estimate which again follows by inspecting the individual integrals and noting that p > d.
The goal of these linearization results is to obtain estimates for a finite element discretization using the framework from [9]. The crucial step now is to show that under some smallness assumptions DF(U ) is an isomorphism from X to Y . with R :W 1 p (Ω) →W 1 p (Ω), RΘ := k(φ)Θ, see [9]. Due to the properties of k(·) and since φ ∈ W 1 p (Ω), R is well defined and an isomorphism.

Finite element discretization and error estimates
Let {T h } h>0 be a quasiuniform, shape regular family of conforming triangulations of Ω with To avoid technical details estimating the mismatch of the triangulation with the exact geometry we (unrealistically) assume that elements on the boundary are curved and match the boundary exactly. Hence Remark 5.1. The quasiuniformity of (T h ) h>0 is required to guarantee the W 1 p stability of the Ritz operator (see below).
To discretizeW 1 p (Ω),W 1 p (Ω) we choose Lagrange elements of polynomial order k ≥ 1. Denote this space by S h = S h (T h ). Other choices, however, are possible and are restricted only by the assumptions in [7,Chapter 8]. Furthermore, for the Navier-Stokes part of the system we choose an inf-sup stable pair of The discrete problem reads:

Error estimates in the norm of X.
Let U = (φ, T, u, p) ∈ X be a solution of F(U ) = 0. To be able to apply the general results from [9] regarding error estimates, we have to show the following.
(1) F : X → Y is differentiable; (2) DF is locally Lipschitz continuous at U ; (5) the following discrete inf-sup condition holds: inf Properties (1)- (3) have been shown in the previous section under some smallness assumption and (4) holds by construction. The remaining point thus is the inf-sup condition (5).
As in the proof of Proposition 4.2 we first consider the reduced operator T . Since we chose a pair of finite element spaces V h × Q h which is inf-sup stable for Navier-Stokes and due to Korn's inequality, for the (u, p) part of T one has: for all w h ∈ V h , r h ∈ Q h and some β 1 > 0.
Next we consider the Θ-part of T . From [9, Thm. 10.1] we infer that for all Θ h ∈ S h and some β 2 > 0. The crucial point in the proof in [9] was the stability of the Ritz-operator In [9] this result was cited from [15], where it was proved for dimension d = 2. A much more general result valid for d = 2 as well as for d = 3 and a variety of finite element spaces can be found in [7].
In order to get an inf-sup estimate for the χ-equation, we rescale the T -equation for F: For λ > 0 define All what have been shown for F and DF remains valid also for F λ and DF λ except that β 2 becomes λβ 2 . Now for the χ-equation we use the same result from [9] as for the Θ-equation and infer for where the last step follows from the boundedness of h(·) and Hölder's inequality.
Putting everything together we arrive at The rest of the inf-sup estimate follows exactly as in the proof of Proposition 4.2: define N λ := DF λ (φ, T, u, p) − T λ . Then we readily have We are now in a state to apply [9, Thm. 7.1]. (Ω), f, g ∈ W j p (Ω).
Moreover, for 0 < h ≤ h 0 the discrete problem (5.2) has a locally unique solution X h U h ∈ B r0 (U ) ⊆ X and the following error estimate holds: Proof: By the above theorem we know that for sufficiently small δ > 0 there is a solution (φ, T, u) ∈ X ∩ W 2 p (Ω) × W 2 p (Ω) × H 2 (Ω). Possibly reducing δ further, the above theorem guarantees the existence of a locally unique discrete solution U h and its quasi- The rest follows by the approximation property of the finite element spaces and the (possibly) higher regularity shown in Corollary 3.5 and Lemma 3.6.

Error estimates in weaker norms.
Error estimates in L p , L 2 can be proved by using rather standard duality techniques.
Recall that B is an isomorphism, iff B * is an isomorphism. As shown in Proposition 4.2 this is for instance the case, if (φ, T, u, p) X is small enough, which we assume hereafter. Now, choose G = (g φ , g T , g u , 0) ∈ L p (Ω)×L p (Ω)×L 2 (Ω)×L 2,0 (Ω) ⊆ X such that g φ 0,p = g T 0,p = g u 0,2 = 1 and Let W ∈ Y be the unique solution of the dual problem (which exists, since B * is an isomorphism) for allŨ ∈ X.
Let W h ∈ Y h . As in [9] we calculate for U, U h the continuous and discrete solution, respectively. The last term on the right side can be treated as follows if U − U h X is sufficiently small and with L the Lipschitz constant of DF. The above estimate may be viewed as a substitute for the orthogonality of the error in the linear case. Using this estimate in Eq. (5.4) and the boundedness of b(·, ·) (which follows from Lemma 4.1) we arrive at It remains to estimate W − W h Y , which will be accomplished by a regularity result for W .
In the above arguments it is of course understood that the corresponding norms are bounded by the right-hand side.
Putting everything together, we arrive at the following error estimate.
Theorem 5.5. Let p > 3 and U = (φ, T, u, p) ∈ X be a solution of the continuous system Eq.
Proof: In view of the definition of the dual solution W from Eq. (5.3) and the estimate Eq. (5.5), take the infimum over all W h ∈ Y h in the latter estimate. The assertion then follows by the regularity of W from Lemma 5.4 and the approximation properties of the finite element spaces.

Computational results
For the computations presented in this section piecewise quadratic finite elements are used for φ and T as well as the P 2 × P 1 -Taylor-Hood element for u and p. In order to get an interesting computational example, where one can see the effect of thermophoresis, we slightly deviate from the set of boundary conditions imposed in the theoretical part. Set Ω =]0, 2[×]0, 1[. For the concentration φ we impose homogeneous Neumann boundary conditions on the whole of ∂Ω and a mean concentration φ m = 0.1. The temperature T is set to T = 1 at the left side wall and T = 0 at the right side wall. On the remaining parts of ∂Ω homogeneous Neumann boundary conditions are imposed. For the velocity u homogeneous no-slip conditions are chosen except for the upper boundary, where a slip condition is enforced, i.e. u 2 = 0 together with vanishing tangential stress µ(φ)(∂ x2 u 1 + ∂ x1 u 2 ) = 0. Note that this condition can be realized as a natural boundary condition for the space {v ∈ H 1 (Ω) | v 2 = 0 on the upper part of ∂Ω, v = 0 else on ∂Ω}.
Let us first consider the case without thermophoretic effects (i.e. φ ≡ φ m ). The flow is driven by buoyancy forces: the liquid heats up at the left lateral wall inducing an upward flow field that turns to the right at the top of the container, transporting warm liquid to the right, cold wall, where it cools down and flows downwards. The cold liquid is flowing back at the bottom of the container to the left hot wall, see Fig. 1, left picture.
Switching on thermophoretic effects, the flow field is strongly enhanced on the upper boundary. This can be understood by inspecting the concentration field, see  In order to quantitatively assess the convergence, the same setting as above, however with different parameters, is used. Since the exact solution is unknown, the computational solution on a very fine grid with nt = 262, 144 triangles is used as reference instead. Starting from a coarse triangulation, the grid is successively refined by two bisection steps each. The corresponding errors and the experimental order of convergence (EOC) are listed in Tab. 1. As expected one gets a convergence order of 3 (although the boundary of the domain is not of class C 3 ).