No pressure? Energy-consistent ROMs for the incompressible Navier-Stokes equations with time-dependent boundary conditions

This work presents a novel reduced-order model (ROM) for the incompressible Navier-Stokes equations with time-dependent boundary conditions. This ROM is velocity-only, i.e. the simulation of the velocity does not require the computation of the pressure, and preserves the structure of the kinetic energy evolution. The key ingredient of the novel ROM is a decomposition of the velocity into a field with homogeneous boundary conditions and a lifting function that satisfies the mass equation with the prescribed inhomogeneous boundary conditions. This decomposition is inspired by the Helmholtz-Hodge decomposition and exhibits orthogonality of the two components. This orthogonality is crucial to preserve the structure of the kinetic energy evolution. To make the evaluation of the lifting function efficient, we propose a novel method that involves an explicit approximation of the boundary conditions with POD modes, while preserving the orthogonality of the velocity decomposition and thus the structure of the kinetic energy evolution. We show that the proposed velocity-only ROM is equivalent to a velocity-pressure ROM, i.e., a ROM that simulates both velocity and pressure. This equivalence can be generalized to other existing velocity-pressure ROMs and reveals valuable insights in their behaviour. Numerical experiments on test cases with inflow-outflow boundary conditions confirm the correctness and efficiency of the new ROM, and the equivalence with the velocity-pressure formulation.


