Reduced order models for finite-volume simulations of turbulent flow around wind-turbine blades.

The computational cost of the design optimisation of wind turbines, as well as the optimisation of the operation and maintenance of offshore wind farms represents a limitation to the application of conventional simulation methods. New techniques such as reduced order modelling (ROM) and technologies such as hybrid-analysis methods and digital twins have increased in popularity due to their ability to deliver numerical results at a significant speed-up with reasonable accuracy. This work presents a hybrid projection-based proper orthogonal decomposition (POD) strategy applied to transient turbulent flow problems. Key feature of this work is the applicability of the methodology to high Reynolds number cases and the stabilisation of pressure in the online phase via the assembly of the so-called pressure Poisson equation. Another significant part of this work is the implementation of an interpolation scheme for the eddy viscosity field in the classical POD-Galerkin strategy. The sampling procedure and the calculation of the reduced operators in the offline phase is carried out using the finite volume method (FVM), OpenFOAM’s libraries specifically, while the construction of the reduced basis and the solution of the online phase is carried out in Python. The capability of the resulting ROM is tested using the two-dimensional flow around a NACA0015 airfoil at 17° angle of attack with Reynolds number of approximately 300 000.


Introduction
The complexity of offshore wind aerodynamics, with varying external flow fields, platform movement effects and upstream wakes requires the study of a large number of cases, making the use of conventional (high-fidelity) methods prohibitive. On the other hand, reduced order modelling (ROM) techniques provide an excellent alternative. ROM techniques aim to build a simplified model that can efficiently approximate the behaviour of the high-fidelity one with reasonable accuracy. A ROM generally consists of two distinct phases. The offline phase describes the process of constructing the simplified model, a typically computationally intensive procedure that is, however, carried out just once, while the online phase describes the solution, or particularisation of the simplified model for given parameters. That could be the solution of an algebraic problem or simply an interpolation, and can often be run in real-time.
This work focuses on what is known as Projection-Based ROM. Such ROM techniques have been applied in several scientific contributions dealing with laminar fluid-dynamics problems [1][2][3][4][5]. In [3] and [4] stabilisation of the online phase via the supremiser enrichment  [6] can enable faster solution of the linear-algebra system in the online stage.
In contrast, applications on turbulent flows are limited due to issues with regard to the stability of the methodology. ROM techniques generally aim to capture the most significant information, i.e. high energy modes. However, smaller scales are responsible for the dissipation of the turbulent kinetic energy. Various strategies have been proposed to create more dissipative ROMs. One option would be the introduction of a turbulence model in the ROM stage. The work in [7] offers an excellent introduction and contains a useful comparison of such strategies applied in POD-Galerkin ROMs. However, the option that is presented in this work is the use of a data-driven ROM, in this case POD-I, for the eddy viscosity as proposed in [8], giving rise to a combined projection and data-driven POD, a hybrid, which allows for application to turbulent problems, while being independent of the turbulence model used in the offline stage.
Aiming to build a methodology usable in industrial environments, we consider the most commonly employed methodologies, i.e. the Finite Volume Method (FVM) [9][10][11] and the Reynolds-Averaged Navier-Stokes (RANS) method for the resolution of turbulent effects. Most offline computations are carried out employing OpenFOAM's open-source libraries. In [12], the FVM solution is mapped onto a conforming finite-element domain before constructing the ROM and solving the online phase. In order to expand this work to turbulent flows, we use the FVM framework for both ROM phases. Moreover the mixed-ROM for turbulent flows proposed in [8] employs the supremizer enrichment. In this work, striving for consistency with the Full-Order Model (FOM), we adapt the pressure Poisson equation to include turbulent effects.
Section 2 introduces the Full-Order Model, viz. the classical incompressible Reynolds-Averaged Navier-Stokes equations. It also briefly describes the Finite Volume Method as implemented in OpenFOAM. Section 3 introduces the affine representation of the approximation and the Reduced Order Model as applied in the FOM of interest. Finally, section 4 uses a numerical experiment to test the viability of the proposed methodology to problems of moderately high Reynolds number. More precisely, we solve the two-dimensional flow around a NACA0015 airfoil at 17 • angle of attack with Reynolds number of approximately 300 000.

