Introduction

Sensitivity analysis (SA) studies how changes in the input of a model affect the output. It is essential for many engineering applications such as uncertainty quantification, optimal design, and answering what if questions. Here we consider systems modelled by partial differential equations (PDEs); the sensitivity variable is the derivative of the state (the model output) with respect to the parameters of interest.

For PDEs, one distinguishes differentiate-then-discretise methods (the state model is formally differentiated with respect to the parameter, then the analytical sensitivity system is discretised) from discretise-then-differentiate methods (the two steps swapped; in general not commutative) [71]. This work focuses on the first class, and in particular on the continuous sensitivity equation method [11], [46], [45], [16], [53].

The main aim is to estimate the variance of the solution of the Navier-Stokes equations when there are uncertain parameters, and to use it to compute confidence intervals — forward uncertainty propagation, part of uncertainty quantification (UQ). Many UQ strategies exist (Monte Carlo [146], polynomial chaos [179], [187], [89], [35], random space partition [1]; a review for fluid dynamics is in [178]). SA-based uncertainty propagation is efficient compared to Monte Carlo, but, being based on Taylor expansions, is intrinsically local [33] — usable only for random variables with small variance. Under that hypothesis, this work provides a first-order estimate of the variance of the solution at a reasonable computational cost.

The physical model

A more detailed description can be found in [54].

State equations

Consider a channel $\Omega$ with walls on the top and bottom and a square-section obstacle at distance $x_D$ from the inflow boundary (the boundaries being $\Gamma_{in}$, $\Gamma_{out}$, $\Gamma_{top}$, $\Gamma_{bottom}$ and $\Gamma_{obst}$).

Channel-with-obstacle domain for the first test case: a parabolic inflow profile on \f$\Gamma_{in}\f$, a square obstacle \f$\Gamma_{obst}\f$ of side \f$d\f$ centred at abscissa \f$x_D\f$, channel height \f$\ell\f$ and length \f$L\f$.

The Navier-Stokes system with its boundary conditions is:

\[\begin{cases} \partial_t \mathbf{u} - \nu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p = \mathbf{f} & \Omega,\ t>0,\\ \nabla\cdot\mathbf{u} = 0 & \Omega,\ t>0,\\ \mathbf{u}(\mathbf{x},0)=0 & \Omega,\ t=0,\\ \mathbf{u} = -g(y)\mathbf{n} & \text{on }\Gamma_{in},\\ \mathbf{u} = 0 & \text{on }\Gamma_w=\Gamma_{obst}\cup\Gamma_{top}\cup\Gamma_{bottom},\\ (\nu\nabla\mathbf{u}-pI)\mathbf{n}=0 & \text{on }\Gamma_{out}, \end{cases} \]

where $\mathbf{u}=(u^x,u^y)^t$ is the velocity, $p$ the pressure, $\mathbf{f}$ the external force and $g(y)$ the prescribed inflow condition (no-slip on the walls, prescribed velocity at the inflow, homogeneous Neumann at the outflow).

Sensitivity equations

Considering $\mathbf{u}=\mathbf{u}(\mathbf{x},t;a)$ a function of a scalar uncertain parameter $a$ and writing a formal Taylor expansion $\mathbf{u}(\mathbf{x},t;a+\delta a)=\sum_{k=0}^\infty \mathbf{u}_k(\mathbf{x},t;a)\,\delta a^k$, the coefficient $\mathbf{u}_k=\frac{d^k}{da^k}\mathbf{u}$ is the $k$-th order sensitivity. Substituting into the state system and factorising by powers of $\delta a$ gives the state system for $k=0$ and the first-order sensitivity equations for $k=1$ (denoting $\mathbf{u}_1=\mathbf{u}_a$; only first order is considered here):

\[\begin{cases} \partial_t \mathbf{u}_a - \nu \Delta \mathbf{u}_a + (\mathbf{u}_a\cdot\nabla)\mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u}_a + \nabla p_a = \mathbf{f}_a & \Omega,\ t>0,\\ \nabla\cdot\mathbf{u}_a = 0 & \Omega,\ t>0,\\ \mathbf{u}_a(\mathbf{x},0)=0 & \Omega,\ t=0,\\ \mathbf{u}_a = -g_a(y)\mathbf{n} & \text{on }\Gamma_{in},\\ \mathbf{u}_a = 0 & \text{on }\Gamma_w,\\ (\nu\nabla\mathbf{u}_a-p_a I)\mathbf{n}=0 & \text{on }\Gamma_{out}. \end{cases} \]

