Eulerian time-stepping schemes for the non-stationary Stokes equations on time-dependent domains

This article is concerned with the discretisation of the Stokes equations on time-dependent domains in an Eulerian coordinate framework. Our work can be seen as an extension of a recent paper by Lehrenfeld and Olshanskii (ESAIM: M2AN 53(2):585–614, 2019), where BDF-type time-stepping schemes are studied for a parabolic equation on moving domains. For space discretisation, a geometrically unfitted finite element discretisation is applied in combination with Nitsche’s method to impose boundary conditions. Physically undefined values of the solution at previous time-steps are extended implicitly by means of so-called ghost penalty stabilisations. We derive a complete a priori error analysis of the discretisation error in space and time, including optimal L2(L2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(L^2)$$\end{document}-norm error bounds for the velocities. Finally, the theoretical results are substantiated with numerical examples.


Introduction
Flows on moving domains (t) ⊂ R d (d = 2, 3) need to be considered in many different applications. Examples include particulate flows or flows around moving objects like biological or mechanical valves, wind turbines or parachutes. Strongly related problems are fluid-structure interactions or multi-phase flows.
There exists a vast literature on time discretisation of the non-stationary Stokes or Navier-Stokes equations on fixed domains, see for example the classical works of Girault and Raviart [30], Baker et al. [2] and Rannacher and Heywood [39], or more recently Bochev et al. [5] and Burman and Fernández [11,12] in the context of stabilised finite element methods. If the computational domain remains unchanged in each time-step, the same spatial discretisation can be used (unless adaptive mesh refinement is considered) and finite difference schemes based on the method of lines can be applied for time discretisation.
In the case of moderate domain movements, these techniques can be transferred to the moving framework by using the Arbitrary Lagrangian Eulerian (ALE) approach [22,24,50]. Here, the idea is to formulate an equivalent system of equations on a fixed reference configurationˆ , for example the initial configuration (0), by means of a time-dependent map T(t) :ˆ → (t). This technique has been used widely for flows on moving domains, see e.g,. [19,22] and fluid-structure interactions [3,29,53]. The analysis of the time discretisation error is then very similar to the fixed framework, as all quantities and equations are formulated on the same reference domainˆ , see e.g. [54]. For a detailed stability analysis of ALE formulations, we refer to Nobile and Formaggia [49] and Boffi and Gastaldi [7].
On the other hand, it is well-known that the ALE method is less practical in the case of large domain deformations [25,53]. This is due to the degeneration of mesh elements both in a finite element and a finite difference context. A re-meshing of the domain (t) becomes necessary. Moreover, topology changes, for example due to contact of particles within the flow or of a particle with an outer wall [18], are not allowed, as the map betweenˆ and (t) can not have the required regularity in this situation.
In such cases an Eulerian formulation of the problem formulated on the moving domains (t) is preferable. This is also the standard coordinate framework for the simulation of multi-phase flows. In the last years a variety of space discretisation techniques have been designed to resolve curved or moving boundaries accurately. Examples include the cut finite element method [16,17,33,35,[44][45][46] within a fictitious domain approach, extended finite elements [20,31,36,47] or locally fitted finite element techniques [27], to name such a few of the approaches.
Much less analysed is a proper time discretisation of the problem. In the case of moving domains, standard time discretisation based on the method of lines is not applicable in a straight-forward way. The reason is that the domain of definition of the variables changes from time step to time step.
As an example consider the backward Euler discretisation of the time derivative within a variational formulation Note that u h (t n−1 ) is only well-defined on (t n−1 ), but is needed on (t n ).
One solution to this dilemma are so-called characteristic-based approaches [38]. Similar time-stepping schemes result when applying the ALE method only locally within one time-step and projecting the system back to the original reference frame after each step [21], or based on Galerkin time discretisations with modified Galerkin spaces [28]. The disadvantage of these approaches is a projection between the domains (t n−1 ) and (t n ) that needs to be computed within each or after a certain number of steps.
Another possibility consists of space-time approaches [37,41], where a d + 1dimensional domain is discretised if (t) ⊂ R d . The computational requirements of these approaches might, however, exceed the available computational resources, in particular within complex three-dimensional applications. Moreover, the implementation of higher-dimensional discretisations and accurate quadrature formulas pose additional challenges.
A simpler approach has been proposed recently in the dissertation of Schott [56] and by Lehrenfeld and Olshanskii [42]. Here, the idea is to define extensions of the solution u h (t n−1 ) from the previous time-step to a domain that spans at least (t n ). On the finite element level these extensions can be incorporated implicitly in the timestepping scheme by using so-called ghost penalty stabilisations [10] to a sufficiently large domain δ (t n−1 ) ⊃ (t n ). These techniques have originally been proposed to extend the coercivity of elliptic bilinear forms from the physical to the computational domain in the context of CutFEM or fictitious domain approaches [10].
While Schott used such an extension explicitly after each time step to define values for u h (t n−1 ) in mesh nodes lying in (t n )\ (t n−1 ), Lehrenfeld and Olshanskii included the extension operator implicitly within each time step by solving a combined discrete system including the extension operator on the larger computational domain δ (t n ). For the latter approach a complete analysis could be given for the corresponding backward Euler time discretisation, showing first-order convergence in time in the spatial energy norm [42]. Moreover, the authors gave hints on how to transfer the argumentation to a backward difference scheme (BDF(2)), which results in second-order convergence. We should also mention that similar time discretisation techniques have been used previously in the context of surface PDEs [43,51], and mixed-dimensional surface-bulk problems [37] on moving domains.
In this work, we apply such an approach to the discretisation of the non-stationary Stokes equations on a moving domain, including a complete analysis of the space and time discretisation errors. Particular problems are related to the approximation of the pressure variable. It is well-known that stability of the pressure is lost in the case of fixed domains, when the discretisation changes from one time-step to another. This can already be observed, when the finite element mesh is refined or coarsened globally at some instant of time, see Besier and Wollner [4] and is due to the fact that the old solution u n−1 h := u h (t n−1 ) is not discrete divergence free with respect to the new mesh. Possible remedies include the use of Stokes or Darcy projections [4,8] to pass u n−1 h to the new mesh. Our analysis will reveal that similar issues hold true for the case of moving domains, even if the same discretisation is used on (t n ) ∩ (t n−1 ). The reason is that u n−1 h is discrete divergence-free with respect to (t n−1 ), but not with respect to (t n ) (div u n−1 h , φ h ) (t n−1 ) = 0, but (div u n−1 h , φ h ) (t n ) = 0 for certain φ h ∈ V h . For space discretisation, we will use the Cut Finite Element framework [35]. The idea is to discretise a larger domain of simple structure in the spirit of the Fictitious Domain approach. The active degrees of freedom consist of all degrees of freedom in mesh elements with non-empty intersection with δ (t n ). Dirichlet boundary conditions are incorporated by means of Nitsche's method [48].
We will consider both the BDF(1)/backward Euler and the BDF(2) variant of the approach. To simplify the presentation of the analysis, we will neglect geometry approximation errors related to the approximation of curved boundaries and, moreover, focus on the BDF(1) variant. The necessary modifications for the BDF(2) variant will be sketched within remarks. Finally, we will use a duality technique to prove an optimal L 2 (L 2 )-norm estimate for the velocities.
The structure of this article is as follows: In Sect. 2 we introduce the equations and sketch how to prove the well-posedness of the system. Then we introduce time and space discretisation in Sect. 3, including the extension operators and assumptions, that will be needed in the stability analysis of Sect. 4 and the error analysis in Sect. 5. Then, we give some three-dimensional numerical results in Sect. 6. We conclude in Sect. 7.

Equations
We consider the non-stationary Stokes equations with homogeneous Dirichlet boundary conditions on a moving domain (t) ⊂ R d , d = 2, 3 for t ∈ I = [0, t fin ] We assume that the domain motion can be described by a W 1,∞ -diffeomorphism with the additional regularity and that the initial domain (0) is piecewise smooth and Lipschitz. In order to formulate the variational formulation we define the spaces and consider the variational formulation: Find u ∈ V I , p ∈ L 0,I such that where We assume that f ∈ L ∞ (I , L(t) d ) a.e. in t ∈ I and u 0 ∈ H 1 ( (0)) d .
Remark 2.1 (Boundary conditions) It might seem unnatural at first sight to use homogeneous Dirichlet boundary conditions for a Stokes problem on a moving domain (t). In fact the assumption that the flow follows the domain motion on ∂ (t) would be a more realistic boundary condition, i.e.
Note, however, that for a sufficiently smooth map T(t) that fulfils div(∂ t T(T −1 )) = 0, one obtains (1) from (7) for u =ũ−∂ t T(T −1 ) and f :=f ). For this reason and in order to simplify the presentation of the error analysis, we will consider homogeneous Dirichlet conditions in the remainder of this article.

Well-posedness
As the spaces in (4) are lacking a tensor product structure, the proof of well-posedness of (5) is more complicated than on a fixed domain. In the case of a fixed domain existence and uniqueness of solutions can be shown under weaker regularity assumptions on the data f and u 0 and the domain (t) in the velocity spacẽ where V * is the dual space to V. Low regularity is, however, not of interest for the present paper, as we will require additional regularity of the solution in the error estimates. On the other hand, working with the space V I under the additional regularity assumptions made above simplifies the proof of well-posedness of (5) significantly.
The well-posedness of the Navier-Stokes problem on time-dependent domains, including an additional nonlinear convective term has in fact been the subject of a number of papers in literature [6,55]. In order to deal with the additional nonlinearity, additional assumptions on the regularity of the domains are typically made. For completeness, we give a proof of the following Lemma (5) in the "Appendix" under the regularity assumptions made above. 0)). For f ∈ L ∞ (I , L(t) d ) and u 0 ∈ H 1 ( (0)) d , Problem (5) has a unique solution u ∈ V I , p ∈ L 0,I . Proof A proof is given in the "Appendix".

