Nonintrusive proper generalised decomposition for parametrised   incompressible flow problems in OpenFOAM

Vasileios Tsiolakis; Matteo Giacomini; Ruben Sevilla; Carsten Othmer,; Antonio Huerta

arXiv:1906.05403·math.NA·February 5, 2020

Nonintrusive proper generalised decomposition for parametrised incompressible flow problems in OpenFOAM

Vasileios Tsiolakis, Matteo Giacomini, Ruben Sevilla, Carsten Othmer,, Antonio Huerta

PDF

TL;DR

This paper introduces a nonintrusive proper generalised decomposition (PGD) method integrated into OpenFOAM, enabling efficient parametric flow simulations with minimal online computational cost, suitable for industrial applications.

Contribution

It presents the first nonintrusive PGD implementation in OpenFOAM, combining reduced order modeling with industrial CFD tools for parametrised incompressible flow problems.

Findings

01

Successful application to 2D and 3D laminar flow benchmarks

02

Efficient parametric interpolation without additional solves

03

Demonstrated industrial relevance with automotive flow control

Abstract

The computational cost of parametric studies currently represents the major limitation to the application of simulation-based engineering techniques in a daily industrial environment. This work presents the first nonintrusive implementation of the proper generalised decomposition (PGD) in OpenFOAM, for the approximation of parametrised laminar incompressible Navier-Stokes equations. The key feature of this approach is the seamless integration of a reduced order model (ROM) in the framework of an industrially validated computational fluid dynamics software. This is of special importance in an industrial environment because in the online phase of the PGD ROM the description of the flow for a specific set of parameters is obtained simply via interpolation of the generalised solution, without the need of any extra solution step. On the one hand, the spatial problems arising from the PGD…

Tables1

Table 1. Table 1: Normalised characteristic lengths of the spatial and parametric discretisations.

$h_{x}$	$8.3 \times 10^{- 2}$	$4 \times 10^{- 2}$	$2 \times 10^{- 2}$	$1 \times 10^{- 2}$	$5 \times 10^{- 3}$	$2.5 \times 10^{- 3}$
$h_{μ}$	$2 \times 10^{- 2}$	$1 \times 10^{- 2}$	$5 \times 10^{- 3}$	$2.5 \times 10^{- 3}$	$1.25 \times 10^{- 3}$	$6.25 \times 10^{- 4}$

Equations66

⎩ ⎨ ⎧ \boldmath \nabla \unboldmath \cdot (\boldmath u \unboldmath \otimes \boldmath u \unboldmath) - \boldmath \nabla \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath) + \boldmath \nabla \unboldmath p \boldmath \nabla \unboldmath \cdot \boldmath u \unboldmath \boldmath u \unboldmath \boldmath n \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath - p I_{d}) = \boldmath s \unboldmath = 0 = \boldmath u \unboldmath_{D} = \boldmath t \unboldmath in Ω, in Ω, on Γ_{D}, on Γ_{N},

⎩ ⎨ ⎧ \boldmath \nabla \unboldmath \cdot (\boldmath u \unboldmath \otimes \boldmath u \unboldmath) - \boldmath \nabla \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath) + \boldmath \nabla \unboldmath p \boldmath \nabla \unboldmath \cdot \boldmath u \unboldmath \boldmath u \unboldmath \boldmath n \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath - p I_{d}) = \boldmath s \unboldmath = 0 = \boldmath u \unboldmath_{D} = \boldmath t \unboldmath in Ω, in Ω, on Γ_{D}, on Γ_{N},

⎩ ⎨ ⎧ \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (\boldmath u \unboldmath \otimes \boldmath u \unboldmath) d V - \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath) d V + \int_{V_{i}} \boldmath \nabla \unboldmath p d V \int_{V_{i}} \boldmath \nabla \unboldmath \cdot \boldmath u \unboldmath d V = \int_{V_{i}} \boldmath s \unboldmath d V, = 0.

⎩ ⎨ ⎧ \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (\boldmath u \unboldmath \otimes \boldmath u \unboldmath) d V - \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath) d V + \int_{V_{i}} \boldmath \nabla \unboldmath p d V \int_{V_{i}} \boldmath \nabla \unboldmath \cdot \boldmath u \unboldmath d V = \int_{V_{i}} \boldmath s \unboldmath d V, = 0.

⎩ ⎨ ⎧ \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (\boldmath u \unboldmath \otimes \boldmath u \unboldmath) d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath p d V d \boldmath I \unboldmath \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot \boldmath u \unboldmath d V d \boldmath I \unboldmath = \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath s \unboldmath d V d \boldmath I \unboldmath, = 0.

⎩ ⎨ ⎧ \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (\boldmath u \unboldmath \otimes \boldmath u \unboldmath) d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (ν \boldmath \nabla \unboldmath \boldmath u \unboldmath) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath p d V d \boldmath I \unboldmath \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot \boldmath u \unboldmath d V d \boldmath I \unboldmath = \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath s \unboldmath d V d \boldmath I \unboldmath, = 0.

{\boldmath u \unboldmath_{_{PGD}}^{n} (\boldmath x \unboldmath, \boldmath μ \unboldmath) p_{_{PGD}}^{n} (\boldmath x \unboldmath, \boldmath μ \unboldmath) = \boldmath u \unboldmath_{_{PGD}}^{n - 1} (\boldmath x \unboldmath, \boldmath μ \unboldmath) + σ_{u}^{n} \boldmath f \unboldmath_{u}^{n} (\boldmath x \unboldmath) ϕ^{n} (\boldmath μ \unboldmath), = p_{_{PGD}}^{n - 1} (\boldmath x \unboldmath, \boldmath μ \unboldmath) + σ_{p}^{n} f_{p}^{n} (\boldmath x \unboldmath) ϕ^{n} (\boldmath μ \unboldmath),

{\boldmath u \unboldmath_{_{PGD}}^{n} (\boldmath x \unboldmath, \boldmath μ \unboldmath) p_{_{PGD}}^{n} (\boldmath x \unboldmath, \boldmath μ \unboldmath) = \boldmath u \unboldmath_{_{PGD}}^{n - 1} (\boldmath x \unboldmath, \boldmath μ \unboldmath) + σ_{u}^{n} \boldmath f \unboldmath_{u}^{n} (\boldmath x \unboldmath) ϕ^{n} (\boldmath μ \unboldmath), = p_{_{PGD}}^{n - 1} (\boldmath x \unboldmath, \boldmath μ \unboldmath) + σ_{p}^{n} f_{p}^{n} (\boldmath x \unboldmath) ϕ^{n} (\boldmath μ \unboldmath),

σ_{u}^{n} := ∥ \boldmath f \unboldmath_{u}^{n} ∥ and σ_{p}^{n} := ∥ f_{p}^{n} ∥,

σ_{u}^{n} := ∥ \boldmath f \unboldmath_{u}^{n} ∥ and σ_{p}^{n} := ∥ f_{p}^{n} ∥,

{\boldmath u \unboldmath_{_{PGD}}^{n} p_{_{PGD}}^{n} = t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} + σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} = t o 0.0 pt p \hss p_{_{PGD}}^{n} + σ_{p}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n} = \boldmath u \unboldmath_{_{PGD}}^{n - 1} + σ_{u}^{n} \boldmath f \unboldmath_{u}^{n} ϕ^{n} + σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}, = p_{_{PGD}}^{n - 1} + σ_{p}^{n} f_{p}^{n} ϕ^{n} + σ_{p}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n},

{\boldmath u \unboldmath_{_{PGD}}^{n} p_{_{PGD}}^{n} = t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} + σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} = t o 0.0 pt p \hss p_{_{PGD}}^{n} + σ_{p}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n} = \boldmath u \unboldmath_{_{PGD}}^{n - 1} + σ_{u}^{n} \boldmath f \unboldmath_{u}^{n} ϕ^{n} + σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}, = p_{_{PGD}}^{n - 1} + σ_{p}^{n} f_{p}^{n} ϕ^{n} + σ_{p}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n},

{δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n} := Δ \boldmath f \unboldmath_{u} ϕ^{n} + \boldmath f \unboldmath_{u}^{n} Δ ϕ + Δ \boldmath f \unboldmath_{u} Δ ϕ, := Δ f_{p} ϕ^{n} + f_{p}^{n} Δ ϕ + Δ f_{p} Δ ϕ .

{δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n} := Δ \boldmath f \unboldmath_{u} ϕ^{n} + \boldmath f \unboldmath_{u}^{n} Δ ϕ + Δ \boldmath f \unboldmath_{u} Δ ϕ, := Δ f_{p} ϕ^{n} + f_{p}^{n} Δ ϕ + Δ f_{p} Δ ϕ .

σ_{u}^{n} := ∥ \boldmath f \unboldmath_{u}^{n} + Δ \boldmath f \unboldmath_{u} ∥ and σ_{p}^{n} := ∥ f_{p}^{n} + Δ f_{p} ∥.

σ_{u}^{n} := ∥ \boldmath f \unboldmath_{u}^{n} + Δ \boldmath f \unboldmath_{u} ∥ and σ_{p}^{n} := ∥ f_{p}^{n} + Δ f_{p} ∥.