These are the Oseen equations [175] (a similar problem, from a numerical point of view, appears in [44], [46] for shape optimisation, where an expansion of the normal $\mathbf{n}$ leads to more complicated boundary conditions).

Note: If $\nu$ is the parameter of interest, the right-hand side of the first equation becomes $\overline{\mathbf{f}}_a:=\mathbf{f}_a+\nu_a\Delta\mathbf{u}$ and the Neumann boundary condition gains the term $\nu_a\nabla\mathbf{u}\,\mathbf{n}$ (not considered in this work).

Non-adiabatic Navier-Stokes equations

A more detailed description can be found in [55].

State equations

On the domain $\Omega=(0,\ell)^2$ with boundary split into $\Gamma_N$ and $\Gamma_D$, the Navier-Stokes equations coupled with the energy conservation law (Boussinesq buoyancy) are:

\[\begin{aligned} \partial_t\mathbf{u}-\nu\Delta\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u}+\frac{1}{\rho}\nabla p &= -\beta(T-T_0)\mathbf{g}, & \Omega,\ t>0,\\ \nabla\cdot\mathbf{u} &= 0, & \Omega,\ t>0,\\ \partial_t T-\frac{\lambda}{\rho_0 c_p}\Delta T+\mathbf{u}\cdot\nabla T &= 0, & \Omega,\ t>0, \end{aligned} \]

with $\mathbf{u}=0$ on $\partial\Omega$, $\mathbf{u}(\mathbf{x},0)=0$, $\nabla T\cdot\mathbf{n}=0$ on $\Gamma_N$ and $T=f(\mathbf{x})$ on $\Gamma_D$, where $\mathbf{g}=(0,-9.81)^t$, $\nu$ the viscosity, $\lambda$ the temperature diffusion parameter, $c_p$ the heat capacity, $\beta$ the thermal expansion coefficient and $\rho$ constant (notation $\lambda'=\frac{\lambda}{\rho c_p}$).

Sensitivity equations

The first-order sensitivity equations (with respect to $a$) are:

\[\begin{aligned} \partial_t\mathbf{u}_a-\nu\Delta\mathbf{u}_a+(\mathbf{u}_a\cdot\nabla)\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u}_a+\frac{1}{\rho}\nabla p_a &= \frac{\rho_a}{\rho^2}\nabla p-\beta_a(T-T_0)\mathbf{g}-\beta(T_a-T_{0,a})\mathbf{g},\\ \nabla\cdot\mathbf{u}_a &= 0,\\ \partial_t T_a-\lambda'\Delta T_a+\mathbf{u}\cdot\nabla T_a &= \lambda'_a\Delta T-\mathbf{u}_a\cdot\nabla T, \end{aligned} \]

with $\mathbf{u}_a=0$ on $\partial\Omega$, $\mathbf{u}_a(\mathbf{x},0)=0$, $\nabla T_a\cdot\mathbf{n}=0$ on $\Gamma_N$, and $T_a=f_a(\mathbf{x})$ on $\Gamma_D$ (the derivative of the Dirichlet temperature condition).

Uncertainty propagation

Here the parameter $a$ is a random variable with known distribution, expected value $\mu_a$ and variance $\sigma_a^2$.

Variance of the model output

For a physical variable $X(\mathbf{x},t;a)$ (a velocity component or the pressure), a Taylor expansion centred at $\mu_a$, $X(\mathbf{x},t;a)=X(\mathbf{x},t;\mu_a)+(a-\mu_a)X_a+o(|a-\mu_a|^2)$ (with $X_a=\partial_a X$ the sensitivity), gives the first-order estimates of the mean and variance:

\[\mu_X(\mathbf{x},t)=E[X]\simeq X(\mathbf{x},t;\mu_a), \qquad \sigma_X^2(\mathbf{x},t)=E[(X-\mu_X)^2]\simeq \sigma_a^2\, X_a^2(\mathbf{x},t;a) \]