Discretisation
For discretisation in time, we split the time interval of interest I = [0, t fin ] into time intervals I n = (t n−1 , t n ] of uniform step size t = t n − t n−1 We follow the work of Lehrenfeld and Olshanskii [42] for parabolic problems on moving domains and use BDF(s) discretisation for s = 1, 2, where s = 1 corresponds to a backward Euler time discretisation. Higher-order BDF formulae are not considered here, due to their lack of A-stability [34]. Following Lehrenfeld and Olshanskii [42] we extend the domain n := (t n ) in each time point t n by a strip of size δ to a domain n δ , which is chosen large enough such that see also the left part of Fig. 1. In particular, we will allow where is the maximum velocity of the boundary movement in normal direction in the Euclidean norm · and c δ > 1 is a constant. If we assume that the domain map T lies in W 1,∞ (I , L ∞ ( (0))) (see Assumption 3.2 below), the lower bound on δ in (9) guarantees (8). The space-time slabs defined by the time discretisation and the space-time domain are denoted by In what follows we denote by c generic positive constants. These are in particular independent of space and time discretisation ( t, N and h) and of domain velocity w max and δ, unless such a dependence is explicitly mentioned.

Space discretisation
Let T n h,δ be a family of (possibly unfitted) quasi-uniform spatial discretisations of n δ into simplices with maximum cell size h. We assume that T n h,δ is based on a common background triangulation T h for all n and may differ only in the elements outside (t n ) that are not present in T k h,δ for k = n. Further, we assume that T n h,δ consists only of elements K with non-empty intersection with n δ , i.e. K ∩ n δ = ∅. The subset of cells with non-empty intersection with is denoted by T n h . An illustration is given in Fig. 1 and F n,ext h,δ , which will be used to define the ghost penalty extensions.
For spatial discretisation, we use continuous equal-order finite elements of degree m ≥ 1 for all variables Note that for the pressure space L n h an extension beyond n h is not required. To deal with the inf-sup stability, we will add a pressure stabilisation term s n h to the variational formulation. In order to simplify the presentation, we will restrict ourselves to the Continuous Interior Penalty method (CIP [14]) in this work, although different pressure stabilisations are possible. We define the CIP pressure stabilisation as The higher derivatives in the boundary elements are necessary to control the derivatives ∇ p h on the extended computational domain n h \ n in the spirit of the ghost penalty stabilisation [10].
We summarise the properties of the pressure stabilisation that will be needed in the following: There exists an operator C n h : V n h ∪ L 2 ( n ) → V n h , such that the following properties are fulfilled for n = 1, . . . , N A suitable projector C n h for the CIP stabilisation is given by the the Oswald or Clément interpolation [14]. For m ≥ 2, we have additionally the consistency property Remark 3.1 (Pressure stabilisation) In general any pressure stabilisation operator that leads to a well-posed discrete problem and that fulfils the assumptions (10)- (14) can be used. The consistency condition s n h ( p, p) = 0 can be relaxed to a weak consistency of order m s > 0 s n h ( p, p) ≤ ch 2m s p 2 H ms ( n ) .
which will limit the spatial convergence order in the error estimates. One possibility is the Brezzi-Pitkäranta stabilisation [9] with order m s = 1. We refer to Burman and Fernández for a review of further possibilities for pressure stabilisation [12].