Introduction
Turbulent flow is an ubiquitous phenomenon in fluid flows such as blood flow in arteries, ocean currents and air flow around wind turbines.Accurate numerical simulation of these turbulent flows requires computationally expensive methods such as Direct Numerical Simulation (DNS) [1] or Large Eddy Simulation (LES) [2].The computational cost of these methods can sometimes be tolerated for a small number of simulations but becomes prohibitive when it comes to tasks that require many simulation runs such as in the case of optimization and uncertainty quantification studies.
To mitigate the computational costs, several techniques have been developed to reduce the complexity of the respective full order model (FOM) resulting in so-called reduced-order models (ROMs).Many of these techniques are projection-based, i.e., they project the high-dimensional discretization of the model obtained from classical discretization methods such as Finite Elements, Finite Volumes or Finite Differences onto a low-dimensional space, the reduced basis space.Methods to obtain such reduced bases include sampling strategies [3] and Proper Orthogonal Decomposition (POD) [4].
The POD, which is also known as Karhunen-Loève expansion and Principal Component Analysis, is optimal in the sense that it extracts the reduced basis from given high-resolution data which minimizes the approximation error with respect to this data in the l 2 norm [4].However, several authors observed that ROMs resulting from the Galerkin projection of such POD reduced bases tend to exhibit instabilities when applied to fluid flow problems [5,6].This observation is often related to energy dissipation at the smallest scales of turbulent flow, which are not included in the POD reduced basis [5].Hence, approaches to address the ROM instability include methods that add such small scale effects either by introducing additional dissipation via closure modelling [7,8] or by using a POD based on a discrete equivalent of the H 1 norm instead of the l 2 norm in order to achieve a better representation of small scales in the reduced basis [9].Alternative approaches use more dissipative ROM bases whose modes are weighted sums of the standard POD modes and randomly chosen higherorder POD modes, which are assumed to represent the small scales [10,11].
Another cause of instability has been identified in the ROM pressure field.This pressure field can violate an inf-sup condition at the ROM level, even if the inf-sup condition of the underlying FOM is satisfied [12,6].When violating the condition, the ROM can exhibit spurious pressure modes that induce instability.Hence, many works have been developed to satisfy this condition on the ROM level [12,6,[13][14][15].At the same time, it has been shown that even for a linear ordinary differential equation (ODE), a stable FOM can result in an unstable POD-Galerkin ROM [16].Consequently, POD-Galerkin ROM instabilities are not exclusively related to the incompressible Navier-Stokes equations.
A promising approach to avoid ROM instabilities, not specific to the incompressible Navier-Stokes equations, is the idea of preserving certain structure of the FOM such that stability is implied.A method widely used for compressible flows [17][18][19][20] preserves the stability of an equilibrium point at the origin by performing POD with respect to energy-based inner products.For Hamiltonian systems, Proper Symplectic Decomposition was proposed as a modification of POD to obtain Galerkin ROMs that preserve symplecticity, so the Hamiltonian is time-invariant and the ROM is stable [21].Alternative symplecticitypreserving bases have been proposed that can be computed more efficiently [22] and recover the FOM Hamiltonian exactly [23].The concept of symplecticity preservation has been extended to model reduction on nonlinear manifolds [24] and non-canonical Hamiltonian systems [25].Furthermore, non-intrusive methods have been developed for canonical and noncanonical Hamiltonian systems that preserve the Hamiltonian structure [26].For mechanical systems, Lagrangian structurepreserving methods have been proposed [27,28] and observed to be stable [28].Stability-preserving methods have also been proposed for port-Hamiltonian systems [29,30] and second-order index-1 descriptor systems [31].For metriplectic systems, a method has been proposed that preserves the metriplectic structure and is observed to be long-term stable [32].For nonlinear conservation laws, stable ROMs have been proposed that preserve a semi-discrete entropy dissipation law [33].
For the incompressible Navier-Stokes equations, the POD-Galerkin ROM proposed in [34] achieves nonlinear stability by preserving the conservation of kinetic energy in the inviscid limit for homogeneous boundary conditions.In this work, this ROM is generalized to arbitrary time-dependent inhomogeneous boundary conditions, while possessing three properties: (i) structure preservation (consistency with the kinetic energy equation), (ii) velocity-only formulation (pressure-free), and (iii) arbitrary inhomogeneous boundary conditions (not separable in space-and time-dependent functions).To our best knowledge such a ROM formulation does not yet exist.
To achieve the first two properties, we propose the use of a specific lifting function.Classically, lifting functions are chosen as a velocity field that satisfies the prescribed inhomogeneous boundary conditions, so that the ROM can be solved with homogeneous boundary conditions [35][36][37][38].Inspired by the Helmholtz-Hodge decomposition [39], we propose a lifting function that is orthogonal to the homogeneous velocity field.This orthogonality allows us to preserve the structure of the kinetic energy equation on the ROM level.
In addition, this choice of lifting function has the benefit that the pressure is eliminated from the formulation, leading to a velocity-only ROM, i.e., the ROM velocity can be simulated without requiring the computation of the pressure.An important benefit of eliminating the pressure is that our velocity-only ROM is not affected by the inf-sup instability mentioned above.We stress that the omission of the pressure is not an approximation as performed in some early works [40,41].As shown in [42], merely neglecting the pressure in ROMs can lead to significant errors.In our ROM, in contrast, the elimination of the pressure is rigorous and exact.In fact, we show that the proposed velocity-only ROM is equivalent to a velocity-pressure ROM (for a specific choice of basis).This equivalence gives us valuable insights in the behaviour of the ROMs.Compared to other formulations that eliminate the pressure such as the streamfunction-vorticity formulation [43], our approach has the advantage of partial decoupling: while streamfunction and vorticity depend on each other, our lifting function does not depend on the remaining velocity component, but only on the prescribed boundary conditions.This partial decoupling has many advantages, for example, time integration errors do not accumulate in the lifting function.
Lastly, to be able to deal with arbitrary inhomogeneous time-dependent boundary conditions, we propose to approximate the boundary conditions by projection onto a POD basis.This explicit approximation of boundary conditions is necessary because standard hyper-reduction of the ROM ODE right-hand side [44,45] would not effect a sufficient efficiency improvement in context of our proposed lifting function.
This article is organized as follows.In Section 2, we introduce the continuous formulation and our finite volume discretization of the incompressible Navier-Stokes equations.In Section 3, we present our novel velocity-only ROM for time-dependent boundary condition.In Section 4, we show that this ROM preserves the structure of the kinetic energy evolution.In Section 5, we show the ROM's equivalence to a velocity-pressure ROM.In Section 6, we test the novel ROM in two test cases involving time-dependent inflow boundary conditions.

The incompressible Navier-Stokes equations
The incompressible Navier-Stokes equations describe the conservation of mass and momentum of a fluid in a domain ⊂ R d over a time interval [0, T ], and can be expressed as with divergence-free initial condition where the vector field u : × [0, T ] → R d describes the velocity, the scalar field p : × [0, T ] → R the pressure, ν is the kinematic viscosity, ρ the (constant) density (which will be absorbed in the pressure in subsequent equations), and f a forcing term (possibly space-and time-dependent).A set of boundary conditions which is appropriate to ensure uniqueness of the solution u (see [46,Section 3.8]) is assumed to be prescribed.

Full order model: finite volume discretization
The Navier-Stokes equations ( 1)-( 2) in d = 2 dimensions are spatially discretized with the second-order finite volume method on a staggered grid proposed in [34] and depicted in Fig. 1.The extension to d = 3 is straightforward, and all the methods described in the subsequent sections also apply to that case.
The uand v-component of velocity, and the pressure p are discretized by respectively, where N p = N x × N y is the number of "pressure finite volumes" and N u , N v ≈ N p are the numbers of "velocity component finite volumes" which depend on the type of boundary conditions.The velocity component vectors are aggregated in one velocity vector The spatial discretization of (1)-( 2) can be summarized as ( Here, h ∈ R N V ×N V is a diagonal matrix with the sizes of the velocity component finite volumes corresponding to the entries of V h on its diagonal.The matrix M h ∈ R N p ×N V describes the mass equation ( 1) such that each row of (4) represents conservation of mass over the finite volume around a pressure point p i, j as depicted in Fig. 1, using the midpoint rule.
The vector y bc (t) ∈ R N bc constitutes the time-dependent parametrization of inhomogeneous boundary conditions.In this work, we only consider time-dependence of inflow boundary conditions.Hence, y bc (t) comprises the prescribed Dirichlet velocity values at finite volume faces on inflow boundaries and N bc is the number of such faces.However, the construction of y bc (t) can be generalized to test cases with other types of time-dependent boundary conditions.Note that N bc scales only by O h −d+1 , while N V and N p scale by O h −d upon decreasing the mesh size h.
The matrix F M ∈ R N p ×N bc maps the vector y bc (t) to the mass equation right-hand side y M (t) ∈ R N p .The term , y bc (t)) contains the discretized convection, diffusion and forcing terms as well as all boundary condition contributions.The discretized forcing term f (t) ∈ R N V is assumed fixed and time-independent, hence neither f (t) nor t appear as an argument of F C D h .The term G h p h is the discretization of the (integrated) pressure gradient ∇ p in (2) and G h ∈ R N V ×N p satisfies the discrete compatibility relation between the divergence and gradient operator for all types of boundary conditions.This is a crucial property of the staggered grid discretization that will be used at several points in the construction of our new energy-consistent ROM approach for time-dependent boundary conditions.To solve the system (4)-( 5), we apply the discrete divergence operator M h to the −1 h -premultiplied momentum equation ( 5) and insert the mass equation ( 4) to get the Poisson equation for the pressure, As a result, when integrating the velocity V h (t) according to the momentum equation ( 5), the pressure p h (t) computed from (8) ensures the velocity to satisfy the mass equation ( 4).Note that this approach requires the invertibility of L h .As described in [47], L h is singular for certain boundary conditions and the pressure p h is only defined up to an additional constant.However, by adding a constraint, we can construct an invertible matrix Lh ∈ R N p ×N p such that the unique solution to the regularized Poisson equation is also a solution to the Poisson equation (8).For boundary conditions for which L h is singular, we therefore compute p h via (9) to simulate the FOM (4)- (5).
To avoid the need to distinguish between cases with L h singular and invertible in our notation, we consider (9) the equation to be solved within the FOM simulation, where we set Lh := L h for boundary conditions for which L h is invertible.In both cases, the inverse of Lh satisfies and In other words, within the image space of M h , L−1 h is symmetric and acts as an inverse to L h .Suitable time discretizations of the system consisting of (5) and (9) are, for example, described in [47,48].

Lifting function inspired by Helmholtz-Hodge decomposition
To construct a ROM for the FOM described in Section 2.2, we generalize the ideas in [34,Section 4] to time-dependent boundary conditions.We decompose the FOM velocity, (12) such that the time-dependent lifting function V inhom (t) satisfies the mass equation, (13) and consequently, the homogeneous component V h,hom (t) satisfies the homogeneous mass equation, The choice of the lifting function V inhom (t) is of pivotal importance.Inspired by the Helmholtz-Hodge decomposition [39], we choose the lifting function as This choice has three decisive properties.First, the lifting function indeed satisfies the mass equation with the inhomogeneous right-hand side, (13).This property implies the homogeneous mass equation ( 14) for the homogeneous component V h,hom (t).We use this homogeneous mass equation to construct a velocity-only ROM in the next section.Second, the lifting function is a linear function of the boundary condition vector y bc (t).We use this property to compute the ROM efficiently for arbitrary time-dependent boundary conditions in Section 3.3.Third, as in the Helmholtz-Hodge decomposition, the components in ( 12) are orthogonal, namely with respect to the This property guarantees the energy consistency of the ROM as described in Section 4.

Velocity-only ROM
To approximate the homogeneous velocity V h,hom (t), we construct a POD basis with respect to the h -inner product computed from snapshots defined by the FOM ( 4)-( 5) at time step t j .Due to the homogeneous mass equation ( 14), the resulting as described in [34].Hence, we can construct the low-dimensional approximation V r,hom (t) = hom a hom (t) To determine the coefficient vector a hom (t) ∈ R R hom , we replace V h (t) by V r (t) in the momentum equation ( 5) and project onto hom , Due to the divergence-gradient duality (7) and property (17), the pressure term vanishes, so the ODE (19) simplifies to d dt This ODE does not contain the pressure p h anymore, so our ROM is velocity-only, i.e., the pressure is not needed to simulate the velocity.This property is very useful for two reasons: First because we avoid additional computational costs, second because we avoid pressure-related inf-sup-instabilities [12].The steps taken so far also hold for other choices of lifting functions, e.g., [36,37].For the next simplification of ODE (20), however, we exploit the specific definition of our lifting function, (15), namely the orthogonality property (16).Due to property (17), this orthogonality property is inherited by the POD basis hom , i.e., The initial condition is given by the projection of the FOM initial condition, a hom (0) = T hom V h (0).

Efficient evaluation of the ODE right-hand side
To simulate our novel ROM, we have to discretize the ODE (22) in time and evaluate the right-hand side.As a naive approach, we could first evaluate the FOM-like right-hand side F C D h ( hom a hom (t) + V inhom (t), y bc (t)) and then perform the projection onto hom , both during the online phase of the ROM simulation.However, the computational costs would depend on the FOM dimensions and hence deteriorate the computational complexity of the ROM in the online phase substantially.
Common ways of avoiding such prohibitive costs are hyper-reduction techniques such as the discrete interpolation method (DEIM) [44] and gappy POD [45].Instead of evaluating the complete FOM-dimensional vector F C D h ( hom a hom (t) + V inhom (t), y bc (t)), these methods only require a small number of entries to construct an approximation of this vector.However, the definition of our lifting function V inhom (t) (15) contains the inverse L−1 h , requiring the solution of a Poisson equation with computational complexity O N 2 p .Consequently, sampling a few entries of V inhom (t) is already prohibitively expensive.
Instead of applying hyper-reduction to the right-hand side of ( 22) as a whole, we therefore propose to approximate only the boundary condition vector y bc (t).In detail, we propose to replace the boundary condition vector y bc (t) in the ODE (22) and in the definition of the lifting function (15) with a linear approximation One option to construct this approximation is to apply hyper-reduction techniques such as DEIM [44].In that case, the basis bc would be computed via a POD of snapshots of y bc (t) during the offline phase.The coefficient vector a bc (t) would be computed based on the evaluation of only R bc entries of y bc (t) with complexity O R 2 bc at each time step during the online phase.However, such a DEIM coefficient vector is only an approximation to the optimal coefficient.Actually, the best approximation to y bc (t) within the space spanned by the basis bc is given by bc T bc y bc (t).Hence, we can improve the accuracy of the approximation ỹbc (t) by computing the coefficient vector via the projection onto bc , a bc (t) = T bc y bc (t), instead of performing DEIM.This computation would be prohibitively expensive when performed during the online phase.However, we can maintain the online efficiency by computing a bc (t j ) for all required time steps t j already in the offline phase.This approach requires that y bc (t) is known a priori, which we assume in this work.If y bc (t) would not be known a priori, we suggest to resort to computing a bc (t) via DEIM [44].
The replacement of the boundary condition vector y bc (t) by the approximation ( 23) introduces an error in the ROM.In particular, the resulting ROM V r (t) = hom a hom (t) + Ṽ inhom (t) does not satisfy the mass equation with respect to the original boundary conditions (18), but with respect to the approximated boundary conditions, The coefficient a hom (t) is defined by the ODE d dt (25) where the approximated lifting function Ṽ inhom (t) is defined by The main advantage of this approximation is the drastic reduction of the number of (regularized) Poisson solves required for computing the lifting function due to the linear dependence on y bc (t) in the definition (15).Indeed, to precompute the matrix F inhom ∈ R N V ×R bc , we need to solve only R bc (regularized) Poisson equations for the R bc columns of the matrix Note that this approximated lifting function satisfies the same orthogonality property as the exact lifting function, (21), We will use this orthogonality property in Section 4 to derive energy consistency.

Exact offline decomposition using matricized third-order tensors
To evaluate the right-hand side of the approximated ODE (25) efficiently, we extend the offline precomputation approach in [34].This approach exploits the low-dimensional structure of the terms hom a hom (t), F inhom a bc (t) and bc a bc (t) as well as the fact that the nonlinearities in F C D h are only second order polynomials.In detail, we formulate these second order polynomials as precomputable "matricized" third order tensors.For example, we can write the term quadratic in hom a hom (t) as C2 ( hom a hom (t) V and the Kronecker product ⊗.Then, we find where hom can be precomputed offline.Similarly, we can precompute operators in R R hom ×R 2 hom and R R hom ×(R hom R bc ) for the term quadratic in bc a bc (t) and the mixed term linear in both, hom a hom (t) and bc a bc (t), respectively.Hence, the resulting method has an online complexity scaling in O R

Pressure recovery
As discussed in Section 3.2, our proposed ROM is velocity-only, so no computation of the pressure is needed to simulate the velocity.However, we can compute the pressure p r (t) of the ROM as a post-processing step by inserting the ROM velocity V r (t) and the boundary condition approximation ỹbc (t) into the (regularized) FOM Poisson equation ( 9) This pressure is consistent in the sense that the (regularized) Poisson equation ( 29) has the same structure as the (regularized) FOM Poisson equation (8).Consequently, the ROM pressure p r (t) converges to the FOM pressure p h (t) if the ROM velocity V r (t) converges to the FOM velocity V h (t) and the boundary condition approximation ỹbc (t) converges to the exact boundary condition vector y bc (t).

Kinetic energy evolution
A distinctive property of our proposed ROM is its kinetic energy consistency, i.e., the ROM kinetic energy evolution has the same structure as the FOM kinetic energy evolution.We define the FOM kinetic energy as , leading to the kinetic energy evolution d dt using the momentum equation ( 5) and the mass equation ( 4).
The key property leading to a kinetic energy evolution with the same structure for our proposed ROM is the orthogonality of V r,hom (t) and Ṽ inhom (t), given by equation (27).Due to this orthogonality the ROM kinetic energy can be decomposed into two terms, The time-derivative of the first term (the energy of the homogeneous component) is given by d dt 1 2 using the momentum equation (25).Furthermore, the time-derivative of the second term (the energy of the approximated lifting function (15)) can be written as where we use the divergence-gradient duality (7), the symmetry of L h and h , the definition of L h (8), and the inverse and Ṽ inhom (t) Here, we have again used the divergence-gradient duality (7), the symmetry of L h and h , the symmetry property of L−1 h (11), and identified the definition of the lifting function (15), Altogether, the kinetic energy evolution of the ROM can be written as This ROM energy evolution has exactly the same structure as the FOM energy evolution (30), which is our definition of consistency.This energy consistency will be useful to transfer stability statements from the FOM to the ROM.Namely, if we assume that a bound on the kinetic energy can be derived from the FOM kinetic energy evolution with respect to certain approximated boundary conditions ỹbc (t), then the same derivation could be used to obtain a bound for the ROM kinetic energy.Note, however, that stability (energy) bounds for the incompressible Navier-Stokes equations including time-dependent boundary conditions are still an open topic [49,50], so this is a suggestion for future research.
As indicated above, the orthogonality property (27) is key to achieve this energy consistency.In Appendix A, we show that a different choice of the lifting function would spoil the energy consistency so that FOM energy bounds could not directly be transferred to the ROM.
To preserve the energy-consistency on the time-discrete level, one could use energy-conserving Runge-Kutta methods [48].However, in our setting of time-dependent boundary conditions, the kinetic energy is not a conserved quantity and energy-conserving methods do not (and should not) yield a constant energy.Instead, we opt for using high-order explicit methods, which are an accurate and computationally attractive alternative, also given the fact that we are using an exact third-order tensor representation of the convective terms.As shown in Appendix B, the difference between explicitly and implicitly integrated ROM velocities is quite small.

Velocity-pressure ROM
Our proposed ROM has the appealing property that it is velocity-only.However, we can interpret this ROM as a velocitypressure ROM, as well.This interpretation provides further insight in the relation between time-dependent boundary conditions and the role of the pressure.
The velocity-pressure ROM that we consider comprises the velocity approximation V r (t) = a(t) ≈ V h (t) with the horthonormal basis ∈ R N V ×R V and the coefficient vector a(t) ∈ R R V , and the pressure approximation p r (t) = b(t) ≈ p h (t) with the basis ∈ R N p ×R p and the coefficient vector b(t) ∈ R R p .In the next section, we will specify the choice of bases such that the resulting velocity-pressure ROM is equivalent to the velocity-only ROM proposed in Section 3. In general, however, the bases can also be chosen differently.
We insert these approximations in the mass equation ( 4) projected by the pressure basis and in the momentum equation (5) projected by the velocity basis , leading to: The initial condition a(t 0 ) = a 0 is defined via the constrained optimization problem of minimizing a 0 − V h (t 0 ) h such that a 0 satisfies the projected mass equation (39) at t = t 0 .Analogously to the FOM described in Section 2, we solve this system by solving the reduced Poisson equation for b(t), Then, due to the pressure term in the ODE (40), a(t) computed from this ODE is guaranteed to satisfy the projected mass equation (39).In contrast to the FOM, the system matrix L r of this reduced Poisson equation is not high-dimensional and sparse, but low-dimensional and dense.Hence, the system is solved (by employing a precomputed LU-factorization) with a computational complexity of O R 2 p at each time step in the online phase.This approach requires the invertibility of the system matrix L r = T M h T G h .We will guarantee this invertibility by choosing the bases and suitably, as explained in the next section.
Analogously to the discussion in Section 3.3 for the velocity-only ROM, computing the right-hand side of the ODE (40) for the exact boundary condition vector y bc (t) would in general be prohibitively expensive.Hence, we have replaced the exact boundary condition vector y bc (t) in ( 39)-( 41) by the lower-dimensional approximation ỹbc (t) := bc a bc (t) as described in Section 3.3.

Equivalence to the proposed velocity-only ROM
The following proposition states that the velocity-pressure ROM (39)-( 40) is equivalent to the velocity-only ROM ( 25)-( 26) if we choose the bases and appropriately.

Proposition 5.1 (Equivalence of velocity-only and velocity-pressure ROM).
If the h -orthonormal velocity basis ∈ R N V ×R V and the pressure basis ∈ R N p ×R p are chosen such that a) the image space of is equivalent to the union of the image spaces of hom in equation (17) and F inhom as defined in equation (26).
b) for all t ∈ [0, T ], there exists a solution to the unprojected mass equation then the velocity of the velocity-pressure ROM, V r,velocity−pressure (t) = a(t), is for all t ∈ [0, T ] equivalent to the velocity of the velocity-only ROM, V r,velocity−only (t) = hom a hom (t) + Ṽ inhom (t) described in Section 3.
Proof.First of all, Assumption (a) implies that there exists an orthonormal matrix ) is an h -orthonormal basis with the same image space as F inhom which can be computed from F inhom , e.g., via a Gram-Schmidt-orthonormalization.We can decompose Then, we write the velocity-pressure ROM velocity as We proceed to show that the velocity components of the two ROMs are equivalent, namely a 1 (t) = a hom (t) for all t ∈ [0, T ], (44) and inhom a 2 (t) = Ṽ inhom (t) for all t ∈ [0, T ]. (45) For this purpose, we note that the decomposition (43) separates the modes in hom , which lie in the kernel of M h , from those in inhom , which lie outside the kernel of M h .Consequently, equations ( 39) and ( 42) simplify to T M h inhom a 2 (t) = T ỹM (t), (46) and M h inhom a 2 (t) = ỹM (t), (47) respectively.Let, for a fixed t, â2 and ā2 be the solutions to (46) and (47), respectively (their existence follows from Assumption (b)).Then, premultiplying (47) by T , we find By construction of inhom , we have rank 48) implies ā2 = â2 .Hence, the solution â2 to the projected mass equation ( 46) also satisfies the unprojected mass equation (47).This solution â2 to the unprojected mass equation is unique.To show this uniqueness, let ǎ2 be a second solution to (47).Then, we find