ν (\boldmath x \unboldmath, \boldmath μ \unboldmath) := ψ (\boldmath μ \unboldmath) D (\boldmath x \unboldmath) = i = 1 \sum n_{ν} ψ_{1, i} (μ_{1}) \dots ψ_{M, i} (μ_{M}) D_{i} (\boldmath x \unboldmath),

ν (\boldmath x \unboldmath, \boldmath μ \unboldmath) := ψ (\boldmath μ \unboldmath) D (\boldmath x \unboldmath) = i = 1 \sum n_{ν} ψ_{1, i} (μ_{1}) \dots ψ_{M, i} (μ_{M}) D_{i} (\boldmath x \unboldmath),

⎩ ⎨ ⎧ \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} ψ \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (D \boldmath \nabla \unboldmath (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n})) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath (σ_{p}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n}) d V d \boldmath I \unboldmath = R_{u}, \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath = R_{p},

⎩ ⎨ ⎧ \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} ψ \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (D \boldmath \nabla \unboldmath (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n})) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath (σ_{p}^{n} δ t o 0.0 pt p \hss p_{_{PGD}}^{n}) d V d \boldmath I \unboldmath = R_{u}, \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (σ_{u}^{n} δ t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath = R_{p},

R_{u} := = R_{u} (\boldmath u \unboldmath_{_{PGD}}^{n - 1}, p_{_{PGD}}^{n - 1}, σ_{u}^{n} \boldmath f \unboldmath_{u}^{n}, σ_{p}^{n} f_{p}^{n}, ϕ^{n}) = R_{u} (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}, t o 0.0 pt p \hss p_{_{PGD}}^{n}) \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath s \unboldmath d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} ψ \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (D \boldmath \nabla \unboldmath t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath t o 0.0 pt p \hss p_{_{PGD}}^{n} d V d \boldmath I \unboldmath,

R_{u} := = R_{u} (\boldmath u \unboldmath_{_{PGD}}^{n - 1}, p_{_{PGD}}^{n - 1}, σ_{u}^{n} \boldmath f \unboldmath_{u}^{n}, σ_{p}^{n} f_{p}^{n}, ϕ^{n}) = R_{u} (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}, t o 0.0 pt p \hss p_{_{PGD}}^{n}) \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath s \unboldmath d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} \otimes t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath + \int_{\boldmath I \unboldmath} ψ \int_{V_{i}} \boldmath \nabla \unboldmath \cdot (D \boldmath \nabla \unboldmath t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) d V d \boldmath I \unboldmath - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath t o 0.0 pt p \hss p_{_{PGD}}^{n} d V d \boldmath I \unboldmath,

R_{p} := = R_{p} (\boldmath u \unboldmath_{_{PGD}}^{n - 1}, σ_{u}^{n} \boldmath f \unboldmath_{u}^{n}, ϕ^{n}) = R_{p} (t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n}) - \int_{\boldmath I \unboldmath} \int_{V_{i}} \boldmath \nabla \unboldmath \cdot t o 0.0 pt \boldmath u \unboldmath \hss \boldmath u \unboldmath_{_{PGD}}^{n} d V d \boldmath I \unboldmath .

(\boldmath v \unboldmath, q) = {(ϕ^{n} δ \boldmath f \unboldmath_{u}, ϕ^{n} δ f_{p}), (σ_{u}^{n} \boldmath f \unboldmath_{u}^{n} δ ϕ, σ_{p}^{n} f_{p}^{n} δ ϕ), for the spatial iteration, for the parametric iteration,

(\boldmath v \unboldmath, q) = {(ϕ^{n} δ \boldmath f \unboldmath_{u}, ϕ^{n} δ f_{p}), (σ_{u}^{n} \boldmath f \unboldmath_{u}^{n} δ ϕ, σ_{p}^{n} f_{p}^{n} δ ϕ), for the spatial iteration, for the parametric iteration,

(\boldmath v \unboldmath, q) = {(ϕ^{n}, ϕ^{n}), (σ_{u}^{n} \boldmath f \unboldmath_{u}^{n}, σ_{p}^{n} f_{p}^{n}), for the spatial iteration, for the parametric iteration .

(\boldmath v \unboldmath, q) = {(ϕ^{n}, ϕ^{n}), (σ_{u}^{n} \boldmath f \unboldmath_{u}^{n}, σ_{p}^{n} f_{p}^{n}), for the spatial iteration, for the parametric iteration .

\left\{\begin{aligned} \int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}{\otimes}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\,dV}{+}\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}\Big{)}\,dV}&\\ -\alpha_{2}\!\!\int_{V_{i}}\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}(D\text{\boldmath$\nabla$\unboldmath}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}))\,dV+\!\alpha_{3}\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}(\sigma_{p}^{n}\varDelta f_{\!p})\,dV}=\,R_{u}^{n},&\\ \alpha_{3}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u})\,dV}=\,R_{p}^{n},&\end{aligned}\right.

\left\{\begin{aligned} \int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}{\otimes}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\,dV}{+}\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}\Big{)}\,dV}&\\ -\alpha_{2}\!\!\int_{V_{i}}\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}(D\text{\boldmath$\nabla$\unboldmath}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}))\,dV+\!\alpha_{3}\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}(\sigma_{p}^{n}\varDelta f_{\!p})\,dV}=\,R_{u}^{n},&\\ \alpha_{3}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u})\,dV}=\,R_{p}^{n},&\end{aligned}\right.

α_{1}^{m} := \int_{\boldmath I \unboldmath} [ϕ^{n}]^{2} ϕ^{m} d \boldmath I \unboldmath, α_{2} := \int_{\boldmath I \unboldmath} [ϕ^{n}]^{2} ψ d \boldmath I \unboldmath, α_{3} := \int_{\boldmath I \unboldmath} [ϕ^{n}]^{2} d \boldmath I \unboldmath .

α_{1}^{m} := \int_{\boldmath I \unboldmath} [ϕ^{n}]^{2} ϕ^{m} d \boldmath I \unboldmath, α_{2} := \int_{\boldmath I \unboldmath} [ϕ^{n}]^{2} ψ d \boldmath I \unboldmath, α_{3} := \int_{\boldmath I \unboldmath} [ϕ^{n}]^{2} d \boldmath I \unboldmath .

\left\{\begin{aligned} \int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}{\otimes}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\,dV}-\alpha_{2}\!\!\int_{V_{i}}\text{\boldmath$\nabla$\unboldmath}{\cdot}(D\text{\boldmath$\nabla$\unboldmath}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}))\,dV\hskip 25.0pt&\\ +\alpha_{3}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}(\sigma_{p}^{n}\varDelta f_{\!p})\,dV}=R_{u}^{n}{-}\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{k-1}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}^{k-1}\Big{)}\,dV}&\\ \alpha_{3}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u})\,dV}=R_{p}^{n},\hskip 198.0pt&\end{aligned}\right.

\left\{\begin{aligned} \int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}{\otimes}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\,dV}-\alpha_{2}\!\!\int_{V_{i}}\text{\boldmath$\nabla$\unboldmath}{\cdot}(D\text{\boldmath$\nabla$\unboldmath}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}))\,dV\hskip 25.0pt&\\ +\alpha_{3}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}(\sigma_{p}^{n}\varDelta f_{\!p})\,dV}=R_{u}^{n}{-}\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{k-1}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u}^{k-1}\Big{)}\,dV}&\\ \alpha_{3}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\varDelta\text{\boldmath$f$\unboldmath}_{\!\!u})\,dV}=R_{p}^{n},\hskip 198.0pt&\end{aligned}\right.

(m = 1 \sum n a_{1}^{m} ϕ^{m} - a_{2} ψ + a_{3}) Δ ϕ = r_{u}^{n} + r_{p}^{n},

(m = 1 \sum n a_{1}^{m} ϕ^{m} - a_{2} ψ + a_{3}) Δ ϕ = r_{u}^{n} + r_{p}^{n},

\!\!\!\!\left\{\begin{aligned} a_{1}^{m}&:=\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\bigl{[}\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\otimes}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m})\bigr{]}\,dV}\\ &\hskip 20.0pt+\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\bigl{[}\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n})\bigr{]}\,dV},\\ a_{2}&:=\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\bigl{[}\text{\boldmath$\nabla$\unboldmath}{\cdot}(D\text{\boldmath$\nabla$\unboldmath}(\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}))\bigr{]}\,dV},\\ a_{3}&:=\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\text{\boldmath$\nabla$\unboldmath}(\sigma_{p}^{n}f_{\!p}^{n})\,dV}+\!\!\int_{V_{i}}{\!\!\!\sigma_{p}^{n}f_{\!p}^{n}\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n})\,dV}.\end{aligned}\right.

\!\!\!\!\left\{\begin{aligned} a_{1}^{m}&:=\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\bigl{[}\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\otimes}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m})\bigr{]}\,dV}\\ &\hskip 20.0pt+\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\bigl{[}\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n})\bigr{]}\,dV},\\ a_{2}&:=\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\bigl{[}\text{\boldmath$\nabla$\unboldmath}{\cdot}(D\text{\boldmath$\nabla$\unboldmath}(\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}))\bigr{]}\,dV},\\ a_{3}&:=\int_{V_{i}}{\!\!\!\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n}{\cdot}\text{\boldmath$\nabla$\unboldmath}(\sigma_{p}^{n}f_{\!p}^{n})\,dV}+\!\!\int_{V_{i}}{\!\!\!\sigma_{p}^{n}f_{\!p}^{n}\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{n}\text{\boldmath$f$\unboldmath}_{\!\!u}^{n})\,dV}.\end{aligned}\right.

