Finite-Volume High-Fidelity Simulation Combined with Finite-Element-Based Reduced-Order Modeling of Incompressible Flow Problems

Abstract

We present a nonintrusive approach for combining high-fidelity simulations using Finite-Volume (FV) methods with Proper Orthogonal Decomposition (POD) and Galerkin Reduced-Order Modeling (ROM) methodology. By nonintrusive we here imply an approach that does not need specific knowledge about the high-fidelity Computational Fluid Dynamics (CFD) solver other than the velocity and pressure results given on an element mesh representing the related discrete interpolation spaces. The key step in the presented approach is the projection of the FV results onto suitable finite-element (FE) spaces and then use of classical POD Galerkin ROM framework. We do a numerical investigation of aerodynamic flow around an airfoil cross-section (NACA64) at low Reynolds number and compare the ROM results obtained with high-fidelity FV-generated snapshots against similar high-fidelity results obtained with FE using Taylor–Hood velocity and pressure spaces. Our results show that we achieve relative errors in the range of 1–10% in both $H_1$-seminorm of the computed velocities and in the $L_2$-norm of the computed pressure with reasonably few ROM modes. Similar accuracy was obtained for lift and drag.

1. Introduction

High-fidelity numerical simulations based on finite-volume (FV) method [1], finite-element (FE) method [2] and finite-difference (FD) method [3] are generally used to study phenomena governed by partial differential equations [4]. In the context of incompressible fluid flow, these equations are called the Navier–Stokes equations and the branch of computational techniques that deals with it called the Computational Fluid Dynamics (CFD) [5]. The solution of these equations for complicated flows around complex geometries such as a turbine blade [6] requires high degrees of freedom and hence are computationally very demanding even on the most advanced computing facilities [7,8]. The industry, on the other hand, has a need for accurate yet computationally light modeling tools for design optimization and real-time predictions. To this end, researchers have been developing reduced-order models. One such development in the field is towards the use of self-excited oscillator models for prediction of the aerodynamic forces (drag and lift) on the bluff bodies [9,10]. Others have exploited reduced basis methods for higher-dimensional systems by extracting the most useful information from the flow field [11,12,13]. In a reduced basis method, high-fidelity numerical simulation results for a set of input parameters are projected onto a significantly lower-dimensional space to construct reduced basis which is then used to build up a reduced-order model [14,15,16]. Such models exhibit efficiency, robustness, reliability, and accuracy.

Usually, the aerodynamic flow-related problems in wind engineering depend heavily on some critical parameters such as angle of attack, wind velocity, airfoil shape, etc. [17]. For faster-repeated sampling, it could be a legitimate alternative to approximate the behavior based on a linear combination of precomputed solutions on a reduced basis for different choices of input parameter space (physical or geometrical) [18]. This required transformation of the non-affine problems in an affine way which—to the best of the knowledge of the authors—have not been addressed before for wind turbines, where the angle of attack undergoes periodic variations over a complete cycle. Concerning order reduction procedures, known methods to develop the basis in the literature are Greedy Method (GM), Centroidal Voronoi Tessellations (CVT) Taylor, Lagrange, Hermite, and Proper Orthogonal Decomposition (POD) [19]. The techniques provide a subspace which optimally represents a dynamical system with fewer degrees of freedom [20,21]. The dominant modes after the decomposition capture the highly energetic coherent structures in a flow field [22,23,24].

To perform POD, input data sets are required and taken from experiments or numerical studies. The reduced basis obtained from POD is ordered by descending energy content and termed as basis functions or modes [25]. The first use of POD in engineering applications was formally introduced in the study by authors [26,27], who employed the technique to study the formation of coherent structures in a flow field. Later, authors of [20,28] used it to solve problems related to boundary layer turbulence. The author of [29] used FE technique to generate the snapshots and analyzed the turbulent structures using POD around a circular cylinder. The author of [30] extended that technique to develop a Galerkin free POD formulation of incompressible flows using snapshots from a DNS solver. The authors of [31] developed a minimum residual projection to build coupled velocity-pressure POD-ROM for incompressible Navier–Stokes equations using the high-fidelity data from the FE approach. The new model enabled to increase the velocity and pressure fields accuracy. The authors of [11] employed POD to develop a Galerkin ROM for the construction of both pressure and velocity spaces. They adopted an enrichment procedure of the velocity space by supremizer solutions to enhance the stability and monitored the fulfillment of an inf-sup Ladyzhenskaya–Babuska–Brezzi (LBB) stability condition at the reduced-order level. Whereas studies related to the affine representation of velocity and pressure fields are important and based on the reduced basis and a transformation between physical and reference geometries. The authors of [32] employed the strategy and developed a ROM for the Stokes problem using a domain decomposition technique. The author of [24] has shown the integration of FVM and POD, whereas authors of [16] demonstrated the applicability of such an approach on a two-dimensional cylinder.

In this paper, the concept of mixing of methods for the development of novel FE-based ROM is demonstrated. We devised mixed method term for the reduced FE solution obtained from the snapshots generated by FV method. Whereas, uniform method is called for the reduced FE solution based on the snapshots created by FE method. The methods are implemented, and their applicability is manifested on a parametrized two-dimensional NACA64 airfoil. For a combination of input physical (velocity) and geometrical (incidence angle) parameters, high-fidelity simulations of FE and FV are conducted using an in-house code and OpenFOAM (4.0) [33] respectively. First, an optimal reduced basis is constructed on the reference domain with the application of POD on the snapshots generated using FE and FV methods. Using the reduced space, a reduced-order model based on FE method is formed in the Nutils [34] framework. The ROM solutions are compared with the corresponding FE, and FV high-fidelity discretization and the differences are illustrated in the energy spectrum, error magnitudes, computational time for the methods and the shape of the modes. Supremizer solutions are added to the velocity space such that non-spurious pressure modes can be developed resulting in a stable solution. The whole procedure is demonstrated with a novel procedure, that transforms wind turbine-related problems from a classical non-affine context to an affine framework thereby motivating the use of precomputed solutions as a linear combination of reduced basis functions.

In this paper, Section 2 describes the systems of governing equations, the discretization technique for the map and the affine transformation. The formulation of reduced-order modeling (ROM) is then presented in Section 3. The section also describes the implementation of our solver and boundary conditions. In Section 4, the results are compared for the two developed models in terms of flow dynamics, spectrum of velocity/pressure, $L_2$ error and the computational efficiencies evaluated in a test case of the NACA64 airfoil. The last section describes the conclusions of the work followed by future recommendations.

2. Theory

This section describes the mathematical formulation of the reduced assembly based on the governing equations, the POD, and the affine representation of the equations.

2.1. Governing Equations

