FINITE ELEMENT APPROXIMATION OF DIELECTROPHORETIC FORCE DRIVEN FLOW PROBLEMS

. In this paper, we propose a full discretization scheme for the instationary thermal-electro-hydrodynamic (TEHD) Boussinesq equations. These equations model the dynamics of a non-isothermal, dielectric fluid under the influence of a dielectrophoretic (DEP) force. Our scheme combines an H 1 - conformal finite element method for spatial discretization with a backward differentiation formula (BDF) for time stepping. The resulting scheme allows for a decoupled solution of the individual parts of this multi-physics system. Moreover, we derive a priori convergence rates that are of first and second order in time, depending on how the individual ingredients of the BDF scheme are chosen and of optimal order in space. In doing so, special care is taken of modeling the DEP force, since its original form is a cubic term. The obtained error estimates are verified by numerical experiments.


Introduction
If a non-isothermal, dielectric fluid is exposed to an alternating, external electric field E, then it experiences the so-called dielectrophoretic (DEP) force. The resulting state of the fluid, given by its velocity u, pressure , temperature and internal electric potential Φ, can be modeled by means of the thermal-electro-hydrodynamic (TEHD) Boussinesq equations in the domain Ω [28], u + (u · ∇)u − ∆u + ∇ = |∇Φ| with fluid density , = ( ) and dielectric constant 0 [28]. In recent years, there has been growing interest in the investigation of DEP-driven flows, since the DEP force can be used to enhance the effectiveness of heat exchanging systems [7]. Moreover, the DEP force can be considered as artificial gravity and thereby allows to simulate Earth mantle convection in physical experiments [8]. While many works on this topic are of experimental nature [10,23,25,27,35], or address linear stability analysis [22,24,37,38,40,45], there are also few contributions on the numerical solution of (1.1)- (1.5). For instance, in [19][20][21]46] the use of finite volumes in combination with low order time stepping schemes is proposed. A spectral approach is presented in [42][43][44] and the finite element method (FEM) is used in [26,35,36].
However, to the best of the authors' knowledge, there is a lack of publications that address the numerical analysis of any kind of discretization method applied to the TEHD Boussinesq equations. We thus propose a full discretization scheme, based on H 1 -conformal FEM and variants of the backward differentiation formula (BDF), and derive corresponding a priori error estimates in this work. The presented analysis is inspired by Tabata and Tagami [39], where error estimates for a FEM/backward Euler discretization of the standard Boussinesq equations with temperature dependent viscosity and thermal diffusion coefficient are derived. Further results on numerical analysis of FEM applied to the standard Boussinesq equations are derived in many works, see for instance [1,2,4,5,[29][30][31]33] and the references therein.
Compared to the natural convection case, the main difficulty that comes with system (1.1)-(1.5) is due to the DEP force term. For the typical W 1,2 -regularity of solutions and Φ of heat equation and Gauss' lawinvolving non-smooth coefficients, Lipschitz domains and mixed boundary conditions -one cannot assume that |∇Φ| 2 ∇ is sufficiently regular to be an element of H −1 (Ω). Therefore, we propose an approximate DEP force term, based on a cut-off operator, to obtain a well-defined variational formulation. The error estimates are derived for this modified DEP term and convergence rates of the form ℎ + with mesh width ℎ and time step are obtained. As usual, the spatial order of convergence depends on the assumed regularity of the exact solution and the approximation properties of the underlying finite element spaces. Concerning the temporal order of velocity, temperature and potential error, we obtain = 1 or = 2 for the one-step and two-step BDF scheme, respectively.
All presented error estimates only hold under a small data condition which originates from the DEP force and which is of the form max with domain dependent constant Ω . Compared to small data conditions in other situations, see the ones for the stationary Boussinesq equations for instance [31], equation (1.6) is rather mild though and its validness can be shown in realistic scenarios, see Section 7 and [12]. The derived convergence rates are validated by numerical experiments and we further discuss in which sense the obtained results also apply to the original DEP term |∇Φ| 2 ∇ .
The proposed time stepping schemes treat all nonlinear terms in a semi-implicit fashion, which allows for solving the incompressible flow equations, heat equation and Gauss' law in a decoupled way. Moreover, the approximation of and the extrapolation of u in the arising convection terms are defined in a general manner. Thus, we are able to cover the aforementioned one-and two-step BDF schemes in a unified analysis.
The outline of this paper is as follows: First, the continuous variational formulation is stated in Section 2. In Section 3 we present the spatial and temporal discretization scheme and show well-posedness of the discrete equations. Corresponding a priori error estimates are derived in Section 4. In Section 5 we address the modeling of the DEP force. Numerical experiments underlining the theoretical convergence rates are presented in Section 6. In Section 7 we discuss the relation between the original and modified DEP term. We conclude this work in Section 8.