σ_{⋄}^{n} < η_{⋄}^{⋆} m = 1 \sum n σ_{⋄}^{m}, for ⋄ = u, p .

σ_{⋄}^{n} < η_{⋄}^{⋆} m = 1 \sum n σ_{⋄}^{m}, for ⋄ = u, p .

\boldmath u \unboldmath (x, y, μ)

\boldmath u \unboldmath (x, y, μ)

p (x, y, μ)

λ (μ) = \frac{1}{2 ν ( μ )} - \frac{1}{( 2 ν ( μ ) ) ^{2}} + 4 π^{2} .

λ (μ) = \frac{1}{2 ν ( μ )} - \frac{1}{( 2 ν ( μ ) ) ^{2}} + 4 π^{2} .

η_{(u, p)} := (\frac{σ _{u}^{n}}{\sum _{m = 1}^{n} σ _{u}^{m}})^{2} + (\frac{σ _{p}^{n}}{\sum _{m = 1}^{n} σ _{p}^{m}})^{2} .

η_{(u, p)} := (\frac{σ _{u}^{n}}{\sum _{m = 1}^{n} σ _{u}^{m}})^{2} + (\frac{σ _{p}^{n}}{\sum _{m = 1}^{n} σ _{p}^{m}})^{2} .

\text{\boldmath$u$\unboldmath}_{\text{jet}}(x,y){=}\!\begin{cases}\!\Bigl{(}{-}1{-}\cos\bigl{(}2\pi y/0.12\bigr{)},0\Bigr{)}&\!\!\text{for $x{=}1$, $y\in[0,0.12]$,}\\ \!\Bigl{(}{-}1{+}\cos\bigl{(}{-}2\pi(y{-}0.88)/0.12\bigr{)},0\Bigr{)}&\!\!\text{for $x{=}1$, $y\in[0.88,1]$,}\\ \!\Bigl{(}{-}1{+}\cos\bigl{(}{-}2\pi(y{-}0.88)/0.12\bigr{)},0\Bigr{)}&\!\!\text{for $x{=}0$, $y\in[0.88,1]$.}\end{cases}

\text{\boldmath$u$\unboldmath}_{\text{jet}}(x,y){=}\!\begin{cases}\!\Bigl{(}{-}1{-}\cos\bigl{(}2\pi y/0.12\bigr{)},0\Bigr{)}&\!\!\text{for $x{=}1$, $y\in[0,0.12]$,}\\ \!\Bigl{(}{-}1{+}\cos\bigl{(}{-}2\pi(y{-}0.88)/0.12\bigr{)},0\Bigr{)}&\!\!\text{for $x{=}1$, $y\in[0.88,1]$,}\\ \!\Bigl{(}{-}1{+}\cos\bigl{(}{-}2\pi(y{-}0.88)/0.12\bigr{)},0\Bigr{)}&\!\!\text{for $x{=}0$, $y\in[0.88,1]$.}\end{cases}

u_{\overset{y}{^}} (\overset{x}{^}, \overset{z}{^}) = 0.0375 (1 - cos (- 2 π \overset{x}{^})) (1 - cos (2 π \overset{z}{^}))

u_{\overset{y}{^}} (\overset{x}{^}, \overset{z}{^}) = 0.0375 (1 - cos (- 2 π \overset{x}{^})) (1 - cos (2 π \overset{z}{^}))

p_{drop} := \frac{1}{A _{in}} i = 1 \sum N_{in} A_{i} p_{i}

p_{drop} := \frac{1}{A _{in}} i = 1 \sum N_{in} A_{i} p_{i}

\displaystyle\begin{aligned} R_{u}^{n}:=&\int_{V_{i}}{\alpha_{4}\text{\boldmath$S$\unboldmath}\,dV}-\sum_{m=1}^{n}\sum_{q=1}^{n}\alpha_{5}^{mq}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{q}\text{\boldmath$f$\unboldmath}_{\!\!u}^{q})\,dV}\\ &+\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}D\text{\boldmath$\nabla$\unboldmath}\Big{(}\sum_{m=1}^{n}\alpha_{6}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\!\Big{)}\,dV}-\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}\Big{(}\sum_{m=1}^{n}\alpha_{7}^{m}\sigma_{p}^{m}f_{\!p}^{m}\Big{)}\,dV},\end{aligned}

\displaystyle\begin{aligned} R_{u}^{n}:=&\int_{V_{i}}{\alpha_{4}\text{\boldmath$S$\unboldmath}\,dV}-\sum_{m=1}^{n}\sum_{q=1}^{n}\alpha_{5}^{mq}\!\!\int_{V_{i}}{\text{\boldmath$\nabla$\unboldmath}{\cdot}(\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}{\otimes}\sigma_{u}^{q}\text{\boldmath$f$\unboldmath}_{\!\!u}^{q})\,dV}\\ &+\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}D\text{\boldmath$\nabla$\unboldmath}\Big{(}\sum_{m=1}^{n}\alpha_{6}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\!\Big{)}\,dV}-\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}\Big{(}\sum_{m=1}^{n}\alpha_{7}^{m}\sigma_{p}^{m}f_{\!p}^{m}\Big{)}\,dV},\end{aligned}

\displaystyle R_{p}^{n}:=-\!\!\int_{V_{i}}{\!\!\text{\boldmath$\nabla$\unboldmath}{\cdot}\Big{(}\sum_{m=1}^{n}\alpha_{7}^{m}\sigma_{u}^{m}\text{\boldmath$f$\unboldmath}_{\!\!u}^{m}\Big{)}\,dV},

⎩ ⎨ ⎧ α_{4} α_{6}^{m} := \int_{\boldmath I \unboldmath} ϕ^{n} η d \boldmath I \unboldmath, := \int_{\boldmath I \unboldmath} ϕ^{n} ϕ^{m} ψ d \boldmath I \unboldmath, α_{5}^{m q} α_{7}^{m} := \int_{\boldmath I \unboldmath} ϕ^{n} ϕ^{m} ϕ^{q} d \boldmath I \unboldmath, := \int_{\boldmath I \unboldmath} ϕ^{n} ϕ^{m} d \boldmath I \unboldmath .

⎩ ⎨ ⎧ α_{4} α_{6}^{m} := \int_{\boldmath I \unboldmath} ϕ^{n} η d \boldmath I \unboldmath, := \int_{\boldmath I \unboldmath} ϕ^{n} ϕ^{m} ψ d \boldmath I \unboldmath, α_{5}^{m q} α_{7}^{m} := \int_{\boldmath I \unboldmath} ϕ^{n} ϕ^{m} ϕ^{q} d \boldmath I \unboldmath, := \int_{\boldmath I \unboldmath} ϕ^{n} ϕ^{m} d \boldmath I \unboldmath .

r_{u}^{n} := a_{4} η + m = 1 \sum n (- q = 1 \sum n a_{5}^{m q} ϕ^{q} + a_{6}^{m} ψ - a_{7}^{m}) ϕ^{m},

r_{u}^{n} := a_{4} η + m = 1 \sum n (- q = 1 \sum n a_{5}^{m q} ϕ^{q} + a_{6}^{m} ψ - a_{7}^{m}) ϕ^{m},

r_{p} := - m = 1 \sum n a_{8}^{m} ϕ^{m},

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Nonintrusive proper generalised decomposition for parametrised incompressible flow problems in OpenFOAM

V. Tsiolakis111Volkswagen AG, Brieffach 011/1777, D-38436, Wolfsburg, Germany ,222Laboratori de Càlcul Numèric (LaCàN), ETS de Ingenieros de Caminos, Canales y Puertos, Universitat Politècnica de Catalunya, Barcelona, Spain ,333Zienkiewicz Centre for Computational Engineering, College of Engineering, Swansea University, Wales, UK

Corresponding author: Matteo Giacomini. E-mail: [email protected] , M. Giacomini222Laboratori de Càlcul Numèric (LaCàN), ETS de Ingenieros de Caminos, Canales y Puertos, Universitat Politècnica de Catalunya, Barcelona, Spain, R. Sevilla333Zienkiewicz Centre for Computational Engineering, College of Engineering, Swansea University, Wales, UK

Corresponding author: Matteo Giacomini. E-mail: [email protected] , C. Othmer111Volkswagen AG, Brieffach 011/1777, D-38436, Wolfsburg, Germany and A. Huerta222Laboratori de Càlcul Numèric (LaCàN), ETS de Ingenieros de Caminos, Canales y Puertos, Universitat Politècnica de Catalunya, Barcelona, Spain

Abstract