and because of Assumption (d), rank(
Since M h inhom has full column rank, we infer â2 − ǎ2 = 0, so the two solutions are equivalent.This uniqueness implies the equivalence of the inhomogeneous velocity components (45) because F inhom has the same image space as inhom and both terms satisfy the unprojected mass equation for all t ∈ [0, T ].
What remains to be shown is the equivalence of the homogeneous velocity components (44).For this purpose, we premultiply the momentum equation (40) by (50) where Q 1 T = T hom and the pressure term T hom G h b(t) disappears due to equation (17).This ODE is equivalent to the velocity-only ODE (25).As both ODEs have the same initial condition, they also have the same solution, so we find (44).
Altogether, we hence find The requirements of the proposition are met, for example, by the choice = [ hom inhom ] and = M h inhom .For this choice, we have Q = I.Note that this choice also guarantees the invertibility of L r .
The proof of the proposition provides an explanation why eliminating the pressure in the velocity-only ROM is possible.As can be seen from ODEs (40) and ( 50), the pressure term does not affect the coefficients of the homogeneous modes (a 1 (t)) but only those of the inhomogeneous modes (a 2 (t)).Since these inhomogeneous coefficients are uniquely determined by the mass equation (39), they do not require simulating ODE (40), so the pressure is superfluous.This observation is in line with the common understanding that the pressure is a (vector of) Lagrange multiplier(s) enforcing the velocity to satisfy the mass equation.This understanding also implies that the pressure does not affect the velocity beyond enforcing the mass equation.Indeed, the error between the approximated pressure p r (t) = b(t) and the FOM pressure p h (t) does not introduce an additional error in the ROM velocity.
In practice, the velocity-pressure formulation is slightly more expensive than the velocity-only formulation because the ODE system is integrated in time for R V = R hom + R inhom instead of only R hom coefficients and additionally the pressure coefficients have to be computed by solving the reduced Poisson equation ( 41) at every time step.Hence, the ROM described in Section 3 remains the method of choice.
Analogously to [51], we can decompose the velocity-pressure ROM velocity into a homogeneous and an inhomogeneous component using a Q R decomposition for arbitrary choices of and provided that L r is invertible.Provided the discretization satisfies the gradient-divergence duality (7), this decomposition can be formulated as a velocity-only ROM.Due to the stability problems of the QR decomposition described in [51], this velocity-only ROM would in general not be useful in practice.

Comparison with standard POD velocity-pressure ROMs
A common choice for the velocity basis in the velocity-pressure ROM defined by ( 39) and ( 40) is the ordinary POD basis.The resulting ROM does not in general meet the requirements of Proposition 5.1 because the POD basis does not necessarily satisfy Assumption (b).Consequently, this standard POD velocity-pressure ROM is in general not equivalent to the proposed velocity-only ROM.In fact, such a standard POD velocity-pressure ROM can suffer from inf-sup instabilities [6] whereas the proposed velocity-only ROM cannot.
On the other hand, the important insight that Proposition 5.1 does give, is that the proposed velocity-only ROM is equivalent to a velocity-pressure ROM, which differs from the standard POD velocity-pressure ROM only in the choice of velocity basis (and possibly pressure basis ).To compare these differing velocity bases, let us define the best approximation error over the snapshots Then, the standard POD bases with R V = R hom and assuming exact boundary conditions ỹbc (t) = y bc (t).The proof can be found in Appendix C.
If the ROM errors are dominated by these best approximation errors, we can consequently expect the accuracy of the proposed velocity-only ROM to lie in between the standard POD velocity-pressure ROMs with R V = R hom + R inhom and with R V = R hom modes, provided that the standard POD velocity-pressure ROMs are stable.
As explained in the previous section, velocity-pressure ROMs with R V = R hom + R inhom modes are slightly more expensive than the velocity-only ROM with R hom and R bc modes.Also, velocity-pressure ROMs with R V = R hom modes are slightly more expensive than the velocity-only ROM due to the need to solve Poisson equations for the pressure.
The accuracy per computation time of the proposed velocity-only ROM is consequently similar to standard POD velocitypressure ROMs.While the latter are known to be prone to instabilities, the velocity-only ROM is energy-consistent and does not suffer from inf-sup instability.Furthermore, while standard POD velocity-pressure ROMs only satisfy a projected mass equation, the velocity-only ROM satisfies the unprojected mass equation with respect to explicitly approximated boundary conditions.

Numerical experiments
In our numerical experiments, we investigate our proposed velocity-only ROM regarding convergence behaviour, mass conservation and energy consistency.Furthermore, we verify the equivalence to the velocity-pressure ROM described in Section 5.2 and analyse the runtime.
with p ∞ = 0 and ν = 10 −2 m 2 s .In contrast to [48, Section 6.4.2], we consider inflow on the left boundary (x = 0) which is not only time-but also space-dependent.For the first testcase varying angle, we consider inflow with constant magnitude 1 but time-and space-dependent angle, As we will see in Section 6.1.3,these boundary conditions can be approximated well by a lower dimensional approximation ỹbc (t) as described in Section 3.3.In contrast, the second testcase exemplifies the extreme and in practice unlikely case that the boundary conditions cannot be reduced effectively.For this second testcase moving mode, we move the parabolic inflow field u parabolic (y) = 1 10 (y − y min )(y max − y) uniformly such that it enters the simulation domain from the bottom at t start and leaves it at the top at t end , This inflow does not have a tangential component, so v b = 0.For both testcases, the Poisson matrix L h is invertible, so Lh = L h in the FOM Poisson equation ( 8) and in the definition of the lifting function (26).
The actuator disk represents a momentum sink and is modelled as a constant force f = 0.25 acting in negative xdirection and we consider the time intervals [0, 4π ] and [0, 20], respectively.

Numerical methods
We discretize the simulation domain with a uniform 200 × 80 staggered grid.For time integration, we use the explicit fourth-order Runge-Kutta method RK4 [52] and discretize for both testcases the respective time interval equidistantly with time step size t = t end 800 .The initial condition is given by V 0 h = V 0 inhom which is computed according to (15) based on the respective boundary conditions.For the moving mode testcase, this definition results in V 0 h = 0 while the initial condition for the varying angle testcase is a nonzero function.
The POD bases hom and bc are constructed from the snapshot matrices and X bc = y bc t 0 . . .y bc t j . . .
where V j h,hom is the time discrete approximation of V h,hom t j with t j = j • t.We compute the coefficient vector a bc (t) = bc T bc y bc (t) during the offline precomputation phase for all required time steps, as described in Section 3.3.
The numbers of ROM modes R hom and R bc can be chosen independently of each other.In the numerical experiments here, however, we choose R hom = R bc = R for convenience only, and consider velocity-only ROMs for different values of R.

Singular value decays
First, we investigate the singular value decays of the snapshot matrices which we construct our ROM bases from.As can be seen in Fig. 3a, only 23 singular values of the boundary condition vector snapshot matrix X bc have a significant magnitude for the varying angle testcase.In fact, only 5 modes are required to capture at least 99% of the snapshot energy [26].
In contrast, Fig. 3b shows that 80 singular values of X bc for the moving mode testcase, which we designed to have irreducible boundary conditions, have a relative magnitude larger than 10 −4 .This number equals the number of grid points at the inflow boundary and the number of nonzero entries in y M , so these entries do not exhibit any lower-dimensional coherent structure.Here, 29 modes are required to capture at least 99% of the snapshot energy.
Similarly, the decays of singular values of the velocity snapshot matrices X hom shown in Fig. 4 differ considerably for large mode numbers.While only 164 singular values are larger than machine precision (when divided by the largest singular value), for the varying angle testcase, 684 singular values are larger than this relative threshold for the moving mode testcase.For small mode numbers, however, both decays are similarly fast, and only 13 and 14 modes capture at least 99% of the snapshot energy, respectively.

Boundary condition approximation accuracy
Before investigating the proposed ROM itself, we investigate the accuracy of the boundary condition approximation for the two test cases.As can be seen in Fig. 5, the best approximations are accurate up to machine precision for R bc ≥ 23 and R bc ≥ 80, respectively.This observation is perfectly in line with Fig. 3 that shows 23 and 80 significant singular values.
For the approximations via DEIM, we observe accuracy up to machine precision for R bc = 23 and R bc = 80, but large errors if the numbers of modes are increased further.These large errors are likely due to the fact that these numbers of modes exceed the rank of the corresponding snapshot matrix, causing high sensitivity to machine precision errors.When using the DEIM approach, it is consequently important to know the rank of the snapshot matrix X bc .In the following, we set R bc = min(R, 23) for the test case varying angle and R bc = R for the test case moving mode.

Convergence
First, we investigate whether the ROM converges to the FOM for increasing number of modes.As can be seen in Fig. 6 the relative velocity error clearly decreases monotonously for increasing number of ROM modes.Consequently, our velocity-only ROM for time-dependent boundary conditions indeed provides a convergent approximation to the FOM for both boundary condition approximation approaches.

Mass conservation
Second, we investigate to which extent our proposed ROM violates the original mass equation (18).As explained in Section 3.3, the proposed ROM satisfies the mass equation exactly, but with respect to the approximated boundary conditions.
The difference between the two is shown in Fig. 7, effectively indicating how well the mass equation right-hand side y M (t) is approximated by ỹM (t).
Indeed, as can be seen in Fig. 7 the mass conservation violation decreases for increasing number of boundary condition approximation modes.For the testcase varying angle, the mass conservation violation is in the order of machine precision for R ≥ 20.In contrast, for the moving mode testcase, R = 80 modes are required to achieve a similar mass conservation accuracy.Even for R = 40 modes, the mass conservation violation is larger than 10 −4 .This observation can directly be explained by the singular value decay depicted in Fig. 3b, which shows that the 80 largest singular values all have a relative magnitude larger than 10 −4 .Note however, that in this case the boundary conditions were chosen on purpose to be highly 'irreducible' with POD and as such this can be seen as a limiting case, which is very unlikely in practice.Furthermore, the violation of the exact mass equation seems not to impair the accuracy of the actual solution V r , as can be observed by comparing Fig. 6a to 6b.

Energy consistency
Third, we investigate whether the structural correspondence of the FOM and ROM kinetic energy evolutions derived in Section 4 is accompanied by small energy errors in practice.
for the two testcases for boundary condition DEIM approximations (solid lines) and best approximations (dashed lines).Indeed, Fig. 8 shows that the relative kinetic energy errors are approximately at the same order of magnitude as the corresponding relative velocity errors depicted in Fig. 6.

Verification of ROM equivalence
Next, we verify the ROM equivalence demonstrated in Section 5.2.For this purpose, we simulate the same testcases as in the previous experiment but now with the velocity-pressure ROM (39)- (40) with the bases = [ hom inhom ] and = M h inhom as described in Section 5.2.Here, we compute inhom via a QR decomposition of F inhom .Hence, the numbers of modes of the velocity-pressure ROM are given by R As can be seen in Fig. 9, the error between the two ROMs ranges for all values of R between 10 −18 and 10 −10 .Hence, the ROM equivalence is confirmed up to machine precision.

Runtime analysis
Fig. 10 shows the computation time of the ROM in comparison with the FOM.For all considered numbers of we observe a significant speed-up of the ROM online phase compared to the FOM simulations.As in [34], this speed-up amounts to approximately two orders of magnitude for moderate numbers of modes.The offline precomputations, on the other hand, have significant costs and even exceed the FOM simulation costs for R = 80.However, when used in practice for a parametric study, the ROM would already pay off for a small number of ROM simulations.

Conclusion
In this paper we have generalized the energy-conserving ROM for the incompressible Navier-Stokes equations proposed by Sanderse [34] to arbitrary time-dependent boundary conditions.Our ROM is velocity-only, which avoids computing the pressure and pressure-related inf-sup instability, and it preserves the structure of the kinetic energy evolution.The ROM is based on a Helmholtz-Hodge-inspired decomposition that separates the velocity into two orthogonal components: i) a velocity that satisfies the mass equation with homogeneous boundary conditions and ii) a lifting function that satisfies the mass equation with the prescribed time-dependent boundary conditions.The homogeneous component is defined by the Galerkin-projected momentum equation.Due to the divergence-free projection basis, the pressure term disappears form the momentum equation resulting in the proposed velocity-only ROM.
To allow efficient simulation of the ROM, we approximate the prescribed boundary conditions by a lower-dimensional approximation via POD.The resulting ROM satisfies the mass equation with respect to these approximated boundary conditions.Our numerical experiments show that the error introduced by this boundary condition approximation does not dominate the velocity error.In these experiments, we use as many modes to approximate the boundary conditions as we use to approximate the homogeneous velocity field.In future research, different choices of these mode numbers should be investigated as well.
We have theoretically shown that the orthogonality of the components leads to a ROM kinetic energy evolution that has the same structure as the FOM kinetic energy evolution.Our numerical experiments confirm that the proposed ROM is able to accurately mimic the kinetic energy evolution.In future work, the structural equivalence of the kinetic energy evolution can be used to derive bounds on the kinetic energy and hence to derive nonlinear stability results.
Furthermore, we have proven the important result that the proposed velocity-only ROM is equivalent to a velocitypressure ROM, when employing a specific choice of basis functions, and confirmed this equivalence in our numerical experiments.This equivalence can function as a key concept to understand the role of the pressure in ROMs for the incompressible Navier-Stokes equations.It gives important insights into the behaviour of the velocity-pressure ROM regarding the choice for the number of velocity, pressure and boundary condition modes.While this ROM equivalence is based on the assumption of the discrete gradient-divergence duality, future research should investigate whether discretizations that do not exhibit this property also allow for a similar equivalence statement.
In numerical experiments, we have shown the convergence of the ROM upon increasing the number of ROM modes and the accuracy of the boundary condition approximation.Instead of performing a hyper-reduction technique for the convective term, we have used an exact offline decomposition.Scaling cubically in the number of modes, this offline decomposition is limited to small numbers of modes.For moderate numbers of modes, the speed-up of the ROM online phase compared to the FOM is approximately two orders of magnitude.If larger numbers of modes are required that prohibit this cubic scaling, we suggest to incorporate energy-conserving hyper-reduction methods [53].
In future work, it should be investigated how the proposed method performs in practically relevant testcases involving higher Reynolds numbers and generalization over parameter spaces, varying boundary and initial conditions as well as longer time intervals.In case of higher Reynolds numbers, it might be useful to combine our approach with methods that improve the representation of dissipative small scales [10,11].For such combinations, we suggest to keep the lifting function unchanged and only adapt the construction of the homogeneous velocity approximation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix C. Best approximation error relations between the proposed velocity-only ROM and a standard POD velocity-pressure ROM
We repeat and prove the statement in Section 5.3 about best approximation errors of standard POD bases and the proposed velocity-only ROM basis.