The steady Navier–Stokes equations can model the flow in a low Reynolds number regime, where the equations are non-linear due to the presence of the convection term [35]. The equations are considered for the cases when $\nu$ is constant. The flow remains incompressible below the Mach number of 0.3. The mathematical description of the Navier-Stoke equations is given as

$$\begin{array}{cc}\hfill -\nu \Delta \mathit{u}+(\mathit{u}·\nabla )\mathit{u}+\nabla p& =\mathit{f}\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \nabla ·\mathit{u}& =0\mathrm{in}\mathsf{\Omega },\hfill \\ \hfill \mathit{u}& =\mathit{g}\mathrm{on}{\mathsf{\Gamma }}_{\mathrm{D}},\hfill \\ \hfill -p\mathit{n}+\nu (\nabla \mathit{u})\mathit{n}& =\mathit{h}\mathrm{on}{\mathsf{\Gamma }}_{\mathrm{N}}.\hfill \end{array}$$

The equation represents distributive force $\mathit{f}$, Dirichlet data $\mathit{g}$ and Neumann fluxes $\mathit{h}$ where $\nu$ is the kinematic viscosity, $\mathit{u}=(u_1,u_2)$ is the velocity field, p, is the pressure and $\mathit{n}$ shows the normal unit vector to ${\mathsf{\Gamma }}_{\mathrm{N}}$. The domain of interest is represented by $\mathsf{\Omega }\subset {\mathbb{R}}^d$, with mixed boundary conditions. Figure 1 represents the boundaries of the domain, where a solid line at the inlet illustrates Dirichlet (${\mathsf{\Gamma }}_{\mathrm{D}}$) and the dotted line represents the Neumann (${\mathsf{\Gamma }}_{\mathrm{N}}$) conditions at the boundary. The mathematical description of boundary is introduced as

$$\partial \mathsf{\Omega }={\mathsf{\Gamma }}_{\mathrm{D}}\cup {\mathsf{\Gamma }}_{\mathrm{N}},{\mathsf{\Gamma }}_{\mathrm{D}}\cap {\mathsf{\Gamma }}_{\mathrm{N}}=\varnothing ,{\mathsf{\Gamma }}_{\mathrm{D}}\ne \varnothing .$$

2.1.1. FE Discretization

The discretization of the Navier–Stokes equations is developed by determining the weak variational formulation of our problem [36,37]. To manage the inhomogeneous Dirichlet boundary conditions, a lifting function is constructed ℓ corresponding to the Dirichlet boundary condition such that $\ell =\mathit{g}$ on ${\mathsf{\Gamma }}_{\mathrm{D}}$, writing $\mathit{u}={\mathit{u}}_0+\ell$ where ${\mathit{u}}_0$ satisfies homogeneous Dirichlet conditions ${\mathit{u}}_0=\mathbf{0}$ on ${\mathsf{\Gamma }}_{\mathrm{D}}$. The governing equations are multiplied with the test functions ($\mathit{w}$, q) and perform integration by parts to obtain the solution onto the physical domain. Hence, the weak form is described as

$$\begin{array}{cc}\hfill \overline{d}({\mathit{u}}_0,\mathit{w})+c({\mathit{u}}_0,{\mathit{u}}_0,\mathit{w})+b(p,\mathit{w})& ={\overline{f}}_1\left(\mathit{w}\right),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill b(q,{\mathit{u}}_0)& =-b(q,\ell ).\hfill \end{array}$$

(2)

the bi- and trilinear forms are defined as

$$\begin{array}{cc}\hfill \overline{d}({\mathit{u}}_0,\mathit{w})& =d({\mathit{u}}_0,\mathit{w})+c({\mathit{u}}_0,\ell ,\mathit{w})+c(\ell ,{\mathit{u}}_0,\mathit{w}),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\overline{f}}_1\left(\mathit{w}\right)& ={\int }_{{\mathsf{\Gamma }}_{\mathrm{N}}}\mathit{h}·\mathit{w}+{\int }_{\mathsf{\Omega }}\mathit{f}·\mathit{w}-d(\ell ,\mathit{w})-c(\ell ,\ell ,\mathit{w}),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill d(\mathit{v},\mathit{w})& =\nu {\int }_{\mathsf{\Omega }}\nabla \mathit{v}:\nabla \mathit{w},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill b(p,\mathit{w})& =-{\int }_{\mathsf{\Omega }}p\nabla ·\mathit{w},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill c(\mathit{z},\mathit{v},\mathit{w})& ={\int }_{\mathsf{\Omega }}(\mathit{z}·\nabla )\mathit{v}·\mathit{w}.\hfill \end{array}$$

(7)

For Galerkin discretization, we then choose finite dimensional subspaces $U^h\subset U$, $P^h\subset P$, $U_0^h=\left\{\mathit{v}\in U^h|\mathit{v}=\mathbf{0}\mathrm{on}{\mathsf{\Gamma }}_{\mathrm{D}}\right\}$ and wish to find $({\mathit{u}}_0,p)\in U_0^h\times P_n$ such that (1) and (2) hold for all $(\mathit{w},q)\in U^h\times P^h$.

Specifically, let $\mathsf{\Omega }={⋓}_i{K}_e$ be a decomposition of the computation domain into quadrilateral pairwise disjoint elements $K_e$, and for each element, define a mapping

$${\pi }_e:K_e->{[0,1]}^2$$

Given a selection of shape functions defined on the reference element ${[0,1]}^2$, we can then define the corresponding shape functions on $K_e$ by pullbacks through ${\pi }_e$. With this in mind, let us define on ${[0,1]}^2$ the bilinears and biquadratics,

$$\begin{array}{cc}\hfill Q_{1,1}& =\mathrm{span}\left\{1,\xi ,\eta ,\xi \eta \right\}\hfill \\ \hfill Q_{2,2}& =\mathrm{span}\left\{1,\xi ,\eta ,\xi \eta ,{\xi }^2,{\eta }^2,{\xi }^2\eta ,\xi {\eta }^2,{\xi }^2{\eta }^2\right\}\hfill \end{array}$$

and define our finite dimensional subspaces as globally $C^0$-continuous biquadratics and bilinears [38], respectively. These are also known as Taylor–Hood [39] type elements.

$$\begin{array}{cc}\hfill U^h& =\left\{\mathit{v}\in H^1{\left(\mathsf{\Omega }\right)}^2\cap C^0{\left(\mathsf{\Omega }\right)}^2|{\pi }_e^{*}\mathit{v}\in Q_{2,2}^2\forall e\right\}\hfill \\ \hfill P^h& =\left\{p\in L^2\left(\mathsf{\Omega }\right)\cap C^0\left(\mathsf{\Omega }\right)|{\pi }_e^{*}p\in Q_{1,1}\forall e\right\}.\hfill \end{array}$$