The computational cost of parametric studies currently represents the major limitation to the application of simulation-based engineering techniques in a daily industrial environment. This work presents the first nonintrusive implementation of the proper generalised decomposition (PGD) in OpenFOAM, for the approximation of parametrised laminar incompressible Navier-Stokes equations. The key feature of this approach is the seamless integration of a reduced order model (ROM) in the framework of an industrially validated computational fluid dynamics software. This is of special importance in an industrial environment because in the online phase of the PGD ROM the description of the flow for a specific set of parameters is obtained simply via interpolation of the generalised solution, without the need of any extra solution step. On the one hand, the spatial problems arising from the PGD separation of the unknowns are treated using the classical solution strategies of OpenFOAM, namely the semi-implicit method for pressure linked equations (SIMPLE) algorithm. On the other hand, the parametric iteration is solved via a collocation approach. The resulting ROM is applied to several benchmark tests of laminar incompressible Navier-Stokes flows, in two and three dimensions, with different parameters affecting the flow features. Eventually, the capability of the proposed strategy to treat industrial problems is verified by applying the methodology to a parametrised flow control in a realistic geometry of interest for the automotive industry.

Keywords: Reduced order models, proper generalised decomposition, finite volume, incompressible laminar Navier-Stokes, pressure Poisson equation, parametrised flows, OpenFOAM, nonintrusiveness

1 Introduction and motivations

Computational fluid dynamics (CFD) is a key component in the current industrial design pipeline. Simulations of incompressible flows are performed on a daily basis to solve different problems both in automotive and aeronautical industries. Owing to its robustness, the most widely spread CFD methodology is the finite volume (FV) method [1, 2, 3, 4, 5, 6]. Using this technique, numerically evaluated quantities of interest (e.g. drag and lift) have proved to match reasonably well experimental results.

Nonetheless, design and optimisation cycles in a production environment require multiple queries of the same problem with boundary conditions, physical properties of the fluid and geometry of the domain varying within a range of values of interest. In this context, parameters act as extra-coordinates of a high-dimensional partial differential equation (PDE). The computational cost of such parametric studies currently represents the major limitation to the application of simulation-based engineering techniques in a daily industrial environment. It is well-known that the computational complexity of approximating the PDEs describing the problems under analysis increases exponentially with the number of parameters considered. In recent years, reduced order models [7], including reduced basis (RB) [8, 9, 10, 11, 12, 13, 14], proper orthogonal decomposition (POD) [15, 16, 17, 18, 19, 20, 21] and hierarchical model reduction (HiMod) [22, 23, 24, 25], have been proposed to reduce the computational burden of parametric analysis for several physical problems, including incompressible flows. The aforementioned techniques rely on an a posteriori reduction based on snapshots computed as solutions of the full-order model for different values of the parameters under analysis. An alternative approach is represented by PGD [26, 27, 28, 29, 30]. This method features an a priori reduction [31, 32, 33], using a separable approximation of the solution, which depends explicitly on the parameters under analysis. In this context, during an offline phase, a reduced basis is constructed with no a priori knowledge of the solution, whereas efficient online evaluations of the generalised solution are performed by simple interpolation in the parametric space. The PGD framework has been first applied to incompressible Navier-Stokes equations in [34] to separate $x$ and $y$ directions in two-dimensional problems and in [35, 36] to separate space and time discretisations of unsteady flows. See also [37, 38, 39, 40, 41, 42] for several applications of PGD to different physical problems.

In the context of flow problems, model reduction techniques based on Galerkin projection have been extensively studied in the literature [43, 44, 45]. In this framework, several strategies have been proposed to construct the trial basis, using POD [46], RB [47] or the empirical interpolation method [48]. Concerning incompressible Navier-Stokes equations, in [49, 50] supremiser stabilisations techniques have been investigated to couple the FV method with POD to solve parametrised turbulent flow problems. Alternative projection methods based on minimisation of the residual of the momentum equation only [51] and on a least-squares Petrov-Galerkin approach [52, 53, 54] have been proposed. More recently, special attention has also been devoted to FV-based structure-preserving ROMs for conservation laws [55].

Another key aspect for the application of simulation-based techniques to industrial problems is the capability of the proposed methods to provide verified and certified results. This problem has been classically treated by equipping numerical methods with reliable and fully-computable a posteriori error estimators using equilibrated fluxes [56, 57, 58] and flux-free approaches [59, 60, 61] to control the error of the solution as well as of quantities of interest [62, 63, 64, 65, 66, 67, 68, 69]. Nonetheless, these approaches require intrusive modifications of existing computational libraries and may not be feasible in the context of commercial software. Hence, although the effort of the academic community in this direction, such solutions have not been successfully and widely integrated in codes utilised by the industry. More recently, to circumvent this issue, great effort has been devoted to nonintrusive implementations in which novel numerical methodologies are externally coupled to existing commercial and open-source software used in industry on a daily basis. Some contributions in this direction have been successfully proposed coupling PGD with Abaqus® for mechanical problems [70] and PGD with SAMTECH® for shape optimisation problems [71]. For flow problems, the coupling of POD and OpenFOAM has been discussed in [72, 50]. The present contribution is the first nonintrusive integration of the PGD framework in OpenFOAM for the solution of parametrised incompressible Navier-Stokes problems in the laminar regime. The resulting algorithm, henceforth referred to as pgdFoam, relies on internal functions and routines of OpenFOAM [73] and exploits the incompressible flow solver simpleFoam for the spatial iteration of the alternating direction scheme.

The rest of this paper is organised as follows. Section 2 recalls the incompressible Navier-Stokes equations and their FV approximation. The parametrised Navier-Stokes equations are introduced in Section 3 as well as their PGD approximation and its nonintrusive implementation in OpenFOAM. Numerical simulations to validate the discussed reduced-order strategy are provided in Section 4, whereas its application to parametrised flow control problems is presented in Section 5. Section 6 summarises the results and two appendices complement the information with some technical details on the formulation and the OpenFOAM spatial solver utilised.

2 Finite volume approximation of the incompressible Navier-Stokes equations

In this section, the steady Navier-Stokes equations for the simulation of incompressible viscous laminar flows in $d$ spatial dimensions are recalled. Let $\Omega\subset\mathbb{R}^{d},\ \partial\Omega{=}\Gamma_{D}\cup\Gamma_{N}$ be an open bounded domain with disjoint Dirichlet, $\Gamma_{D}$ , and Neumann, $\Gamma_{N}$ , boundaries. The flow problem under analysis consists of computing the velocity field $u$ and the pressure $p$ such that

[TABLE]

where the first equation describes the balance of momentum and the second one the conservation of mass. In Equation (1), $s$ represents a volumetric source term, $\nu{>}0$ is the dynamic viscosity and $\mathbf{I}_{d}$ is the $d{\times}d$ identity matrix. On the Dirichlet boundary $\Gamma_{D}$ , the value $\text{\boldmath$ u $\unboldmath}_{D}$ of the velocity is imposed, whereas on $\Gamma_{N}$ the pseudo-traction $t$ is applied. From the modelling point of view, inlet surfaces and physical walls are described as Dirichlet boundaries with an imposed entering velocity profile and a homogeneous datum, respectively, whereas outlet surfaces feature homogeneous Neumann boundary conditions. For the sake of simplicity and without loss of generality, $\Gamma_{N}$ is henceforth assumed to be an outlet boundary, that is a null $t$ is considered.

2.1 A cell-centred finite volume approximation using OpenFOAM

In this section, the formulation of a FV scheme for the incompressible Navier-Stokes equations is briefly recalled to introduce the notation needed for the high-dimensional parametrised problem of Section 3. The domain $\Omega$ is partitioned in $N$ nonoverlapping cells $V_{i},\,i{=}1,\ldots,N$ such that $\Omega{:=}\bigcup_{i=1}^{N}V_{i}$ and $V_{i}{\cap}V_{j}{=}\emptyset,\,\text{for}\,i{\neq}j$ . The FV discretisation is constructed starting from the integral formulation of Equation (1), namely find $(\text{\boldmath$ u $\unboldmath},p)$ , constant on each cell $V_{i},\,i{=}1,\dots,N$ , such that $\text{\boldmath$ u $\unboldmath}{=}\text{\boldmath$ u $\unboldmath}_{D}\,\text{on}\,\Gamma_{D}$ and it holds

[TABLE]

OpenFOAM implements a cell-centred FV rationale in which piecewise constant approximations are sought for velocity and pressure in each cell of the computational mesh and the degrees of freedom of the discretised problem are located at the centroid of each finite volume. Employing Gauss’s theorem, the integrals in Equation (2) are rewritten in terms of fluxes over the boundaries of the cells and approximated using classical central differencing schemes [4, 5]. Moreover, to handle the nonlinearity in the convection term, OpenFOAM considers a relaxation approach introducing a fictitious time variable. The resulting solution strategy relies on the SIMPLE algorithm which belongs to the family of fractional-step projection methods [74, 75]. A brief description of this method is provided in B.

3 Nonintrusive proper generalised decomposition for parametrised flow problems