The Reynolds-averaged Navier-Stokes equations
For the simulation of the incompressible flow by means of RANS equations, the coupled variables (u, p) are decomposed into their mean flow component, (u, p), from now on (u, p), and their respective perturbations, (u , p ). Given an open-bounded domain Ω with boundaries such that ∂Ω = Γ W ∪ Γ In ∪ Γ Out , where the three distinct boundaries describe the walls, inlet and outlet surfaces, respectively, the transient RANS equations for the mean flow components read Where ν, ν t are the physical and turbulent viscosity, respectively, f , u d and u 0 are the external body forces, the imposed velocity datum on the inlet and the prescribed velocity initial condition. We assume that Γ Out corresponds to a free traction outlet with, accordingly, homogeneous traction. It should be noted that u · n should be non-negative on Γ Out for stability. The eddy viscosity is computed using a turbulence model of choice. However, in accordance with the fact 3 that the proposed methodology is independent of the closure model, no turbulence model is specified in (1).

The finite volume method
Given a domain, physical properties, boundary conditions and initial conditions, the Navier-Stokes equations are approximated using the finite-volume method as implemented in OpenFOAM [13]. Here, for the sake of completeness, the finite-volume rationale is briefly recalled. For more information the reader may refer to [13]. The computational domain Ω is where Ω h defines the discretised spatial domain, dependent on the characteristic mesh size h. The Navier-Stokes equations (1) are integrated over each cell where (u, p) are conceived of as piecewise constant functions, represented by degrees of freedom (DOFs) that are located at cell centres.
The system is solved using the PIMPLE algorithm, an iterative solver that combines the functionality of the SIMPLE [14] and PISO [15] algorithms. Employing Gauss's theorem, the integrals of (2) containing spatial differential operators are rewritten in terms of fluxes over the respective cell faces. The transition from cell centres to cell faces is done via linear interpolation and the approximation of the spatial operators is typically carried out using central differencing schemes. Moreover, the convection term in (2) is linearised using a half-step velocity value. The temporal term is discretised using the implicit Euler scheme, which is the default scheme in early stage simulations. For a more detailed view of the discretisation schemes as implemented in OpenFOAM, the interested reader is referred to [13]. Remark 1. The implicit Euler scheme is the default approach for initial stage simulations in view of its stability and ease of use. However, it is important to note that is is only first order accurate. Other options such as the Crank-Nicolson and 2 nd order backward differencing scheme generally provide better accuracy. Because the scope of this paper pertains to the ROM strategy, more than the underlying high-fidelity model for the Navier-Stokes equations, we will also adopt this baseline time-discretisation scheme.

The Reduced Order Model
Using reduced order modelling, we aim to decompose the transient solution vector (u , p) into a set of (deterministic) functions that each capture part of the kinetic energy. More precisely, using POD an approximation of said coupled variables is sought such that where M is the number of terms required to capture enough of the kinetic energy, u r i , p r i are spatial functions (or modes) for velocity and pressure and α i , β i are their respective time coefficients. The fundamental assumption behind all model reduction techniques, including the one presented in this work, is that M N without essentially compromising the accuracy of the approximation. Effectively, the full-order solution manifold can be accurately approximated by a reduced space whose dimension is much lower, thus allowing the rapid, even real-time, evaluation of complex systems. There exist several techniques for the evaluation of the reduced basis such as Proper Orthogonal Decomposition (POD), Reduced Basis (RB) methods and Proper Generalised Decomposition (PGD). Interested readers are referred to the above methodologies and some applications to computational fluid dynamics in Refs. [16][17][18][19][20][21][22].

Construction of the reduced basis via POD-Galerkin
The spatial functions in (3) are computed using the standard Proper Orthogonal Decomposition strategy. POD proceeds by sampling the FOM in time and generating an ensemble of solutions (snapshots) of (1). The following description uses the vector field of velocity, but it stands similarly for the scalar field of pressure: where N is the number of snapshots. At this stage a covariance matrix is constructed where · , · corresponds to the L 2 inner product.
, the M basis functions are extracted as: where is, by construction, orthonormal and optimal (see [23], Proposition 6.1). Moreover, the approximative power of the reduced model depends on the square of the relative amplitude of the truncated modes as follows: where ε > 0 is the average lower error bound of interest. Following the same procedure for pressure, we end up with the following two orthonormal bases where it should be noted that the number of spatial modes for velocity and pressure (M u , M p ) is not necessarily equal. An approach to assembling a system in the online phase that satisfies the inf-sup condition is to retain fewer modes for pressure than for velocity (M u > M p ). However, in our experience using the proposed methodology in the numerical experiment below, retaining the same amount of modes for both pressure and velocity does not seem to cause instability. Section 3.2 goes into further detail on pressure stabilisation. At this stage, all that remains is for the time coefficients to be calculated. Considering, for the sake of simplicity, the laminar problem and assuming, without loss of generality, vanishing external body forces and substituting (3) in (1) (8), which gives rise to the following algebraic system: where are the vectors of unknown time coefficients for the time step t and the rest of the terms are defined as follows: The above matrices are expressions of V r and Q r and can be pre-computed in the offline stage. It should be stressed that the quadratically non-linear convective form is represented by a trilinear operator whose size increases cubically in the number of basis functions, M u . Therefore, its computation and storage might become cumbersome for large systems.