Proposition C.1 (Best approximation errors of velocity-only ROM and standard POD velocity-pressure ROMs).
Let, as in (52), (C.2) with R inhom , hom ∈ R N V ×R hom and inhom ∈ R N V ×R inhom as used in the proof of Proposition 5.1, and the assumption that ỹbc (t) = y bc (t) for all t ∈ [0, T ].
Proof.We use the well-known optimality property of POD bases [4] that the h -orthonormal POD basis ∈ R N V ×R V minimizes ε(φ) over all h -orthonormal φ ∈ R N V ×R V .Since

( 21 )
Due to this orthogonality, the left-hand side of(20) simplifies to d dt a hom , so the ODE for the ROM coefficient a hom (t) becomes d dt

Fig. 2 .
Fig. 2. Sketch of testcase setup with actuator disk and time-and space-dependent boundary conditions.

Fig. 3 .
Fig. 3. Singular values of the snapshot matrices X bc for the two testcases divided by largest singular value.

Fig. 4 .
Fig. 4. Singular values of the snapshot matrices X hom for the two testcases.

Fig. 5 .
Fig. 5. Relative boundary condition vector approximation errorỹbc (t j )− y bc (t j ) h average j y bc (t j ) h for the two test cases for boundary condition DEIM approximations (solid

Fig. 7 .
Fig. 7. Mass conservation violation Mh V j r − y M (t j ) 2for the two testcases for boundary condition DEIM approximations (solid lines) and best approxima-

Fig. 8 .
Fig. 8. Relative kinetic energy errorK j h −K j r average j K j hfor both testcases for boundary condition DEIM approximations (solid lines) and best approximations

Fig. 9 .
Fig. 9. ROM equivalence error V j r,velocity−only − V j r,velocity−pressure h

Fig. B. 11 .Fig. B. 12 .
Fig. B.11. h -norm of velocity difference, V j r,explict − V j r,implicit h , of ROMs integrated in time with explicit RK4 and with implicit midpoint for the two testcases with boundary condition best approximation and R hom = R bc = R.

1 )
the best approximation error over a set of snapshotsV j h K j=1.Then, the h -orthonormal POD bases over the set of snapshotsV j h K j=1 with R V = R hom and R V = R hom + R inhom modes, [R hom ] and [R hom +R inhom ] , satisfy relation (53), Staggered grid.Velocity values are located on the faces of the pressure finite volumes.