Continuous problem
In this section, we state a continuous variational formulation for the TEHD Boussinesq equations (1.1)-(1.5). As previously noted, a variational formulation of (1.1) which involves a DEP force of the form |∇Φ| 2 ∇ would imply some problematic consequences. Using W 1,2 (Ω) ˓→ L 6 (Ω), the regularity conditions Φ ∈ W 1, (Ω), ∈ W 1, (Ω) with 1 = 1 + 2 + 1 6 (which is only possible to hold for ≥ 12 5 ) are required to infer, by means of Hölder's inequality, that ⟨F dep ( , Φ), v⟩ H −1 := (|∇Φ| 2 ∇ , v) is well-defined. Similarly, a Lipschitz property of F dep (needed in the proposed error analysis) relies on the norms ‖·‖ 1, and ‖·‖ 1, . However, the natural norm in this analysis for bounding , Φ is ‖·‖ 1,2 and the conversion between these norms, e.g. by means of inverse estimates, would introduce negative powers of ℎ, which could not be compensated. For this reason, we replace F dep by a general body force F for now, which is supposed to model the DEP term F dep , and also incorporates additional terms such as gravitation and other source terms. Below, we impose conditions on F that are sufficient to perform the error analysis, see Assumption 3.13, and propose a suitable realization for F as an approximation to the original DEP force F dep and possible additional terms, see Section 5. At this point, we only require the following assumption.
Remark 2.6. f v , , denote arbitrary source terms, that are independent of the unknown solution and boundary conditions. The strong formulation (1.1)-(1.5) corresponds to Problem 2.4 for f v = = = 0.
Remark 2.7. If is additionally in 2 and the boundary lifting is assumed to be in 2 as well, then Green's formula can be applied in the usual way on the term ( + , ); first with ∈ C ∞ (Ω) to show that + solves (1.3); second with ∈ ∞ (Ω) to obtain ∇( + ) · n = 0 on Γ . Analogous results hold for Φ, Φ .
Remark 2.8. In [12], existence and stability of solutions of the instationary TEHD equations in solenoidal form is shown. These solutions, however, are of lower temporal regularity, u ∈ L 1 (0, ; V * ), ∈ L 1 (0, ; Θ * ), and the instationary equations only hold in a distributional sense, i.e. in C ∞ ((0, )) * . Moreover, the underlying assumptions on F in [12] differ from those imposed in Section 3. For this existence result it is assumed that weak convergence of sequences ( , Φ ) ⇀ ( , Φ) in H 1 ×H 1 should imply pointwise convergence F( , Φ ) → F( , Φ) in U * , which is not satisfied by the DEP model proposed for the numerical analysis in this work. Corresponding modelizations of F are based on linearization w.r.t. Φ around a given reference potential. To the best of the authors' knowledge, this is the only work on well-posedness of the TEHD equations.

Discrete problem
In this section, we propose a full discretization of Problem 2.4, leading to the discrete Problem 3.12. Afterwards, unique existence and stability of discrete solutions is shown.
In the upcoming error analysis, we make use of the well-known Stokes projection, e.g. [15], to bound the approximation error of velocity and pressure. In this way, some terms of the instationary error equation naturally cancel out.
The Stokes Projection is defined as Moreover, there is , independent of ℎ and , such that for all (u, ) If the Stokes problem is regular, i.e. the mapping (u, ) ↦ → − ∆u + ∇ is an isomorphism from Proof. See Lemma 4.42 and 4.43 in [18] for the stability result and ℎ convergence, and Theorem II.1.9 in [14] for the improved L 2 convergence rates.
Remark 3.5. According to Theorem I.5.4 and Remark I.5.6 in [14], the Stokes problem is regular if Γ ∈ 2 or if = 2 and Ω is a convex polygon.
The full discretization scheme proposed in the next section will lead to a sequence of stationary, spatially discretized problems of the following form. Problem 3.6 (Linearized stationary TEHD equations). Let > 0, source terms f v ∈ U * , ∈ Θ * , ∈ Υ * and functionsū ∈ U,¯,˜∈ Θ be given.

Full discretization
To set up the temporal discretization, let 0 = 0 < 1 < . . . < = denote an equidistant partition of the time interval [0, ], i.e. = with := and ∈ N. The approximation of and the extrapolation in the nonlinear terms are described by multi-step difference operators. Definition 3.9 ( -step difference operator). Let denote a vector space, ∈ N and c ∈ R +1 . The -step difference operator c , applied to a sequence ( ) ⊂ , is defined as c is called explicit, if 0 = 0, and 1 denotes the 1-step operator determined by c = (1, −1).
Assumption 3.10 (Difference operators). Let -step difference operators c , d , f , g be given with (iii) There exists a constant¯≥ 0 such that for each sequence ( ) ≥0 ⊂ with ∈ {L 2 (Ω), L 2 (Ω)}, there are sequences for all ≥ and only depends on for ≤ .
The following lemma, see e.g. [9,39], gives two possible choices, for which condition (iii) of Assumption 3.10 is satisfied.
(ii) If d = ( 3 2 , −2, 1 2 ), then Assumption 3.10 (iii) holds with The full discretization of Problem 2.4 is now given by the following sequence of stationary problems.