Variational formulation
To cope with the evolving geometry from one time-step to another, we extend the velocity variable u n h to n δ , which will be needed in the following time-step, by using so-called ghost penalty terms g n h . We will describe different possibilities to define g n h in the next subsection. For k = s, . . . , n we define the following time-stepping scheme: where D (s) t is an approximation of the time derivative ∂ t by the BDF(s) backward difference formula, i.e.

The bilinear form A h is defined by
It includes the Stokes part and Nitsche terms to weakly impose the Dirichlet boundary conditions In (18) the last term can be seen as a penalty term to weakly impose the homogeneous Dirichlet condition for the velocities. The first term on the right-hand side makes the variational formulation consistent (in space). Finally, the second term, which vanishes for u n h = 0, yields a formulation, which is symmetric for the velocities, but skewsymmetric for the pressure. The skew-symmetry in the pressure variable leads to a stable variational formulation, as the pressure terms cancel out by diagonal testing (v n h = u n h , q n h = p n h ), see for example [13]. The parameters γ D , γ p and γ g are positive constants.
To include the initial condition, we set u 0 h := π 1 h Eu 0 , where π n h denotes the L 2 ( n )projection onto T n h and E denotes an L 2 -stable extension operator, which is introduced in the next section. Summing over k = 1, . . . , n in time, the complete system reads for s = 1 In order to simplify the presentation, we will assume that the integrals in (19) are evaluated exactly. If the integrals are only roughly approximated, for example due to a discrete level-set function φ n h ∈ V n h which is only an approximation of a continuous function φ n , an additional geometry approximation error needs to be considered. We refer to the work of Lehrenfeld and Olshanskii [42], where these additional error contributions have been analysed in detail for parabolic problems on moving domains. An advantage of the CutFEM methodology compared to standard finite elements is that besides the geometry approximation no additional discretisation errors related to the approximation of curved boundaries within the finite element spaces need to be considered.
To initialise the BDF(2) scheme the value u 1 h needs to be computed with sufficient accuracy before the first full BDF(2) step can be made. We will comment on the specific requirements and on different possibilities below in Remark 5.6.

Extension operators
Due to the evolution of the domain, we will frequently need to extend variables defined on smaller domains to larger ones. Therefore, we will use W k, p -stable extension operators E n : n → n δ to extend functions u(t n ) ∈ W k, p ( n ). We make the following assumption for the regularity of the domains (t) and the domain movement, depending on the polynomial degree m of the finite element spaces.

Assumption 3.2
We assume that the boundary of the initial domain (0) is piecewise smooth and Lipschitz, and that the domain motion T(t) is a W 1,∞diffeomorphism for each t and smooth in the sense that T ∈ L ∞ (I , W m+1,∞ ( (0)))∩ W 1,∞ (I , W m,∞ ( (0))).
If Assumption 3.2 is fulfilled for m ∈ N, suitable extension operators E n : n → n δ exist with the properties For a proof of (20) we refer to Stein [57], Theorem 6 in Chapter VI. The estimate (21) has been shown in [42], Lemma 3.3. The estimate (22) follows by the same argumentation. In order to alleviate the notation we will in the following skip the operator E n frequently and denote the extension also by u(t n ).