valid only for small variances $\sigma_a^2$, but obtained with just one state simulation and one sensitivity simulation — a minimal cost compared to Monte Carlo. With several uncertain parameters, one state simulation plus one sensitivity simulation per parameter is needed [56] (the problem being easily parallelisable).

Confidence intervals

SA does not provide the distribution of the output, so the Chebyshev inequality [83] $P(|X-\mu_X|\geq\lambda)\leq\frac{\sigma_X^2}{\lambda^2}$ is used. Imposing the confidence level $\alpha=\frac{\sigma_X^2}{\lambda^2}$ gives $\lambda=\frac{\sigma_X}{\sqrt{\alpha}}$, hence

\[CI_X=\left[\mu_X-\frac{\sigma_X}{\sqrt{\alpha}},\ \mu_X+\frac{\sigma_X}{\sqrt{\alpha}}\right] \]

(standard choices $\alpha=0.05$ or $0.01$).

Sensitivity analysis using the first-order polynomial chaos method

Further details can be found in [128], [129], [127], [130].

SA methods are classified as local or global. Global SA (variance-based/Sobol, Morris, Monte Carlo) measures sensitivity across the whole parameter space; local SA (Taylor expansion, Polynomial Chaos Method, adjoint) examines sensitivity around a specific parameter set. Here the intrusive Polynomial Chaos Method (PCM) is considered.

In PCM the uncertain variables are decomposed on a basis of complete orthogonal polynomials. The Polynomial Chaos Expansion (PCE) of a random variable $Y(\mathbf{x},t;a)$ is $Y(\mathbf{x},t;a)=\sum_{i=0}^{n}Y_i(\mathbf{x},t)\psi_i(a)$, where the $Y_i$ are deterministic PC coefficients and the $\psi_i$ are random polynomials of degree $i$, orthogonal with respect to $\langle\psi_i,\psi_j\rangle\equiv\int\zeta(a)\psi_i(a)\psi_j(a)\,da=\langle\psi_i,\psi_i\rangle\delta_{ij}$, with $\zeta(a)$ the weighting function. Only first order ( $n=1$) is considered; taking $\zeta(a)$ as the PDF of $a$ gives $\psi_0(a)=1$, $\psi_1(a)=\frac{a-\mu}{\sigma}$. First-order PCE is computationally efficient and sufficient for small input uncertainties (higher-order or non-intrusive methods would be needed for larger uncertainties).

Application to the non-adiabatic Navier-Stokes equations

The velocity, pressure and temperature are expressed by their PCEs (e.g. $\mathbf{v}=\mathbf{v}_0+\frac{a-\mu}{\sigma}\mathbf{v}_1$) and inserted into the non-adiabatic system; multiplying by $\psi_0$ and $\psi_1$ and taking the inner product gives two coupled systems for $i=0$ and $i=1$. The $i=0$ momentum, mass and energy equations couple $\mathbf{u}_0,\mathbf{u}_1$ through cross terms (e.g. $\mathbf{u}_0\cdot\nabla\mathbf{u}_0+\mathbf{u}_1\cdot\nabla\mathbf{u}_1$), and similarly for $i=1$ (e.g. $\mathbf{u}_0\cdot\nabla\mathbf{u}_1+\mathbf{u}_1\cdot\nabla\mathbf{u}_0$).

To make them solvable, a decoupling is performed via the change of variables $\mathbf{u}_1=\sigma\tilde{\mathbf{u}}_1$, $p_1=\sigma\tilde p_1$, $T_1=\sigma\tilde T_1$ (and the data), and $\mathbf{u}_0=\tilde{\mathbf{u}}_0+\sigma\tilde{\mathbf{u}}_1$, etc., yielding a triangular system. The state system (for $\tilde{\mathbf{u}}_0,\tilde p_0,\tilde T_0$) is the non-adiabatic Navier-Stokes system itself:

\[\begin{aligned} \partial_t\tilde{\mathbf{u}}_0-\nu\Delta\tilde{\mathbf{u}}_0+(\tilde{\mathbf{u}}_0\cdot\nabla)\tilde{\mathbf{u}}_0+\rho'\nabla\tilde p_0 &= -\beta(\tilde T_0-\tau)\mathbf{g},\\ \nabla\cdot\tilde{\mathbf{u}}_0 &= 0,\\ \partial_t\tilde T_0-\lambda'\Delta\tilde T_0+\tilde{\mathbf{u}}_0\cdot\nabla\tilde T_0 &= 0, \end{aligned} \]