Problem 3.12 (Discretized instationary TEHD equations). Let ≥ 1 and initial conditions
Assume that the -step difference operators c , d , f , g satisfy Assumption 3.10 and that the source terms are given by f := f ( ), where f is defined as in Problem 2.4 for ∈ {v, , }.
Due to Assumption 3.10 (i), the nonlinear terms are treated semi-implicitly. In this way, the resulting system for each time step is of the same form as Problem 3.6 and can be solved by solving a series of linear systems. Moreover, the corresponding proof of a priori error estimates does not require a time step restriction when applying Grönwall's inequality, see Lemma A.1. Assumption 3.10 (ii) is a consistency condition, that implies f = g = for the (not time-dependent) boundary liftings ∈ { , Φ }. The stability condition (iii) is used to obtain the L 2 -norm of the discretization error in the proof of Theorem 4.1. It remains to specify requirements on the general body force F, that are sufficient for proving a priori error estimates. In particular, F should satisfy a Lipschitz condition w.r.t. both arguments. Concerning the second argument Φ, local Lipschitz continuity is sufficient, since an ∞ (H 1 ) bound on the discrete potential can be derived. For the same reason, the Lipschitz moduli w.r.t. the first argument may depend on ‖Φ‖ 1,2 . (i) There are non-decreasing functions Well-posedness of the discrete Problem 3.12 is provided by Lemma 3.14.
Lemma 3.14 (Existence, uniqueness and stability of discrete solutions). Let Assumptions 2.1-3.13 hold and assume that Here, 0 , denote the constants of Friedrich's inequality w.r.t. H 1 0 (Ω) and H 1 (Ω), respectively. Further, f is determined by f . Moreover, suppose that the starting values are chosen in such a way, that their respective norms can be bounded uniformly in ℎ:
Remark 3.17. The need for the small data condition (3.11) imposed by Lemma 3.14 stems from the occurence of ‖∇ ‖ when bounding the U * -norm of the body force term F, see (3.21). For this reason, the right-hand side in Assumption 3.13 (i) is split into H 1 and L 2 terms, where only the first one contributes to (3.11).

Error estimate
The main result of this work, Theorem 4.1, provides an a priori error estimate for the numerical scheme defined by Problem 3.12. The derived upper bounds depend on spatial approximation errors, which can be bounded in terms of ℎ, errors induced by the initial conditions, and approximation properties of the temporal difference operators. At this point, no regularity on temporal derivatives of the exact solution is imposed and, consequently, no convergence rates w.r.t. are derived. This is done as a separate step in form of Corollaries 4.7 and 4.8 in Section 4.2.

Main a priori estimate
Before stating Theorem 4.1, some notation has to be introduced. The exact solution at time is denoted by (u , , , Φ ) and corresponding approximations in the underlying finite element spaces are obtained by using the Stokes projection Π and projections provided by Assumption 3.1, The discretization errors between these projections and the discrete solution of Problem 3.12 are defined as When applying the estimate stated in Assumption 3.10 (iii) for = d u,ℎ and = ,ℎ , respectively, one further obtains real sequences ( ) ≥ −1 , ( ) ≥ −1 , ( ) ≥ , ( ) ≥ such that the following estimates hold for all ≥ , The error will be bounded in terms of the subsequent expressions. The first set of expressions measures the error contributions by the initial conditions: Error contributions that are purely driven by the temporal discretization are given by In Section 4.2, ( ) is bounded in terms of . Contributions by both temporal and spatial discretization are denoted by with constants , * that depend on the spatial regularity of the exact solution and the approximation quality of the finite element spaces. Finally, the term which involves the unknown velocity discretization error will be bounded by other terms and is used to bound the pressure discretization error. The actual error bounds are now composed by the previously defined terms: The following theorem states that the errors of velocity and temperature, contributed by the spatial discretization and measured in L 2 (H 1 ), ∞ (L 2 )-norm, converge with rate ℎ for fixed time step size . The analogous result in L 2 (H 1 ) holds for the potential error. Concerning the pressure error, the factor − 1 2 enters the corresponding upper bound ℰ through the term d ( , ℎ), thus yielding a reduced temporal convergence rate compared to velocity, temperature and potential. Here, the term denotes a bound on the ∞ (H 1 )-norm of discrete velocity (4.57) and one can show that it stays bounded for ℎ, → 0, supposed that the exact solution is sufficiently regular and the initial conditions are chosen appropriately, see Corollaries 4.7, 4.8 and the subsequent remarks.
Theorem 4.1 (Error of numerical scheme). In the framework of Problem 3.12, let a solution of Problem 2.4 be given that satisfies the following additional regularity requirements for some ≥ 1: Moreover, the exact solution should be pointwise well-defined for all ∈ [0, ].
Suppose that Assumption 3.1 holds with¯=¯= . Moreover, let * = 1 if the underlying Stokes problem is regular and * = 0, otherwise. Let further the assumptions of Lemma 3.14 hold; in particular, assume that F satisfies Assumption 3.13 and the small data conditions denoting the respective constants from Friedrich's inequality w.r.t.
with Π given by Assumption 3.1 and Φ given by Lemma 3.14. Then, The underlying principle of the corresponding proof is inspired by comparable results on the error of fully discretized flow problems, see e.g. [39] for the Boussinesq problem with temperature-dependent coefficients.
Proof. The error of the respective variables is decomposed into approximation and discretization contributions:

Error equation
Subtracting the discrete equations, Problem 3.12, from the continuous ones, Problem 2.4, at each time , with residuals By definition of (w ℎ , ℎ ) and the Stokes projection, there holds ⟨ 2 , v ℎ ⟩ = 0 and ⟨ 5 , ℎ ⟩ = 0 and for residual In the following, the introduced constants { , } do not depend on the discrete solution, the discretization parameters ℎ, and the diffusion constants , . However, they may include norms of the exact solution.

Finite element approximation error
The finite element approximation errors can be bounded by means of Assumption 3.1, Lemma 3.4, and the imposed regularity of the exact solution: By using the linearity of Π , Π Θ , Π ϒ , Assumption 3.1 and Lemma 3.4, one obtains for an arbitrary -step difference operator h , Final error estimation for velocity, temperature and potential Using (4.44)-(4.49) with h = d and the discrete stability result, Lemma 3.14, the following estimate follows: (4.51) Using (4.50), the right-hand side in (4.43) is equal to ℰ u, (ℎ, , ) 2 , up to a multiplicative constant that is independent of ℎ and . In this way, the estimates (4.55) follow. Similarly, the square root of the right-hand side in (4.51) can be estimated from above by ℰ Φ (ℎ, , ) (under the use of (4.52)). Thus, is obtained. The stated assertions on the complete error of velocity, temperature and potential now follow by the decomposition ·,ℎ = ·,ℎ + ·,ℎ , (4.52), (4.53), (4.56), estimation (4.44)-(4.47) for ·,ℎ and the triangle inequality.

Bound on temporal variation of estimation error
It remains to derive an upper bound for d (ℎ, ). Setting v ℎ = d d u,ℎ in (4.17), using (w ℎ , v ℎ ) = (u ℎ , v ℎ ) = 0 for all v ℎ , ≥ 0, and using Young's inequality yields Summing this inequality from = , . . . , and using the previously derived bounds on , (4.59)-(4.61), the approximation estimates (4.44)-(4.49) and the error estimates for velocity and temperature, (4.52), (4.53), yields Remark 4.2. We required that 2 4, < 2˜˜|f | −2 2 −1 . However, this assumption could be relaxed to 2 4, < 4˜˜|f | −2 2 −1 , to the expense of a slightly more technical presentation. (1) In particular, they depend on the typically small physical constants , in the same way and on the diameter of the domain by virtue of Friedrich's inequality. The DEP model presented in the next section is globally Lipschitz continuous w.r.t. with a constant independent of Φ. In this case, F . Remark 4.4. A time step size restriction when applying Grönwall's inequality could be avoided, since the right-hand side in (4.42) only contains terms ‖d u,ℎ ‖, ‖ ,ℎ ‖ for < . For this reason, it is assumed that c , f , g are explicit difference operators.
One drawback of the presented proof lies in the exponential factor that arises due to the application of Grönwall's inequality. This exponent is proportional to negative powers of viscosity and thermal diffusion and may thus take very large values for typical fluids. Strictly speaking, one obtains an exponent exp with The term −3 + −2 −1 comes from the estimation of the convection terms˜v,˜, see (4.40), (4.31) and (4.32). In [34], an approach is shown that avoids the explicit occurrence of in the exponential factor when deriving error estimates for the instationary incompressible Navier-Stokes equations. This is achieved by the use of H( )conforming elements that yield exactly divergence free discrete velocities in combination with ∇u ∈ L 1 (0, ; ∞ ) and a W 1,∞ -stability assumption on the Stokes projection. Conditions for this stability are given in [15], for instance. The remaining factor −1 arises due to the body force and is multiplied by the constant Thus, large values of −1 could be compensated if the domain has a small diameter, i.e. small 0 , and if the spatial variations of temperature and potential are rather low.