Ghost penalty extension
The discrete quantities are extended implicitly by adding so-called ghost penalty terms to the variational formulation. We will consider three variants for the ghost penalty stabilisation, and refer to [15,32] for a more abstract approach on how to design suitable ghost penalties for a PDE problem at hand. The first "classical" variant [10,15,45] is to penalise jumps of derivatives over element edges This variant has the advantage that it is fully consistent, i.e. it vanishes for u ∈ H m+1 ( n δ ) d , which implies ∂ k n u | e = 0 for k ≤ m. A disadvantage is that higher derivatives need to be computed for polynomial degrees m > 1.
To define two further variants, let us introduce the notation K e,1 and K e,2 for the two cells surrounding a face e ∈ F n,g h,δ , such that e = K e,1 ∩ K e,2 .
We denote the union of both cells by w e := K e,1 ∪ K e,2 and use the L 2 -projection π w e : L 2 ( n δ ) → P m (w e ), which is defined by We define the "projection variant" of the ghost penalty stabilisation [10] g n,proj h The last equality is a direct consequence of the definition of the L 2 -projection. The third variant, which has first been used in [52], uses canonical extensions of polynomials to the neighbouring cell instead of the projection π w e u. Let us therefore denote the polynomials that define a function u ∈ V n h in a cell K e,i by u e,i = u| K e,i . We use the same notation for the canonical extension to the neighbouring cell, such that u e,i ∈ P m (w e ). Using this notation, we define the so-called "direct method" of the ghost penalty stabilisation For the analysis, we extend the definition of the stabilisation to functions u, v ∈ L 2 ( n δ ). Here, we set u e,i := π K e,i u| K e,i for i = 1, 2, where π K e,i denotes the L 2projection to P m (K e,i ) and extend this polynomial canonically to the neighbouring cell. In contrast to the classical variant, g n,proj h and g n,dir h are only weakly consistent, i.e. they fulfil the estimate We will summarise the properties of these stabilisation terms, that we will need below, in the following lemma. Therefore, we assume that from each cell K ∈ T n h,δ with K ∩ n = ∅, there exists a path of cells K i , i = 1, . . . , m, such that two subsequent cells share one common face e = K i ∩ K i+1 , and the final element lies in the interior of n , i.e. K m ⊂ n . In addition the path shall fulfil the following properties. Let K be the maximum number of cells needed in the path among all cells K ∈ T n h,δ . We assume that where the second inequality follows from (9). Moreover, we assume that the number of cases in which a specific interior element K m ⊂ n is used as a final element among all the paths is bounded independently of t and h. These assumptions are reasonable, as one can choose for example the final elements by a projection of distance δ towards the interior. For a detailed justification, we refer to Lehrenfeld and Olshanskii [42], Remark 5.2.

Lemma 3.3 For v h ∈ V n h and the three variants g n
Proof The first four properties have been proven for the three possibilities introduced above by Lehrenfeld and Olshanskii [42]. The last inequality in (24) follows similarly.

Properties of the bilinear form
We start with a continuity result for the combined bilinear form including the Nitsche terms in the functional spaces.

Proof We apply integration by parts in (17)
For the Nitsche terms standard estimates result in The estimate for (A n S + a n D )(v, q; u, p) can be shown in exactly the same way by inverting the role of test and trial functions.
Next, we show continuity and coercivity of the discrete bilinear form. To this end, we introduce the triple norm

Lemma 3.5 (Coercivity and Continuity in the discrete setting) For the bilinear form
A h defined in (16) and u h ∈ V n h and p h ∈ L n h , it holds for γ D sufficiently large Proof To show coercivity (26), we note that To estimate the term −2(u n h , ∂ n u n h ) ∂ n , we apply a Cauchy-Schwarz and Young's inequality for > 0, followed by an inverse inequality on n for γ D sufficiently large. Concerning continuity, we estimate For the Nitsche terms, we have using inverse inequalities and Lemma 3.3 Finally, Lemma 3.3 and the assumption (10) for the pressure stabilisation yield Moreover, we have the following modified inf-sup condition for the discrete spaces.
Proof We follow Burman and Hansbo [14] and define v n Such a solution exists, see Temam [58], and fulfils v n p H 1 ( n ) ≤ c. We introduce an L 2 -stable interpolation i n h v n p (for example the Clément interpolation) to get We apply integration by parts in the first term and use that v n p vanishes in ∂ n The statement follows by noting that The well-posedness of the discrete system (15) for sufficiently large γ p , γ g , γ D and given u n−1 h (and u n−2 h for BDF (2)) follows by standard arguments, see for example [14].

Stability analysis
In order to simplify the analysis, we restrict ourselves in this and the next section to the case s = 1 first, i.e., the backward Euler variant of the time discretisation and comment on the case s = 2 in remarks. In order to abbreviate the notation, we write for the space-time Bochner norms where m ∈ Z and H 0 ( (t)) := L 2 ( (t)).
We start with a preliminary result concerning the extension of discrete functions to n δ . For with constants c 1 (w max ) := 1/2 + cs 2 w 2 max , c 2 (w max ) := cw 2 max h 2 + 1 and c > 0.
Proof These results follow similarly to Lemmas 3.4 and 5.3 in [42]. Nevertheless, we give here a sketch of the proof due to the importance of the Lemma in the following estimates. We define We apply a multiplicative trace inequality and Young's inequality for with a constant c 0 depending on the curvature of ∂ n . Integration over r ∈ (0, δ) yields (34).
Using (9) and choosing = for h < 1 with the constants c 1 (w max ), c 2 (w max ) given in the statement. The inequality (35) follows by combining (37) Now we are ready to show a stability result for the discrete formulation (15).
Under the regularity assumptions stated above, it holds for n ≥ 1 that with c 1 (w max ) given in Lemma 4.1 and u 0 h := π 1 h Eu 0 .
Proof Testing (15) with v h = 2 tu n h , q h = 2 t p n h , using the coercivity (26) and the relation We bring the term u n−1 h 2 n to n−1 by using Lemma 4.1 Inserting (42) into (41) we have for γ g ≥ c 2 (w max )K and γ D sufficiently large. For n = 1, we have instead of (41) In both cases (n ≥ 1) we use the Cauchy-Schwarz and Young's inequality for the first term on the right-hand side to get Summing over k = 0, . . . , n in (43) and using the L 2 -stability of the extension of the initial value yields Application of a discrete Gronwall lemma yields the statement. (2)) For the BDF(2) variant, we get the stability estimate (39) with