2.1.2. FV Discretization

The Navier–Stokes equations in space are discretized using the FV approximation technique. The PDEs are discretized over the control volume using the integral form approach. To avoid additional errors due to the differences between the domain subdivisions into finite cells/elements, the same number of non-overlapping cells are used in both methods. The discretized equations for momentum and continuity are represented for each term using the Gauss theorem which converts the volume integrals into surface integrals [40]. The integral form of the momentum equation over each control volume is given by

$${\int }_V^{}(\mathit{u}·\nabla )\mathit{u}dV-\nu {\int }_V^{}\Delta \mathit{u}dV+{\int }_V^{}\nabla pdV=0$$

(8)

Here, the discretization of each term in the Navier–Stokes is discussed explicitly. The convective term, using the linear variation of $\mathit{u}$ and applying the generalized form of Gauss theorem is discretized

$${\int }_V^{}(\mathit{u}·\nabla )\mathit{u}dV={\int }_V^{}\nabla ·\left(\mathit{u}\mathit{u}\right)dV={\int }_S^{}d{\mathit{S}}_{\mathit{f}}·\left(\mathit{u}\mathit{u}\right)\approx \sum _f^{}{\mathit{S}}_{\mathit{f}}·{\mathit{u}}_{\mathit{f}}{\mathit{u}}_{\mathit{f}}=\sum _f^{}F{\mathit{u}}_{\mathit{f}}$$

(9)

the ${\mathit{u}}_f$ represents the velocity at the center of the face of the control volume. where $S_f$ is the surface area of the face and F is the face flux. The F represents the mass flux over the surface of each face constituting the control volume. It is formed from the product of ${\mathit{S}}_{\mathit{f}}·{\mathit{u}}_{\mathit{f}}$.

The Laplacian, continuity, and pressure terms are discretized as follows

$${\int }_V^{}\nabla ·(\nu \nabla \mathit{u})dV={\int }_S^{}d{\mathit{S}}_{\mathit{f}}·(\nu \nabla \mathit{u})\approx \sum _f^{}\nu {\mathit{S}}_{\mathit{f}}·{(\nabla \mathit{u})}_f=\nu \sum _f^{}F$$

(10)

$${\int }_V^{}\nabla ·\mathit{u}dV={\int }_S^{}d\mathit{S}·\left(\mathit{u}\right)=\sum _f^{}{\mathit{S}}_{\mathit{f}}·{\mathit{u}}_{\mathit{f}}=\sum _f^{}F$$

(11)

$${\int }_V^{}\nabla pdV={\int }_S^{}d\mathit{S}·\left(p\right)=\sum _f^{}{\mathit{S}}_{\mathit{f}}·p$$

(12)

p represents the values at the center of each control face of the control volume. The integrals appearing in the equations above are linearized employing the first/second-order upwind, blended difference, and the central differencing schemes available in OpenFOAM (4.0) [33]. The continuity equation is solved taking the divergence of momentum and solving as pressure Poisson equation. For further explanation, the reader is referred to [33]. The schemes used for the solution of numerical examples are described in Section 4.1.

2.2. Projection of FV Results to FE Space

Earlier work on combining the FV high-fidelity solvers with POD Galerkin ROM is presented in [16,24]. They have adopted a consistent approach and included the contribution of the face flux field (face area vector multiplied by the face flux) in addition to the velocity and pressure fields into the POD Galerkin ROM. This is done to ensure that the method employed to construct the ROM remain consistent with the Full-order Model (FOM). However, the drawback with this is that it introduces an extra field variable and that it requires access to the high-fidelity face flux field, which is not always easily accessible.

Herein, we explore a nonintrusive approach, i.e., no specific knowledge is needed about the high-fidelity CFD solver other than the velocity and pressure snapshots given on an element mesh representing the related discrete interpolation spaces. The key step in our approach is that the FV high-fidelity velocity and pressure fields are projected onto suitable FE spaces and then the classical POD Galerkin ROM algorithm is applied. Unfortunately, the projection from FV-spaces to FE spaces introduces an inconsistency between the FV-FOM and the resulting ROM, i.e., the projected high-fidelity FV velocities and pressures do not fulfill Equations (1) and (2). The significance of this inconsistency error relative to the error introduced by ROM is an important issue in our numerical investigation presented below.

We have used OpenFOAM in this study to produce the high-fidelity FV results and the computed velocities and pressures are transferred to the Taylor–Hood FE space in a two-stage process:

• Cell centers to cell nodes interpolation:
  In the first stage, the velocity and pressure fields that are calculated by OpenFOAM at the center of the control volumes are interpolated onto the control volume vertices using Inverse Distance Weighting (IDW) interpolation before they are stored at VTK-files. Thus, the stored velocity and pressure results on the VTK-files are implicitly assumed to be interpolated as bilinear fields.
• Projection to taylor hood FE space:
  In the second stage, the velocity and pressure results stored on the VTK-files are projected onto the Taylor–Hood FE space. Herein, we have used $L_2$-projection for both velocities and pressure.

To minimize the inconsistency error, we generated first the FE-mesh with characteristic element size $2h$ and then uniformly refined every element into four new elements with characteristic element size h. As illustrated in Figure 2a this was done in such a way that the geometry of the domains ${\mathsf{\Omega }}^h$ and ${\mathsf{\Omega }}^{2h}$ is exactly the same for the two different meshes.

The Taylor–Hood FE space is characterized by biquadratic interpolation of velocities and bilinear interpolation of the pressure, whereas the OpenFOAM results stored on the VTK-files are assumed interpolated bi-linearly for both variables. We have here chosen to project the FV results computed on a mesh with characteristic size h onto an FE-mesh of characteristic size $2h$ as this gives the same number of degrees of freedoms (DOFs) for the FV and FE velocity fields. Furthermore, the original FV results are piecewise (cell-centered) constant, and therefore considered to be less accurate than Taylor–Hood FE of the same mesh size.

2.3. Parametric Dependency

We consider here the case where the problem defined by (1) and (2) is depends on several parameters. Denote by $\mathcal{P}$ the parameter space and by $\mathit{\mu }$ any given element of $\mathcal{P}$. For flow problems we typically have the following classes of parameters: physical parameters, boundary condition parameters and geometric parameters. We shall denote the linear and bilinear forms dependence on the parameters $\mathit{\mu }$ with notation such as $d(u,w;\mathit{\mu })$, and similarly for other forms.

Physical parameters are very often material parameters describing the continuum at hand e.g., the viscosity $\nu$ of an incompressible fluid. Boundary condition parameters are related to the Dirichlet and Neumann conditions, e.g., g and h in the governing equations of the Navier–Stokes equation. In particular, non-homogeneous Dirichlet parameters g (that introduce the lifting function l) needs to be handled properly [41].