Temporal convergence rates
We now consider two concrete realizations of the temporal discretization by setting the difference operators introduced in Problem 3.12. The resulting schemes are given by a 1-step BDF (similar to semi-implicit Euler) and a 2-step BDF method.  , c = f = g := (0, 2, −1). Both schemes are widely used for discretizing viscous flow problems, see e.g. [39] for the use of BDF1 for the Boussinesq equations and [9,13] for the use of BDF2 for the incompressible Navier-Stokes equations.
It is easy to verify that difference operators given by Definition 4.5 and 4.6 satisfy Assumption 3.10 (i) and (ii). The stability condition (iii) is stated by Lemma 3.11. For both schemes, the -dependent error contributions in Theorem 4.1, ( ), (ℎ, ) and d (ℎ, ), are bounded from above in order to derive convergence rates w.r.t. by employing Taylor series, see Corollaries 4.7 and 4.8. In doing so, one has to impose temporal regularity conditions on the solution which imply non-local compatibility conditions on the initial data which may be uncheckable in practice [16]. Among others, these critical conditions include the finiteness of the following quantities [16]: ‖u‖ ∞;H 3 ‖ u‖ ∞;H 1 , ‖ u‖ 2;H 2 , ‖ u‖ 2;L 2 . In this work, however, we do not address this issue and impose temporal regularity that is sufficient to obtain optimal convergence rates. We refer to [12] for a more detailled consideration on this topic, where reduced convergence rates for weaker regularity conditions are derived.
Corollary 4.7 (Convergence rates of BDF1 scheme). Let the assumptions of Theorem 4.1 hold and suppose that the exact solution satisfies the following, additional regularity conditions: Let the difference operators in Problem 3.12 be given according to Definition 4.5 and ℎ, < 1. Then, Proof. Application of Lemma A.2 (i), (iii), with = 0, yields Moreover, by Lemma A.2 (i) with = 0, Thus, Further, by Lemma 3.11 (i) and (4.55), Let the difference operators in Problem 3.12 be given according to Definition 4.6 and ℎ, < 1. Then, Proof. Application of Lemma A.2 (ii), (iv) with = 0 yields and by Lemma A.2 (v) with = 0, = 1, Plugging these estimates into the definitions of the indiviudal error terms yields Remark 4.9. For the BDF1 scheme, there holds = 0 according to Lemma 3.11. Thus, the terms −1 , −1 in the error bound 1 ( ) vanish. In the BDF2 case, we have = ⃦ ⃦ 2 − −1 ⃦ ⃦ and therefore with an analogous estimate holding for | −1 |. In summary, if the initial conditions (u ℎ , ℎ ) =0 , are chosen appropriately, for instance by applying the projection operators of Assumption 3.1 to the continuous initial conditions (u 0 , 0 ) and by computing (u ℎ , ℎ ) =1 by means of BDF1 with time step size 2 , then the corresponding error contributions ( ) are of order (ℎ ) (BDF1) and (ℎ + 2 ) (BDF2). Corollaries 4.7, 4.8 and Remarks 4.9, 4.10 show that the proposed BDF1 and BDF2 methods are indeed of first and second order w.r.t. the temporal discretization. To be precise, the error of velocity, temperature and potential, as defined by Theorem 4.1, is of order (ℎ + ) with = 1 for BDF1 and = 2 for BDF2. Regarding the pressure error, a factor 1 2 is lost, similar to the estimates derived in [39].

Modeling of DEP force
We now propose a modification of the standard DEP and buoyancy force term that satisfies Assumption 3.13. The critical term F dep = |∇Φ| 2 ∇ contains the squared magnitude of the electric field |∇Φ| 2 which is only in L 1 (Ω), since the natural space for the solution Φ of Gauss' law is H 1 (Ω). In order to increase the integrability of this term, we thus replace it by |m ∘ ∇Φ| 2 , where m denotes a suitable cut-off function. In this way, ∞ (Ω) integrability is obtained and Assumption 3.13 is satisfied, as shown in the subsequent Definition 5.1 and Lemma 5.2.
An example of a suitable cut-off function m is given by the metric projection of a vector w ∈ R onto the ball {x ∈ R : |x| ≤ }.

Definition 5.3 (Metric projection). For
> 0 the metric projection is defined by The following Lemma 5.4 states that the previously defined cut-off function satisfies the required properties as stated in Definition 5.1. Proof. The stated assertion ‖m ‖ 0,∞ = is clear and the Lipschitz continuity follows from the fact, that m can be written as m (w) = w 0 where w 0 is the unique solution of the minimization min r∈ |w − r|, with the compact and convex set = {x ∈ R : |x| ≤ }, see e.g. [3].
Concerning the choice of , note that Lemma 3.14 and Theorem 4.1 require a small data condition of the form