and the first-order sensitivity system (for $\tilde{\mathbf{u}}_1,\tilde p_1,\tilde T_1$), solvable once the state is known, is:

\[\begin{aligned} \partial_t\tilde{\mathbf{u}}_1-\nu\Delta\tilde{\mathbf{u}}_1+(\tilde{\mathbf{u}}_0\cdot\nabla)\tilde{\mathbf{u}}_1+(\tilde{\mathbf{u}}_1\cdot\nabla)\tilde{\mathbf{u}}_0+2\sigma(\tilde{\mathbf{u}}_1\cdot\nabla)\tilde{\mathbf{u}}_1+\rho'\nabla\tilde p_1 &= -\beta\tilde T_1\mathbf{g},\\ \nabla\cdot\tilde{\mathbf{u}}_1 &= 0,\\ \partial_t\tilde T_1-\lambda'\Delta\tilde T_1+\tilde{\mathbf{u}}_0\cdot\nabla\tilde T_1+\tilde{\mathbf{u}}_1\cdot\nabla\tilde T_0+2\sigma\tilde{\mathbf{u}}_1\cdot\nabla\tilde T_1 &= 0, \end{aligned} \]

with the same boundary/initial conditions as the state, and $\tilde T_1=f_1(\mathbf{x})$ on $\Gamma_D$. The parameters $\nu$, $\beta$, $\lambda'$, $\tau$ can also be considered uncertain, in which case they too are expressed by their PCEs (adding $\tilde\nu_1$, $\tilde\beta_1$, etc. terms); see [130].

Confidence intervals

The PCE gives the mean and variance directly: $\mu_Y(\mathbf{x},t)=Y_0=\tilde Y_0+\sigma\tilde Y_1$ and $\sigma_Y(\mathbf{x},t)=Y_1^2=\sigma^2\tilde Y_1^2$. The Chebyshev inequality then defines the confidence interval as above. This strategy requires only $M+1$ simulations (one state, $M$ sensitivities for $M$ uncertain parameters), the sensitivity simulations being independent and thus easily parallelisable.

Optimisation of a Navier-Stokes problem: adjoint method

The adjoint method efficiently computes gradients of objective functions with respect to many design/control variables [68]. For the incompressible steady-state Navier-Stokes equations with Dirichlet conditions ( $-\nu\Delta\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u}+\nabla p=\mathbf{f}$, $\nabla\cdot\mathbf{u}=0$ in $\Omega$, $\mathbf{u}=\mathbf{u}_D$ on $\partial\Omega$) and an objective function $J(\alpha;\mathbf{u})=\int_\Omega j_\Omega\,dx+\int_{\partial\Omega}j_{\partial\Omega}\,ds$, the adjoint variables $\mathbf{v}$, $q$ (Lagrange multipliers) satisfy the adjoint system [61] :

\[\begin{cases} -\nu\Delta\mathbf{v}+(\mathbf{v}\cdot\nabla^T)\mathbf{u}-(\mathbf{u}\cdot\nabla)\mathbf{v}+\nabla q = -\frac{\partial j_\Omega}{\partial\mathbf{u}}+\nabla\cdot\left(\frac{\partial j_\Omega}{\partial\nabla\mathbf{u}}\right) & \text{in }\Omega,\\ \nabla\cdot\mathbf{v}=0 & \text{in }\Omega,\\ \mathbf{v}=\mathbf{0} & \text{on }\partial\Omega. \end{cases} \]

The gradient is then $\frac{dJ}{d\alpha}=\frac{\partial J}{\partial\alpha}+\langle\mathbf{v},\frac{\partial\mathcal{R}}{\partial\alpha}\rangle$, where $\mathcal{R}$ is the residual of the state equations. The cost of computing the gradient with respect to all parameters is essentially that of a single flow solve, and the gradients are analytical (more accurate and stable than finite differences).

Note: This adjoint formulation is specific to the steady incompressible Navier-Stokes equations, assuming a pressure-independent objective function and Dirichlet boundary conditions only; more general problems require an extended formulation.