Geometric parameters imply that the computational domain $\mathsf{\Omega }=\mathsf{\Omega }\left(\mathit{\mu }\right)$ varies. Let $\widehat{\mathsf{\Omega }}$ be a reference domain which herein is described by the physical domain of $\mathsf{\Omega }({\mu }_2=\phi =0)$ (i.e., angle of attack equal to zero) and a corresponding FE and FV discretization of $\widehat{\mathsf{\Omega }}$. A one-to-one mapping is then introduced to map the physical domain to the reference geometry and is defined by

$${\pi }_{\mathit{\mu }}:\widehat{\mathsf{\Omega }}\left(\mathit{\mu }\right)\to \mathsf{\Omega }.$$

Suitable functions are created in the physical domain $\mathsf{\Omega }\left(\phi \right)$ by applying the pullback on the transformation which maps the solution from the reference geometry onto the physical geometry and is given by (see [41] for more details)

$$\begin{array}{c}\hfill \mathit{v}={\pi }_{\mathit{\mu }}^{*}\widehat{\mathit{v}},p={\pi }_{\mathit{\mu }}^{*}\widehat{p}.\end{array}$$

(13)

The geometry variation affects the linear and bilinear forms. Hence for the convective form c, the mapping takes the form

$$c({\pi }_{\phi }^{*}\widehat{\mathit{u}},{\pi }_{\phi }^{*}\widehat{\mathit{v}},{\pi }_{\phi }^{*}\widehat{\mathit{w}};\phi )={\int }_{\widehat{\mathsf{\Omega }}}(\widehat{\mathit{u}}·{\mathit{J}}^{-⊺}\nabla )\widehat{\mathit{v}}·\widehat{\mathit{w}}\approx \sum _{i=0}^{2n}{\phi }^i{\int }_{\widehat{\mathsf{\Omega }}}(\widehat{\mathit{u}}·{\mathit{B}}_i^{(-)}\nabla )\widehat{\mathit{v}}·\widehat{\mathit{w}}.$$

and for the diffusive part d it reads

$$\begin{array}{cc}\hfill d({\pi }_{\phi }^{*}\widehat{\mathit{u}},{\pi }_{\phi }^{*}\widehat{\mathit{w}};\phi )& =\nu {\int }_{\widehat{\mathsf{\Omega }}}({\mathit{J}}^{-⊺}\nabla )\widehat{\mathit{u}}:({\mathit{J}}^{-⊺}\nabla )\widehat{\mathit{w}}=\nu {\int }_{\widehat{\mathsf{\Omega }}}\nabla \widehat{\mathit{u}}:({\mathit{J}}^{-1}{\mathit{J}}^{-⊺}\nabla )\widehat{\mathit{w}},\hfill \end{array}$$

and finally, the incompressibility constraints b becomes

$$b({\pi }_{\phi }^{*}\widehat{p},{\pi }_{\phi }^{*}\widehat{\mathit{w}};\phi )={\int }_{\widehat{\mathsf{\Omega }}}\widehat{p}({\mathit{J}}^{-⊺}\nabla )·\widehat{\mathit{w}}={\int }_{\widehat{\mathsf{\Omega }}}\widehat{p}{\mathit{J}}^{-⊺}:\nabla \widehat{\mathit{w}},$$

To find the detail derivation of the above ROM equations and the explanation of expressions ${\mathit{B}}_i^{(-)}$ the reader is referred to the article by Fonn et al. [41]

2.4. ROM Using POD

The notion of ROM is introduced to obtain an accurate solution of the Navier–Stokes equations over a discrete solution space whose dimensions far exceeds the natural dimensions of the model itself [42]. In other words, the practical dimensions of the space

$$\left\{({\widehat{u}}_0\left(\mathit{\mu }\right),p\left(\mathit{\mu }\right))|\mathit{\mu }\in \mathcal{P}\right\}$$

(14)

is considerably more manageable than the dimension of any discretization necessary to achieve good accuracy with conventional FE/FV methods. We will solve the steady Navier–Stokes problem and proceed by sampling and generate the ensemble of solutions which corresponds to different snapshots to produce an approximation of the space of Equation (14). The POD is employed to reduce the dimension of the space such that the solution can be described by an optimal combination of the most energy rich modes [21,42].

First the ensemble is developed by performing the high-fidelity simulations at each state chosen by a combination of input parameters with FE and FV methods (${\mathit{\mu }}_i=({\phi }_i,{\mathit{g}}_i$))