Numerical experiments
In this section we consider a 2D test case to underline the theoretically derived error estimates of Section 4. > 0, respectively. All physical parameters are listed in Figure 1. We further assume that the permittivity depends linearly on the temperature [28] and set for¯= 10 2 , ∈ (0,¯− 1 ) to meet the requirements of Assumption 2. Thus, for = per the exact solution satisfies the regularity conditions for both, BDF1 and BDF2, assumed by Corollaries 4.7 and 4.8, respectively. On the other hand, for = alg the regularity conditions for BDF1 are still satisfied, whereas the conditions for BDF2 do not hold. The initial condition is set to the exact solution at = 0 and for BDF2, the approximate solution for the first time step is set to the exact solution at = . Thus, errors induced by the initial conditions, given by ( ), can be neglected.
The spatial discretization is based on a mesh consisting of quadrilaterals with finite element spaces given by Taylor-Hood elements for velocity and pressure (Q 2 /Q 1 ), and continuous quadratic Lagrange elements for temperature and potential (Q 2 ). The arising linear systems are solved by a preconditioned GMRES method [32]. We use a block Jacobi method as preconditioner whereat the inverse of each diagonal block is approximated by an incomplete LU factorization. The implementation is based on the open-source, general purpose FEM package HiFlow 3 [11]. Concerning the modified DEP force F , the cut-off threshold should be small enough such that the requirements imposed by Lemma 3.14 and Theorem 4.1 are satisfied. As pointed out in Section 5, these conditions are of the form with domain dependent constant As 0 , denote the respective constants in Friedrich's inequalities for H 1 0 , H 1 on a square domain, they can be bounded as 0 ≤ 1 and ≤ 1, see [12] for instance. Concerning the time stepping contributions, we have In summary, equations (6.5), (6.6) and = 0 2 imply that ≤ * with is sufficient for (6.4) to hold. On the other hand, the exact potential solution is chosen such that Using max ∈[0,1] | (2) ( )| ≤ 2, there holds Φ := max ∈[0,1] ‖∇Φ * ( )‖ 0,∞ ≤ 10 . Setting ≥ 10 thus implies that the modified DEP force reduces to the original version for the exact potential. Therefore, as long as ≤ 5.5 × 10 5 , it is possible to define such that both F ( * , Φ * ) = F ( * , Φ * ) and < * do hold. For the first test series we set = = 10 3 ≪ * , i.e. the small data condition is satisfied and F ( * , Φ * ) ̸ = F ( * , Φ * ). Tables 1 (BDF1) and 2 (BDF2) list the corresponding errors and convergence rates  Table 1. BDF1, = per , = = 10 3 , ℎ = 1.25 .
Error for decreasing and ℎ ∝ in case of full temporal regularity of the exact solution. Since ℎ ∝ , the ℎ terms in the error estimates given by Corollaries 4.7, 4.8 are transformed to . We obtain = 2 due to the spatial regularity of the exact solution and due to the polynomial order of the underlying finite element space. Thus, the temporal contributions in the derived error estimates dominate the spatial contributions ℎ ∝ and the observed errors should resemble the temporal convergence order. This is indeed the case for velocity, temperature and potential whose corresponding error rates are inline with Corollaries 4.7, 4.8, see Tables 1 and 2.  Tables 3 and 4 illustrate the errors for = alg . As pointed out by (6.3), the exact solution does not satisfy the regularity requirements of Corollary 4.8 any more, whereas the weaker assumptions supposed by Corollary 4.7 do still hold. In practice, one can observe = 1 for the BDF1 scheme, in accordance with Corollary 4.7 and = 1.6 for the velocity error when BDF2 is employed. Thus, the regularity condition that is supposed for second order convergence appears not only sufficient, but also necessary. Since the condition u * ∈ L 2 (0, ; H −1 ) in Corollary 4.8 is violated, but not u * ∈ L 2 (0, ; H +1 ), we conclude that the overall convergence is reduced due to the purely temporal contributions ( ), but not due to the mixed contributions (ℎ, ).
Next, we consider the case when all conditions imposed by Theorem 4.1 and Corollaries 4.7, 4.8 are satisfied, except of the small data condition. To this end, we set = 10 7 and = 10 12 . With this choice, equation (6.4) is violated since ≫ * . Moreover, the modified DEP force F reduces to the original version F , because the Euclidean norm of the potential gradients of the exact and discrete solution (both are of order (10 )) are way below . Tables 5 and 6 list the corresponding errors. Apparently, the respective convergence orders are almost identical to those obtained for = = 10 3 . These results indicate that the small data condition imposed by Theorem 4.1 might not be a necessary condition for the stated error estimates to hold.
In all considered cases, the pressure error converges faster than predicted by the theory. In fact, it converges with the same order as the velocity error in the high temporal regularity case = per . In the low regularity case = alg the pressure convergence rate is slightly higher (BDF1) and lower (BDF2) than its velocity counterpart. Overall, one cannot observe a rate that is reduced by 1 2 . In the proof of Theorem 4.1, the factor − 1 2 stems from estimating the term d (ℎ, ), (4.63), which basically measures the temporal variation of the discretization error (d u,ℎ ) . The similar quantity (ℎ, ), which measures the temporal variation of the exact solution (u ) , does Table 3. BDF1, = alg , = = 10 3 , ℎ = 1.25 .
Error not cause the occurrence of − 1 2 , since the temporal regularity of u could be exploited. From a theoretical point of view, this temporal regularity could not be derived for (d u,ℎ ) . In the underlying test problem, however, this discretization error might have a similar temporal smoothness as the exact solution, resulting in comparable -convergence rates for pressure and velocity.  Finally, Figures 2 and 3 visualize the long-term error evolution for the BDF1 and BDF2 scheme, respectively. To this end, the normalized velocity error is plotted over time for the different temporal functions = per (black curve) and = alg (red curve). Note that (1) is just the temporal factor in the exact solution. Thus, ( ) equals the relative velocity error, up to a time-independent constant. Apparently, the actual error does not exhibit the exponential growth that enters the estimates by means of Gröwall's inequality, see (4.43). Instead, ( ) stabilizes on a fairly constant level. However, the lack of temporal regularity of alg at = 0 can be observed through high relative errors in the beginning of the simulation.