The Pressure Poisson Equation
In the FVM we use a Chorin-type splitting, i.e. the so-called pressure Poisson equation for the stabilisation of the full-order model. However, the reduced basis V r , Q r does not necessarily inherit the stability of the FOM solution. Therefore, on the reduced-order level, we must separately ensure that it is. Aiming to stay as consistent as possible to the underlying fullorder model, for the calculation of the pressure time coefficients, β, it is chosen to also assemble the pressure Poisson equation, resulting in a well-posed system. The pressure Poisson equation is effectively a reformulation of the incompressibility constraint. Initially proposed by Chorin [24] and Temam [25], it is an approach that has been studied extensively in literature [26][27][28]. Taking the divergence of the momentum equation and exploiting the divergence free velocity modes, one can derive the following equation where the Neumann boundary condition is determined by projecting the momentum equation onto the boundary faces. Substituting the POD approximation into (11) and carrying out a Galerkin projection onto the reduced basis spaces, we get the following algebraic problem for the pressure coefficients β j : where the coefficient vector α(t) is known from the solution of the momentum equation, (9a).
The new matrices are defined as follows: Finally, for the laminar Navier-Stokes equations, the online stage can be summarised as solving the algebraic system consisting of the momentum equation (9a), and the pressure Poisson equation (12): The system is solved by means of a segregated iterative solver. Initially, the momentum equation is solved for α(t), followed by the solution of the pressure Poisson equation for β(t). For the solution of the momentum equation, the term associated to the trilinear form is linearised, B r ik = B r ijk α l−1 j , where the superscript l − 1 describes the previous iteration.

Hybrid POD-Galerkin strategy for turbulent flows
Similar to the full-order model, when tackling flow problems of high Reynolds number, it is required to resolve turbulent effects. In the RANS approach a diffusive term is introduced in the momentum equation and an auxiliary turbulence model is solved for the eddy viscosity.
In this work a POD-I element is considered. Similar to (3), an approximation of the eddy viscosity is sought as follows: where M ν is the number of modes, ν r t i are the spatial basis functions and γ i are their time coefficients. Following the strategy described in Sec. 3.1, we construct the reduced basis Y r and evaluate the time dependence matrix via projection of the solution matrix, , on the reduced basis The j th row of Ξ contains the coefficients γ i of the ROM approximation (15) at instance t j . In the online phase, for any valid particularisation of time, the problem consists of calculating the time coefficients for eddy viscosity, γ i (t). In the scope of this work a piecewise linear interpolation is employed, as follows: where χ j (t − t j ) are piecewise linear interpolation functions centred around t j and supported on (t j−1 , t j+1 ) with range [0 , 1] and g j i are weights. Owing to the fact that we use a piecewise linear interpolation g j i = γ i (t j ). The turbulent contribution in (1) is expressed in the diffusive term. Focusing on the turbulent part, substituting the POD approximations gives rise to the following term Carrying out a Galerkin projection on the reduced basis V r , we derive the algebraic form where CT r ijk = u r i , ∇ · (ν r t j ∇ s u r k ) . The enriched momentum equation (9a) that includes the turbulent term reads: Similar contributions are derived for the pressure Poisson equation and its boundary conditions. Solving (21) and the corresponding enriched Poisson equation for pressure in the online phase returns the time coeficients (α, β) for turbulent flow problems. For a more detailed description of the proposed strategy, the interested reader is referred to [29].