Remark 4.3 (BDF
To this end, one uses the relation instead of (40).

Stability estimate for the pressure
We show the following stability estimates for the L 2 -and H 1 -semi-norm of pressure.
To this end, we extend ∇ p n h by zero to n δ,h \ n h , using the same notation for the extended function. We insert ±C n h ∇ p n h , where C n h : L 2 ( n δ,h ) d → V n h is the interpolation operator used in (10)- (13), and integrate by parts For the first term, we have by means of (12) and Young's inequality The last term in (49) will be absorbed into the left-hand side of (48). For the second term on the right-hand side of (48), we use that (u n h , p n h ) solves the discrete system To estimate the first term on the right-hand side of (50), we use the Cauchy-Schwarz inequality and (12) to get Similarly, we get for the last term on the right-hand side of (50) For the second-term on the right-hand side of (50), we use an inverse inequality on n h and (12) For the Nitsche term a n D , we have as in (29) −h 2 a n In the last step (12) has been used. Note that the boundary term on the right-hand side will cancel out with the third term in (48). For the ghost penalty we have by means of an inverse inequality and (12) Altogether, (50)- (51) and (53) In the last step we have applied Young's inequality. Combination of (48), (49) an (55) yields (47). To show (46) we start using the modified inf-sup condition (Lemma 3.6) By (15), we have To estimate the right-hand side of (57), we use the continuity of the bilinear form A n h (27) and the Cauchy-Schwarz inequality In the last step, we have used that v n For s = 2, we have Proof We start by proving (59) for s = 1. To this end, we distinguish between the cases t ≥ h 2 and t < h 2 . In the first case, we note that, by (47) that follows from the triangle inequality and (35). The estimate (60) follows by a similar argumentation by distinguishing between the cases t ≶ h. is not discrete divergence-free with respect to n (div u n−1 h , ξ n h ) n = 0 for certain ξ n h ∈ L n h . Moreover, the domain mismatch n−1 = n causes additional problems in the transfer of the term |||u n−1 h ||| h,n = |||u n−1 h ||| h,n−1 from one time level to the previous one. In the error analysis developed in the following section, we will therefore use the H 1 -stability results in Corollary 4.5 for the pressure variable.

Error analysis
The energy error analysis for the velocities follows largely the argumentation of Lehrenfeld and Olshanskii [42] and is based on Galerkin orthogonality and the stability result of Theorem 4.2. We write u n := u(t n ), p n := p(t n ) and introduce the notation e n u := u n − u n h , η n u := u n − I n h u n , ξ n h,u := I n h u n − u n h , e n p := p n − p n h , η n p := p n − i n h p n , ξ n h, p := i n h p n − p n h for n ≥ 1, where I n h denotes the standard Lagrangian nodal interpolation to T n h,δ and i n h a generalised L 2 -stable interpolation (for example the Clément interpolation) to T n h . Moreover, we set This is possible, as u 0 h cancels out in the summed space-time system (19). The following estimates for the interpolation errors are well-known We will again make use of the extension operators E n introduced in Sect. 3.3. For better readability, we will sometimes skip the operators E n assuming that quantities that would be undefined on the domains of integration are extended smoothly.
For the error analysis, we assume that the solution (u, p) to (5) lies in H m ( (t))) for m ≥ 1. Then, we can incorporate the Nitsche terms in the variational formulation on the continuous level and see that (u, p) is the solution to whereṼ (t) := H 1 ( (t)) d .

Energy error
As a starting point for the error estimation, we subtract (15) from (65) to obtain the orthogonality relation for n ≥ s with the consistency error E n c (v h , q h ). Note that this relation holds in particular also for n = s, as we have defined e 0 u = 0. We have used a different splitting in the pressure stabilisation s n h compared to the other terms, in order to include the case p ∈ H 1 ( ) (m = 1), where s n h ( p n , q h ) would not be well-defined.
We further split (66) into interpolation and discrete error parts where the interpolation error is defined by We will apply the stability result of Theorem 4.2 to (67), which will be the basis of the error estimate. For better readability, we will restrict restrict ourselves again to the case s = 1 first. Let us first estimate the consistency and interpolation errors.
Proof For the first part of the consistency error, we have using integration by parts and a Cauchy-Schwarz inequality in time The extension operator E n is needed, as the integration domain in the left-hand side of (69) includes parts, that lie outside the physical domain Q n . For the ghost penalty part, we have with Lemma 3.3 and the H m+1 -stability of the extension (20) Concerning the pressure stabilisation, we note that for p n ∈ H 1 ( n ) the term s n h ( p n , p n ) is not well-defined. For this reason we distinguish between the cases m = 1 and m ≥ 2. In the first case, we estimate using (10) and the H 1 -stability of the interpolation For m ≥ 2, we insert ± p n and use (10), (14) and the interpolation error estimate (64) Proof We estimate the interpolation error (68) term by term. For the first term we use that we can exchange time derivative and interpolation operator ∂ t I h u n = I h ∂ t u(t n ) We note again that the integration domain in the first norm on the right-hand side includes parts, that might lie outside the physical domain Q n . By means of (21) we conclude For the second term in (68), we use Lemma 3.4 Finally, we get for the ghost penalty part from (24) and the H m+1 -stability of the extension Now, we are ready to show an error estimate for the velocities.
be the discrete solution of (15) for s = 1 and (u, p) the continuous solution of (5). Further, let γ g ≥ c 2 (w max )K with c 2 (w max ) defined in Lemma 4.1, γ D , γ p sufficiently large and t ≥ ch 2 for some c > 0. Under the assumptions stated in Sect. 3, including Assumption 3.2, it holds for the error e k where e 0 u := 0 and c 1 (w max ) is defined in Lemma 4.1. Proof As in the stability proof (Theorem 4.2, (43)), we obtain from (67) We multiply (73) by > 0 and add it to (72). Due to the assumption t ≥ ch 2 the first three terms on the right-hand side of (73) can be absorbed into the left-hand side of (72) for sufficiently small We sum over k = 1, . . . , n and apply a discrete Gronwall lemma to find Finally, the interpolation estimates (62)-(64) and the argumentation used in (70) yield Addition of (76) and (77) proves the statement.