Discussion on the constructed DEP force
Summing up the previous sections, we derived a discrete well-posedness result and error estimates for the TEHD Boussinesq equations, provided a DEP force F satisfying Assumption 3.13. Since these conditions are not fulfilled by the original force F , we constructed a family of modified forces {F , > 0}, based on a cut-off function applied to the gradient of the potential, see Section 5. With this choice, the theoretical results of Section 3 and 4 do hold for F = F , as long the cut-off threshold satisfies Note that (7.1) implies both previously used small data conditions, (3.11) and (4.16), when Ω is defined as in (6.5).
The domain-dependent constant Ω depends on the constants in Friedrich's inequality and time-stepping parameters and can be explicitly computed or estimated. This is done in Section 6 for an hypercube and below for a vertical annulus that is commonly used in physical experiments. In these cases, the upper bound can be calculated and one can draw the following conclusions on the error between the discrete solution obtained for F = F and the exact solution for F = F if the initially stated small data condition (1.6) does hold: Further suppose that the physical parameters and discretization fulfill Assumptions 2.2-3.10 and that the starting values are chosen to fulfill condition (3.15). Let the discrete solution (u ℎ , ℎ , ℎ , Φ ℎ ) be obtained by solving Problem 3.12 with cut-off DEP force F = F . Then, the error estimates stated by Corollary 4.7 (case I) and Corollary 4.8 (case II) remain valid for (u, , , Φ) = (u * , * , * , Φ * ).
The stated assertion now follows by the fact that F satisfies Assumption 3.13 and the validness of the small data conditions supposed by Lemma 3.14 and Theorem 4.1 through the definition of and Ω , see (6.5).
Corollary 7.1 creates the link between the original DEP formulation and the modified version, under the assumption that the exact electric field E = −∇Φ satisfies (7.2).
Summing up, condition (3.11) and (4.16) have their origin in the fact that there is a force term in the momentum equation, that is an element of U * and depends on ∇ . Bounding this term from above yields the product ‖∇ ‖ 1,2 ‖∇u‖ 1,2 before applying Young's inequality. After this step, the small data conditions are needed to ensure positive factors when collecting all ‖∇ ‖ 2 1,2 , ‖∇u‖ 2 1,2 terms on the left-hand side. Thus, (3.11) and (4.16) are induced by the applied proof technique. On the other hand, condition (7.2) stems from the fact that F also depends on |∇Φ| 2 , which cannot be handled properly if Φ is only in 1 . We introduced the cut-off based modification F to circumvent this issue. With the chosen construction, equations (3.11) and (4.16) are summarized by condition (7.1) which can actually always be met (if is made sufficiently small). Also by the way F is constructed, the additional condition (7.2) onto the exact potential implies that the modified force term reduces to the original one. The higher the critical value , the more likely it is that we find a that satisfies (7.1) and F ( , Φ) = F ( , Φ) for the exact solution; thus solving the original, unmodified equations.
We now present experiments for a realistic scenario, for which (7.2) is indeed fulfilled (for the numerical solution). To this end, we consider the DEP-driven flow inside the vertical annulus Ω := {( 1 , 2 , 3 ) : < 2 1 + 2 2 < , 3 ∈ (0, )} of height and inner / outer radii 0 < < . The top and bottom endplates Γ = Ω ∩ { 3 = 0}, Γ = Ω ∩ { 3 = } are supposed to be thermally and electrically insulating, i.e. there holds ∇ · n = ∇Φ · n = 0 with n denoting the unit outward normal. We further apply a temperature and potential difference between the inner and outer wall Γ = Ω ∩ { 2 1 + 2 2 = }, Γ = Ω ∩ { 2 1 + 2 2 = } by imposing Dirichlet boundary conditions for some temperature difference ≥ 0 and potential difference Φ ≥ 0. Suitable boundary liftings , Φ can easily be obtained by using cylinder coordinates and linear interpolation in the radius direction. Here, we will vary the potential difference Φ in the interval [700 V, 7000 V] and the initial conditions are given by u 0 , 0 as stationary solution of the TEHD equations with = 0 (i.e. the standard Boussinesq equations for natural convection) and Φ 0 is defined as solution of (1.4) with ( ) = = (i.e. the standard Gauss' law with constant permittivity, evaluated at room temperature). The considered fluid is the silicone oil AK5 with properties listed in Table 4.
The proposed scenario is often used in physical experiments to investigate the influence of DEP force onto heat transfer enhancement between inner and outer wall, see [35] and the references therein, for instance. Corresponding experiments are conducted under laboratory conditions (with standard gravitational acceleration g), as well as under microgravity conditions obtained during parabolic flights (g ≈ 0), see e.g. [27].
Here, we present numerical experiments for the laboratory case as done in chapter 6.3 of the first author's Ph.D. thesis [12]. In this scenario, the fluid flow undergoes a rapid change from the initial unicellular, natural convection-like, pattern to a final, stationary state which is made up by multiple, axially aligned vortices. See Figure 5 for a prototypical flow evolution.
The flow behaviour can be best understood by looking at the strength of axial vorticity plotted over time for different values of Φ, see Figure 6. The previously described transition from unicellular pattern to multiple vortices manifests itself in an exponential increase of axial vorticity from 0 to a certain, fixed final value. Both, exponential growth rate and final value, are monotonically increasing w.r.t. Φ. However, one can also observe   that for values of Φ below some critical threshold (here for Φ ∈ {2800 V, 3500 V}), axial vorticity stays around 0 and the flow pattern does not significantly deviate from the natural convection base state. This behaviour can also be observed experimentally [35] and can be explained by linear stability analysis [45].
Coming back to the small data condition, we note that the corresponding domain dependent constant Ω can be calculated in a similar way as in Section 6, where the constants of Friedrich's inequality for the annular shaped cylinder can be estimated from above by 0 , see Lemma A.107 in [12]. Using (6.5), (7.1), (7.3) the critical value can now be calculated: = {︃ 2.5 × 10 6 for BDF1 1.41 × 10 6 for BDF2. (7.4) In Figure 7 we plot minimal and maximal values of the discrete electrical field magnitude |∇Φ ℎ | over different applied potential differences Φ. The critical values are marked as dotted lines for both time stepping schemes. Apparently, the discrete potential satisfies (7.2) in case of BDF1 for all considered potential differences, whereas (7.2) is satisfied for BDF2 as long as Φ <≈ 5000 V. We thus showed that our crucial small data assumption (7.2) is fulfilled by the numerically obtained solution in an application relevant scenario over a parameter range of physical interest (note that the previously described transition of flow pattern can be observed for Φ ≥ 4200 V).
We finally provide some data that underlines the derived discrete energy estimates of Lemma 3.14. Figures 8 and 9 plot kinetic energy and energy dissipation over time for BDF1 and BDF2 and different time step sizes. The spatial resolution is kept fixed and Φ = 7000 V. Except of the case BDF2, = 0.1, only small differences between the different energy courses can be observed during the transition phase from unicellular to multi vortex flow pattern. In particular, the respective final states hardly differ and -independent energy bounds are observed, as predicted by Lemma 3.14. Only in case of BDF2, a time step size of = 0.1 appears to be too coarse to obtain an accurate solution, while the energy still stays within reasonable bounds.
Recall that we set F = F * , with * = 2.4 × 10 6 < (BDF1) for the presented simulation data, i.e. this choice for F fulfills the small data conditions (3.11) and (4.16) which are imposed by the stability and error  estimate results when BDF1 is used. In case of BDF2, however, F * does not satisfy the required assumptions.
It remains an open question, whether this is the reason for the observed behaviour for = 0.1.
Summing up, we have presented an application relevant scenario, for which the numerical potential indeed fulfills the required small data condition over a parameter range of physical interest. Further, we want to point out that if the sensitivity of the potential solution w.r.t. temperature is rather low for a given dielectric fluid (here = (10 −3 )), then the potential obtained by solving the TEHD equations only modestly deviates from the one obtained by setting = ( ) = . In case of constant , however, it is often possible to derive analytical expressions for Φ if the considered geometry takes the form of typical condensators which are often considered when investigating DEP driven flow, e.g. spherical gaps [46] and infinitely long cylindrical gaps [45]. For such an analytical form of Φ, one can easily determine a parameter range, over which the small data condition (7.2) holds. Regarding the previously investigated scenario for instance, it was shown in [12], that the relative difference between Φ ℎ and the analytical solution Φ for Gauss' law with constant in an infinitely long cylinder is only of order (10 −3 ) at each time instance, measured in ‖·‖ 1,2 norm. This ideal potential Φ fulfills (7.2) for Φ < 8620 V (BDF1) and Φ < 4860 V (BDF2), respectively.

Conclusion
In this paper, we proposed an H 1 -conformal FEM/BDF discretization for the instationary TEHD Boussinesq equations and showed well-posedness of the discrete equations. Moreover, we derived a priori error estimates in a general framework, that allows some freedom in the selection of the temporal discretization components. Based on this framework we defined a semi-implicit 1-and 2-step method and showed that the resulting schemes are first and second order accurate in time, respectively. In doing so, the main difficulty lied in the modeling of the cubic DEP term. Here, we proposed a modification of the original term by means of a cut-off operator. The theoretically obtained convergence rates for velocity, temperature and potential could be verified by numerical experiments in a 2D test case, whereas the pressure error converged faster than predicted. The experiments also demonstrated that a reduced convergence rate is obtained if the exact solutions lacks temporal regularity.
The presented error analysis relies on a small data condition that involves the maximal value of the electric field. The validness of this condition was shown for the presented test case and for a realistic scenario. However, the conducted experiments indicate that this might be not a necessary condition and further research has to be done to relax or even remove this condition by a more sophisticated analysis or a modification of the spatial discretization scheme. We further showed that the derived error estimates do also hold for the original DEP term, if the norm of the exact electric field is below a computable threshold.

=0
, then the restriction < 1 is not necessary and can be set to 1.