Consider now the case in which the user-prescribed data in Equation (1), i.e. the viscosity coefficient, the source term and the boundary conditions, depend on a set of parameters $\text{\boldmath$ \mu $\unboldmath}\in\text{\boldmath$ \mathcal{I} $\unboldmath}\subset\mathbb{R}^{M}$ , with $M$ being the number of parameters. More presicely, the set $\mathcal{I}$ describing the range of admissible parameters can be defined as the Cartesian product of the domains of the $M$ parameters, namely, $\text{\boldmath$ \mathcal{I} $\unboldmath}{:=}\mathcal{I}_{1}\times\mathcal{I}_{2}\times\dotsb\times\mathcal{I}_{M}$ with $\mu_{i}\in\mathcal{I}_{i}$ for $i{=}1,\dots,M$ . Within this context, $\mu$ is treated as a set of additional independent variables, or parametric coordinates, instead of problem parameters. For the purpose of discretisation, each interval $\mathcal{I}_{i}$ is subdivided in $N_{\mu}$ subintervals. The unknown pair $(\text{\boldmath$ u $\unboldmath},p)$ is thus sought in a high-dimensional space described by the independent variables $(\text{\boldmath$ x $\unboldmath},\text{\boldmath$ \mu $\unboldmath})\in\Omega{\times}\text{\boldmath$ \mathcal{I} $\unboldmath}$ and fulfils the following parametrised Navier-Stokes equations on each cell $V_{i},\,i{=}1,\ldots,N$

[TABLE]

In the following sections, the rationale for the construction of a separated solution of the parametrised Navier-Stokes equations is recalled and the proposed nonintrusive implementation of the alternating direction scheme in OpenFOAM is presented.

3.1 The proper generalised decomposition rationale

PGD constructs an approximation $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n})$ of the solution $(\text{\boldmath$ u $\unboldmath},p)$ of Equation (3) in terms of a sum of $n$ separable functions, or modes. Each mode is the product of functions depending solely on one of the arguments $\text{\boldmath$ x $\unboldmath},\mu_{1},\dots,\mu_{M}$ . For the sake of readability and without loss of generality, only the spatial coordinates $x$ and the parametric ones $\mu$ are henceforth separated.

Following [41], the so-called single parameter approximation is detailed. That is, for each mode, a unique scalar parametric function $\phi(\text{\boldmath$ \mu $\unboldmath})$ is considered for all the variables and the resulting separated form of the unknowns is

[TABLE]

where the superindex $n$ denotes the, a priori unknown, number of terms in the PGD expansion and the positive scalar coefficients $\sigma_{u}^{n}$ and $\sigma_{p}^{n}$ represent the amplitude of the $n$ -th mode for velocity and pressure, respectively. These coefficients are obtained normalising the modal functions, namely

[TABLE]

with $\|\phi^{n}\|=1$ . Appropriate user-defined norms on the spatial and parametric domains need to be introduced for each function. For all the simulations in Section 4 and 5, the $\mathcal{L}_{2}$ norm has been considered for normalisation.

Remark 1.

The normalisation coefficients play a critical role in checking the convergence of the PGD algorithm and may be used as quantitative stopping criterion in the PGD enrichment procedure described in Section 3.2.

For a discussion on alternative formulations of the separation in Equation (4), involving both scalar and vector-valued parametric functions, the interested reader is referred to [41]. Henceforth and except in case of ambiguity, the dependence of the modes on $x$ and $\mu$ is omitted.

Considering a linearised approach to compute each new mode, Equation (4) can be rewritten as the following predictor-corrector single parameter approximation

[TABLE]

where $\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n}{:=}\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n-1}+\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n}\phi^{n}$ and $\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}p_{{}_{\texttt{PGD}}}^{\>n}{:=}p_{{}_{\texttt{PGD}}}^{n-1}+\sigma_{p}^{n}f_{\!p}^{n}\phi^{n}$ account for the $n{-}1$ previously computed terms and a prediction of the current mode. More precisely, $(\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n}\phi^{n},\sigma_{p}^{n}f_{\!p}^{n}\phi^{n})$ play the role of predictors in the computation of the $n$ -th mode, whereas $(\sigma_{u}^{n}\delta\!\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n},\sigma_{p}^{n}\delta\!\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}p_{{}_{\texttt{PGD}}}^{\>n})$ are the corresponding correctors featuring the variations $\varDelta$ in the spatial and parametric functions, namely

[TABLE]

Note that the last term in Equation (6) represents a high-order variation which is henceforth neglected. As for the classical single parameter approximation, $\sigma_{u}^{n}$ and $\sigma_{p}^{n}$ represent the amplitudes of the $n$ -th velocity and pressure modes. That is, setting $\|\phi^{n}+\varDelta\phi\|=1$ , they are defined as

[TABLE]

3.2 Predictor-corrector alternating direction scheme

In order to compute $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n})$ in Equation (5), a greedy algorithm is implemented. The first PGD mode $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{0},p_{{}_{\texttt{PGD}}}^{0})$ is arbitrarily chosen to fulfil the Dirichlet boundary conditions of the problem and the $n$ -th mode is successively computed assuming that term $n{-}1$ is available [29, 30]. Some variations of this strategy based on Arnoldi-type iterations have been investigated in [76, 77]. In this section, the alternating direction scheme used to compute the PGD modes is described. A key assumption for the application of this method is the separability of the data. For the sake of simplicity and without any loss of generality, the separated form of the viscosity coefficient, see e.g. [38], is reported

[TABLE]

and analogous separations are considered for all the parametric data in the problem under analysis.

By plugging (5) into (3) and gathering the unknown increments $(\sigma_{u}^{n}\delta\!\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n},\sigma_{p}^{n}\delta\!\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{p}} $}\hss}}p_{{}_{\texttt{PGD}}}^{\>n})$ on the left-hand side while leaving on the right-hand side the residuals computed using the previous modes $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n-1},p_{{}_{\texttt{PGD}}}^{n-1})$ and the predictions $(\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n}\phi^{n},\sigma_{p}^{n}f_{\!p}^{n}\phi^{n})$ of the current one, the following equations are obtained

[TABLE]

where the residuals are defined as

[TABLE]

As classical in ROMs [78, 79], an affine dependence of the forms in (8), (9) and (10) on the parameters is required to construct the PGD approximation. The spatial (respectively, parametric) component of each mode is thus computed by restricting Equation (8) to the tangent manifold associated with the spatial (respectively, parametric) coordinate. Following from Equation (6) and setting a fixed value for the parametric function $\phi^{n}$ , the pair $(\sigma_{u}^{n}\varDelta\text{\boldmath$ f $\unboldmath}_{\!\!u},\sigma_{p}^{n}\varDelta f_{\!p})$ is determined by solving a purely spatial PDE. Recall that the PGD alternating direction scheme handles homogeneous Dirichlet boundary conditions at each iteration of the spatial solver [30], whereas inhomogeneous data are treated by the first arbitrary PGD mode $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{0},p_{{}_{\texttt{PGD}}}^{0})$ introduced above. In a similar fashion, the increment $\varDelta\phi$ is computed as the solution of an algebraic system of equations in the parameter $\mu$ while the spatial functions $(\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n},\sigma_{p}^{n}f_{\!p}^{n})$ are considered known.

Note that at each iteration of the alternating direction scheme, $\mathchoice{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \displaystyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \textstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}{\hbox to0.0pt{\raisebox{-0.86108pt}{$ \scriptscriptstyle\kern 0.0pt\widetilde{\phantom{\text{\boldmath $u$ \unboldmath}}} $}\hss}}\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n}$ is known and may be expressed in separated form as $\sum_{m=1}^{n}\sigma_{u}^{m}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{m}\phi^{m}$ . Thus, exploiting the separated structure of the unknowns and the affine parametric decomposition of the involved integral forms, the numerical complexity of the high-dimensional PDE is reduced by alternatively solving for the spatial and the parametric unknowns, as detailed in the next subsections.

Remark 2.

By restricting Equation (8) to the tangent manifold in the spatial (respectively, parametric) direction, the integral forms are multiplied by $\phi^{n}$ (respectively, $(\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n},\sigma_{p}^{n}f_{\!p}^{n})$ ). This is equivalent to the projection of the high-dimensional PDE to the tangent manifold discussed for PGD in the context of finite element approximations [41]. More precisely, being $\text{\boldmath$ v $\unboldmath}(\text{\boldmath$ x $\unboldmath},\text{\boldmath$ \mu $\unboldmath})$ (respectively, $q(\text{\boldmath$ x $\unboldmath},\text{\boldmath$ \mu $\unboldmath})$ ) the test function in the finite element weak form of the momentum (respectively, continuity) equation, the projection on the tangent manifold leads to

[TABLE]

where $(\delta\!\text{\boldmath$ f $\unboldmath}_{\!\!u},\delta\!f_{\!p})$ and $\delta\!\phi$ are test functions depending solely on the spatial and parametric variables, respectively. In the framework of FV discretisations, these test functions are set equal to $1$ to retrieve the classical integral formulation of the PDE under analysis. The corresponding restriction to the tangent manifold thus leads to the following functions multiplying the integral forms in Equations (8), (9) and (10)

[TABLE]

3.2.1 The spatial iteration

First, the parametric function $\phi^{n}$ is fixed and the increments $(\sigma_{u}^{n}\varDelta\text{\boldmath$ f $\unboldmath}_{\!\!u},\sigma_{p}^{n}\varDelta f_{\!p})$ are determined by solving a spatial PDE. More precisely, restricting Equation (8) to the tangent manifold in the spatial direction and neglecting the high-order terms, a pair $(\sigma_{u}^{n}\varDelta\text{\boldmath$ f $\unboldmath}_{\!\!u},\sigma_{p}^{n}\varDelta f_{\!p})$ , constant element-by-element, is sought such that in each cell $V_{i},\,i{=}1,\ldots,N$ it holds

[TABLE]