Remark 5.4 (Optimality)
which is of second order in time t, if we assume that the initial error is bounded by The initialisation will be discussed in the following remark. The main modifications in the proof concern the approximation of the time derivative in Lemmas 5.1 and 5.2. In (69) we estimate see [12,34]. In order to estimate the analogue of (70), we use Then the argumentation used in (70) can be applied to both terms on the right-hand side of (79). (2)) To initialise the BDF(2) scheme, the function u 1 h needs to be computed with sufficient accuracy. The simplest possibility is to use one BDF(1) step by solving

Remark 5.6 (Initialisation of BDF
. Similar to the proof of Theorem 5.3, the error after one BDF(1) step can be estimated by where in the last step a Sobolev inequality has been applied in time to show ∂ 2 t u 2 Q 1 ≤ t ∂ 2 t u 2 ∞,0,I 1 ≤ c t u 2 H 3 (I ,L 2 ( n )) .

L 2 (L 2 )-norm error of pressure
The energy estimate in Theorem 5.3 includes an optimal bound for the H 1 -norm of the pressure. To show an optimal bound in the L 2 -norm seems to be non-trivial, due to the fact that u n−1 h is not discrete divergence-free with respect to n and V n h , see the discussion in Sect. 4.1. We show here only a sub-optimal bound for s = 1. An optimal estimate is subject to future work.
Proof We use the modified inf-sup condition for the discrete part ξ n h, p = i n h p n − p n h and standard interpolation estimates The second term on the right-hand side is bounded by the energy estimate. For the first term, we use Galerkin orthogonality (66), followed by Cauchy-Schwarz and Poincaré inequalities After summation in (80), we obtain Using the standard interpolation estimate η k , we see that (81) holds for ξ k h, p replaced by e k p . Finally, Theorem 5.3 yields the statement. Unfortunately, the factor 1 t in front of the first term on the right-hand side of (81) leads to a loss of t −1/2 in the final estimate. (2)) For s = 2 we can only control t 3 D (2) t e n u 2

Remark 5.8 (BDF
n for s = 1), which leads to a further loss of t −1 in the above estimate:

L 2 (L 2 )-norm error of velocity
To obtain an optimal bound for the velocity error in the L 2 -norm, we introduce a dual problem. The argumentation of Burman and Fernández [12], that does not require a dual problem, but is based on a Stokes projection P h (u, p) of the continuous solution, can not be transferred in a straight-forward way to the case of moving domains, as it requires an estimate for the time derivative ∂ t (u − P u h u). Time derivative and Stokes projection do, however, not commute in the case of moving domains, as P u h u(t) depends on the domain (t). For this reason an estimate for the time derivative is nontrivial. We focus again on the case s = 1 first and remark on how to transfer the argumentation to the case s > 1 afterwards. The argumentation will be based on a semi-discretised (in time) dual problem. Before we introduce the dual problem, let us note that the semi-discretised primal problem is given by : Find (u k , p k where E k denotes the smooth extension operator to k δ introduced in Sect. 3.3. The corresponding semi-discretised dual problem, which will be needed in the following, reads: Note that the Dirichlet conditions are imposed strongly in this formulation and the bilinear form A k S does not include the Nitsche terms. We start by showing the well-posedness of the problem (83).

Lemma 5.10
Let s = 1, e k u ∈ L 2 ( k ) for k = 1, . . . , n and assume Assumption 3.2. The semi-discrete dual problem (83) defines unique solutions (z k u , z k p ) n k=1 with regularity z k u ∈ H 2 ( k ), z k p ∈ H 1 ( k ). Moreover, the following regularity estimates are valid, where S k δ := k δ \ k and D (1) Proof By testing (83) withφ l u = δ kl φ k u ,φ l p = δ kl φ k p , l = 1, . . . , n, where δ kl is the Kronecker delta, we observe that the system splits into n separate time steps, where each step corresponds to a stationary Stokes system with an additional L 2 -term coming from the discretisation of the time derivative. For k < n we have and for k = n As the corresponding reduced problems are coercive in the velocity space V 0 (t k ) (cf. Sect. 2.1), existence and uniqueness of solutions z k u ∈ H 1 0 ( k ), z k p ∈ L 2 0 ( k ) follow inductively by standard arguments for k = n, . . . , 1, see e.g. Temam [58], Section I.2.
To show the regularity estimates (84) and (85), let us re-formulate the problems (86) and (87) in the following way: For k < n we have and for k = n A n S (φ n u , φ n p ; z n u , z n p ) = e n u , φ n u n − If we can prove that F k lies in the dual space [L 2 ( k ) d ] * , Proposition I.2.2 in Temam's book [58] guarantees the regularity estimate We need to show that the right-hand side is bounded. Splitting the first integral on the right-hand side into an integral over k and S k δ , we have for k < n and thus, For k = n, we obtain The boundedness of F k follows by induction for k = n, . . . , 1 and by using the stability of the extension operator E k . Combination of (89) and (90), resp. (91), yield the regularity estimates (86) and (87).
Next, we derive a stability estimate for the semi-discretised dual problem (83). We remark that a stability estimate for the continuous dual problem, including the first time derivative ∂ t z, could be obtained as well. This is however not enough to bound the consistency error of the time derivative in a sufficient way for an optimal L 2 -norm error estimate.

Lemma 5.11 Let the assumptions made in Sect. 3 be valid, including Assumption 3.2.
For sufficiently small t < ξ, where ξ depends only on c δ , w max and the domains k , k = 1, . . . , n, the solution (z k u , z k p ) n k=1 to the semi-discretised dual problem (83) for s = 1 fulfils the stability estimate Proof We show a stability estimate for the first derivatives ∇z k u first. For better readability we will in the following skip the extension operators E k and denote the extension E k z k+1 u still by z k+1 u and similarly for other variables. Diagonal testing in (83) As z k u vanishes on ∂ k , a Poincaré-like estimate gives in combination with (9) and the stability of the extension operator where c p denotes a constant depending on the domain k and c δ > 1 is the constant in (9). Using Young's inequality, this implies for t ≤ (2c 2 We obtain For the right-hand side in (92), we apply the Cauchy-Schwarz, a Poincaré and Young's inequality to get For k < n Lemma 5.10 gives us We estimate the term on S k For sufficiently small t < 1 cc 2 p c 3 δ w 3 max we obtain from (103)-(105) For the third term in (102), we use a telescope argument To bring the last term to k+1 , we estimate using (34) For the boundary term in (102), we use Green's theorem on S k For the second term on the right-hand side in (107) we use (34) twice with < 1, followed by (9), the stability of the extensions and Young's inequality For the last term in (107) we obtain as in (104) In the last inequality, we have used (9) and Young's inequality. Together, (107)-(108) yield the estimate To estimate the pressure term in (102), we obtain as in (108) To summarise we have shown that For k = n we obtain from (100) tested with φ n u = z n u and φ n p = z n p that Summation in (109) over k = 1, . . . , n−1 and addition of (110) and (98) multiplied by a factor of 3 yields for < 1 12c 0 Using (96) we can estimate the last term by which completes the proof. Now we are ready to prove an error estimate for the L 2 (L 2 )-norm of the velocities. First, we note that, due to the regularity proven in Lemma 5.10, the solution (z k u , z k p ) n k=1 of (83) is also the unique solution to the Nitsche formulation: The main difference in the proof is that the energy norm estimate does not give a bound for t D (2) t e k u k , see Remark 5.5. We have using (34) The L 2 -term on the right-hand side can then be absorbed into the left-hand side of (116) to obtain (117).

Numerical example
To substantiate the theoretical findings, we present numerical results for polynomial degrees m = 1, 2 and BDF formulas of order s = 1, 2. The results have been obtained using the CutFEM library [16], which is based on FeNiCS [1]. We consider flow through a 3-dimensional rectangular channel with a moving upper and lower wall in the time interval is I = [0, 2]. The moving domain is given by Due to the simple polygonal structure of the domain (t), the integrals in (19) are evaluated exactly within the CutFEM library [16] and we can expect higher-order convergence in space for m ≥ 2.
The data f and u D is chosen in such a way that the manufactured solution  To study the convergence orders in space and time, we show values for four different time-step and four different mesh sizes in Table 1. For P 1 finite elements, the finest mesh contains approximately 143.000 degrees of freedom. We observe that the temporal error is barely visible in the L 2 (L 2 )-norm and L 2 (H 1 )-semi-norm of velocities, as the spatial error is dominant. The spatial component of the velocities converges as expected by the theory (Theorems 5.3, 5.12) with orders 2 and 1. On the other hand, the temporal error shows up clearly in the pressure norms. To compute an estimated order of convergence (eoc), let us assume that the overall error can be separated into a temporal and a spatial component To estimate for instance the temporal order of convergence eoc t , we fit the three parameters g h , c t and eoc t of the function for a fixed mesh size h ∈ { 1 2 , 1 4 , 1 8 , 1 16 } against the computed values. This is done by means of a least-squares fit using gnuplot [40]. The values for g h and eoc t in the first row are for example computed by fitting the previous values in the same row (i.e. those obtained with h = 1 2 for different time-step sizes). A spatial order of convergence eoc h is estimated similarly using the values for a fixed time-step size t ∈ {0.4, 0.2, 0.1, 0.05}, i.e. those in the same column. For the pressure norms the estimated temporal order of convergence is very close to 1 in both the L 2 -and the H 1 semi-norm. This is expected for the L 2 (H 1 )-semi-norm by Theorem 5.3, but better than proven in Lemma 5.7 for the L 2 (L 2 )-norm. The spatial component of the error converges much faster than expected with eoc h around 2 for both norms (compared to O(1), which has been shown for the H 1 -semi-norm, and O(h) for the L 2 -norm). This might be due to superconvergence effects, as frequently observed for CIP stabilisations (see e.g. [26]), and possibly due to the sub-optimality of the pressure estimates.
The convergence orders of both pressure norms are very similar, especially for larger h and t. Here it seems that due to the superconvergence of the L 2 (H 1 )-semi-norm the simple Poincaré estimate e k p k ≤ c P ∇e k p k is optimal for the L 2 (L 2 )-norm. Only for smaller t and h, the convergence of the L 2 (L 2 )-norm seems to be slightly faster compared to the L 2 (H 1 )-semi-norm.

P 2 -BDF(1)
In order to increase the visibility of the temporal error component, we increase the order of the spatial discretisation first. In Table 2 we show results for P 2 finite elements and BDF(1) (m = 2, s = 1) on three different mesh levels. For P 2 the finest mesh level has again around 143.000 degrees of freedom, which is similar to P 1 elements on the next-finer mesh level. Again the spatial error is dominant in the velocity norms on coarser meshes and shows convergence orders of approximately 3 in the L 2 (L 2 )-norm and 2 in the L 2 (H 1 )-semi-norm, as shown in Theorems 5.12 and 5.3. In contrast to P 1 elements, the temporal error is however visible on the finest mesh level, where eoc t is close to 1, as expected.
In the L 2 (L 2 )-norm of pressure, the temporal error is dominant and shows again a convergence order of O( t). Due to the dominance of the temporal component, it is less clear to deduce the spatial error contribution. From the values and the eoc h it seems to converge again faster as predicted. Concerning the L 2 (H 1 )-norm of pressure, the assumption (118) that the spatial and temporal error are separated, which was assumed in order to compute eoc t and eoc h , is not valid, as the extrapolated values g h and g t do not or converge only very slowly towards zero. For this reason, the computed convergence orders eoc t and eoc h are not meaningful in this case. This does not contradict the theory, as Theorem 5.3 guarantees only the bound n k=1 t ∇e k

P 2 -BDF(2)
Finally, we show results for m = 2 and s = 2 in Tables 3 and 4 . In Table 3, we use an extension of the analytically given solution u(x, t) to t < 0 for initialisation, i.e. we use the starting values u 0 = 0 and u −1 := u(− t) in the first time step. Due to the (expected) second-order convergence in time, the temporal error is barely visible in the velocity norms on the finer mesh levels, in contrast to the results for BDF (1). The estimated order of convergence of the spatial component lies slightly below the orders 3 and 2 in the L 2 (L 2 )-norm and L 2 (H 1 )-semi-norm, respectively, that have been shown analytically. In the L 2 (L 2 )-norm of pressure both temporal and spatial errors are visible. Both eoc h and eoc t are around 2, which has been shown in Sect. 5.1.1 for the spatial part. For the temporal part only a reduced order of convergence of O( t) has been shown theoretically. This bound seems not to be sharp in the numerical example studied here. In the L 2 (H 1 )-semi-norm of pressure the spatial error is dominant, which is in contrast to the BDF(1) results. However, the assumption (118) that the error allows for a separation into spatial and temporal error components is again not valid, which makes the computed values of eoc t and eoc h meaningless.
In Table 4 we show results, where -instead of analytical values-one BDF(1) step has been used for initialisation, according to the discussion in Remark 5.6. We see that for large t the errors are slightly larger, due to the additional initial error. In fact the velocity norm errors are still relatively close, in particular for smaller t, while a stronger impact is visible in the pressure norms. These deviations get, however, smaller for t → 0. Moreover and most importantly, the estimated orders of convergence are very similar in the velocity norms and lie still significantly above the theoretical predictions in the pressure norms. This confirms numerically that the initialisation with BDF(1) is indeed sufficient to preserve the theoretically predicted convergence orders. The experimental orders of convergence (eoc) have been computed as in Table 1

Conclusion
We have derived a detailed a priori error analysis for two Eulerian time-stepping schemes based on backward difference formulas applied to the non-stationary Stokes equations on time-dependent domains. Following Schott [56] and Lehrenfeld and Olshanskii [42] discrete quantities are extended implicitly by means of ghost penalty terms to a larger domain, which is needed in the following step of the time-stepping scheme.
In particular, we have shown optimal-order error estimates for the L 2 (H 1 )-seminorm and the L 2 (L 2 )-norm error for the velocities. The main difficulties herein consisted in the transfer of quantities between domains n and n−1 at different The experimental orders of convergence (eoc) have been computed as in Table 1 time-steps and in the estimation of the pressure error. Optimal L 2 (H 1 )-norm errors for the pressure can be derived under the inverse CFL conditions t ≥ ch 2 for the CIP pressure stabilisation and BDF(1) ( t ≥ ch for BDF(2)), or unconditionally, when the Brezzi-Pitkäranta pressure stabilisation is used. Fortunately, these estimates are sufficient to show optimal bounds for the velocities in both the L 2 (H 1 )and the L 2 (L 2 )-norms. All these estimates are in good agreement with the numerical results presented. For the L 2 (L 2 )-norm error of the pressure, we have shown suboptimal bounds in terms of the time step t. The derivation of optimal bounds seems to be non-trivial and The experimental orders of convergence (eoc) have been computed as in Table 1 needs to be investigated in future work. Moreover, it would be interesting to further investigate if the exponential growth in the stability and error estimates can indeed be observed in numerical computations, for example by considering more complex domain motions. Further directions of research are the application of the approach to the non-linear Navier-Stokes equations, multi-phase flows and fluid-structure interactions, as well as the investigation of different time-stepping schemes, such as Crank-Nicolson or the fractional-step θ scheme within the framework presented and investigated in the present work.
We will show the well-posedness of (122) by a Galerkin argumentation. A basis ŵ j j∈N of the time-dependent spaceV 0 (t) is given by the inverse Piola transform of an L 2 -orthonormal basis φ j j∈N of the spaceV 0 (0) Under the given regularity assumptions on the domain movement T , the basis functions lie in W 1,∞ (I , H 1 ( (0)) d ).