$$\mathsf{\Phi }={\left\{{\varphi }_i=({\widehat{\mathit{u}}}_0^i({\phi }_i,{\mathit{g}}_i),({\widehat{p}}^i({\phi }_i,{\mathit{g}}_i)\right\}}_{i=1}^N$$

A correlation matrix is then constructed from the ensemble for velocity and pressure respectively,

$$C_{ij}=a({\varphi }_i,{\varphi }_j)$$

Given eigenpairs $({\lambda }_i,{\mathit{v}}^i)$ of $\mathit{C}$, one may then pick an M such that

$$\begin{array}{c}\hfill \sum _{i=1}^M{\lambda }_i\ge \left(1-{ϵ}^2\right)\sum _{i=1}^N{\lambda }_i\end{array}$$

(15)

where $ϵ>0$ is a suitable error tolerance. It is hoped that $M\ll N$, as M will be the dimension of the reduced space (it can be interpreted as the practical dimension of the space given by Equation (14) Error $ϵ$, as measured in the norm ${\left|\right|·\left|\right|}_a=\sqrt{\left(a\right(·,·\left)\right)}$. The reduced basis functions are then given by

$${\psi }_j=\frac{1}{\sqrt{{\lambda }_j}}\sum _{i=1}^N{v}_i^j{\varphi }_i,1\le j\le M$$

(16)

It can be noted that the ${\psi }_j$ are orthonormal in the $a(·,·)$ inner product [43]. In practice, this construction is performed separately on the velocity and pressure fields. Thus, first we choose a velocity correlation function

$$a_u((\widehat{{\mathit{u}}_0},\widehat{p}),(\widehat{{\mathit{v}}_0},\widehat{q}))={(\nabla \widehat{{\mathit{u}}_0},\nabla \widehat{{\mathit{v}}_0})}_{L_2\left(\mathsf{\Omega }\right)}$$

Secondly, we choose the pressure correlation function

$$a_p((\widehat{{\mathit{u}}_0},\widehat{p}),(\widehat{{\mathit{v}}_0},\widehat{q}))={(\widehat{p_0},\widehat{q_0})}_{L_2\left(\mathsf{\Omega }\right)}$$

The velocity and pressure basis are constructed as

$${\widehat{\mathit{u}}}_j^{\mathrm{R}}=\frac{1}{\sqrt{{\lambda }_j}}\sum _{i=1}^N{v}_i^{\left(u\right),j}{\widehat{\mathit{u}}}_0^i,{\widehat{p}}_j^{\mathrm{R}}=\frac{1}{\sqrt{{\lambda }_j}}\sum _{i=1}^N{v}_i^{\left(p\right),j}{\widehat{p}}_i,1\le j\le M$$

2.5. Mixed and Uniform Methods

Classical reduced basis methods generally use the same high-fidelity method both for constructing the ensemble and creating the reduced basis method itself. In other words, given that the ensemble solutions $\mathit{u}\left({\mathit{\mu }}_i\right)$ satisfy the equation

$$\mathit{A}\left({\mathit{\mu }}_i\right)\mathit{u}\left({\mathit{\mu }}_i\right)=\mathit{b}\left({\mathit{\mu }}_i\right),$$

where $\mathit{A},\mathit{b}$ denote a parameter-dependent system matrix and right-hand side, respectively, the reduced system is given as

$${\mathit{A}}^{\mathrm{R}}\left(\mathit{\mu }\right)={\mathit{V}}^{⊺}\mathit{A}\left(\mathit{\mu }\right)\mathit{V},{\mathit{b}}^{\mathrm{R}}\left(\mathit{\mu }\right)={\mathit{V}}^{⊺}\mathit{b}\left(\mathit{\mu }\right),$$

the “tall” matrix $\mathit{V}$ is generated from the relation between high-fidelity basis and the reduced-order basis as given in Equation (16)

$${\left({\mathit{V}}^T\right)}_{ij}=\frac{v_i^j}{\sqrt{{\lambda }_j}}$$

In the following, the above method will be denoted as uniform, to distinguish it from the mixed method. By the latter it is meant that a model is obtained by generating ensemble solutions with a given high-fidelity method $\left(1\right)$,

$${\mathit{A}}_{\left(1\right)}\left({\mathit{\mu }}_i\right)\mathit{u}\left({\mathit{\mu }}_i\right)={\mathit{b}}_{\left(1\right)}\left({\mathit{\mu }}_i\right),$$

and then using a different high-fidelity method $\left(2\right)$ as a basis for the reduction, i.e.

$${\mathit{A}}^{\mathrm{R}}\left(\mathit{\mu }\right)={\mathit{V}}^{⊺}{\mathit{A}}_{\left(2\right)}\left(\mathit{\mu }\right)\mathit{V},{\mathit{b}}^{\mathrm{R}}\left(\mathit{\mu }\right)={\mathit{V}}^{⊺}{\mathit{b}}_{\left(2\right)}\left(\mathit{\mu }\right).$$

This is motivated by the idea of enable nonintrusive ROM. So long as the two high-fidelity methods purport to solve the same physical model (and that this physical model accurately reflects reality), one may reasonably hope that mixing methods in this way will be sufficiently accurate.

In the following, the two high-fidelity methods used in the present study are the FE and FV methods described in Section 2.1.1 and Section 2.1.2 respectively. By the “Uniform” method, it is meant that the reduced model is based fully on the FE method, and by the “Mixed” method, it is meant that the snapshots (and therefore the basis) were generated with the FV method, but the method underlying the reduction was the FE method. As such, the difference between the two reduced models is entirely in their bases.

3. Development of the Solver

ROM Solver

The reduce space after the POD stage is obtained. The solution for the reduced-order model requires solving several systems of the form

$$\left(\begin{array}{cc}{\mathit{A}}_{vv}& {\mathit{B}}_{vp}\\ {\mathit{B}}_{pv}\end{array}\right)\left(\begin{array}{c}\mathit{v}\\ \mathit{p}\end{array}\right)=\left(\begin{array}{c}{\mathit{b}}_v\\ {\mathit{b}}_p\end{array}\right)$$

The matrix ${\mathit{B}}_{vp}$ is observed to be numerically rank-deficient, i.e., it has a nontrivial kernel, therefore leading to unstable pressure solutions. To remedy this, the velocity space is artificially enriched with supremizers basis functions whose purpose is not to achieve greater approximative power but rather to keep the system well-conditioned [11]. Therefore, for each pressure snapshot ${\widehat{p}}_i$, represented by the high-fidelity coefficient vector $\mathit{p}$, at a parameter ${\mathit{\mu }}_i$, to solve the system

$${\mathit{M}}_V\left({\mathit{\mu }}_i\right)\mathit{v}={\mathit{B}}^{⊺}\left({\mathit{\mu }}_i\right)\mathit{p}$$

for the velocity coefficient vector $\mathit{v}$, which then yields the supremizer function at that particular parameter-value. Here, ${\mathit{M}}_V$ is the $H^1$-norm mass matrix, and $\mathit{B}$ is the matrix associated with the $b(·,·;\mathit{\mu })$ bilinear form.

The supremizers are then reduced using the POD as well, just like the regular velocity solutions, into a separate and independent reduced space, yielding the system

$$\left(\begin{array}{ccc}{\mathit{A}}_{vv}& {\mathit{A}}_{vs}& {\mathit{B}}_{vp}\\ {\mathit{A}}_{sv}& {\mathit{A}}_{ss}& {\mathit{B}}_{sp}\\ {\mathit{B}}_{pv}& {\mathit{B}}_{ps}\end{array}\right)\left(\begin{array}{c}\mathit{v}\\ \mathit{s}\\ \mathit{p}\end{array}\right)=\left(\begin{array}{c}{\mathit{b}}_v\\ {\mathit{b}}_s\\ {\mathit{b}}_p\end{array}\right)$$

(17)

In practice, the sizes of each basis are kept equal, so with M degrees of freedom each for velocity, supremizers and pressure, one obtains a $3M$-sized system [11].

4. Results and Discussion

The aim of the numerical investigation is to see how well the mixed method compares with the uniform method. In particular, we want to study the significance of the inconsistency error compared to the error introduced by choosing a reduced-order basis. Our aim regarding accuracy is to achieve a difference between the high-fidelity simulation and the reduced-order simulation in the range of 1–10% relative error in suitable norms for the quantity of interest.

In the following, we first describe in detail the high-fidelity simulation setup and the choice of parametric values for the snapshot creation. To introduce the study of the obtained ROM results we display the most energy-rich modes and corresponding spectrum of the reduced-order basis. Then we focus on the relative error convergence plots of $H_1$-seminorm of the computed ROM-velocities and the $L_2$-norm of the computed ROM-pressures. Similarly, we display the relative error convergence plot of the computed ROM lift and drag. We then report the parametric dependence of the difference in the FV and FE high-fidelity solutions. In the end, we report the observed speedup.

4.1. High-Fidelity Simulation Setup

The developed ROM is benchmarked on a two-dimensional NACA64, which is the airfoil located at the outer most section of NREL 5MW wind turbine [44,45]. For testing our methodology, the Reynolds number range of the simulation is designed to remain in the laminar range (2–20). Under the current operating conditions, the flow is creeping and produces insignificant changes while flowing over the airfoil. The present aim of the numerical examples is to establish the accuracy of the developed methods, rather than studying the details of flow dynamics across the airfoil. To simulate the flow field, a cylindrical domain is constructed with a diameter of $30\mathrm{m}$ with the airfoil chord length of $1\mathrm{m}$. The discretized mesh for the two-dimensional case is shown in Figure 3. In particular, inlet velocity is chosen $2 \mathrm{m}/\mathrm{s}<\mathit{g}<20 \mathrm{m}/\mathrm{s}$ as the physical parameters and angle of attack $-25{}^{\circ }<\phi <25{}^{\circ }$ as the geometric parameter as shown in Figure 4. The viscosity is set at $\nu$ = 1 m${}^2$/s, with strongly enforced no-slip conditions at the airfoil boundary.

To handle variation in geometry due to variation in the angle of attack $\phi$ we introduce the following geometrical mapping written in reference variables as a function of the angle of attack and radius r

$$\begin{array}{c}\hfill {\pi }_{\phi }^{-1}\widehat{\mathit{r}}=\mathit{R}\left(\phi \theta \left(r\right)\right)\widehat{\mathit{r}},\end{array}$$

(18)

where $r=\parallel \widehat{\mathit{r}}\parallel$, $\widehat{\mathit{r}}=\left(\widehat{x},\widehat{y}\right)$ and $\mathit{R}$ is the rotation matrix

$$\begin{array}{c}\hfill \mathit{R}\left(a\right)=\left[\begin{array}{cc}\cos a& -\sin a\\ \sin a& \cos a\end{array}\right].\end{array}$$

(19)

The function $\theta$ is chosen as sufficiently smooth function satisfying the properties

$$\begin{array}{ccc}\hfill \theta \left(r\right)& =1\hfill & \hfill r<r_{\min },\end{array}$$

(20)

$$\begin{array}{ccc}\hfill \theta \left(r\right)& =0\hfill & \hfill r>r_{\max }.\end{array}$$

(21)

Figure 5 describes the variation of r with the changes in the incidence angle. One possible choice is with $\overline{r}=(r-r_{\min )}/(r_{\max }-r_{\min })$

$$\begin{array}{c}\hfill \theta \left(r\right)=\left\{\begin{array}{cc}1,\hfill & r<r_{\min },\hfill \\ 0,\hfill & r>r_{\max },\hfill \\ {(1-\overline{r})}^3(3\overline{r}+1),\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

(22)

Noting that the airfoil is within the region $r<r_{\min }$, and the external boundaries of the domain are entirely contained within the region $r>r_{\max }$, this produces a suitable parameter-dependent domain $\mathsf{\Omega }\left(\phi \right)$ as shown in Figure 3. Here we have set $r_{\min }$ = 1 m and $r_{\max }$ = 10 m.

The solution strategy is divided into offline and online stages, where the basis is precomputed in the offline stage such that fast parametric estimates can be made in the online stage. An ensemble of high-fidelity solution is created with different choices of the parameters $(\mathit{g},\phi )$, see Section 4.2. For the FV solutions, a mesh is used that is twice as dense as for FE as described in Section 2.2. Three meshes are shown in Figure 3 has been used for providing FV high-fidelity velocity and pressure fields. Based on the computed results we choose the ${\mathbb{G}}_2$ that is characterized by having the following number of elements $N_{\theta }\times N_r=80\times 36=2880N_{\mathrm{el}}$.

The FV and FE-based solvers are created in OpenFOAM (4.0) [33] and Nutils (open source numerical tool in Python [34]) respectively. To provide an accurate estimate of initial guess over the discretize points for higher speedups and rapid convergence rates, the solution for the Stokes equation is first obtained. Later the convective terms are introduced to develop the standard system of Navier–Stokes. The FE-based solver employs Newton iteration to obtain the solution of velocity and pressure, which is calculated until the desired value of the residuals is achieved (iteration tolerance is ${10}^{-10}$ in velocity $H^1$-seminorm). For the FV solution, the face values are determined from the cell centers using an interpolation scheme. Linear upwind (${\varphi }_f={\varphi }_N+0.5({\varphi }_P-{\varphi }_N)$) scheme is employed to discretize the velocity in the convective term of the momentum equation, where ${\varphi }_f$ represents the values at the face of the control volume, and ${\varphi }_N$ and ${\varphi }_P$ shows the values at the neighboring upwind and downwind cell centers respectively. The pressure and diffusion terms are discretized using linear schemes (${\varphi }_f=f_x{\varphi }_N+(1-f_x){\varphi }_N$), where $f_x$ is the interpolation factor defined as the ratio of distances between the face and the downwind cell center to the distance between the upwind and downwind cell centers [33]. For the solution of FV algebraic equations, the divergence of momentum equation produces an elliptic equation for the modified pressure. The resulting equation along with the momentum is solved in a segregated manner using the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm. The converged solution of the field variables ($\mathit{u},p$) is located at the cell centers. IDW method ($z={\sum }_{i=1}^n(z_i/d_i)/{\sum }_{i=1}^n(1/d_i)$, where d is the distance between the cell nodes and cell centers, and z is the corresponding cell center value) is used to interpolate the values from the cell centers to the cell nodes. This procedure provides a well-defined solution for both FE and FV methods at the same locations (see Figure 2), which eventually ease the projection of the solution on the reduced basis at the later stage. The FE and FV ensembles are created on a desktop computer with Intel Core i7-4930MX 3.00 GHz, 8 logical cores, 32GB RAM with 4MB of cache memory.

4.2. Snapshots Creation

To generate datasets of high-fidelity simulations, $10\times 11=110$ snapshots for the pressure and velocity fields over the range of parametric values ($\mathit{g},\phi$) are considered as shown in the Figure 4. The Reynolds number values depict a laminar flow regime and vary between 2–20, thus snapshots are only considered for the steady state solutions of the Navier–Stokes. To establish the fact that the flow field remains consistent under steady and unsteady conditions, one simulation is formulated under each condition. The results are plotted in Figure 6 which manifests that the two simulations converge with minimal differences. Thus, steady solutions over the parametric space are employed for the simulations conducted afterward. Each subsequent snapshot contains the velocity and pressure field variables over the domain. The ensemble is generated for two sets of case definitions employing FV and FE methods. Using these ensembles, two ROM models are developed, i.e., FE-based ROM with snapshots generated from FV (mixed ROM) and FE-based ROM with snapshots generated from FE (uniform ROM). The detailed comparison of each ROM is presented in the upcoming section.

4.3. Spatial Development of Modes/Reduced Basis

To gain insight into the type of eigenmodes we get using POD we display the first six reduced-order basis function for velocities in Figure 7 and Figure 8 and pressures in Figure 9 and Figure 10 obtained with the uniform and mixed ROM methods, respectively. The geometrical variation is highlighted in the modes with vortex fronts of various spatial directions. The POD produces in general optimal linear subspace of modes by acting as an energy filter which identifies the large coherent structures and passes the less coherent structures to the lower modes [46,47]. The figures represent the modes sorted in a hierarchical manner of high to low energetic states of the system with its spatial evolution. All the sorted modes are orthogonal to each other (in respective $H_1$-seminorm for velocity and $L_2$-norm for pressure) and exhibit a distinct spatial pattern. The change in the spatial distribution of the small coherent structures in the modes is considered to be due to the redistribution of the energy between scales of varying length.

The plots of mixed and uniform reduced methods demonstrate that the first few velocity modes show a mean behavior of the flow field. The first mode, in particular, carries the highest portion of energy and exhibits a mean flow pattern of the velocity field. The second mode exhibits a wake region behind the airfoil; however, the mean behavior is still reflected in it. The third mode interestingly manifests two distinct vortex regions behind the airfoil [48]. These dipole vortex structures turn into tripole structure in the downstream direction of the fourth mode. In the subsequent modes, an addition of one vortex structure is seen in all modes. Increasing vortex structures show the limited coherence of energy in the flow field. The unique large coherent structures are only present in the first few modes, and the lower modes characterize the vorticity structure present inside the flow field. The vorticity mixing and the eddies in the wake are observed to collapse and became less pronounced. The initial modes represent the global behavior (large scale coherent motion), while the following modes have predominantly shown the local response (such as the representation of small eddies in the wake region). For instance, the scales depicted in the corresponding modes principally describes small scales with a subsequent increase in the number of vortices in the wake. Interestingly, the modes from both high-fidelity input methods have shown similar behavior, which demonstrates that irrespective of the input method used, the six first reduced basis functions remain nearly the same.

4.4. Spectrum and Related Error

The decay of eigenvalues of the mixed and uniform methods is depicted for the sixty-first modes in the Figure 11a whereas the accumulated energy for the twenty-first modes is shown in Figure 11b and tabulated in

Table 1. NACA64: Accumulated energy for the 20 first velocity and pressure modes.

	Mixed-Velocity	Uniform-Velocity	Mixed-Pressure	Uniform-Pressure
DoFs	${\sum }_{\mathit{i}=\mathbf{1}}^{\mathbf{20}}{\mathit{\lambda }}_{\mathit{i}}$	${\sum }_{\mathit{i}=\mathbf{1}}^{\mathbf{20}}{\mathit{\lambda }}_{\mathit{i}}$	${\sum }_{\mathit{i}=\mathbf{1}}^{\mathbf{20}}{\mathit{\lambda }}_{\mathit{i}}$	${\sum }_{\mathit{i}=\mathbf{1}}^{\mathbf{20}}{\mathit{\lambda }}_{\mathit{i}}$
1	0.5855145451	0.5764554385	0.7334786512	0.7413536071
2	0.9320673270	0.9276187838	0.9947750576	0.9948088692
3	0.9819303885	0.9801400555	0.9993205450	0.9993294807
4	0.9900400584	0.9890380980	0.9998532791	0.9998502524
5	0.9954221934	0.9947546689	0.9999682922	0.9999676828
6	0.9975491011	0.9971518139	0.9999827134	0.9999820036
7	0.9990781402	0.9988937420	0.9999945064	0.9999940851
8	0.9993870184	0.9992443210	0.9999970640	0.9999968682
9	0.9996459298	0.9995575066	0.9999984249	0.9999982012
10	0.9997908626	0.9997329171	0.9999995168	0.9999994521
11	0.9998805961	0.9998491443	0.9999996861	0.9999996457
12	0.9999385575	0.9999185397	0.9999998352	0.9999998220
13	0.9999560220	0.9999403728	0.9999998957	0.9999998900
14	0.9999725492	0.9999616808	0.9999999427	0.9999999450
15	0.9999810290	0.9999733205	0.9999999688	0.9999999681
16	0.9999892611	0.9999845707	0.9999999772	0.9999999772
17	0.9999926070	0.9999891758	0.9999999843	0.9999999857
18	0.9999957257	0.9999935479	0.9999999881	0.9999999904
19	0.9999969535	0.9999953122	0.9999999918	0.9999999940
20	0.9999979067	0.9999968180	0.9999999945	0.9999999967

. It highlights the number of basis functions required for a suitable construction of ROM for each method. The decay of energy for velocity is similar for the uniform and mixed method up to the fifty-first modes, whereas, for the scalar pressure, the two methods differ after the forty-first modes.

From Figure 11b we observe that the accumulated energy in the pressure field converges more rapidly towards 1 (i.e., 100%) than correspondingly for the velocity. Table 1 verify this observation. We may expect that the error in $H_1$-seminorm for velocity and $L_2$-norm is proportional to $ϵ$ if the accumulated energy (measured in the same norms) by the reduced-order modes in the corresponding fields are equal to $1-{ϵ}^2$ [14]. However, Figure 12 shows that the error for pressure does not comply very well with this assertion, only the velocity obtained using the uniform method comply reasonably well. Notice furthermore that the error obtained for the mixed method using nine and more modes is decreasing very slowly (saturation).

4.5. Accuracy of the Mixed and Uniform Methods

In Figure 13 the relative $H_1$-seminorm for velocity and relative $L_2$-error for pressure versus the number of reduced-order modes (equivalent to the number of DOFs in the reduced-order model) are displayed. The relative $H_1$-seminorm error in the velocity decay rapidly towards $5\%$ at 5 modes for both the uniform and the mixed methods. The convergence rate for the uniform method is unchanged for a higher number of modes, whereas the mixed method just slowly reduces to $3\%$ at 20 modes. We observe similar behavior for the relative $L_2$-error in the pressure, but here we see a drastic change in convergence rate starting after using 8 modes (with $6\%$ relative error) and then very slow convergence down to $5\%$ relative error using 20 modes.

In Figure 14 we have reported the relative error in aerodynamic lift and drag versus DOFs obtained by integrating the tractions derived from the computed reduced-order model solutions. The computed drag is quite similar for the uniform and mixed method, whereas the error in the computed lift is nearly one order lower for the uniform than the mixed method. We see that the relative error in both lift and drag is about $3\%$ for the mixed method when using 20 modes.

To explain the saturation of the results obtained using the higher number of modes in the mixed method we report in Figure 15 the relative errors between the POD modes obtained from mixed and uniform methods. We clearly see that for velocities the relative error in each POD mode increases steadily up to 18 modes and then flatten outs, whereas for the pressure it flattens out around 8 modes. As the accumulated energy converge rapidly towards unity for the first 10 modes for the velocity and the first 5 modes for pressure the high relative difference in the POD modes higher than 18 and 8 modes for velocity and pressure respectively, results in saturation at acceptable level (i.e., around 3–$5\%$ relative error in the investigated quantities).

As mentioned in Section 2.2 we must expect an inconsistency error in the mixed method. Even though that the inconsistency error and the error due to the use of a reduced basis are not necessarily orthogonal to each other (in the proper norms) we may expect that relative error of a quantity computed by the presented mixed method to be bounded by the relative difference in the corresponding high-fidelity FE and FV quantity:

$$ϵ⪆\frac{1}{\left|\mathcal{P}\right|}{\int }_{\mathcal{P}}\frac{\parallel \mathrm{FE}\left(\mathit{\mu }\right)-\mathrm{FV}\left(\mathit{\mu }\right)\parallel }{\parallel \mathrm{FE}\left(\mathit{\mu }\right)\parallel }$$

(23)

While the error of the “mixed” methods is bounded, the “uniform” methods continue to converge. However, it must be noted that the errors reported are those between the reduced-order models and the high-fidelity counterparts, which are not necessarily indicative of the true error of the method (as compared to the analytical solution of the PDE). Indeed, one may claim that there is no desire to create reduced models that approximate the high-fidelity model better than that model approximates reality. Since we have that

$$\parallel \mathrm{True}\left(\mathit{\mu }\right)-\mathrm{ROM}\left(\mathit{\mu }\right)\parallel \le \parallel \mathrm{True}\left(\mathit{\mu }\right)-\mathrm{HiFi}\left(\mathit{\mu }\right)\parallel +\parallel \mathrm{HiFi}\left(\mathit{\mu }\right)-\mathrm{ROM}\left(\mathit{\mu }\right)\parallel$$

one should aim for the compromise

$$\parallel \mathrm{HiFi}\left(\mathit{\mu }\right)-\mathrm{ROM}\left(\mathit{\mu }\right)\parallel \approx \parallel \mathrm{True}\left(\mathit{\mu }\right)-\mathrm{HiFi}\left(\mathit{\mu }\right)\parallel$$

Assuming that the two methods given the same number of DOFs are comparable in accuracy, one can expect the following a posteriori error estimate to roughly hold.

$$\parallel {\mathrm{ROM}}_{\mathrm{Uniform}}\left(\mathit{\mu }\right)-{\mathrm{ROM}}_{\mathrm{Mixed}}\left(\mathit{\mu }\right)\parallel \le 2\parallel \mathrm{True}\left(\mathit{\mu }\right)-\mathrm{HiFi}\left(\mathit{\mu }\right)\parallel .$$

Given these considerations, it can be seen that the error bound given by Equation (23) does not impose undue problems for “mixed” methods, and that it is only a natural consequence of the reduction process.

In Figure 16, the parametric dependency of the difference between the snapshots generated by the two methods is investigated. It is observed that the differences are not particularly dependent on $\phi$, but that higher differences are reported for slower flows compared to rapid ones. The same behavior is observed in both velocity and pressure norms. An unsymmetrical shape of the profile is observed for the positive and negative values of $\phi$, which is due to the unsymmetrical shape of the airfoil used in the analysis.

4.6. Computational Speedup

The

Table 2. NACA64: Illustration for the speed gains for mixed and uniform methods. $Speedup$ is, time taken for high-fidelity divided by the time taken for the reduced solution on average. $Relativeerror$ is absolute $H_1$ seminorm error divided by $H_1$ seminorm of the velocity of the reference/high-fidelity solution.

	♯ DoFs	Speedup	Relative Error (Velocity)	Relative Error (Pressure)
High-fidelity	110	1	0	0
Uniform method	5	15,370	$1.0\times {10}^{-1}$	$1.1\times {10}^{-1}$
	10	4122	$3.3\times {10}^{-2}$	$4.9\times {10}^{-2}$
	15	1972	$1.0\times {10}^{-2}$	$1.5\times {10}^{-2}$
	20	1011	$3.1\times {10}^{-3}$	$3.7\times {10}^{-3}$
Mixed method	5	25,981	$1.1\times {10}^{-1}$	$1.2\times {10}^{-1}$
	10	6902	$4.8\times {10}^{-2}$	$8.8\times {10}^{-2}$
	15	2936	$3.5\times {10}^{-2}$	$7.2\times {10}^{-2}$
	20	1764	$3.2\times {10}^{-2}$	$4.2\times {10}^{-2}$

illustrates speedups achieved from ROM using the different number of modes (i.e., DOFs). By employing 5 DOF’s, massive speedup of 25,981 is observed for the mixed method, whereas uniform methods show a speedup of 15,370. However, the accuracy is about 10% in relative $H_1$-seminorm of the velocity, about 10% in relative $L_2$-norm for the pressure and more than 40% relative error in drag and lift. Thus, only 5 reduced-order modes are not acceptable in this case.

If we use 15 modes we observe a speedup of approximately 2000 for the uniform method and close to 3000 for the mixed method. The corresponding relative errors for the four quantities investigated herein are all approximately 5%, which is within the aim for our desired level of accuracy.

5. Conclusions

In this paper, a nonintrusive approach for combining high-fidelity simulations using FV methods with POD and Galerkin Reduced-Order Modeling (ROM) methodology is explored. The key ingredient is to project the FV velocity and pressure results onto suitable FE spaces and then use the classical POD Galerkin ROM framework.

The main finding from our study is that the use of FV high-fidelity snapshots and FE-based ROM, herein denoted mixed method, compute results with an error measured in relative $H_1$-norm for the velocity and relative $L_2$-norm for the pressure as well as relative error in lift and drag that are all about 5%. The corresponding computational speedup of the mixed method to achieve this level of accuracy was 3000 compared to the FV high-fidelity simulation.

The inherent inconsistency in the mixed method due to the projection of FV results onto the FE Taylor–Hood spaces shows up in higher-order modes but does not preclude the aim of reaching accuracy in the range of 1–10% relative error of the typical quantity of interests for laminar flow problems studied herein.

A natural extension of this work is to develop a mixed method for high Reynolds flow problems, and that will be pursued in future studies.

References