where $R_{u}^{n}$ and $R_{p}^{n}$ are the spatial residuals associated with the discretisation of the momentum and mass equations, respectively, and each coefficient $\alpha_{i},\,i{=}1,\ldots,3$ depends solely on the parametric function $\phi^{n}$ and on the data of the problem

[TABLE]

Note that given the separable form of (9)-(10), an efficient implementation of the right-hand side of the spatial iteration may be devised and the corresponding FV discretisation is obtained. A detailed description of the residuals acting as linear functionals on the right-hand side of Equation (11) is provided in A.

The terms in Equation (11) feature a structure similar to the original incompressible Navier-Stokes problem in the spatial domain $\Omega$ , see Equation (2). The discretisation is thus performed using the cell-centred FV method implemented in OpenFOAM, see Section 2.1. The main difference is represented by the first two integrals on the left-hand side of the momentum equation in (11).

On the one hand, the first integral is a relaxation of the classical nonlinear convection term in Navier-Stokes equations, based on the last computed approximation $\sum_{m=1}^{n}\!\alpha_{1}^{m}\sigma_{u}^{m}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{m}$ of the unknown velocity field. From a practical point of view, this treatment of the convection term is equivalent to the one performed by the SIMPLE algorithm which solves a linearised version of the Navier-Stokes equations by introducing a fictitious time variable $t^{k}$ and substituting the unknown convection field with its approximation at time $t^{k-1}$ , see B.

On the other hand, the second integral does not appear in the classical Navier-Stokes equations and is thus not treated by the SIMPLE algorithm. In order to preserve the nonintrusiveness of the discussed PGD approach, the standard solution strategy implemented in OpenFOAM for such problem, namely simpleFoam, is applied. Hence, a relaxation is introduced in the SIMPLE iterations and this term is handled explicitly as part of the right-hand side of the momentum equation, leading to

[TABLE]

where the index $k{-}1$ is now associated with the lastly computed increment $\sigma_{u}^{k-1}\varDelta\text{\boldmath$ f $\unboldmath}_{\!\!u}^{k-1}$ in the SIMPLE algorithm, see B.

3.2.2 The parametric iteration

In the parametric step, the value of the previously computed spatial functions is fixed $(\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n},f_{\!p}^{n}){\leftarrow}(\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n}{+}\varDelta\text{\boldmath$ f $\unboldmath}_{\!\!u},\sigma_{p}^{n}f_{\!p}^{n}{+}\varDelta f_{\!p})$ and the parametric increment $\varDelta\phi$ acts as unknown. Within the single parameter approximation rationale, a unique scalar function depending on $\mu$ is sought. Following the strategy described for the spatial iteration, the high-order terms are neglected in the restriction of Equation (8) to the parametric direction of the tangent manifold and $\varDelta\phi$ is computed by solving the following algebraic equation

[TABLE]

where $r_{u}^{n}$ and $r_{p}^{n}$ are the parametric residuals associated with the discretisation of the momentum and mass equations, respectively, and each coefficient $a_{i},\,i{=}1,\ldots,3$ depends solely on the spatial functions $(\sigma_{u}^{n}\text{\boldmath$ f $\unboldmath}_{\!\!u}^{n},\sigma_{p}^{n}f_{\!p}^{n})$ and on the data of the problem, namely

[TABLE]

The unknown $\varDelta\phi$ is discretised at the nodes of the parametric domain $\mathcal{I}$ and the resulting algebraic equation is solved via a collocation method. Similarly to the spatial iteration, the separable form of (9)-(10) is exploited to perform computationally efficient pointwise evaluations of the residuals at the nodes of $\mathcal{I}$ . The complete derivation of the separated form of the right-hand side is detailed in A.

3.3 A nonintrusive implementation of the proper generalised decomposition in OpenFOAM

In order to construct an efficient PGD strategy applicable to engineering problems of interest for the industry, a critical aspect is its nonintrusiveness with respect to the OpenFOAM solving procedure simpleFoam. As discussed in Section 3.2, inhomogeneous Dirichlet boundary conditions are treated by means of a spatial mode computed using the full-order solver, whereas the corresponding parametric mode is set equal to $1$ (Algorithm 1 - Step 1). Then, the enrichment process is started and at each iteration of the alternating direction scheme a spatial mode is computed using simpleFoam (Algorithm 1 - Steps 7 to 10) and a linear system is solved to determine the corresponding parametric term (Algorithm 1 - Steps 11 to 14). The alternating direction iterations stop when the computed corrections $\varDelta\text{\boldmath$ f $\unboldmath}_{\!\!\diamond}$ , $\varDelta\phi$ are negligible with respect to the amplitudes $\sigma_{\diamond}^{n}$ , $\sigma_{\phi}$ of the current mode for $\diamond=u,p$ and the residuals $\varepsilon_{\circ}^{r}$ are sufficiently small for $\circ=u,p,\phi$ (Algorithm 1 - Steps 6 and 15). The global enrichment strategy ends when the amplitude of the current mode $\sigma_{\diamond}^{n}$ is negligible with respect to the first one $\sigma_{\diamond}^{1}$ for $\diamond=u,p$ (Algorithm 1 - Step 3). The complete flowchart of pgdFoam is displayed in Figure 1.

Remark 3.

Alternative criterions may be considered to stop the greedy algorithm, e.g. when the magnitude of the last mode normalised with respect to the sum of the amplitudes of all the computed terms is lower than a user-defined tolerance $\eta_{\diamond}^{\star}$ , namely

[TABLE]

4 Numerical validation

In this section, numerical examples are presented to validate the proposed methodology. First, a test case with known analytical solution is considered to verify the optimal convergence rate of the high-dimensional FV approximation of the velocity and pressure fields, measured in the $\mathcal{L}_{2}(\Omega{\times}\text{\boldmath$ \mathcal{I} $\unboldmath})$ norm, for a parametrised viscosity coefficient. In this context, special emphasis is given to the additional error introduced by the PGD, highlighting the range of applicability of the discussed reduced-order strategy in terms of expected accuracy of the parametric solution. Moreover, a classical benchmark test for incompressible flow solvers, namely the nonleaky lid-driven cavity, is studied parametrising the imposed velocity of the lid in a range of values of the Reynolds number spanning from $1,\!000$ to $4,\!000$ .

4.1 Kovasznay flow with parametrised viscosity

Consider the Kovasznay flow [80] for a parametrised viscosity $\nu(\mu){=}\mu$ . The analytical solution is

[TABLE]

where the constant $C$ is determined by fixing a reference value for the pressure field in one point of the domain, whereas $\lambda$ is a function of the parametrised viscosity and changes when the Reynolds number is modified, namely,

[TABLE]

The parameter $\mu$ is sought in the space $\mathcal{I}{=}[5{\times}10^{-3},10^{-2}]$ , which is discretised with uniform intervals. The corresponding values of the Reynolds number span from $100$ to $200$ . The spatial domain $\Omega{=}[-1,1]^{2}$ is discretised with a family of Cartesian meshes of quadrilateral cells. The characteristic lengths $h_{x}$ and $h_{\mu}$ of the spatial and parametric discretisations, respectively, are provided in Table 1.

A convergence study under uniform mesh refinement is performed for the linearised Navier-Stokes equations using the meshes described in Table 1. In this context, a convective field $a$ given by the analytical expression of the Kovasznay velocity is introduced in Equation (1) and the convective term $\text{\boldmath$ \nabla $\unboldmath}{\cdot}(\text{\boldmath$ u $\unboldmath}\otimes\text{\boldmath$ u $\unboldmath})$ is replaced by $\text{\boldmath$ \nabla $\unboldmath}{\cdot}(\text{\boldmath$ u $\unboldmath}\otimes\text{\boldmath$ a $\unboldmath})$ . As detailed in Section 3.2, an affine separation of the data is required to run PGD. Thus, the convective field $a$ is separated a priori considering the first four terms of the Taylor expansion of $e^{\lambda x}$ in the analytical form of the velocity, see Equation (16). For $\mu{=}10^{-2}$ , the relative $\mathcal{L}_{2}(\Omega)$ error of the resulting separated velocity field with respect to the exact one is $4.3{\times}10^{-3}$ and, consequently, a target error of $10^{-2}$ in the spatial discretisation is considered for the following convergence study. Moreover, the Dirichlet boundary datum $\text{\boldmath$ u $\unboldmath}_{D}$ requires five modes to be described in a separated form.

The $\mathcal{L}_{2}(\Omega{\times}\text{\boldmath$ \mathcal{I} $\unboldmath})$ error between the PGD approximation $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n})$ computed using fifteen modes and the high-dimensional analytical solution $(\text{\boldmath$ u $\unboldmath},p)$ as a function of the characteristic mesh size $h_{x}$ is displayed in Figure 2.

The optimal first-order convergence rate for pressure and second-order one for velocity are obtained.

To run the iterative procedure pgdFoam described in Algorithm 1, a stopping criterion $\eta_{(u,p)}\leq 10^{-5}$ is considered, where $\eta_{(u,p)}$ accounts for the relative amplitude of both the velocity and pressure modes, namely

[TABLE]

In Figure 3(a), the evolution of the amplitude $\eta_{(u,p)}$ , $\eta_{u}$ and $\eta_{p}$ is displayed for the finest mesh described in Table 1. After ten computed modes, the stopping criterion is fulfilled and the PGD enrichment stops.