Numerical experiment
In this section the proposed methodology is tested on the two-dimensional NACA 0015 airfoil suspended in a mesh. Aiming to capture the transient flow characteristics the airfoil is rotated to 17 • angle of attack, leading to an unsteady flow topology. The Reynolds number is chosen at Re ≈ 320 000, so that it is comparable to experiments run in smaller wind tunnels.
The C-domain is discretised using 240 000 cells. The inlet extends 20c, where c = 1 m is the chord length of the airfoil, while the outlet is 30c downstream and the total height of the domain is 40c. The freestream velocity datum (u d = 5 m/s) is imposed on the inlet, top and bottom boundaries, while a free traction is imposed on the outlet. Finally, a no-slip condition is imposed on the airfoil walls. After the flow is developed (t ≈ 0.5 s), we observe an oscillatory behaviour, which repeats itself with approximate period 2.5 s. Owing to the observed periodicity, a time range from t 0 = 1 s to T = 3.5 s is examined. The time step is set to 0.005 s.
The results of the Hybrid POD-Galerkin strategy will be compared to the full-order results, measuring the relative L 2 -error for both velocity and pressure over time. The aforementioned errors are evaluated on a control box that covers the airfoil and extends 1 chord length upstream, 5 chord lengths downstream and 1 chord length both above and below the airfoil (2 chord lengths total height). The other quantities of interest for this case are the lift and drag coefficients C l and C d , respectively.

Offline phase
We start by collecting a set of full-order snapshots, storing a solution on equal intervals. From the resulting full-order solution matrix for velocity, pressure and eddy viscosity we construct the spatial basis functions and pre-compute the reduced operators. The spatial modes (Fig.1) represent significant flow structures.
The number of modes to be retained, M u , M p and M t , is determined by the relative amplitude of the truncated terms (Fig.2). Using Eq. (7) and considering a lower average error bound of 1%, we require approximately 30, 60 and 90 modes for pressure, velocity and eddy viscosity, respectively.

Online results
The online phase consists of calculating the time-coefficients such that the solution fields can be approximated using (3). We define the time-range and time-step of interest. Note that for the  proposed strategy, the time-range needs to be contained in that of the FOM. Figure 3 displays a comparison between the online approximation of velocity and the full-order solution using OpenFOAM at two time-steps. It is observed that the proposed methodology can capture the drastic changes in the flow topology over time, which is quantified in Fig. 4 where the relative L 2 -error for velocity and pressure over time is plotted. It should be noted that the validation points include a combination of sample points and time instances that were not used in the construction of the ROM. The average error is less than 5% for pressure and less than 2% for velocity, which is in line with the expected lower average error bound of 1%. Similar  Figure 4. Plot of the relative velocity and pressure error over the examined time range. The dash-lines correspond to the average error, less than 5% for pressure and less than 2% for velocity. observations can be made for the approximated force coefficients over time (Fig.5), where the maximum relative error is approximately 3.5% for both C l and C d . The approximative power of the ROM is tested for various bounds of the projection error ε, i.e. various values of M u , M p and M t . Figure 6 shows the mean error of the velocity and pressure approximations as a function of the square of the the projection error, ε 2 . Both the velocity and pressure approximation achieve the expected convergence rate until ε = 10 −2 . We observe stagnation of the error at approximately 1%, which is assumed to be a limitation derived from the mesh size. We also observe that the approximation error never reaches exactly the error bound, as is expected due to truncation errors introduced by numerical schemes. That is more prevalent for pressure, which also suffers from the induced velocity error and stagnates at approximately 4% error.  Figure 6. Plot of the mean relative L 2 -error of the velocity and pressure approximation as a function of the square of the projection error bound, ε 2 .

Conclusion
The main goal of this work was to outline a POD strategy that can be employed in turbulent flow problems using the finite volume method, namely OpenFOAM's libraries. The proposed approach is a hybrid combining interpolation-and projection-based POD introducing a nonintrusive reduced-order model of the eddy viscosity in the POD-Galerkin approach, allowing a more stream-lined use in industrial applications.
Initially, we formulated the fundamentals for the construction of a ROM for incompressible transient Navier-Stokes equations. Aiming to stay as consistent as possible with the FOM setting, we derived the pressure Poisson equation for the online evaluation of the reduced pressure approximation. The resulting ROM is then modified to account for turbulent effects, while the reduced model of the eddy viscosity employs a linear interpolation scheme.
The potential of the proposed strategy is finally tested using the flow around a NACA airfoil at a moderately high Reynolds number (≈ 320 000) and an angle of attack of 17 • . The proposed ROM strategy has proven to be able to approximate the high-fidelity results with good accuracy. The average L 2 -errors are close to 1% and 5% for velocity and pressure, respectively. It showed quantitatively great ability in closely following significant temporal effects of the flow topology, as well as reliability in the approximation of industrially relevant quantities of interest with maximum relative errors well below 5%.