Numerical Analysis of Unsteady Implicitly Constituted Incompressible Fluids: Three-Field Formulation
NNUMERICAL ANALYSIS OF UNSTEADY IMPLICITLYCONSTITUTED INCOMPRESSIBLE FLUIDS: THREE-FIELDFORMULATION ∗ P. E. FARRELL † , P. A. GAZCA-OROZCO ‡ , AND
E. SÜLI § Abstract.
In the classical theory of fluid mechanics a linear relationship between the shear stressand the symmetric velocity gradient tensor is often assumed. Even when a nonlinear relationship isassumed, it is typically formulated in terms of an explicit relation. Implicit constitutive models pro-vide a theoretical framework that generalises this, allowing for general implicit constitutive relations.Since it is generally not possible to solve explicitly for the shear stress in the constitutive relation, anatural approach is to include the shear stress as a fundamental unknown in the formulation of theproblem. In this work we present a mixed formulation with this feature, discuss its solvability andapproximation using mixed finite element methods, and explore the convergence of the numericalapproximations to a weak solution of the model.
Key words.
Implicitly constituted models, non–Newtonian fluids, finite element method
AMS subject classifications.
1. Implicitly constituted models.
In the classical theory of continuum me-chanics the balance laws of momentum, mass, and energy do not determine completelythe behaviour of a system. Additional information that captures the specific prop-erties of the material to be studied is needed; this is what is commonly known as a constitutive relation . The constitutive law usually expresses the stress tensor in termsof other kinematical quantities (e.g. the symmetric velocity gradient) and, even if itis nonlinear, it is typically formulated by means of an explicit relationship. It hasbeen known for some time that in many cases explicit constitutive relations are notadequate when modeling materials with viscoelastic or inelastic responses (see e.g.[51, 52]), which has led to the introduction of many ad-hoc models that try to fitthe experimental data. Implicitly constituted models, introduced in [51], provide atheoretical framework that not only serves to justify these ad-hoc models, but alsogeneralises them. The physical justification of these types of models, including a studyof their thermodynamical consistency, is available and will not be discussed here; theinterested reader is referred to [53, 52, 54].If a fluid occupies part of a space represented by a simply-connected open set Ω ⊂ R d , where d ∈ { , } , then the evolution of the system during a given timeinterval [0 , T ) , for T > , is determined by the usual equations of balance of mass,momentum, angular momentum and energy, which in Eulerian coordinates take the ∗ Submitted to the editors 2019.
Funding:
This research is supported by the Engineering and Physical Sciences Research Council[grant numbers EP/K030930/1, EP/R029423/1], and by the EPSRC Centre For Doctoral Trainingin Partial Differential Equations: Analysis and Applications [grant number EP/L015811/1]. Thesecond author was partially supported by CONACyT (Scholarship 438269). † Mathematical Institute, University of Oxford, UK ([email protected]). ‡ Mathematical Institute, University of Oxford, UK ([email protected]). § Mathematical Institute, University of Oxford, UK ([email protected]).1 a r X i v : . [ m a t h . NA ] D ec P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI form: ∂ρ∂t + div ( ρ u ) = 0 ,∂ ( ρ u ) ∂t + div ( ρ u ⊗ u ) = div T + ρ f , (1.1) T = T T ,∂ ( ρe ) ∂t + div ( ρe u ) = div ( T u − q ) . Here: • u : [0 , T ) × Ω → R d is the velocity field; • ρ : [0 , T ) × Ω → R is the density; • T : (0 , T ) × Ω → R d × d is the Cauchy stress; • e : [0 , T ) × Ω → R is the internal energy; • q : (0 , T ) × Ω → R d is the heat flux.The constitutive law relates the Cauchy stress (or some other appropriate measureof the stress) and the heat flux to other kinematical variables such as the shear strain,temperature, etc. In the following we will assume that the material is incompressible,homogeneous and undergoes an isothermal process. This implies that the energyequation decouples from the system and that the Cauchy stress can be split in twocomponents:(1.2) T = − p I + S , where I is the identity matrix, p : (0 , T ) × Ω → R is the pressure (mean normal stress),and S : (0 , T ) × Ω → R d × d sym is the shear stress (hereafter referred only as “stress”). Inthis work we will consider constitutive relations of the form(1.3) G ( · , S , D ( u )) = , where G : Q × R d × d sym × R d × d sym → R d × d sym and D ( u ) := ( ∇ u + ( ∇ u ) T ) is the symmetricvelocity gradient; here Q is used to denote the parabolic cylinder (0 , T ) × Ω . Theprecise assumptions on this implicit function will be stated in the next section.For a rigorous mathematical analysis of models of implicitly constituted fluids thereader is referred to [13, 14]. Existence of weak solutions for problems of this typewas obtained in [13] and [14] for the steady and unsteady cases, respectively. Someextensions include [15, 46, 50], where additional physical responses are incorporatedinto the system.As for the numerical analysis of these systems, very few results have been pub-lished so far. In [21] the convergence of a finite element discretisation to a weaksolution of the problem was proved for the steady case, and the corresponding a-posteriori analysis was carried out in [43]. More recently, this approach was extendedto the time-dependent case in [61]. Also, several finite element discretisations werecompared computationally in [41] for problems with Bingham and stress-power-law-like rheology.Numerical methods for the incompressible Navier–Stokes equations are usuallybased on a velocity-pressure formulation, and extensive studies have been carried outover the years in relation to this (see e.g. [33, 10]). Such a formulation is possible,because in the case of a Newtonian fluid the explicit constitutive relation S = 2 µ D ( u ) allows one to eliminate the deviatoric stress S from the momentum equation. In -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS
2. Preliminaries.2.1. Function spaces.
Throughout this work we will assume that Ω ⊂ R d ,with d ∈ { , } , is a bounded Lipschitz polygonal domain (unless otherwise stated),and use standard notation for Lebesgue, Sobolev and Bochner–Sobolev spaces (e.g. ( W k,r (Ω) , (cid:107) · (cid:107) W k,r (Ω) ) and ( L q (0 , T ; W n,r (Ω)) , (cid:107) · (cid:107) L q (0 ,T ; W n,r (Ω)) ) ). We will define W k,r (Ω) for r ∈ [1 , ∞ ) as the closure of the space of smooth functions with compactsupport C ∞ (Ω) with respect to the norm (cid:107) · (cid:107) W k,r (Ω) and we will denote the dualspace of W ,r (Ω) by W − ,r (cid:48) (Ω) . Here r (cid:48) is used to denote the Hölder conjugate of r ,i.e. the number defined by the relation /r + 1 /r (cid:48) = 1 . The duality pairing will bewritten in the usual way using brackets (cid:104)· , ·(cid:105) . The space of traces on the boundary offunctions in W ,r (Ω) will be denoted by W /r (cid:48) ,r ( ∂ Ω) .If X is a Banach space, C w ([0 , T ]; X ) will be used to denote the space of continuousfunctions in time with respect to the weak topology of X . For r ∈ [1 , ∞ ) we also definethe following useful subspaces: L r (Ω) := (cid:26) q ∈ L r (Ω) : (cid:90) Ω q = 0 (cid:27) ,L div (Ω) d := { v ∈ C ∞ (Ω) d : div v = 0 } (cid:107)·(cid:107) L ,W ,r , div (Ω) d := { v ∈ C ∞ (Ω) d : div v = 0 } (cid:107)·(cid:107) W ,r (Ω) ,L r tr ( Q ) d × d := { τ ∈ L r ( Q ) d × d : tr ( τ ) = 0 } ,L r sym ( Q ) d × d := { τ ∈ L r ( Q ) d × d : τ T = τ } . In the definition of the space L r tr ( Q ) d × d above, tr ( τ ) denotes the usual matrixtrace of the d × d matrix function τ . In the various estimates the letter c will denote ageneric positive constant whose exact value could change from line to line, wheneverthe explicit dependence on the parameters is not important. The following embeddings will be useful whenderiving various estimates. Assume that the Banach spaces ( W , W , W ) form aninterpolation triple in the sense that (cid:107) v (cid:107) W ≤ c (cid:107) v (cid:107) λW (cid:107) v (cid:107) − λW , for some λ ∈ (0 , , P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI and W (cid:44) → W (cid:44) → W . Then (cf. [56]) L r (0 , T ; W ) ∩ L ∞ (0 , T ; W ) (cid:44) → L r/λ (0 , T ; W ) ,for r ∈ [1 , ∞ ) and(2.1) (cid:107) v (cid:107) L r/λ (0 ,T ; W ) ≤ c (cid:107) v (cid:107) − λL ∞ (0 ,T ; W ) (cid:107) v (cid:107) λL r (0 ,T ; W ) . An example of an interpolation triple that can be combined with this result is givenby the Gagliardo–Nirenberg inequality, which states that for given p, r ∈ [1 , ∞ ) , thereis a constant c p,r > such that [20]:(2.2) (cid:107) v (cid:107) L s (Ω) ≤ c p,r (cid:107)∇ v (cid:107) λL r (Ω) (cid:107) v (cid:107) − λL p (Ω) ∀ v ∈ W ,r (Ω) ∩ L p (Ω) , provided that s ∈ [1 , ∞ ) and λ ∈ (0 , satisfy λ = p − s d − r + p . A particularly useful example can be obtained if we assume that r> dd +2 and take p = 2 and λ = dd +2 : (2.3) (cid:107) v (cid:107) L r ( d +2) d ( Q ) ≤ c (cid:107)∇ v (cid:107) λL r ( Q ) (cid:107) v (cid:107) − λL ∞ (0 ,T ; L (Ω)) ∀ v ∈ L r (0 , T ; W ,r (Ω)) ∩ L ∞ (0 , T ; L (Ω)) . In this work we will use Simon’scompactness lemma (see [60]) instead of the usual Aubin–Lions lemma to extractconvergent subsequences when taking the discretisation limit in the time–dependentproblem. Assume that X and H are Banach spaces such that the compact embedding X (cid:44) → (cid:44) → H holds. Simon’s lemma states that if U ⊂ L p (0 , T ; H ) , for some p ∈ [1 , ∞ ) ,and it satisfies: • U is bounded in L loc (0 , T ; X ) ; • (cid:82) T − (cid:15) (cid:107) v ( t + (cid:15), · ) − v ( t, · ) (cid:107) pH → , as (cid:15) → , uniformly for v ∈ U ;then U is relatively compact in L p (0 , T ; H ) .Let X and V be reflexive Banach spaces such that X (cid:44) → V densely and let V ∗ bethe dual space of V . The following continuity properties (see [56]) will be importantwhen identifying the initial condition: v ∈ L (0 , T ; V ∗ ) , ∂ t v ∈ L (0 , T ; V ∗ ) = ⇒ v ∈ C ([0 , T ]; V ∗ ) , (2.4) v ∈ L ∞ (0 , T ; X ) ∩ C w ([0 , T ]; V ) = ⇒ v ∈ C w ([0 , T ]; X ) . (2.5) In the mathe-matical analysis of these systems it is more convenient to work not with the function G , but with its graph A , which is introduced in the usual way:(2.6) ( D , S ) ∈ A ( · ) ⇐⇒ G ( · , S , D ) = . We will assume that A is a maximal monotone r -graph for some r > , which meansthat the following properties hold for almost every z ∈ Q :(A1) [ A includes the origin ] ( , ) ∈ A ( z ) .(A2) [ A is a monotone graph ] For every ( D , S ) , ( D , S ) ∈ A ( z ) , ( S − S ) : ( D − D ) ≥ . -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS A is maximal monotone ] If ( D , S ) ∈ R d × d sym × R d × d sym is such that (ˆ S − S ) : ( ˆ D − D ) ≥ for all ( ˆ D , ˆ S ) ∈ A ( z ) , then ( D , S ) ∈ A ( z ) .(A4) [ A is an r -graph ] There is a non-negative function m ∈ L ( Q ) and a constant c > such that S : D ≥ − m + c ( | D | r + | S | r (cid:48) ) for all ( D , S ) ∈ A ( z ) . (A5) [ Measurability ] The set-valued map z (cid:55)→ A ( z ) is L ( Q ) – ( B ( R d × d sym ⊗ R d × d sym )) measurable; here L ( Q ) denotes the family of Lebesgue measurable subsets of Q and B ( R d × d sym ) is the family of Borel subsets of R d × d sym .(A6) [ Compatibility ] For any ( D , S ) ∈ A ( z ) we have thattr ( D ) = 0 ⇐⇒ tr ( S ) = 0 . Assumption (A6) was not included in the original works [13, 14, 21], but it is neededfor consistency with the physical property that S is traceless if and only if the velocityfield is divergence-free (see the discussion in [62]). A very important consequence ofAssumption (A5) (see [62]) is the existence of a measurable function (usually calleda selection ) D : Q × R d × dsym → R d × dsym such that ( D ( z, σ ) , σ ) ∈ A ( z ) for all σ ∈ R d × dsym .In the existence results it will be useful to approximate the selection using smoothfunctions. To that end, let us define the mollification:(2.7) D k ( · , σ ) := (cid:90) R d × d sym D ( · , σ − τ ) ρ k ( τ ) d τ , where ρ k ( τ ) = k d ρ ( k τ ) , k ∈ N , and ρ ∈ C ∞ ( R d × d sym ) is a mollification kernel. It ispossible to check (see e.g. [62]) that this mollification satisfies analogous monotonicityand coercivity properties to those of the selection D , i.e. we have that • For every τ , τ ∈ R d × d sym and for almost every z ∈ Q the monotonicity condi-tion(2.8) ( D k ( z, τ ) − D k ( z, τ )) : ( τ − τ ) ≥ holds. • There is a constant C ∗ > and a nonnegative function g ∈ L ( Q ) such thatfor all k ∈ N , for every τ ∈ R d × d , and for almost every z ∈ Q we have(2.9) τ : D k ( z, τ ) ≥ − g ( z ) + C ∗ ( | τ | r (cid:48) + | D k ( z, τ ) | r ) . • For any sequence { S k } k ∈ N bounded in L r (cid:48) ( Q ) d × d , we have for arbitrary B ∈ R d × d sym and φ ∈ C ∞ ( Q ) with φ ≥ :(2.10) lim inf k →∞ (cid:90) Q ( D k ( · , S k ) − D ( · , B )) : ( S k − B ) φ ( · ) ≥ . It is important to remark that (2.8), (2.9) and (2.10) are the essential properties; theexplicit form (2.7) of the approximation to the selection is not very important. Thereare other ways to achieve the same result; for instance piecewise affine interpolation ora generalised Yosida approximation could also be used (see [61, 62]). The following isa localized version of Minty’s lemma that will aid in the identification of the implicitconstitutive relation (for a proof see [12]).
P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI
Lemma
Let A be a maximal monotone r -graph satisfying (A1)–(A4) for some r > . Suppose that { D n } n ∈ N and { S n } n ∈ N are sequences of functions defined on ameasurable set ˆ Q ⊂ Q , such that: ( D n ( · ) , S n ( · )) ∈ A ( · ) a.e. in ˆ Q, D n (cid:42) D , weakly in L r ( ˆ Q ) d × d , S n (cid:42) S , weakly in L r (cid:48) ( ˆ Q ) d × d , lim sup n →∞ (cid:90) ˆ Q S n : D n ≤ (cid:90) ˆ Q S : D . Then, ( D ( · ) , S ( · )) ∈ A ( · ) a.e. in ˆ Q. The goal of this work is to prove convergence of a three-field finite element approxi-mation of the following system:(2.11) ∂ t u − div ( S − u ⊗ u ) + ∇ p = f in (0 , T ) × Ω , div u = 0 in (0 , T ) × Ω , ( D ( u ) , S ) ∈ A ( · ) a.e. in (0 , T ) × Ω , u = on (0 , T ) × ∂ Ω , u (0 , · ) = u ( · ) in Ω , where A ( · ) satisfies (A1)–(A6). The next section introduces the notation and toolsthat will be useful in the analysis of the discrete problem. In this section, the notation and as-sumptions regarding the finite element approximation will be presented. Essentiallythe same arguments would work for any method based on a Galerkin approxima-tion, but here we will focus only on finite element methods. Consider a family oftriangulations {T n } n ∈ N of Ω satisfying the following assumptions: • (Affine equivalence). Given n ∈ N and an element K ∈ T n , there is an affineinvertible mapping F K : K → ˆ K , where ˆ K is the closed standard referencesimplex in R d . • (Shape-regularity). There is a constant c τ , independent of n , such that h K ≤ c τ ρ K for every K ∈ T n , n ∈ N , where h K := diam ( K ) and ρ K is the diameter of the largest inscribed ball. • The mesh size h n := max K ∈T n h K tends to zero as n → ∞ .Define the conforming finite element spaces associated with the triangulation T n : V n := (cid:110) v ∈ W , ∞ (Ω) d : v | K ◦ F − K ∈ ˆ P V , K ∈ T n , v | ∂ Ω = 0 (cid:111) ,M n := (cid:110) q ∈ L ∞ (Ω) : q | K ◦ F − K ∈ ˆ P M , K ∈ T n (cid:111) , Σ n := (cid:110) σ ∈ L ∞ (Ω) d × d : σ | K ◦ F − K ∈ ˆ P S , K ∈ T n (cid:111) , where ˆ P V ⊂ W , ∞ ( ˆ K ) d , ˆ P M ⊂ L ∞ ( ˆ K ) and ˆ P S ⊂ L ∞ ( ˆ K ) d × d are finite-dimensionalpolynomial subspaces on the reference simplex ˆ K . Each of these spaces will be as-sumed to have a finite and locally supported basis. As in the continuous case, it will -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS r > : M n := M n ∩ L r (cid:48) (Ω) , Σ n tr := Σ n ∩ L r tr (Ω) d × d , Σ n sym := Σ n ∩ L r sym (Ω) d × d ,V n div := (cid:26) v ∈ V n : (cid:90) Ω q div v = 0 , ∀ q ∈ M n (cid:27) , Σ n div ( f ) := (cid:26) σ ∈ Σ n sym : (cid:90) Ω σ : D ( v ) = (cid:104) f , v (cid:105) , ∀ v ∈ V n div (cid:27) . Assumption
For every s ∈ [1 , ∞ ) we have that inf v ∈ V n (cid:107) v − v (cid:107) W ,s (Ω) → as n → ∞ ∀ v ∈ W ,s (Ω) d , inf q ∈ M n (cid:107) q − q (cid:107) L s (Ω) → as n → ∞ ∀ q ∈ L s (Ω) , inf σ ∈ Σ n (cid:107) σ − σ (cid:107) L s (Ω) → as n → ∞ ∀ σ ∈ L s (Ω) d × d . Assumption Π n Σ ). For each n ∈ N there is a linear projector Π n Σ : L sym (Ω) d × d → Σ n sym such that: • (Preservation of divergence). For every σ ∈ L (Ω) d × d we have (cid:90) Ω σ : D ( v ) = (cid:90) Ω Π n Σ ( σ ) : D ( v ) ∀ v ∈ V n div . • ( L s –stability). For every s ∈ (1 , ∞ ) there is a constant c > , independent of n ,such that: (cid:107) Π n Σ σ (cid:107) L s (Ω) ≤ c (cid:107) σ (cid:107) L s (Ω) ∀ σ ∈ L s sym (Ω) d × d . Assumption Π nV ). For each n ∈ N there is a linear projector Π nV : W , (Ω) d → V n such that the following properties hold: • (Preservation of divergence). For every v ∈ W , (Ω) d we have (cid:90) Ω q div v = (cid:90) Ω q div (Π nV v ) ∀ q ∈ M n . • ( W ,s –stability). For every s ∈ (1 , ∞ ) there is a constant c > , independent of n , such that: (cid:107) Π nV v (cid:107) W ,s (Ω) ≤ c (cid:107) v (cid:107) W ,s (Ω) ∀ v ∈ W ,s (Ω) d . Assumption Π nM ). For each n ∈ N there is a linear projector Π nM : L (Ω) → M n such that for all s ∈ (1 , ∞ ) there is a constant c > , independentof n , such that: (cid:107) Π nM q (cid:107) L s (Ω) ≤ c (cid:107) q (cid:107) L s (Ω) ∀ q ∈ L s (Ω) . It is not difficult to show that the approximability and stability properties imply thatfor s ∈ [1 , ∞ ) we have: (cid:107) σ − Π n Σ σ (cid:107) L s (Ω) → as n → ∞ ∀ σ ∈ L s sym (Ω) d × d , (cid:107) v − Π nV v (cid:107) W ,s (Ω) → as n → ∞ ∀ v ∈ W ,s (Ω) d , (2.12) (cid:107) q − Π nM q (cid:107) L s (Ω) → as n → ∞ ∀ q ∈ L s (Ω) . P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI
Remark s ∈ (1 , ∞ ) , of two positive constants β s , γ s > , independent of n ,such that the following discrete inf-sup conditions hold: inf q ∈ M n sup v ∈ V n (cid:82) Ω q div v (cid:107) v (cid:107) W ,s (Ω) (cid:107) q (cid:107) L s (cid:48) (Ω) ≥ β s , (2.13) inf v ∈ V n div sup τ ∈ Σ n sym (cid:82) Ω τ : D ( v ) (cid:107) τ (cid:107) L s (cid:48) (Ω) (cid:107) v (cid:107) W ,s (Ω) ≥ γ s . (2.14) Example P – P element and the Taylor–Hood element P k – P k − for k ≥ d (see [5, 8, 21, 34, 18]). In addition to stability, the Scott–Vogelius elementalso satisfies the property that the discretely divergence-free velocities are pointwisedivergence-free (the stability can be guaranteed by assuming for example that themesh has been barycentrically refined, see [59]); another example of a velocity-pressurepair with this property is given by the Guzmán–Neilan element [37, 36]. To satisfyAssumption 2.5, one could use the Clément interpolant [17].Sometimes it is easier to prove the inf-sup condition directly. For example, if thespace of discrete stresses consists of discontinuous P k polynomials (with k ≥ ): Σ n = { σ ∈ L ∞ (Ω) d × d : σ | K ∈ P k ( K ) d × d , for all K ∈ T n } , and we have that D ( V n ) ⊂ Σ n (e.g. we could take the Taylor–Hood element P k +1 – P k for the velocity and the pressure), then the inf-sup condition follows from the factthat for s ∈ (1 , ∞ ) there is a constant c > , independent of h , such that for any σ ∈ Σ n there is τ ∈ Σ n such that [58]: (cid:90) Ω τ : σ = (cid:107) σ (cid:107) sL s (Ω) and (cid:107) τ (cid:107) L s (cid:48) (Ω) ≤ c (cid:107) σ (cid:107) s − L s (Ω) . In case a continuous piecewise polynomial approximation of the stress is preferred,one could use the conforming Crouzeix–Raviart element for the discrete velocity andpressure and the following space for the stress [57] : Σ n = { σ ∈ C (Ω) d × d : σ | K ∈ ( P ( K ) ⊕ B ) d × d , for all K ∈ T n } , where B := span { λ λ λ , λ λ λ , λ λ λ } , and { λ j } j =1 are barycentric coordinates on K . Remark V n div ⊂ W ,r , div (Ω) d , and D ( V n ) ⊂ Σ n , then the stress-velocity inf-supcondition also holds for the subspace of traceless stresses. Consequently, fewer degreesof freedom are needed to compute the stress unknowns. In this section we will describe the notation thatwill be used when performing the time discretisation of the problem. Let { τ m } m ∈ N be a sequence of time steps such that T /τ m ∈ N and τ m → , as m → ∞ . For each m ∈ N we define the equidistant grid: { t mj } T/τ m j =0 , t j = t mj := jτ m . -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS Q ji := ( t i , t j ) × Ω , where ≤ i ≤ j ≤ T /τ m . Also, given a set of functions { v j } T/τ m j =0 belonging to a Banach space X ,we can define the piecewise constant interpolant v ∈ L ∞ (0 , T ; X ) as:(2.15) v ( t ) := v j , t ∈ ( t j − , t j ] , j ∈ { , . . . , T /τ m } , and the piecewise linear interpolant ˜ v ∈ C ([0 , T ]; X ) as:(2.16) ˜ v ( t ) := t − t j − τ m v j + t j − tτ m v j − , t ∈ [ t j − , t j ] , j ∈ { , . . . , T /τ m } . For a given function g ∈ L p (0 , T ; X ) , with p ∈ [1 , ∞ ) , we define the time averages:(2.17) g j ( · ) := 1 τ m (cid:90) t j t j − g ( t, · ) d t, j ∈ { , . . . , T /τ m } . Then the piecewise constant interpolant g defined by (2.15) satisfies [56]:(2.18) (cid:107) g (cid:107) L p (0 ,T ; X ) ≤ (cid:107) g (cid:107) L p (0 ,T ; X ) , and(2.19) g → g strongly in L p (0 , T ; X ) , as m → ∞ .
3. Weak formulation.
In this section we will present a weak formulation forthe problem (2.11), where now we assume that f ∈ L r (cid:48) (0 , T ; W − ,r (cid:48) (Ω) d ) , u ∈ L div (Ω) d and the graph A satisfies the assumptions (A1)–(A6) for some r > dd +2 .Similarly to previous works on the analysis of implicitly constituted fluids, a Lipschitztruncation technique will be required when proving that the limit of the sequenceof approximate solutions satisfies the constitutive relation. The theory of Lipschitztruncation for time-dependent problems is not as well developed as in the steady case;here it will be necessary to work locally and the equation plays a vital role (severalversions of parabolic Lipschitz truncation have appeared in the literature, see e.g.[22, 14, 9, 23]). Since the pressure will not be present in the weak formulation, it willbe more convenient to use the construction developed in [9] because it preserves thesolenoidality of the velocity. The following lemma states the main properties of thissolenoidal Lipschitz truncation. Lemma
Let p ∈ (1 , ∞ ) , σ ∈ (1 , min( p, p (cid:48) )) and let Q = I × B ⊂ R × R be a parabolic cylinder, where I is an open interval and B is an open ball.Denote by αQ , where α > , the α -scaled version of Q keeping the barycenter thesame. Suppose { e l } l ∈ N is a sequence of divergence-free functions that is uniformlybounded in L ∞ ( I ; L σ ( B ) d ) and converges to zero weakly in L p ( I ; W ,p ( B ) d ) andstrongly in L σ ( Q ) d . Let { G l } l ∈ N and { G l } l ∈ N be sequences that converge to zeroweakly in L p (cid:48) ( Q ) d × d and strongly in L σ ( Q ) d × d , respectively. Define G l := G l + G l and suppose that, for any l ∈ N , the equation (3.1) (cid:90) Q ∂ t e l · w = (cid:90) Q G l : ∇ w ∀ w ∈ C ∞ , div ( Q ) d . is satisfied. Then there is a number j ∈ N , a sequence { λ l,j } l,j ∈ N with j ≤ λ l,j ≤ j +1 − , a sequence of functions { e l,j } l,j ∈ N ⊂ L ( Q ) d , a sequence of open sets B λ l,j ⊂ Q , for l, j ∈ N , and a function ζ ∈ C ∞ ( Q ) with Q ≤ ζ ≤ Q with thefollowing properties: P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI e l,j ∈ L q ( I ; W ,q , div ( B ) d ) for any q ∈ [1 , ∞ ) and supp ( e l,j ) ⊂ Q , for any j ≥ j and any l ∈ N ; e l,j = e j on Q \ B λ l,j , for any j ≥ j and any l ∈ N ; There is a constant c > such that lim sup l →∞ λ pl,j |B λ l,j | ≤ c − j , for any j ≥ j ; For j ≥ j fixed, we have as l → ∞ : e l,j → , strongly in L ∞ ( Q ) d , ∇ e l,j (cid:42) , weakly in L q ( Q ) d × d , ∀ q ∈ [1 , ∞ ); There is a constant c > such that: lim sup l →∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Q G l : ∇ e l,j (cid:12)(cid:12)(cid:12)(cid:12) ≤ c − j , for any j ≥ j ; There is a constant c > such that for any H ∈ L p (cid:48) ( Q ) d × d : lim sup l →∞ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Q ( G l + H ) : ∇ e l,j ζ B cλl,j (cid:12)(cid:12)(cid:12)(cid:12) ≤ c − j/p , for any j ≥ j . Before we presentthe weak formulation, let us define ˇ r := min (cid:26) r ( d + 2)2 d , r (cid:48) (cid:27) . The weak formulation for (2.11) then reads as follows.
Formulation ˇ A. Find functions S ∈ L r (cid:48) sym ( Q ) d × d ∩ L r (cid:48) tr ( Q ) d × d , u ∈ L r (0 , T ; W ,r , div (Ω) d ) ∩ L ∞ (0 , T ; L div (Ω) d ) ,∂ t u ∈ L ˇ r (0 , T ; ( W , ˇ r (cid:48) , div (Ω) d ) ∗ ) , such that (cid:104) ∂ t u , v (cid:105) + (cid:90) Ω ( S − u ⊗ u ) : D ( v ) = (cid:104) f , v (cid:105) ∀ v ∈ W , ˇ r (cid:48) , div (Ω) d , a.e. t ∈ (0 , T ) , ( D ( u ) , S ) ∈ A ( · ) , a.e. in (0 , T ) × Ω , ess lim t → + (cid:107) u ( t, · ) − u ( · ) (cid:107) L (Ω) = 0 . Remark r = 2 ) the pressure is only a distribution in time, when workingwith a no-slip boundary condition (see e.g. [31]). An integrable pressure can beobtained if Navier’s slip boundary condition is used instead [14], but in this work wewill confine ourselves to the more common no-slip boundary condition. -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS Remark u ∈ C ([0 , T ]; ( W , ˇ r (cid:48) , div (Ω) d ) ∗ ) (cid:44) → C w ([0 , T ]; ( W , ˇ r (cid:48) , div (Ω) d ) ∗ ) , and since ˇ r ≤ r (cid:48) we also know that L div (Ω) d (cid:44) → ( W , ˇ r (cid:48) , div (Ω) d ) ∗ . Combined with (2.5)this yields u ∈ C w ([0 , T ]; L div (Ω) d ) and hence the initial condition only makes sensea priori in this weaker sense. However, for this problem it will be proved that it alsoholds in the stronger sense described above.For a given time step τ m and j ∈ { , . . . , T /τ m } , let f j ∈ W − ,r (cid:48) (Ω) d and D kj : Ω × R d × d → R d × d be the time averages associated with f and D k , respec-tively (recall (2.17)). The time derivative will be discretised using an implicit Eulerscheme; higher order time stepping techniques might not be more advantageous herebecause higher regularity in time of weak solutions to the problem is not guaranteeda priori. The discrete formulation of the problem can now be introduced. Formulation ˇ A k , n , m , l . For j ∈ { , . . . , T /τ m } , find functions S k,n,m,lj ∈ Σ n sym and u k,n,m,lj ∈ V n div such that: (cid:90) Ω ( D kj ( · , S k,n,m,lj ) − D ( u k,n,m,lj )) : τ = 0 ∀ τ ∈ Σ n sym , τ m (cid:90) Ω ( u k,n,m,lj − u k,n,m,lj − ) · v + 1 l (cid:90) Ω | u k,n,m,lj | r (cid:48) − u k,n,m,lj · v + (cid:90) Ω ( S k,n,m,lj : D ( v ) + B ( u k,n,m,lj , u k,n,m,lj , v )) = (cid:104) f j , v (cid:105) ∀ v ∈ V n div , u k,n,m,l = P n div u . Here P n div : L (Ω) d → V n div is simply the L –projection defined through(3.2) (cid:90) Ω P n div v · w = (cid:90) Ω v · w ∀ w ∈ V n div . The form B is meant to represent the convective term and is defined for functions u , v , w ∈ C ∞ (Ω) d as: B ( u , v , w ) := − (cid:90) Ω u ⊗ v : D ( w ) , if V n div ⊂ W ,r , div (Ω) d , (cid:90) Ω u ⊗ w : D ( v ) − u ⊗ v : D ( w ) , otherwise . This definition guarantees that B ( v , v , v ) = 0 for every v for which this expressionis well defined, regardless of whether v is pointwise divergence-free or not, which isvery useful when obtaining a priori estimates; it reduces to the usual weak form ofthe convective term whenever the velocities are exactly divergence-free. It is nownecessary to check that B can be continuously extended to the spaces involving time.By standard function space interpolation, we have that for almost every t ∈ (0 , T ) : (cid:90) Ω | u ( t, · ) ⊗ v ( t, · ) : D ( w ( t, · )) | ≤ (cid:107) u ( t, · ) (cid:107) L r (Ω) (cid:107) v ( t, · ) (cid:107) L r (Ω) (cid:107) D ( w ( t, · )) (cid:107) L ˇ r (cid:48) (Ω) ≤ (cid:107) u ( t, · ) (cid:107) L r ( d +2) d (Ω) (cid:107) v ( t, · ) (cid:107) L r ( d +2) d (Ω) (cid:107) D ( w ( t, · )) (cid:107) L ˇ r (cid:48) (Ω) ≤ c (cid:107) u ( t, · ) (cid:107) W ,r (Ω) (cid:107) v ( t, · ) (cid:107) W ,r (Ω) (cid:107) w ( t, · ) (cid:107) W , ˇ r (cid:48) (Ω) . P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI
As in the steady case (cf. [21]), a more restrictive condition is needed in order tobound the additional term in B whenever the elements are not exactly divergence-free. Namely, if we assume that r ≥ d +1) d +2 (this is the analogue of the condition r ≥ dd +1 in the steady case) then there is a q ∈ (1 , ∞ ] such that r + dr ( d +2) + q = 1 ,and therefore (cid:90) Ω | u ( t, · ) ⊗ w ( t, · ) : D ( v ( t, · )) | ≤ (cid:107) u ( t, · ) (cid:107) L r ( d +2) d (Ω) (cid:107) D ( v ( t, · )) (cid:107) L r (Ω) (cid:107) w ( t, · ) (cid:107) L q (Ω) ≤ c (cid:107) u ( t, · ) (cid:107) W ,r (Ω) (cid:107) v ( t, · ) (cid:107) W ,r (Ω) (cid:107) w ( t, · ) (cid:107) W , ˇ r (cid:48) (Ω) . On the other hand, using Hölder’s inequality we can also obtain the estimate (cid:107)B ( u , v , w ) (cid:107) L (0 ,T ) ≤ (cid:107) u (cid:107) L r (cid:48) ( Q ) (cid:107) v (cid:107) L r (cid:48) ( Q ) (cid:107) w (cid:107) L r (0 ,T ; W ,r (Ω)) + (cid:107) u (cid:107) L r (cid:48) ( Q ) (cid:107) w (cid:107) L r (cid:48) ( Q ) (cid:107) v (cid:107) L r (0 ,T ; W ,r (Ω)) , which means that if the L r (cid:48) ( Q ) d norm of u is finite, then the additional restriction r ≥ d +1) d +2 is not needed. Moreover, this would also imply that the velocity is anadmissible test function, which is useful in the convergence analysis. This motivatesthe introduction of the penalty term in Formulation ˇ A k , n , m , l . Remark ˇ A k , n , m , l does not contain the pressure, in practicethe incompressibility condition is enforced through the addition of a Lagrange mul-tiplier p k,n,m,lj ∈ M n , which could be thought of as the pressure in the system (thereason for the omission of the pressure in the analysis is explained in Remark 3.2). Forthis reason it is necessary to consider additional assumptions that guarantee inf-supstability of the spaces V n and M n (see Assumptions 2.4 and 2.5). In case the problemdoes have an integrable pressure p , then it is expected that the sequence of discretepressures converges to it in L ( Q ) . Remark S : Q × R d × dsym → R d × dsym such that ( τ , S ( z, τ )) ∈ A ( z ) for all τ ∈ R d × dsym , and some modelscan be written more naturally with a selection of this form; the same analysis asthe one presented in this work can be applied to that situation. In fact, in practiceit is not necessary to find a selection in order to perform the computations, i.e. inthe simulations it is possible to work directly with the implicit function G . Whenperforming the analysis though, the function G is not appropriate because manydifferent expressions could lead to the same constitutive relation, but have differentmathematical properties. Remark H ( div ; Ω) , because for the unsteady problem we do not have at our disposal resultsthat guarantee the integrability of div S .In the next theorem, convergence of the sequence of discrete solutions to a weaksolution of the problem is proved. Since the ideas and arguments contained in theproof are similar to the ones presented in the previous sections and follow a similarapproach to [61], we will not include here all the details of the calculations unlessthere is a significant difference. Theorem
Assume that r > dd +2 , let { Σ n , V n , M n } n ∈ N be a family of finiteelement spaces satisfying Assumptions 2.2–2.4. Then for k, n, m, l ∈ N there exists asequence { ( S k,n,m,lj , u k,n,m,lj ) } T/τ m j =1 of solutions of Formulation ˇ A k , n , m , l , and a couple ( S , u ) ∈ L r (cid:48) sym ( Q ) d × d ∩ L r (cid:48) tr ( Q ) d × d × L r (0 , T ; W ,r , div (Ω) d ) ∩ L ∞ (0 , T ; L div (Ω) d ) such -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS that the corresponding time interpolants (recall (2.15) and (2.16) ) u k,n,m,l , ˜ u k,n,m,l and S k,n,m,l satisfy (up to a subsequence): S k,n,m,l (cid:42) S weakly in L r (cid:48) ( Q ) d × d , u k,n,m,l (cid:42) u weakly in L r (0 , T ; W ,r (Ω) d ) , (3.3) u k,n,m,l , ˜ u k,n,m,l ∗ (cid:42) u weakly* in L ∞ (0 , T ; L (Ω) d ) , and ( S , u ) solves Formulation ˇ A , with the limits taken in the order k → ∞ , ( n, m ) →∞ and l → ∞ .Proof. The idea of the proof is common in the analysis of nonlinear PDE: weobtain a priori estimates and use compactness arguments to pass to the limit inthe equation. In order to prove the existence of solutions of Formulation ˇ A k , n , m , l ,we need to check that given ( S k,n,m,lj − , u k,n,m,lj − ) , we can find ( S k,n,m,lj , u k,n,m,lj ) , for j ∈ { , . . . , T /τ m } . Testing the equation with ( S k,n,m,lj , u k,n,m,lj ) , we see that: (3.4) (cid:90) Ω D k ( · , S k,n,m,lj ) : S k,n,m,lj + 1 l (cid:107) u k,n,m,lj (cid:107) r (cid:48) L r (cid:48) (Ω) ≤ (cid:104) f , u k,n,m,lj (cid:105) + 1 τ m (cid:90) Ω u k,n,m,lj − · u k,n,m,lj . On the other hand, since all norms are equivalent in a finite-dimensional normed linearspace, there is a constant C n > such that:(3.5) (cid:107) v (cid:107) W ,r (Ω) ≤ C n (cid:107) v (cid:107) L r (cid:48) (Ω) ∀ v ∈ V n div . The constant C n may blow up as n → ∞ , but since n is fixed for now this does notpose a problem. Now, recalling (2.9) and combining (3.4) and (3.5) with a standardcorollary of Brouwer’s Fixed Point Theorem (cf. [33]) we obtain the existence of so-lutions of Formulation ˇ A k , n , m , l . In the first time step (i.e. j = 1 ), it is essential to usethe fact that the projection P n div is stable:(3.6) (cid:107) P n div u (cid:107) L (Ω) ≤ (cid:107) u (cid:107) L (Ω) . The estimate (3.5) suffices to guarantee the existence of discrete solutions, but inorder to pass to the limit n → ∞ , an estimate that does not degenerate as n → ∞ is required. This uniform estimate is a consequence of the discrete inf-sup condition(2.14):(3.7) γ r (cid:107) u k,n,m,lj (cid:107) W ,r (Ω) ≤ (cid:107) D k ( · , S k,n,m,lj +1 ) (cid:107) L r (Ω) . Therefore, the following a priori estimate holds: sup j ∈{ ,...,T/τ m } (cid:107) u k,n,m,lj (cid:107) L (Ω) + T/τ m (cid:88) j =1 (cid:107) u k,n,m,lj − u k,n,m,lj − (cid:107) L (Ω) + τ m T/τ m (cid:88) j =1 (cid:107) S k,n,m,lj (cid:107) L r (cid:48) (Ω) + τ m T/τ m (cid:88) j =1 (cid:107) u k,n,m,lj (cid:107) rW ,r (Ω) (3.8) + T/τ m (cid:88) j =1 (cid:107) D k ( · , · , S k,n,m,lj ) (cid:107) L r ( Q jj − ) + τ m l T/τ m (cid:88) j =1 (cid:107) u k,n,m,lj (cid:107) r (cid:48) L r (cid:48) (Ω) ≤ c, where c is a positive constant that depends on the data; in particular, c is indepen-dent of k, n, m and l . Let u k,n,m,l ∈ L ∞ (0 , T ; V n div ) and ˜ u k,n,m,l ∈ C ([0 , T ]; V n div ) P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI be the piecewise constant and piecewise linear interpolants defined by the sequence { u k,n,m,lj } T/τ m j =1 (see (2.15) and (2.16)) and let S k,n,m,l ∈ L ∞ (0 , T ; Σ n sym ) be the piece-wise constant interpolant defined by the sequence { S k,n,m,lj } T/τ m j =1 . Furthermore, definealso the piecewise constant interpolants: f ( t, · ) := f j ( · ) , D k ( t, · , · ) := D kj ( · , · ) , t ∈ ( t j − , t j ] , j ∈ { , . . . , T /τ m } Then the discrete formulation can be rewritten as: (cid:90) Ω ( D k ( t, · , S k,n,m,l ) − D ( u k,n,m,l )) : τ = 0 ∀ τ ∈ Σ n sym , (cid:90) Ω ∂ t ˜ u k,n,m,l · v + 1 l (cid:90) Ω | u k,n,m,l | r (cid:48) − u k,n,m,l · v + (cid:90) Ω ( S k,n,m,l : D ( v ) + B ( u k,n,m,l , u k,n,m,l , v )) = (cid:104) f , v (cid:105) ∀ v ∈ V n div , ˜ u k,n,m,l (0 , · ) = P n div u ( · ) . The a priori estimate (3.8) can in turn be written as: (cid:107) u k,n,m,l (cid:107) L ∞ (0 ,T ; L (Ω)) + τ m (cid:107) ∂ t ˜ u k,n,m,l (cid:107) L ( Q ) + (cid:107) S k,n,m,l (cid:107) r (cid:48) L r (cid:48) ( Q ) (3.9) + (cid:107) u k,n,m,l (cid:107) rL r (0 ,T ; W ,r (Ω)) + (cid:107) D k ( · , · , S k,n,m,l ) (cid:107) rL r ( Q ) + 1 l (cid:107) u k,n,m,l (cid:107) r (cid:48) L r (cid:48) ( Q ) ≤ c. Using the equivalence of norms in finite-dimensional spaces we also obtain (cid:107) ∂ t ˜ u k,n,m,l (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ c ( n ) (cid:107) ∂ t ˜ u k,n,m,l (cid:107) L ( Q ) , and together with the a priori estimate this implies that(3.10) (cid:107) ˜ u k,n,m,l (cid:107) W , ∞ (0 ,T ; L (Ω)) ≤ c ( n, m ) . Therefore, up to subsequences, as k → ∞ we have: u k,n,m,l → u n,m,l strongly in L ∞ (0 , T ; L (Ω) d ) , ˜ u k,n,m,l → ˜ u n,m,l strongly in W , ∞ (0 , T ; L (Ω) d ) , u k,n,m,l → u n,m,l strongly in L r (cid:48) ( Q ) d , u k,n,m,l → u n,m,l strongly in L r (0 , T ; W ,r (Ω) d ) , S k,n,m,l → S n,m,l strongly in L r (cid:48) ( Q ) d × d , D k ( · , · , S k,n,m,l ) (cid:42) D n,m,l weakly in L r ( Q ) d × d , D k ( · , · , S k,n,m,l ) (cid:42) D n,m,l weakly in L r ( Q ) d × d , D kj ( · , S k,n,m,lj ) (cid:42) D n,m,lj weakly in L r (Ω) d × d , for j ∈ { , . . . , T /τ m } . Since the function D kj is simply an average in time, the uniqueness of the weak limitimplies that(3.11) D n,m,lj ( · ) = 1 τ m (cid:90) t j t j − D n,m,l ( t, · ) d t, j ∈ { , . . . , T /τ m } , and that D n,m,l is the piecewise constant interpolant determined by the sequence { D n,m,lj } T/τ m j =1 . Moreover, since the convergence of the velocity and stress sequences -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS k → ∞ and thus we obtain (cid:90) Ω ( D n,m,l − D ( u n,m,l )) : τ = 0 ∀ τ ∈ Σ n sym , (cid:90) Ω ∂ t ˜ u n,m,l · v + 1 l (cid:90) Ω | u n,m,l | r (cid:48) − u n,m,l · v + (cid:90) Ω ( S n,m,l : D ( v ) + B ( u n,m,l , u n,m,l , v )) = (cid:104) f , v (cid:105) ∀ v ∈ V n div . It is also clear that the initial condition ˜ u n,m,l (0 , · ) = P n div u ( · ) holds, since the expres-sion on the right-hand side is independent of k . The identification of the constitutiverelation can be carried out using (2.10) in exactly the same manner as in [61], whichmeans that (the strong convergence is again essential):(3.12) ( D n,m,l , S n,m,l ) ∈ A ( · ) , a.e. in (0 , T ) × Ω . The next step is to take the limit in both the time and space discretisations simul-taneously. The weak lower semicontinuity of the norms and the estimate (3.9) implythat: (cid:107) u n,m,l (cid:107) L ∞ (0 ,T ; L (Ω)) + τ m (cid:107) ∂ t ˜ u n,m,l (cid:107) L ( Q ) + (cid:107) S n,m,l (cid:107) r (cid:48) L r (cid:48) ( Q ) (3.13) + (cid:107) u n,m,l (cid:107) rL r (0 ,T ; W ,r (Ω)) + (cid:107) D n,m,l (cid:107) rL r ( Q ) + 1 l (cid:107) u n,m,l (cid:107) r (cid:48) L r (cid:48) ( Q ) ≤ c, and(3.14) (cid:107) ˜ u n,m,l (cid:107) L ∞ (0 ,T ; L (Ω)) = (cid:107) u n,m,l (cid:107) L ∞ (0 ,T ; L (Ω)) ≤ c, where c is a constant, independent of n, m and l . Consequently, there exist (notrelabelled) subsequences such that, as n, m → ∞ : u n,m,l ∗ (cid:42) u l weakly* in L ∞ (0 , T ; L (Ω) d ) , ˜ u n,m,l ∗ (cid:42) u l weakly* in L ∞ (0 , T ; L (Ω) d ) , u n,m,l (cid:42) u l weakly in L r (0 , T ; W ,r (Ω) d ) , S n,m,l (cid:42) S l weakly in L r (cid:48) ( Q ) d × d , D n,m,l (cid:42) D l weakly in L r ( Q ) d × d , D n,m,l (cid:42) D l weakly in L r ( Q ) d × d , l (cid:90) Q | u n,m,l | r (cid:48) − u n,m,l (cid:42) l (cid:90) Q | u l | r (cid:48) − u n,m,l weakly in L (2 r (cid:48) ) (cid:48) ( Q ) d . At this point it is a standard step to use the Aubin–Lions lemma to obtain strongconvergence of subsequences. However, following [61], we will instead use Simon’scompactness lemma; this choice is made to avoid the need for stability estimates of P n div in Sobolev norms, which would require additional assumptions on the mesh. Toapply this lemma, it will be more convenient to work with the modified interpolant: ˆ u n,m,l ( t, · ) := u n,m,l ( · ) , if t ∈ [0 , t ) , ˜ u n,m,l ( t, · ) , if t ∈ [ t , T ] . P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI
Let (cid:15) > be such that s + (cid:15) < T and let v ∈ V n div . Then, using the definition of ˆ u n,m,l we have (cid:90) Ω ( ˆ u n,m,l ( s + (cid:15), x ) − ˆ u n,m,l ( s + (cid:15), x )) · v ( x ) d x = (cid:90) s + (cid:15) max( s,τ m ) (cid:90) Ω ∂ t ˆ u n,m,l ( t, x ) · v ( x ) d x d t = (cid:90) s + (cid:15) max( s,τ m ) (cid:90) Ω ∂ t ˜ u n,m,l ( t, x ) · v ( x ) d x d t = (cid:90) s + (cid:15) max( s,τ m ) (cid:18) − l (cid:90) Ω | u n,m,l ( t, x ) | r (cid:48) − u n,m,l ( t, x ) · v ( x ) d x − (cid:90) Ω ( S n,m,l ( t, x ) : D ( v ( x )) + B ( u n,m,l ( t, x ) , u n,m,l ( t, x ) , v ( x ))) d x + (cid:104) f ( t ) , v (cid:105) (cid:19) d t ≤ c ( l ) (cid:32)(cid:90) s + (cid:15) max( s,τ m ) (cid:107) v (cid:107) rW ,r (Ω) d t (cid:33) /r + (cid:32)(cid:90) s + (cid:15) max( s,τ m ) (cid:107) v (cid:107) r (cid:48) L r (cid:48) (Ω) d t (cid:33) / r (cid:48) ≤ c ( l )( (cid:15) /r + (cid:15) / r (cid:48) ) (cid:16) (cid:107) v (cid:107) W ,r (Ω) + (cid:107) v (cid:107) L r (cid:48) (Ω) (cid:17) . Choosing v = ˆ u n,m,l ( s + (cid:15), · ) − ˆ u n,m,l ( s, · ) we conclude that (cid:90) T − (cid:15) (cid:107) ˆ u n,m,l ( s + (cid:15), · ) − ˆ u n,m,l ( s, · ) (cid:107) L (Ω) d s → , as (cid:15) → . On the other hand, the a priori estimates imply that ˆ u n,m,l is bounded (uniformly in n, m ∈ N ) in L ( Q ) d and L (0 , T ; W ,r (Ω) d ) . Moreover, since r > dd +2 , the embedding W ,r (Ω) d (cid:44) → L (Ω) d is compact and thus Simon’s compactness lemma guarantees thestrong convergence:(3.15) ˆ u n,m,l → u l strongly in L ( Q ) d . Since the interpolants converge to the same limit as τ m → , using standard functionspace interpolation (and recalling (2.3)) we also obtain that, as n, m → ∞ : ˜ u n,m,l → u l strongly in L p (0 , T ; L (Ω) d ) , (3.16) u n,m,l → u l strongly in L p (0 , T ; L (Ω) d ) ∩ L q ( Q ) , (3.17)for p ∈ [1 , ∞ ) and q ∈ [1 , max(2 r (cid:48) , q ( d +2) d )) .Now, using the property (2.12), we can check that u l is actually divergence-free: (3.18) (cid:90) T (cid:90) Ω φ Π nM q div u n,m,l → (cid:90) T (cid:90) Ω φ q div u l ∀ q ∈ L r (cid:48) (Ω) , φ ∈ C ∞ (0 , T ) . Furthermore, (2.12) also yields convergence of the initial condition, as n, m → ∞ :(3.19) ˜ u n,m,l (0 , · ) = P n div u → u strongly in L (Ω) d . The functions D l and D l can easily be identified using the property (2.19) and thedefinition of the piecewise constant interpolant (3.11). Indeed, for an arbitrary σ ∈ C ∞ ( Q ) we have, as n, m → ∞ :(3.20) (cid:90) T (cid:90) Ω D n,m,l : σ = (cid:90) T (cid:90) Ω D n,m,l : σ → (cid:90) T (cid:90) Ω D l : σ . -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS D l = D l .Combining all these properties and using an analogous computation to (3.18) itis possible to prove that the limiting functions are a solution of the following problem: (cid:90) T (cid:90) Ω ( D l − D ( u l )) : τ ϕ = 0 ∀ τ ∈ C ∞ , sym (Ω) d × d , ϕ ∈ C ∞ (0 , T ) , − (cid:90) T (cid:90) Ω u l · v ∂ t ϕ − (cid:90) Ω u · v ϕ (0) + (cid:90) T (cid:90) Ω ( S l − u l ⊗ u l ) : D ( v ) ϕ + 1 l (cid:90) T (cid:90) Ω | u l | r (cid:48) − u l · v ϕ = (cid:90) T (cid:104) f , v (cid:105) ϕ ∀ v ∈ C ∞ , div (Ω) d , ϕ ∈ C ∞ ( − T, T ) . From the equation above and the estimate (2.3) we then see that the distributionaltime derivative belongs to the spaces: ∂ t u l ∈ L min( r (cid:48) , (2 r (cid:48) ) (cid:48) ) (0 , T ; ( W ,r , div (Ω) d ∩ L r (cid:48) (Ω) d ) ∗ ) , (3.21) ∂ t u l ∈ L min(ˇ r, (2 r (cid:48) ) (cid:48) ) (0 , T ; ( W , ˇ r (cid:48) , div (Ω) d ) ∗ ) . (3.22)It is important to note that (3.22) holds uniformly in l ∈ N , while (3.21) does not.Now, observe that W ,r , div (Ω) d ∩ L r (cid:48) (Ω) d (cid:44) → L div (Ω) d (cid:44) → ( L div (Ω) d ) ∗ (cid:44) → ( W ,r , div (Ω) d ∩ L r (cid:48) (Ω) d ) ∗ . Combining this with (2.4), (2.5), and the fact that u l ∈ L ∞ (0 , T ; L div (Ω) d ) guaranteesthat u l ∈ C w ([0 , T ] , L div (Ω) d ) . Let v ∈ C ∞ , div (Ω) d and ϕ ∈ C ∞ ( − T, T ) be such that ϕ (0) = 1 ; then the following equality holds:(3.23) (cid:90) T (cid:90) Ω ∂ t ( u l ϕ ) · v = − (cid:90) Ω u l (0 , · ) · v ϕ (0) . On the other hand, using the equation we also have that: (3.24) (cid:90) T (cid:90) Ω ∂ t ( u l ϕ ) · v = (cid:90) T (cid:90) Ω ∂ t u l · v ϕ + (cid:90) T (cid:90) Ω u l · v ∂ t ϕ = − (cid:90) Ω u · v ϕ (0) . Comparing (3.23) and (3.24) we conclude that u l (0 , · ) = u ( · ) . This proves that theinitial condition is attained in the weak sense expected a priori from the embeddings;however, in this case the stronger condition(3.25) ess lim t → + (cid:107) u l ( t, · ) − u ( · ) (cid:107) L (Ω) = 0 holds. To see this, note that (3.16) guarantees that, up to a subsequence, ˜ u n,m,l ( t, · ) → ˜ u l ( t, · ) in L (Ω) d for almost every t ∈ [0 , T ] , and therefore (cid:107) u l ( t, · ) − u ( · ) (cid:107) L (Ω) = lim sup n,m →∞ (cid:107) ˜ u n,m,l ( t, · ) − ˜ u n,m,l (0 , · ) (cid:107) L (Ω) = lim sup n,m →∞ (cid:16) (cid:107) ˜ u n,m,l ( t, · ) (cid:107) L (Ω) − (cid:107) ˜ u n,m,l (0 , · ) (cid:107) L (Ω) +2 (cid:90) Ω ( ˜ u n,m,l (0 , · ) − ˜ u n,m,l ( t, · )) · ˜ u n,m,l (0 , · ) (cid:19) ≤ lim sup n,m →∞ (cid:18)(cid:90) t (cid:104) f , u n,m,l (cid:105) + 2 (cid:90) Ω ( ˜ u n,m,l (0 , · ) − ˜ u n,m,l ( t, · )) · ˜ u n,m,l (0 , · ) (cid:19) ≤ (cid:90) t (cid:104) f , u l (cid:105) + 2 (cid:90) Ω ( u l (0 , · ) − u l ( t, · )) · u l (0 , · ) , P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI for almost every t ∈ [0 , T ] . Observe also that the monotonicity of the constitutiverelation was used to obtain the next to last inequality. Taking the limit t → + thenyields (3.25).The identification of the constitutive relation, i.e. proving that ( D l , S l ) ∈ A ( · ) almost everywhere, can be carried out with the help of Lemma 2.1. In order to applythe lemma, the only thing that remains to be proved, since we already know that ( D n,m,l , S n,m,l ) ∈ A ( · ) almost everywhere, is that:(3.26) lim sup n,m →∞ (cid:90) t (cid:90) Ω S n,m,l : D n,m,l ≤ (cid:90) t (cid:90) Ω S l : D l , for almost every t ∈ [0 , T ] ; then taking t → T we obtain the result in the wholedomain Q . The proof of this fact is essentially the same as in [61] and we will notreproduce it here. Moreover, the following energy identity holds: (3.27) (cid:107) u l ( t, · ) (cid:107) L (Ω) + (cid:90) t (cid:90) Ω S l : D ( u l ) + 1 l (cid:90) t (cid:107) u l (cid:107) r (cid:48) L r (cid:48) (Ω) = (cid:90) t (cid:104) f , u l (cid:105) + (cid:107) u (cid:107) L (Ω) , In time-dependent problems obtaining an energy identity of this kind is not alwayspossible; in this case the energy equality (3.27) can be proved, since the velocity is anadmissible test function in space thanks to the fact that its L r (cid:48) norm is under control(some mollification is needed to overcome the low integrability in time, see [62, 44]).Now, (3.13) and the weak and weak* lower semicontinuity of the norms implythat (3.28) (cid:107) u l (cid:107) L ∞ (0 ,T ; L (Ω)) + (cid:107) S l (cid:107) r (cid:48) L r (cid:48) ( Q ) + (cid:107) u l (cid:107) rL r (0 ,T ; W ,r (Ω)) + (cid:107) D l (cid:107) rL r ( Q ) + 1 l (cid:107) u l (cid:107) r (cid:48) L r (cid:48) ( Q ) ≤ c, where c is a constant independent of l . From this we see that, up to subsequences, as l → ∞ : u l ∗ (cid:42) u weakly* in L ∞ (0 , T ; L (Ω) d ) , u l (cid:42) u weakly in L r (0 , T ; W ,r (Ω) d ) , S l (cid:42) S weakly in L r (cid:48) ( Q ) d × d , (3.29) D l (cid:42) D weakly in L r ( Q ) d × d , l (cid:90) Q | u l | r (cid:48) − u l → strongly in L ( Q ) d . Furthermore, since ˇ r ≤ r (cid:48) and r > dd +2 , the embedding W , ˇ r (cid:48) , div (Ω) d (cid:44) → L div (Ω) d iscompact and hence by the Aubin–Lions lemma (taking into account (3.22)) we havethe strong convergence:(3.30) u l → u strongly in L r (0 , T ; L div (Ω) d ) . With the convergence properties (3.29) and (3.30) it is then possible to pass to thelimit and prove that the limiting functions satisfy: (cid:90) Ω ( D − D ( u )) : τ = 0 ∀ τ ∈ C ∞ , sym (Ω) d × d , a.e. t ∈ (0 , T ) , (cid:104) ∂ t u , v (cid:105) + (cid:90) Ω ( S − u ⊗ u ) : D ( v ) = (cid:104) f , v (cid:105) ∀ v ∈ C ∞ , div (Ω) d , a.e. t ∈ (0 , T ) . -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS ess lim t → + (cid:107) u ( t, · ) − u ( · ) (cid:107) L (Ω) = 0 . Moreover, since the penalty term vanishes in the limit l → ∞ , we can improve theintegrability in time:(3.32) ∂ t u l ∈ L ˇ r (0 , T ; ( W , ˇ r (cid:48) , div (Ω) d ) ∗ ) . To show that ( D , S ) ∈ A ( · ) , Lemma 2.1 will once again be employed. The maindifficulty at this stage, just like in the previous works [21, 61], is that the velocity isno longer an admissible test function (and therefore we do not have an energy equalitysimilar to (3.27)). The idea is now to work with Lipschitz truncations of the error e l := u l − u ; it should be noted however that in the present case we need to verify anumber of additional hypotheses before Lemma 3.1 can be applied.Note that equation (3.1) in Lemma 3.1 is written in divergence form. We thenneed to make a preliminary step and write the penalty term in this form (see [61]).Let B ⊂⊂ Ω be an arbitrary ball compactly contained in Ω and let q ∈ [1 , (2 r (cid:48) ) (cid:48) ) .Then from the standard theory of elliptic operators we know that for almost every t ∈ [0 , T ] there is a unique g l ( t, · ) ∈ W ,q ( B ) d ∩ W ,q ( B ) such that: (cid:90) B ∇ g l ( t, · ) : ∇ v = 1 l (cid:90) B | u l ( t, · ) | r (cid:48) − u l ( t, · ) · v ∀ v ∈ C ∞ , div (Ω) d , (cid:107) g l ( t, · ) (cid:107) W ,q ( B ) ≤ c (cid:13)(cid:13)(cid:13)(cid:13) l | u l ( t, · ) | r (cid:48) − u l ( t, · ) (cid:13)(cid:13)(cid:13)(cid:13) L q ( B ) . This means in particular (by (3.29) and standard function space interpolation) thatfor a fixed time interval I ⊂⊂ (0 , T ) we have:(3.33) g l → strongly in L q ( I ; W ,q ( B ) d ) , ∀ q ∈ [1 , (2 r (cid:48) ) (cid:48) ) . Defining Q := I × B and G l := S l − S , G l := u l ⊗ u l − u ⊗ u − ∇ g l , we readily see that the error e l satisfies the equation(3.34) (cid:90) Q ∂ t e l · w = (cid:90) Q ( G l + G l ) : ∇ w ∀ w ∈ C ∞ , div ( Q ) d . Additionally, as a consequence of (3.29), (3.33) and (3.30) we also have that for any q ∈ [1 , min(ˇ r, (2 r (cid:48) ) (cid:48) ) , the sequence u l is bounded in L ∞ ( I ; W ,q ( Q ) d ) and that: G l (cid:42) weakly in L r (cid:48) ( Q ) d × d , G l → strongly in L q ( Q ) d × d , u l → u strongly in L q ( Q ) d . Consequently, the assumptions of Lemma 3.1 are satisfied. It now suffices to provefor an arbitrary θ ∈ (0 , that(3.35) lim sup l →∞ (cid:90) Q [( D ( u l ) − D ( · , S )) : ( S l − S )] θ ≤ , P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI
Once this has been shown, Chacon’s biting lemma and Vitali’s convergence theoremwill imply, together with Lemma 2.1, that ( D , S ) ∈ A ( · ) almost everywhere in Q (see the details e.g. in [14]). From here then the result follows by observing that Q can be covered by a union of such cylinders (e.g. by using a Whitney covering).In order to prove (3.35), first let B λ l,j ⊂ Ω be the family of open sets and let { e l,j } l,j ∈ N be the sequence of Lipschitz truncations described in Lemma 3.1. If wedefine(3.36) H l ( · ) := ( D ( u l ) − D ( · , S )) : ( S l − S ) ∈ L ( Q ) , then we have by Hölder’s inequality that (cid:90) Q | H l | θ ≤ | Q | − θ (cid:32)(cid:90) Q \B λl,j H l (cid:33) θ + |B λ l,j | − θ (cid:32)(cid:90) Q H l (cid:33) θ . The second term on the right-hand side can be dealt with easily, since H l is boundeduniformly in L ( Q ) thanks to the a priori estimate (3.28), and the properties describedin Lemma 3.1 imply that(3.37) lim sup l →∞ |B λ l,j | − θ ≤ lim sup l →∞ | λ rl,j B λ l,j | − θ ≤ c − j (1 − θ ) , for j ≥ j , where c is a positive constant. For the first term, observe that (cid:90) Q \B λl,j H l = (cid:90) Q H l ζ B cλl,j = (cid:90) Q D ( e l ) : ( S l − S ) ζ B cλl,j + (cid:90) Q \B λl,j ( D ( u ) − D ( · , S )) : ( S l − S ) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Q D ( e l,j ) : G l ζ B cλl,j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Q ( D ( u ) − D ( · , S )) : ( S l − S ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) B λl,j ( D ( u ) − D ( · , S )) : ( S l − S ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , where ζ ∈ C ∞ , div ( Q ) is the function introduced in Lemma 3.1. Taking lim sup l →∞ the assertion follows by taking j → ∞ . In particular, we used for the first termLemma 3.1 part 6, with H = , for the second term the weak convergence of S l andfor the third term the fact that { S l } l ∈ N is bounded, together with (3.37). To concludethe proof, note that the fact that u is divergence-free and Assumption (A6) implythat tr ( S ) = 0 , and so S ∈ L r (cid:48) sym (Ω) d × d ∩ L r (cid:48) tr (Ω) d × d . Remark ˇ A k , n , m , l is a four-step approximation in which the in-dices k, n, m, l refer to the approximation of the graph by smooth functions, the finiteelement discretisation, the discretisation in time, and the penalty term, respectively.The same approach can be used to define a 3-field formulation for the steady prob-lem and the unsteady problem without convection and the proof remains valid withsome simplifications; for instance, for the steady system without convective term,only the indices k and n are needed. Furthermore, in those cases the convergence ofthe sequence of discrete pressures can be guaranteed in the corresponding Lebesguespaces. -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS Remark (cid:107) u k,n,m,lj (cid:107) W ,r (Ω) cannot be deduced from Formulation ˇ A k , n , m , l by simplytesting with the solution. An alternative approach could be to include in the equationan additional diffusion term of the form: k (cid:90) Ω | D ( u k,n,m,lj ) | r − D ( u k,n,m,lj ) : D ( v ) , which would be completely acceptable if we only cared about the existence of weaksolutions, but is undesirable from the point of view of the computation of the finiteelement approximations, since it introduces an additional nonlinearity in the discreteproblem. Remark k → ∞ , ( n, m ) → ∞ and l → ∞ were taken successively. In contrast to the steady case considered in [21],here it is not known whether we can take the limits at once. The result is likely tohold as well, but the proof would require a discrete version of the parabolic Lipschitztruncation, which is not available at the moment. Remark ( D , S , u , p ) . The only additional assumption needed in that case would bean inf-sup condition of the form:(3.38) inf σ ∈ Σ n div ( ) sup τ ∈ Σ n sym (cid:82) Ω σ : τ (cid:107) σ (cid:107) L s (cid:48) (Ω) (cid:107) τ (cid:107) L s (Ω) ≥ δ s , where δ s > is independent of n .
4. Numerical experiments.
According to the analysis carried out in the previ-ous section, the addition of the penalty term is necessary when r ∈ ( dd +2 , d +2 d +2 ] . How-ever, in the examples we observed that the method converges regardless of whetherthe penalty term is present or not. This could be an indication that the requirementto include this penalty term is only a technical obstruction and that there might be adifferent approach to showing convergence of the numerical method that could avoidits inclusion in the numerical method. On the other hand, it could also be the casethat exact solutions with more severe singularities than the ones considered in our nu-merical experiments are needed to demonstrate pathological behaviour. In any case,it appears that in most applications the penalty term can be safely omitted and forthis reason it is not discussed in the numerical examples below. The framework presentedin this work is so broad that in general it is not possible to guarantee uniqueness ofsolutions; in particular it is not clear how error estimates could be obtained. However,as this computational example will show, the discrete formulations presented hereappear to recover the expected orders of convergence in the cases where these ordersare known.In the first part of this numerical experiment we solved the steady problem with-out convection with the Carreau constitutive law (as stated in Remark 3.8, the same3-field approximation can be applied in this setting):(4.1) S ( D ) := 2 ν (cid:0) ε + | D | (cid:1) r − D , P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI where r ≥ and ε, ν > . This is one of the most common non-Newtonian models thatpresent a power-law structure (note that for r = 2 we recover the Newtonian model),and has the advantage that it is not singular at the origin (i.e. when D = ), unlikethe usual power-law constitutive relation. Observe that the constitutive relation issmooth, and therefore only the limit n → ∞ is needed in the results from the previoussection. The problem was solved on the unit square Ω = (0 , with a Dirichletboundary condition for the velocity defined so as to match the value of the exactsolution, which was chosen as:(4.2) u ( x ) = | x | a − ( x , − x ) T , p ( x ) = | x | b , where a, b are parameters used to control the smoothness of the solutions. Define theauxiliary function F := R d × d → R d × d sym as:(4.3) F ( B ) := ( ε + | B sym | ) r − B sym , where B sym := ( B + B T ) . In [5, 38] it was proved for systems of the form (4.1) thatif F ( D ( u )) ∈ W , (Ω) d × d and p ∈ W ,r (cid:48) (Ω) then the following error estimates hold: (cid:107) F ( D ( u )) − F ( D ( u n )) (cid:107) L (Ω) ≤ ch min { , r (cid:48) } n , (cid:107) p − p n (cid:107) L r (cid:48) (Ω) ≤ ch min { r (cid:48) , r (cid:48) } n . In our case, the conditions F ( D ( u )) ∈ W , (Ω) d × d and p ∈ W ,r (cid:48) (Ω) amount torequiring that a > and b > r − . These parameters were then chosen to be a = 1 . and b = r − . in order to be close to the regularity threshold. Wediscretised this problem with the Scott–Vogelius element for the velocity and pressureand discontinuous piecewise polynomials for the stress variables: Σ n = { σ ∈ L ∞ (Ω) d × d : σ | K ∈ P k ( K ) d × d , for all K ∈ T n } ,V n = { w ∈ W ,r (Ω) d : w | ∂ Ω = u , w | K ∈ P k +1 ( K ) d for all K ∈ T n } ,M n = { q ∈ L ∞ (Ω) : q | k ∈ P k ( K ) for all K ∈ T n } . The problem was solved using firedrake [55] with ν = 0 . , ε = 10 − and k = 1 on abarycentrically refined mesh (obtained using gmsh [32]) to guarantee inf-sup stability.The discretised nonlinear problems were linearised using Newton’s method with the L line search algorithm of PETSc [3, 11]; the Newton solver was deemed to haveconverged when the Euclidean norm of the residual fell below × − . The linearsystems were solved with a sparse direct solver from the umfpack library [19]. In theimplementation, the uniqueness of the pressure was recovered not by using a zeromean condition but rather by orthogonalising against the nullspace of constants. Theexperimental orders of convergence in the different norms are shown in Tables 1 and 2(note that the tables do not contain the values of the numerical error, but rather theorder of convergence corresponding to the norm indicated in each column).From Tables 1 and 2 it can be seen that the algorithm recovers the expectedorders of convergence. In the case of the stress we obtain the same order as for thepressure, which seems natural from the point of view of the equation. In [38] it isclaimed that for r < the order of convergence for the velocity should be equal to 1;in our numerical simulations the experimental order of convergence seems to approach r , which is slightly larger than 1. This difference may be due to the fact that in [38] -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS r = 1 . . h n (cid:107) F ( D ( u )) (cid:107) L (Ω) (cid:107) u (cid:107) W ,r (Ω) (cid:107) p (cid:107) L r (cid:48) (Ω) (cid:107) S (cid:107) L r (cid:48) (Ω) r = 1 . . h n (cid:107) F ( D ( u )) (cid:107) L (Ω) (cid:107) u (cid:107) W ,r (Ω) (cid:107) p (cid:107) L r (cid:48) (Ω) (cid:107) S (cid:107) L r (cid:48) (Ω) u ( t, x ) = t | x | a − ( x , − x ) T , p ( t, x ) = t | x | b . In [25], the following error estimate for the approximation of time-dependent systemsof this form, but without convection, was obtained for r ∈ [ dd +2 , ∞ ) : (cid:107) u − u n,m (cid:107) L ∞ (0 ,T ; L (Ω)) + (cid:107) F ( D ( u )) − F ( D ( u n,m )) (cid:107) L ( Q ) ≤ c (cid:16) τ m + h min { , r } n (cid:17) , assuming that u ∈ W ,r , div (Ω) d and that the following additional regularity propertiesof the solution and the data hold: (cid:107)∇ F ( D ( u )) (cid:107) L (Ω) + (cid:107)∇ S ( D ( u )) (cid:107) L (Ω) ≤ c, (cid:107) u (cid:107) W , (0 ,T ; L (Ω)) + (cid:107) u (cid:107) L (0 ,T ; W , (Ω)) + (cid:107) F ( D ( u )) (cid:107) L (0 ,T ; W , (Ω)) ≤ c. The same order of convergence was obtained in [6] for r ∈ ( , in 3D for a semi-implicit discretisation of the unsteady system with convection assuming that u ∈ W , , div (Ω) d , div S ( D ( u )) ∈ L (Ω) d and that the slightly different regularity assump-4 P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI tions hold: (cid:107) ∂ t u (cid:107) L ∞ (0 ,T ; L (Ω)) + (cid:107) F ( D ( u )) (cid:107) W , ( Q ) + (cid:107) F ( D ( u )) (cid:107) L r − / (2 − r )) (0 ,T ; W , ) ≤ c. The problem was solved until the final time T = 0 . with the same parameters asabove; observe that this choice of parameters guarantees that the required regularityproperties are satisfied. Table 3 shows the experimental order of convergence for r = 1 . . The order of convergence for the natural norm (cid:107) F ( D ( u )) (cid:107) L ( Q ) agrees withTable 3: Experimental order of convergence for the full problem with r = 1 . . h n τ m (cid:107) F ( D ( u )) (cid:107) L ( Q ) (cid:107) u (cid:107) L ∞ (0 ,T ; L (Ω)) In this section we will considerthe classical lid–driven cavity problem with the non–standard constitutive relation:(4.4) D = δ s S | S | + ν S , if | D | ≥ δ s , S = 0 , if | D | < δ s , if ( x − ) + ( y − ) ≤ ( ) , D = ν S , otherwise , where ν > is the viscosity and δ s ≥ . This is an example of an activated fluid thatin the middle of the domain transitions between a Newtonian fluid (i.e. Navier–Stokes)and an inviscid fluid (i.e. Euler) depending on the magnitude of the symmetric velocitygradient (for a more thorough discussion of activated fluids see [7]). It is analogousto the Bingham constitutive equation for a viscoplastic fluid, but with the roles ofthe stress and symmetric velocity gradient interchanged; the fact that we can swapthe roles of the stress and the symmetric velocity gradient in constitutive relationswithout any problem is a significant advantage of the framework presented here.The problem was solved on the unit square Ω = (0 , with the rest state as theinitial condition and with the following boundary conditions: ∂ Ω = (0 , × { } , ∂ Ω := ∂ Ω \ ∂ Ω , u = on (0 , T ) × ∂ Ω , u = ( x (1 − x ) y , T on (0 , T ) × ∂ Ω . Although (4.4) has a complicated form, there is a continuous (in D ) selection -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS S = S ( x, y, D ) := ν (cid:16) | D | − δ s B / (1 / ( x, y ) (cid:17) + D | D | , if | D | (cid:54) = 0 , , if | D | = 0 . While the selection stated in (4.5) is already continuous in D , Newton’s methodrequires Fréchet-differentiability of S with respect to D and the constitutive law is notsmooth when | ( x − , y − ) | < ; therefore some regularisation was required for thepurpose of applying Newton’s method (an alternative would have been to use a non-smooth generalisation such as a semismooth Newton method). For this problem wechose a Papanastasiou-like regularisation (cf. [48]); the Papanastasiou regularisationhas been successfully applied to several problems with Bingham rheology [16, 24, 47].The regularised constitutive relation reads:(4.6) D = 12 ν (cid:18) δ s (1 − exp( − M | S | )) | S | + 1 (cid:19) S for ( x − ) + ( y − ) ≤ ( ) , where M > is the regularisation parameter (as M → ∞ we recover the constitutiverelation (4.4), see Figure 1); note that this is not related to the regularisation (2.7),which has the goal of turning the measurable selection into a continuous function.For the velocity and pressure we used Scott–Vogelius elements and discontinuouspiecewise polynomials were used for the stress (cf. subsection 4.1); the problem wasimplemented in firedrake with k = 1 , ν = , using the same parameters for thelinear and nonlinear solvers described in the previous section, and continuation wasemployed to reach the values M = 200 and δ s = 2 . ; more precisely, the problem wasinitially solved with M = 100 and δ s = 0 and that solution was used as the Newtonguess for the problem with M + 1 and δ s + 0 . , repeating the procedure until thedesired values were reached. The time step was chosen as τ m = 5 × − and thealgorithm was applied until the L norm of the difference of solutions at subsequenttime steps was less than × − . . . . . . . | D | . . . . . . . . | S | M = ∞ M = 120 M = 50 M = 15 M = 5 Fig. 1: Regularised constitutive relation for different values of M and δ s = 2 .Note that when the ‘yield strain’ parameter δ s vanishes, we recover the usualNavier–Stokes system. On the other end, if δ s is taken to be very large this could6 P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI be taken as an approximation of the incompressible Euler system in the center of thesquare; notice how in Figure 2 the fluid picks up more speed in the middle of thedomain when δ s > due to the absence of viscosity. This could be an attractiveapproach to simulating the effects of boundary layers, because it is backed up by arigorous convergence result; near the boundary the fluid could behave in a Newtonianway and far away δ s could be taken arbitrarily large so as to make the effects of theviscosity negligible. This is just one of the possibilities that are yet to be exploredwithin this framework of implicitly constituted fluids and mixed formulations and willbe studied in more depth in future work.Fig. 2: Streamlines of the steady state for the problem with δ s = 2 . (left) and theNewtonian problem (right).Figure 3 shows the magnitudes of S and D along the line x = 0 . for the steadystate of the non-Newtonian problem; it can be clearly seen that the stress is negligiblysmall for low values of the symmetric velocity gradient in the center of the square andit then suddenly becomes proportional to it. This transition is not the sharpest inthe figure because the regularisation parameter M was not taken sufficiently large,but in the limit this would recover the non-smooth relation. In a sense this is similarto solving a Navier–Stokes problem with high Reynolds number, so for high valuesof M some stabilisation would be required in order to solve this systems efficiently(even more so if the Newtonian fluid outside of the activation region also has a highReynolds number); this will be the subject of future research. The flow betweentwo parallel plates induced by the movement at constant speed of one of the platesreceives the name of (plane) Couette flow. It is one of the few examples of a configu-ration that allows us to find an exact solution for the steady Navier–Stokes equationsand it is well known that this solution has a linear profile. In this numerical ex-periment we will take the Couette flow as the initial condition and investigate thebehaviour of the system when the plates stop moving. Physically it is expected thatthe viscosity and no–slip boundary condition will slow down the flow until it finallystops; it can be seen in [49] that in the Newtonian case the flow does reach the reststate, albeit in infinite time. -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS S and D at x = 0 . for the problem with δ s = 2 . .In this section we will solve system (2.11) with the Bingham constitutive relation: S = τ y D | D | + 2 ν D , if | S | ≥ τ y , D = 0 , if | S | < τ y , where ν > is the viscosity and τ y ≥ is called the yield stress. This is the mostcommon model for a viscoplastic fluid, which is a material that for low stresses (i.e.with a magnitude below the yield stress τ y ) behaves like a solid and like a Newtonianfluid otherwise. Interestingly, viscoplastic fluids in the configuration described abovereach the rest state in a finite time and there are theoretical upper bounds for theso called cessation time (see [35, 42]), which makes this a good problem to test thenumerical algorithm. Just as in the previous section, for this problem there is also acontinuous selection available:(4.7) D = D ( S ) := ν ( | S | − τ y ) + S | S | , if | S | (cid:54) = 0 , , if | S | = 0 . For this experiment we again applied the Papanastasiou regularisation to the non-smooth constitutive relation, in order to be able to apply Newton’s method. Afternondimensionalisation this regularised constitutive law takes the form (compare with(4.6)):(4.8) S ( D ) = (cid:18) Bn | D | (1 − exp( − M | D | )) + 1 (cid:19) D , where Bn = τ y LνU is the Bingham number (here U and L are a characteristic velocityand length of the problem, respectively), and M > is the regularisation parameter(as M → ∞ we recover the non–smooth relation; compare with Figure 1). The prob-8 P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI lem was solved on the unit square
Ω = (0 , with the following boundary conditions: ∂ Ω = { } × (0 , ∪ { } × (0 , , ∂ Ω := (0 , × { } ∪ (0 , × { } , u = on (0 , T ) × ∂ Ω , u τ = 0 on (0 , T ) × ∂ Ω , − p + Sn · n = 0 , on (0 , T ) × ∂ Ω , where u τ denotes the component of the velocity tangent to the boundary and n is theunit vector normal to the boundary. The initial condition was taken as a standardCouette flow: u (0 , x ) = (1 − x , T . For the velocity and pressure we used Taylor–Hood elements and discontinuous piece-wise polynomials for the stress. This problem was implemented in FEniCS [45] usingthe same parameters for the nonlinear and linear solvers described in the previoussection, with k = 1 and a timestep τ m between × − and × − for the differ-ent values of the Bingham number. We quantify the change in the flow through thevolumetric flow rate (observe that it is constant in x ): Q ( t ) := (cid:90) (1 , · u ( t, x ) d x , whose evolution in time is shown in Figure 4 for different values of the Binghamnumber. An exponential decay of the flow rate is observed in Figure 4, while forpositive values of the Bingham number this decay is much faster; these results agreewith the ones reported in [42, 16]. In [16] the problem was solved by integrating aone-dimensional equation for u ; the framework presented here recovers the resultsobtained there but at the same time has the advantage that it can be applied to amuch broader class of problems and geometries. Q ( t ) Bn = 0.0Bn = 2.0Bn = 4.0
Fig. 4: Evolution of the volumetric flow rate. -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS
5. Conclusions.
In this work we presented a 3-field finite element formulationfor the numerical approximation of unsteady implicitly constituted incompressiblefluids and identified the necessary conditions that guarantee the convergence of thesequence of numerical approximations to a solution of the continuous problem. Al-though the convergence analysis was written in terms of a selection D , the finiteelement formulation presented here can be used in practice with a fully implicit rela-tion; this is in contrast to the works [21, 61], where the algorithms relied on findingan approximate constitutive law expressing the stress S k in terms of the symmetricvelocity gradient D k , which, while always theoretically possible, is not practical formany models. We also presented numerical experiments that showcase the variety ofmodels that the framework of implicitly constituted models can incorporate. The second author is grateful for useful discussionswith J. Málek, T. Tscherpel, and J. Blechta.
REFERENCES[1]
D. N. Arnold, F. Brezzi, and J. Douglas , PEERS: A new mixed finite element for planeelasticity , Japan J. Appl. Math., 1 (1984), pp. 347–367.[2]
D. N. Arnold and R. Winther , Mixed finite elements for elasticity , Numer. Math., 92(2002), pp. 401–419.[3]
S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman,L. Dalcin, V. Eijhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C.McInnes, K. Rupp, S. Smith, B. F. Zampini, H. Zhang, and H. Zhang , Petsc usersmanual
M. A. Behr, L. P. Franca, and T. E. Tezduyar , Stabilized finite element methods forthe velocity-pressure-stress formulation of incompressible flows , Comput. Methods Appl.Mech. Eng., 104 (1993), pp. 31–48.[5]
L. Belenki, L. Berselli, L. Diening, and M. Růžička , On the finite element approximationof p-Stokes systems , SIAM J. Numer. Anal., 50 (2012), pp. 373–397, https://doi.org/10.1137/10080436X.[6]
L. Berselli, L. Diening, and M. Růžička , Optimal error estimate for semi-implicit space-time discretization for the equations describing incompressible generalized Newtonian flu-ids , IMA J. Numer. Anal., 35 (2015), pp. 680–697.[7]
J. Blechta, J. Málek, and K. R. Rajagopal , On the classification of incompressible fluidsand a mathematical analysis of the equations that govern their motion , ArXiv Preprint:1902.04853, (2019), https://arxiv.org/abs/1902.04853.[8]
D. Boffi, F. Brezzi, and M. Fortin , Mixed Finite Element Methods and Applications ,Springer, 2013.[9]
D. Breit, L. Diening, and S. Schwarzacher , Solenoidal Lipschitz truncation for parabolicPDEs , Math. Models Methods Appl. Sci., 23 (2013), pp. 2671–2700.[10]
F. Brezzi and M. Fortin , Mixed and Hybrid Finite Element Methods , Springer Ser. Comput.Math., 1991.[11]
P. R. Brune, B. F. Knepley, B. F. Smith, and X. Tu , Composing scalable nonlinearalgebraic solvers , SIAM Rev., 57 (2015), pp. 535–565.[12]
M. Bulíček, P. Gwiazda, J. Malek, K. R. Rajagopal, and A. Świerczewska-Gwiazda , On flows of fluids described by an implicit constitutive equation characterized by a maximalmonotone graph , vol. 402 of London Math. Soc. Lecture Note Ser, Cambridge Univ. Press:Cambridge, 2012.[13]
M. Bulíček, P. Gwiazda, J. Málek, and A. Świerczewska-Gwiazda , On steady flowsof incompressible fluids with implicit power-law-like rheology , Adv. Calc. Var., 2 (2009),pp. 109–136.[14]
M. Bulíček, P. Gwiazda, J. Málek, and A. Świerczewska-Gwiazda , On unsteady flowsof implicitly constituted incompressible fluids , SIAM J. Math. Anal., 44 (2012), pp. 2756–2801, https://doi.org/10.1137/110830289.[15]
M. Bulíček, J. Málek, V. Průša, and E. Süli , PDE analysis of a class of thermodynami-cally compatible viscoelastic rate-type fluids with stress-diffusion , in Mathematical Analysisin Fluid Mechanics: Selected Recent Results, vol. 710 of AMS Contemporary Mathematics, P. E. FARRELL, P. A. GAZCA-OROZCO, AND E. SÜLI2018, pp. 25–51.[16]
M. Chatzimina, G. C. Georgiou, I. Argyropaidas, E. Mitsoulis, and R. R. Huilgol , Cessation of Couette and Poiseuille flows of a Bingham plastic and finite stopping times ,J. Non-Newtonian Fluid Mech., 129 (2005), pp. 117–127.[17]
P. Clément , Approximation by finite element functions using local regularization , RAIRO,Anal. Numér., R2, (1975), pp. 77–84.[18]
M. Crouzeix and P. A. Raviart , Conforming and nonconforming finite element methodsfor solving the stationary Stokes equations I , ESAIM: M2AN, (1973), pp. 33–75.[19]
T. A. Davis , Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method ,ACM Trans. Math. Softw., 30 (2004), pp. 196–199.[20]
E. DiBenedetto , Degenerate Parabolic Equations , Springer Ver., 1993.[21]
L. Diening, D. Kreuzer, and E. Süli , Finite element approximation of steady flows ofincompressible fluids with implicit power-law-like rheology , SIAM J. Numer. Anal., 51(2013), pp. 984–1015, https://doi.org/10.1137/120873133.[22]
L. Diening, M. Ružička, and J. Wolf , Existence of weak solutions for unsteady motions ofgeneralized Newtonian fluids , Ann. Scuola Norm. Sup. Pisa Cl. Sci., IX (2010), pp. 1–46.[23]
L. Diening, S. Schwarzacher, V. Stroffolini, and A. Verde , Parabolic Lipschitz trun-cation and caloric approximation , Calc. Var., 56 (2017).[24]
Y. Dimakopoulos and J. Tsamopoulos , Transient displacement of a viscoplastic material byair in straight and suddenly constricted tubes , J. Non-Newtonian Fluid Mech., 112 (2003),pp. 43–75.[25]
S. Eckstein and M. Růžička , On the full space-time discretization of the generalized Stokesequations: The Dirichlet case , SIAM J. Numer. Anal., 56 (2018), pp. 2234–2261, https://doi.org/10.1137/16M1099741.[26]
V. J. Ervin, J. S. Howell, and I. Stanculescu , A dual-mixed approximation method for athree-field model of a nonlinear generalized Stokes problem , Comput. Methods Appl. Mech.Eng., 197 (2008), pp. 2886–2900.[27]
M. Farhloul and M. Fortin , New mixed finite element for the Stokes and elasticity problems ,SIAM J. Numer. Anal., 30 (1993), pp. 971–990, https://doi.org/10.1137/0730051.[28]
M. Farhloul and H. Manouzi , Analysis of non-singular solutions of a mixed Navier-Stokesformulation , Comput. Methods Appl. Mech. Eng., 129 (1996), pp. 115–131.[29]
M. Farhloul, S. Nicaise, and L. Paquet , A refined mixed finite-element method for thestationary Navier-Stokes equations with mixed boundary conditions , IMA J. Numer. Anal.,28 (2008), pp. 25–45.[30]
M. Farhloul, S. Nicaise, and L. Paquet , A priori and a posteriori error estimations forthe dual mixed finite element method of the Navier-Stokes problem , Numer. Meth. Part.Differ. Equat., 25 (2009), pp. 843–869.[31]
G. P. Galdi , An Introduction to the Mathematical Theory of the Navier-Stokes Equations:Steady State Problems , Springer, Second Edition ed., 2011.[32]
C. Geuzaine and J. F. Remacle , Gmsh: a three-dimensional finite element mesh generatorwith built-in pre- and post-processing facilities , Int. J. Numer. Meth. Eng., 79 (2009),pp. 1309–1331.[33]
V. Girault and P. A. Raviart , Finite Element Methods for Navier-Stokes Equations: The-ory and Algorithms , Springer Verlag, 1986.[34]
V. Girault and L. R. Scott , A quasi-local interpolation operator preserving the discretedivergence , Calcolo, (2003), pp. 1–19.[35]
R. Glowinski , Numerical Methods for Nonlinear Variational Problems , Springer–Verlag, NewYork, 1984.[36]
J. Guzmán and M. Neilan , Conforming and divergence-free Stokes elements in three dimen-sions , IMA J. Numer. Anal., (2014), pp. 1489–1508.[37]
J. Guzmán and M. Neilan , Conforming and divergence-free Stokes elements on generaltriangular meshes , Math. Comput., 83 (2014), pp. 15–36.[38]
A. Hirn , Approximation of the p-Stokes equations with equal-order finite elements , J. Math.Fluid Mech, 15 (2013), pp. 65–88.[39]
J. S. Howell , Dual-mixed finite element approximation of Stokes and nonlinear Stokes prob-lems using trace-free velocity gradients , J. Comput. Appl. Math, 231 (2009), pp. 780–792.[40]
J. S. Howell and N. J. Walkington , Dual-mixed finite element methods for the Navier-Stokes equations , ESAIM: M2AN, 47 (2013), pp. 789–805.[41]
J. Hron, J. Málek, J. Stebel, and K. Touška , A novel view on computations of steadyflows of Bingham fluids using implicit constitutive relations , Project MORE Preprint,(2017).[42]
R. R. Huilgol, B. Mena, and J. M. Piau , Finite stopping time problems and rheometry of -FIELD FEM FOR UNSTEADY IMPLICITLY CONSTITUTED FLUIDS Bingham fluids , J. Non-Newtonian Fluid Mech., 102 (2002), pp. 97–107.[43]
C. Kreuzer and E. Süli , Adaptive finite element approximation of steady flows of incom-pressible fluids with implicit power-law-like rheology , ESAIM: M2AN, 50 (2016), pp. 1333–1369.[44]
J. L. Lions , Quelques Méthodes De Résolution Des Problèmes Aux Limites Non Linéaires ,Dunod, Paris, 1969.[45]
A. Logg, K. A. Mardal, and G. N. Wells , FEniCS : Automated Solution of Differ-ential Equations by the Finite Element Method , Springer, 2011, https://doi.org/10.1007/978-3-642-23099-8.[46]
E. Maringová and J. Žabenský , On a Navier–Stokes–Fourier-like system capturing tran-sitions between viscous and inviscid fluid regimes and between no-slip and perfect-slipboundary conditions , Nonlinear Anal. Real World Appl., 41 (2018), pp. 152–178.[47]
E. Mitsoulis and R. R. Huilgol , Entry flows of Bingham plastics in expansions , J. Non-Newtonian Fluid Mech., 122 (2004), pp. 45–54.[48]
T. C. Papanastasiou , Flows of materials with yield , J. Rheol, 31 (1987), pp. 385–404.[49]
T. C. Papanastasiou, G. Georgiou, and A. Alexandrou , Viscous Fluid Flow , CRC Press,Boca Raton, 1999.[50]
V. Průša and K. R. Rajagopal , A new class of models to describe the response of electrorhe-ological and other field dependent fluids , Adv. Struct. Mater., 89 (2018), pp. 655–673.[51]
K. R. Rajagopal , On implicit constitutive theories , Appl. Math., 48 (2003), pp. 279–319.[52]
K. R. Rajagopal , On implicit constitutive theories for fluids , J. Fluid Mech., 550 (2006),pp. 243–249.[53]
K. R. Rajagopal , The elasticity of elasticity , Z. Angew. Math. Phys, 5807 (2007), pp. 309–317.[54]
K. R. Rajagopal and A. R. Srinivasa , On the thermodynamics of fluids defined by implicitconstitutive relations , Z. Angew. Math. Phys, 59 (2008), pp. 715–729.[55]
F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. Mcrae,G.-T. Bercea, G. R. Markall, and P. H. J. Kelly , Firedrake: automating the finiteelement method by composing abstractions , ACM Trans. Math. Softw., 43 (2016).[56]
T. Roubiček , Nonlinear Partial Differential Equations with Applications , Birkhäuser, secondedition ed., 2013.[57]
V. Ruas , An optimal three-field finite element approximation of the Stokes system with con-tinuous extra stresses , Japan J. Appl. Math., (1994), pp. 113–130.[58]
D. Sandri , A posteriori estimators for mixed finite element approximations of a fluid obeyingthe power law , Comput. Methods Appl. Mech. Eng., (1998), pp. 329–340.[59]
L. R. Scott and M. Vogelius , Norm estimates for a maximal right inverse of the divergenceoperator in spaces of piecewise polynomials , Math. Modelling Numer. Anal., 19 (1985),pp. 111–143.[60]
J. Simon , Compact sets in the space L p (0 , T ; B ) , Ann. Mat. Pura Appl., 4 (1987), pp. 65–96.[61] E. Süli and T. Tscherpel , Fully discrete finite element approximation of unsteady flows ofimplicitly constituted incompressible fluids , IMA J. Numer. Anal., (2018), pp. 1–49.[62]