As previously mentioned, five terms are required to describe the Dirichlet boundary conditions in a separated form. Henceforth, only the computed modes, starting from the sixth term of the PGD approximation are displayed. In Figures 3(b), the first six normalised computed parametric modes are displayed. The corresponding computed spatial modes for pressure and velocity are presented in Figure 5 and 5, respectively.

The PGD approximation $(\text{\boldmath$ u $\unboldmath}_{{}_{\texttt{PGD}}}^{n},p_{{}_{\texttt{PGD}}}^{n})$ using $n{=}6,8,15$ , that is with $1$ , $3$ and $10$ computed modes, respectively, is compared to the analytical solution for the case of $\text{Re}{=}200$ , in Figure 6.

4.2 Two-dimensional cavity with parametrised lid velocity

In this section, the classical benchmark problem of the nonleaky lid-driven cavity is studied [81]. The unitary square $\Omega{=}[0,1]^{2}$ is considered as spatial domain and homogeneous Dirichlet boundary conditions are imposed on the lateral and bottom walls. On the top wall, a velocity $\text{\boldmath$ u $\unboldmath}_{\text{lid}}(\text{\boldmath$ x $\unboldmath},\mu){=}400\mu\text{\boldmath$ u $\unboldmath}_{\text{lid}}(\text{\boldmath$ x $\unboldmath})$ is enforced, where the parameter $\mu\in[0.25,1]$ acts as a scaling factor of the maximum velocity of the lid, whereas $\text{\boldmath$ u $\unboldmath}_{\text{lid}}(\text{\boldmath$ x $\unboldmath})$ is a velocity profile featuring two ramps on the top-left and top-right corners of the domain to account for the change between null and maximum velocity. As classical in the literature treating the lid-driven cavity example, for $x\in[0,0.06]$ and $x\in[0.94,1]$ , the horizontal component of the lid velocity changes linearly from [math] to $400\mu$ and vice versa. The dynamic viscosity is set to $\nu{=}$ 0.1\text{,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$ and the values considered for the Reynolds number span from $1,\!000$ to $4,\!000$ .

The nonlinear term of the Navier-Stokes equations is now treated as described in Section 3.2. The mode handling the boundary conditions is obtained as a full-order solution of the Navier-Stokes equations using the simpleFoam algorithm for a lid velocity computed using $\mu{=}1$ , that is for a maximum horizontal velocity of $400\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ . The corresponding parametric boundary condition mode is set to be linearly evolving from $\mu{=}0.25$ to $\mu{=}1$ , that is $\phi(\mu){=}\mu$ .

Following the rationale described in the previous section, two different stopping criterions are considered for the PGD enrichment strategy, namely $\eta_{(u,p)}\leq 10^{-3}$ and $\eta_{(u,p)}\leq 10^{-4}$ . Figure 7(a) displays the relative amplitude of the computed modes. Note that the first stopping point is achieved after seven computed modes, whereas seventeen terms are required to fulfil the lower tolerance.

The corresponding computed parametric modes are presented on Figure 7(b). It is worth noting that all the computed parametric modes are close or equal to [math] for $\mu{=}1$ . This is due to the fact that the boundary conditions of the problem are imposed by means of a full-order solution computed for the maximum value of $\mu$ in the parametric space. Hence, the case of $\mu{=}1$ is accurately described by the PGD approximation using solely the mode obtained via simpleFoam.

Now, the online evaluations of the PGD approximation of the velocity and pressure fields for different values of the parameter $\mu$ are compared to the full-order solutions computed using simpleFoam. The corresponding relative $\mathcal{L}_{2}(\Omega)$ errors are presented in Figure 8 as a function of the number of modes utilised in the PGD approximation.

The first three computed modes, for which $\eta_{(u,p)}\leq 10^{-2}$ , provide a good approximation of both velocity and pressure and limited corrections are introduced by the following modes until the stopping criterion of $10^{-3}$ is fulfilled at the vertical dotted line. For the case $\mu{=}1$ , a small error of the order of $10^{-4}$ appears starting from the fifth computed mode, i.e. $n{=}6$ . This is due to the fact that the boundary condition mode already captures all the features of the flow, being a full-order solution of the Navier-Stokes equations as previously mentioned.

A qualitative comparison of the reduced-order and full-order solutions of the parametrised lid-driven cavity problem is displayed in Figure 9.

Using the first seven computed modes, the PGD approximations for $\mu{=}0.25$ , $\mu{=}0.625$ and $\mu{=}1$ are presented as long as their corresponding simulations obtained using simpleFoam. The cases under analysis correspond to a maximum horizontal velocity of the lid of $100$ , $250$ and $400\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ , respectively. It is worth noting that pgdFoam is able to capture the topological changes of the flow with great accuracy, managing to identify location and size of the vortices, as well as their appearance and disappearance according to the values of the Reynolds number considered in the analysis.

5 Application to parametrised flow control problems

Dynamically controlling the features of a flow is a challenging problem with several high-impact applications including, e.g., drag minimisation, stall control and aerodynamic noise reduction [82, 83, 84]. A major bottleneck to the design of flow control devices is represented by the large number of simulations involved in the tuning of the control loop. In this section, the potential of the described nonintrusive PGD implementation in OpenFOAM is demonstrated for parametrised flow control problems. Two- and three-dimensional internal flows with blowing jets are studied. Specifically, a parametric study involving the peak velocity of the jets as extra-coordinate of the problem is considered to test the proposed PGD methodology.

5.1 Lid-driven cavity with parametrised jet velocity

Consider the nonleaky lid-driven cavity problem introduced in Section 4.2. The lid velocity is defined with two linear ramps, increasing from [math] to $10\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ on the top-left corner and decreasing correspondingly on the top-right one. Three jets of size $0.12\text{\,}\mathrm{m}$ are introduced on the vertical walls, two on the right wall and one on the left, respectively. The parametrised velocity of the jets is $\text{\boldmath$ u $\unboldmath}_{\text{jet}}(\text{\boldmath$ x $\unboldmath},\mu){=}\mu\text{\boldmath$ u $\unboldmath}_{\text{jet}}(\text{\boldmath$ x $\unboldmath})$ , where the maximum velocity is controlled by the parameter $\mu\in[0,1]$ and the profile $\text{\boldmath$ u $\unboldmath}_{\text{jet}}(\text{\boldmath$ x $\unboldmath})$ is defined as

[TABLE]

An outlet boundary is added on the left vertical wall for $y\in[0,0.12]$ and a free-traction condition is enforced. The dynamic viscosity is set to $\nu{=}$ 0.01\text{,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$, therefore the corresponding Reynolds number is $\text{Re}{=}1,\!000$ .

The boundary conditions of the problem are enforced through two modes computed as full-order solutions via simpleFoam as shown in Figure 10: the first one, for $\mu{=}0$ , corresponds to lid velocity of $10\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ and inactive jets; the second one, for $\mu{=}1$ , is associated with the maximum velocity of the jets and a zero velocity of the lid. The corresponding parametric modes for the boundary conditions are set to $\phi(\mu){=}1$ and $\phi(\mu){=}\mu$ , respectively.

Following the rationale previously discussed, the PGD enrichment process is stopped when $\eta_{(u,p)}\leq 10^{-4}$ . In Figure 11, the generalised solution computed by the PGD is particularised for several values of the parameter under analysis and compared with the corresponding full-order solutions provided by simpleFoam.

The flows for $\mu{=}0.1$ , $\mu{=}0.3$ , $\mu{=}0.5$ , $\mu{=}0.7$ , $\mu{=}0.8$ and $\mu{=}1$ are displayed, covering a wide range of flow regimes in the cavity. It is worth noting that the discussed reduced-order strategy is able to capture the topological changes in the flow features and accurately reproduce the appearance and disappearance of vortices in different regions of the domain.

The accuracy of the PGD approximation with respect to the full-order solution is also verified by computing the relative $\mathcal{L}_{2}(\Omega)$ error of the spatial discretisation while enriching the modal description of the solution. Specifically, Figure 12 shows that using seven computed modes all approximations present relative errors lower than $10^{-2}$ .

5.2 S-Bend with flow control driven by a jet

In this section, the proposed PGD methodology is applied to a flow control problem using a three-dimensional geometry of industrial interest. The model of a heating, ventilation and air conditioning (HVAC) duct section provided by Volkswagen AG is shown in Figure 13.

A jet is introduced on the red patch, at the first bend of the duct. The velocity profile of the jet is a sinusoidal function defined on the reference planar square $[0,1]^{2}$ as

[TABLE]

and pointing in the direction $\hat{y}$ orthogonal to the plane $(\hat{x},\hat{z})$ . The parametrisation is constructed as a scaling of the jet velocity from $u_{y}{=}$ -0.015\text{,}\mathrm{m}\mathrm{/}\mathrm{s} $, i.e. blowing, to suction with $u_{y}{=}$0.15\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ . A single parameter $\mu$ is introduced and the parametric domain considered for the analysis is $\mathcal{I}{=}[-0.1,1]$ . Note that this problem is particularly challenging due to the change of sign in the interval of parametric values considered leading to different physical phenomenons. The remaining boundary conditions feature homogeneous velocity on all the lateral walls, a parabolic velocity profile with mean value $\text{\boldmath$ u $\unboldmath}{=}(0.83,0,0)$ \text{,}\mathrm{m}\mathrm{/}\mathrm{s} $on the inlet and a free-traction on the outlet. The dynamic viscosity is set to $\nu{=}1.588{\times}10^{-5}\,$\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ and the corresponding value of the Reynolds number is $\text{Re}{=}280$ . The quantity of interest in this problem is the pressure drop computed along the duct.

As previously done for the lid-driven cavity with jets, two modes to account for the boundary conditions are computed using simpleFoam. The first mode is a full-order solution corresponding to the case of inactive jet and given inlet parabolic profile; the second one, is obtained setting a zero inlet velocity and a jet of maximum velocity $u_{y}{=}$ 0.15\text{,}\mathrm{m}\mathrm{/}\mathrm{s}$$. The corresponding parametric modes are $\phi(\mu){=}1$ and $\phi(\mu){=}\mu$ , respectively.

Setting a tolerance of $10^{-3}$ , pgdFoam computes three modes before fulfilling the stopping criterion for $\eta_{(u,p)}$ , see Equation (17), as displayed in Figure 14(a).

The PGD approximation obtained using three computed modes is compared with the full-order solutions given by simpleFoam for the values $\mu{=}-0.1$ , $\mu{=}0.45$ and $\mu{=}1$ of the parameter under analysis. In Figure 14(b), the relative $\mathcal{L}_{2}(\Omega)$ error for these configurations is reported. The numerical experiments confirm that an accuracy of $10^{-2}$ is achieved using one computed mode additionnally to the two terms accounting for the boundary conditions. It is worth noting that the first computed mode is two orders of magnitude more relevant than the following ones (Fig. 14(a)). Thus, after one computed mode, the additional terms only introduce limited corrections to the existing PGD approximation.

A qualitative comparison of the pressure and velocity fields computed using the PGD solution particularised for different values of the parameter $\mu$ and the corresponding full-order discretisations is presented in Figures 16 and 16.

As mentioned at the beginning of this section, the quantity of engineering interest in the analysis of this problem is the pressure drop computed along the duct. The weighted average pressure drop is defined as

[TABLE]

where $A_{\text{in}}$ is the area of the inlet surface, $N_{\text{in}}$ the number of faces $S_{i},\,i{=}1,\ldots,N_{\text{in}}$ on the inlet patch and $p_{i}$ , $A_{i}$ are the pressure and area on the face $S_{i}$ , respectively. For $\mu{=}-0.1$ , $\mu{=}0.1$ , $\mu{=}0.45$ , $\mu{=}0.55$ and $\mu{=}1$ , the pressure drop is evaluated as a particularisation of the generalised PGD solution and using the full-order solver simpleFoam. Figure 17(a) presents the convergence history of the error in the pressure drop as a function of the number of modes in the PGD approximation.

It is straightforward to observe again that using the modes accounting for the boundary conditions and one computed mode is sufficient to capture the flow features of a wide range of parameters. Moreover, by comparing the pressure drop with respect to the maximum velocity of the jet for different configurations with the corresponding values provided by the full-order solver, the capability of the discussed reduced-order strategy to accurately capture the evolution of a quantity of interest throughout the range of values of the parameter $\mu$ is confirmed (Fig. 17(b)).

6 Concluding remarks

A nonintrusive PGD implementation in OpenFOAM has been proposed in the context of parametrised incompressible laminar flows. The main novelty of such approach is represented by the seamless exploitation of OpenFOAM native SIMPLE solver, making the resulting reduced-order strategy suitable for application in a daily industrial environment. The pgdFoam algorithm relies on the industrially-validated solver simpleFoam to compute the spatial modes of the solution, whereas the parametric ones are determined via the solution of a linear system of algebraic equations.

The developed strategy has been validated using a manufactured solution to verify the optimal order of convergence of the PGD-FV approximation and a classical benchmark test case in the literature of CFD techniques for incompressible flows. Moreover, the potential of the proposed PGD approach to rapidly and accurately simulate incompressible flows for different sets of user-defined parameters has been tested in the context of flow control problems. The pgdFoam algorithm has been applied both to an academic test case and an industrial one with a 3D geometrical model provided by Volkswagen AG.

The proposed PGD methodology has proved to be able to compute an accurate reduced basis for the problems under analysis with no a priori knowledge of the expected solutions. Moreover, it has shown robustness when dealing with a large range of values of the parameters, accuracy in capturing significant topological changes in the flow features and reliability in evaluating quantities of engineering interest, with an extremely reduced computing time.

Acknowledgements

This work was partially supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions (Grant agreement No. 675919) that financed the Ph.D. fellowship of the first author. The second, third and last author were also supported by the Spanish Ministry of Economy and Competitiveness (Grant agreement No. DPI2017-85139-C2-2-R). The second and last authors are grateful for the financial support provided by the Generalitat de Catalunya (Grant agreement No. 2017-SGR-1278).

Appendix A Separated representation of the residuals

Consider a separable expression of the source term $\text{\boldmath$ s $\unboldmath}(\text{\boldmath$ x $\unboldmath},\text{\boldmath$ \mu $\unboldmath}){:=}\eta(\text{\boldmath$ \mu $\unboldmath})\text{\boldmath$ S $\unboldmath}(\text{\boldmath$ x $\unboldmath})$ . For the spatial iteration, the residuals in separated form read as

[TABLE]

where the following expressions for the coefficients are devised

[TABLE]

For the parametric iteration, the separated expression of the residuals is

[TABLE]

where the coefficients depend solely on the spatial modes, namely

[TABLE]

Appendix B simpleFoam: the semi-implicit method for pressure linked equations in OpenFOAM

In OpenFOAM, the steady Navier-Stokes equations are approximated by means of an iterative procedure, namely simpleFoam. This algorithm implements the SIMPLE method proposed in [74]. SIMPLE is a fractional-step Chorin-Temam projection method [85] that has been extensively studied in the literature [86, 87]. First, an intermediate velocity $\text{\boldmath$ u $\unboldmath}^{k}$ is computed starting from the momentum equation and neglecting the contribution of pressure, see Equation (25a); second, the step involving the incompressibility constraint is rewritten in terms of a Poisson equation for the pressure $p$ , see Equation (25b); eventually, a correction is applied to the intermediate velocity field to determine the final value $u$ in Equation (25c). Special attention is required to impose the correct set of boundary conditions in each step of the algorithm [88].

[TABLE]

Note that the algorithm in Equation (25) may also be rewritten in the framework of an algebraic splitting method [89]. For a complete introduction to the subject, interested readers are referred to [75].

Bibliography89

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] R. J. Le Veque, Finite volume methods for hyperbolic problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2002.
2[2] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics, 3rd Edition, Springer-Verlag, Berlin, 2009, a practical introduction.
3[3] K. W. Morton, T. Sonar, Finite volume methods for hyperbolic conservation laws, Acta Numer. 16 (2007) 155–238.
4[4] T. Barth, R. Herbin, M. Ohlberger, Finite Volume Methods: Foundation and Analysis, in: Encyclopedia of Computational Mechanics Second Edition, American Cancer Society, 2017, pp. 1–60.
5[5] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, Handbook of Numerical Analysis 7 (2000) 713 – 1018, solution of Equation in ℝ n superscript ℝ 𝑛 \mathbb{R}^{n} (Part 3), Techniques of Scientific Computing (Part 3).
6[6] R. Sevilla, M. Giacomini, A. Huerta, A face-centred finite volume method for second-order elliptic problems, Int. J. Numer. Methods Eng. 115 (8) (2018) 986–1014.
7[7] F. Chinesta, A. Huerta, G. Rozza, K. Willcox, Model Reduction Methods, in: E. Stein, R. de Borst, T. J. R. Hughes (Eds.), Encyclopedia of Computational Mechanics Second Edition, Vol. Part 1 Solids and Structures, John Wiley & Sons, Ltd., Chichester, 2017, Ch. 3, pp. 1–36.
8[8] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations, C. R. Acad. Sci. Ser. I-Math. 339 (9) (2004) 667 – 672.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Nonintrusive proper generalised decomposition for parametrised incompressible flow problems in OpenFOAM

Abstract

1 Introduction and motivations

2 Finite volume approximation of the incompressible Navier-Stokes equations

2.1 A cell-centred finite volume approximation using OpenFOAM

3 Nonintrusive proper generalised decomposition for parametrised flow problems

3.1 The proper generalised decomposition rationale

Remark 1**.**

3.2 Predictor-corrector alternating direction scheme

Remark 2**.**

3.2.1 The spatial iteration

3.2.2 The parametric iteration

3.3 A nonintrusive implementation of the proper generalised decomposition in OpenFOAM

Remark 3**.**

4 Numerical validation

4.1 Kovasznay flow with parametrised viscosity

4.2 Two-dimensional cavity with parametrised lid velocity

5 Application to parametrised flow control problems

5.1 Lid-driven cavity with parametrised jet velocity

5.2 S-Bend with flow control driven by a jet

6 Concluding remarks

Acknowledgements

Appendix A Separated representation of the residuals

Appendix B simpleFoam: the semi-implicit method for pressure linked equations in OpenFOAM

Remark 1.

Remark 2.

Remark 3.