The PPE formulation: Difference between revisions

From openpipeflow.org
Jump to navigation Jump to search
(Created page with "{{latexPreamble}} The incompressibility condition may be 'replaced' by specifying an equation for the pressure. Taking the divergence of the Navier--Stokes equation, we may ...")
 
mNo edit summary
 
(20 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{latexPreamble}}
{{latexPreamble}}
The incompressibility condition may be 'replaced' by specifying an equation for the pressure - the Pressure-Poisson Equation (PPE).  Taking the divergence of the Navier--Stokes equation leads to the system


The incompressibility condition may be 'replaced' by specifying an equation for the pressure.  Taking the divergence of the Navier--Stokes equation, we may write
$\label{eq:NSeqs}
(1) \qquad
\left\{\begin{array}{rcl}
\partial_t \vec{u} & = & L \,\vec{u} + \vec{N} - \bnabla p \, , \\
\nabla^2 p & = & \bnabla\cdot\vec{N} ,
\end{array}\right.$


$\label{eq:NSeqs}
where $L$ is a linear operator and $\vec{N}$ represents nonlinear terms.  The no-slip boundary conditions are $\vec{u}=\vec{0}$,
        \left\{\begin{array}{rcl}
    \partial_t \vec{u} & = & L \,\vec{u} + \vec{N} - \bnabla p \, , \\
                \nabla^2 p & = & \bnabla\cdot\vec{N} ,
        \end{array}\right.$
where $L$ is a linear operator and $\vec{N}$ represents nonlinear terms.  The boundary conditions are $\vec{u}=\vec{0}$,
and we retain that $\bnabla\cdot\vec{u}=0$ must be satisfied everywhere, i.e. also on the boundary.  There is no boundary condition explicitly on the pressure.
and we retain that $\bnabla\cdot\vec{u}=0$ must be satisfied everywhere, i.e. also on the boundary.  There is no boundary condition explicitly on the pressure.
(See [[Equations_and_parameters#Boundary_conditions|Boundary_conditions]].)


== PPE-formulation with correct boundary conditions ==
== PPE-formulation with correct boundary conditions ==
Line 15: Line 18:
Consider the case where the spatial discretisation splits the Navier--Stokes equations into a set one-dimensional problems, here into problems for radially-dependent Fourier modes.
Consider the case where the spatial discretisation splits the Navier--Stokes equations into a set one-dimensional problems, here into problems for radially-dependent Fourier modes.


Let $\vec{u}$ denote the velocity, $\vec{N}$ denote nonlinear terms, and $\mat{X}$ and $\mat{Y}$ be matrices associated with implicit timestepping of the viscous terms for a particular mode.
Let $\vec{u}$ denote the velocity, $\vec{N}$ denote nonlinear terms, and $\mat{X}$ and $\mat{Y}$ be matrices associated with implicit timestepping of the viscous terms for a particular Fourier mode.
The time-discretised Navier–Stokes equations for this mode may be written in the form
The time-discretised Navier–Stokes equations for this mode may be written in the form


$\label{eq:NSdisc}
$\label{eq:NSdisc}
        \left\{\begin{array}{rcl}
(2) \qquad
    \mat{X}\, \vec{u}^{q+1} & = & \mat{Y}\, \vec{u}^q  
\left\{\begin{array}{rcl}
                + \vec{N}^{q+\frac{1}{2}} - \bnabla p \, , \\
\mat{X}\, \vec{u}^{q+1} & = & \mat{Y}\, \vec{u}^q  
                \nabla^2 p & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q  
+ \vec{N}^{q+\frac{1}{2}} - \bnabla p \, , \\
                + \vec{N}^{q+\frac{1}{2}}) ,
\nabla^2 p & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q  
        \end{array}\right.$
+ \vec{N}^{q+\frac{1}{2}}) ,
\end{array}\right.$


where $q$ denotes time $t_q$, which is sixth order in $r$ for $\vec{u}^{q+1}$ and second order for $p$, where the solenoidal condition is not explicitly imposed. Symmetry conditions provide the conditions at the axis. The difficulty is in imposing the remaining four — this system should be inverted, in principle, simultaneously for $p$ and $\vec{u}^{q+1}$ with boundary conditions $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$ (Rempfer 2006). In practice it would be preferable to invert for $p$ first then for $\vec{u}^{q+1}$, but the boundary conditions to not involve $p$ directly.
where $q$ denotes time $t_q$, which is sixth order in $r$ for $\vec{u}^{q+1}$ and second order for $p$, where the solenoidal condition is not explicitly imposed. Symmetry provides the conditions at the axis. The difficulty is in imposing the remaining four boundary conditions — this system should be inverted, in principle, simultaneously for $p$ and $\vec{u}^{q+1}$ with boundary conditions $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$ (Rempfer 2006). In practice it would be preferable to invert for $p$ first then for $\vec{u}^{q+1}$, but the boundary conditions to not involve $p$ directly.


Note that the $\mat{Y}\,\vec{u}^q$ term has been included in the right-hand side of the pressure-Poisson equation, the divergence of which should be small. Assume that pressure boundary condition is known: the right-hand side of the Navier–Stokes equation is then projected onto the space of solenoidal functions though $p$ and hence after inversion, $\vec{u}^{q+1}$ will be solenoidal.
Note that the $\mat{Y}\,\vec{u}^q$ term has been included in the right-hand side of the pressure-Poisson equation, the divergence of which should be small. Assume that pressure boundary condition is known: the right-hand side of the Navier–Stokes equation is then projected onto the space of solenoidal functions though $p$ and hence after inversion, $\vec{u}^{q+1}$ will be solenoidal.
Line 33: Line 37:


$\label{eq:NSbulk}
$\label{eq:NSbulk}
        \left\{\begin{array}{rcl}
(3) \qquad
    \mat{X}\, \bar{\vec{u}} & = & \mat{Y}\, \vec{u}^q  
\left\{\begin{array}{rcl}
                + \vec{N}^{q+\frac{1}{2}} - \bnabla \bar{p} \, , \\
\mat{X}\, \bar{\vec{u}} & = & \mat{Y}\, \vec{u}^q  
                \nabla^2 \bar{p} & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q  
+ \vec{N}^{q+\frac{1}{2}} - \bnabla \bar{p} \, , \\
                + \vec{N}^{q+\frac{1}{2}}) ,
\nabla^2 \bar{p} & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q  
        \end{array}\right.$
+ \vec{N}^{q+\frac{1}{2}}) ,
\end{array}\right.$


with boundary conditions $\bar{\vec{u}}=\vec{0}$ and $\partial_{r}\bar{p}=0$. Introduce the following systems:
with boundary conditions $\bar{\vec{u}}=\vec{0}$ and $\partial_{r}\bar{p}=0$. Introduce the following systems:


$\label{eq:NSp0}
$\label{eq:NSp0}
        \left\{\begin{array}{rcl}
(4) \qquad
          \vec{u}' & = & -\bnabla p' \, , \\
\left\{\begin{array}{rcl}
                \nabla^2 p' & = & 0,
\mat{X}\,\vec{u}' & = & -\bnabla p' \, , \\
        \end{array}\right.$
\nabla^2 p' & = & 0,
\end{array}\right.$


with boundary condition $\partial_{r}p'=1$ on $r=R$, and
with boundary conditions $\vec{u}'=\vec{0}$ and $\partial_{r}p'=1$ on $r=R$, and


$\label{eq:NSu0}
$\label{eq:NSu0}
        \left\{\begin{array}{rcl}
(5) \qquad
    \mat{X}\, \vec{u}' & = & \vec{0},
\left\{\begin{array}{rcl}
        \end{array}\right.$
\mat{X}\, \vec{u}' & = & \vec{0},
\end{array}\right.$


with boundary conditions $u'_+=1$, $u'_-=1$, $u'_z=\mathrm{i}$ on $r=R$ (see [[Equations_and_parameters#Decoupling_the_equations|Decoupling_the_equations]]). The system \ref{eq:NSp0} provides a linearly independent function $\vec{u}'_1$ that may be added to $\bar{\vec{u}}$ without affecting the right-hand side in \ref{eq:NSbulk}, but altering (to correct) the boundary condition applied. ($\mat{X}$ is dropped on the left-hand side of \ref{eq:NSp0} since $\bnabla^2(\bnabla p')=\bnabla(\nabla^2 p')-\bnabla \wedge\bnabla \wedge\bnabla p'=0$.) Similarly the system \ref{eq:NSu0} provides a further three functions $\vec{u}'_j$ (where $j=2,3,4$; the +,- and $z$ components may be considered separately). The superposition
with boundary conditions $u'_+=1$, $u'_-=1$, $u'_z=\mathrm{i}$ on $r=R$ (see [[Equations_and_parameters#Decoupling_the_equations|Decoupling_the_equations]]). The system $(4)$ provides a linearly independent function $\vec{u}'_1$ that may be added to $\bar{\vec{u}}$ without affecting the right-hand side in $(3)$, but altering (to correct) the boundary condition applied.  
<!--
($\mat{X}$ has been dropped on the left-hand side of $(4)$ since $\mat{X}$ represents a combination of the identity and the Laplace operator, but $\bnabla^2(\bnabla p')=\bnabla(\nabla^2 p')-\bnabla \wedge\bnabla \wedge\bnabla p'=0$, i.e., only the identity is left behind.)
-->
Similarly the system $(5)$ provides a further three functions $\vec{u}'_j$ (where $j=2,3,4$; the +,- and $z$ components may be considered separately). The superposition


$\label{eq:usuperpos}
$\label{eq:usuperpos}
  \vec{u}^{q+1} = \bar{\vec{u}} + \sum_{j=1}^4 a_j\, \vec{u}'_j \,$
(6) \qquad
\vec{u}^{q+1} = \bar{\vec{u}} + \sum_{j=1}^4 a_j\, \vec{u}'_j \,$


may be formed in order to satisfy the four boundary conditions, $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$. Substituting \ref{eq:usuperpos} into the boundary conditions, they may be written
may be formed in order to satisfy the four boundary conditions, $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$. Let $\vec{g}(\vec{u})$ be a 4-vector composed of these boundary conditions, such that they are satisfied when $\vec{g}(\vec{u})=\vec{0}$.  Substituting the superposition $(6)$ into the boundary conditions, they may be written  


$\label{eq:velBCs}
$\label{eq:velBCs}
  \mat{A}\,\vec{a} = -\vec{g}(\bar{\vec{u}}) ,$
(7) \qquad
\mat{A}\,\vec{a} = -\vec{g}(\bar{\vec{u}}) ,$


where $\mat{A}=\mat{A}(\vec{g}(\vec{u}'))$ is a 4$\times$4 matrix. The appropriate coefficients required to satisfy the boundary conditions are thereby recovered from solution of this small system for $\vec{a}$.  If needed, the pressure may be calculated  
where $\mat{A}=\mat{A}(\vec{g}(\vec{u}'))$ is a 4$\times$4 matrix and the 4-vector $\vec{a}$ is composed of the $a_j$. Thus, the appropriate coefficients required to satisfy the boundary conditions are recovered from solution of this small system for $\vec{a}$.   
using the same $a_j$:
 
The error in the boundary conditions $g_j(\vec{u}^{q+1})$ using the influence-matrix technique is at the level of the machine epsilon, typically order 1e-16.
The functions $u'_j(r)$, the matrix $\mat{A}$ and its inverse may all be precomputed. The boundary conditions for $\vec{u}'$ have been chosen so that that $u'_\pm$ are pure real, $u'_z$ is pure imaginary, and $\mat{A}$ is real. For each timestep, this application of the influence matrix technique requires only evaluation of the deviation from the boundary condition, multiplication by a 4$\times$4 real matrix, and the addition of only two functions to each component of $\vec{u}$, each either pure real or pure imaginary. Compared to the evaluation of nonlinear terms, the computational overhead is negligible.
 
If needed, the pressure may be calculated  
using the same $a_1$ as in $(6)$ using the $p'$ from $(4)$:


$\label{eq:pressure}   
$\label{eq:pressure}   
  p^{q+1} = \bar{p} + \sum_{j=1}^4 a_j\, p'_j \, .$
(8a) \qquad
\tilde{p} = \bar{p} + a_1\, p' \, \qquad$ (for each Fourier mode, then)


$
(8b) \qquad p = \tilde{p} - \frac{1}{2}\vec{u}\cdot\vec{u}\, \qquad$ (in physical space).


The error in the boundary conditions $g_j(\vec{u}^{q+1})$ using the influence-matrix technique is at the level of the machine epsilon, typically order 1e-18.
The adjustment in $(8b)$ arises when the Navier-Stokes equations are in rotational form (see [[Equations_and_parameters]]).
The functions $u'_j(r)$, the matrix $\mat{A}$ and its inverse may all be precomputed. The boundary conditions for $\vec{u}'$ have been chosen so that that $u'_\pm$ are pure real, $u'_z$ is pure imaginary, and $\mat{A}$ is real. For each timestep, this application of the influence matrix technique requires only evaluation of the deviation from the boundary condition, multiplication by a 4$\times$4 real matrix, and the addition of only two functions to each component of $\vec{u}$, each either pure real or pure imaginary. Compared to the evaluation of nonlinear terms, the computational overhead is negligible.

Latest revision as of 08:12, 10 March 2020

$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $ The incompressibility condition may be 'replaced' by specifying an equation for the pressure - the Pressure-Poisson Equation (PPE). Taking the divergence of the Navier--Stokes equation leads to the system

$\label{eq:NSeqs} (1) \qquad \left\{\begin{array}{rcl} \partial_t \vec{u} & = & L \,\vec{u} + \vec{N} - \bnabla p \, , \\ \nabla^2 p & = & \bnabla\cdot\vec{N} , \end{array}\right.$

where $L$ is a linear operator and $\vec{N}$ represents nonlinear terms. The no-slip boundary conditions are $\vec{u}=\vec{0}$, and we retain that $\bnabla\cdot\vec{u}=0$ must be satisfied everywhere, i.e. also on the boundary. There is no boundary condition explicitly on the pressure. (See Boundary_conditions.)


PPE-formulation with correct boundary conditions

Consider the case where the spatial discretisation splits the Navier--Stokes equations into a set one-dimensional problems, here into problems for radially-dependent Fourier modes.

Let $\vec{u}$ denote the velocity, $\vec{N}$ denote nonlinear terms, and $\mat{X}$ and $\mat{Y}$ be matrices associated with implicit timestepping of the viscous terms for a particular Fourier mode. The time-discretised Navier–Stokes equations for this mode may be written in the form

$\label{eq:NSdisc} (2) \qquad \left\{\begin{array}{rcl} \mat{X}\, \vec{u}^{q+1} & = & \mat{Y}\, \vec{u}^q + \vec{N}^{q+\frac{1}{2}} - \bnabla p \, , \\ \nabla^2 p & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q + \vec{N}^{q+\frac{1}{2}}) , \end{array}\right.$

where $q$ denotes time $t_q$, which is sixth order in $r$ for $\vec{u}^{q+1}$ and second order for $p$, where the solenoidal condition is not explicitly imposed. Symmetry provides the conditions at the axis. The difficulty is in imposing the remaining four boundary conditions — this system should be inverted, in principle, simultaneously for $p$ and $\vec{u}^{q+1}$ with boundary conditions $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$ (Rempfer 2006). In practice it would be preferable to invert for $p$ first then for $\vec{u}^{q+1}$, but the boundary conditions to not involve $p$ directly.

Note that the $\mat{Y}\,\vec{u}^q$ term has been included in the right-hand side of the pressure-Poisson equation, the divergence of which should be small. Assume that pressure boundary condition is known: the right-hand side of the Navier–Stokes equation is then projected onto the space of solenoidal functions though $p$ and hence after inversion, $\vec{u}^{q+1}$ will be solenoidal.

Consider the ‘bulk’ solution, $\{\bar{\vec{u}},\bar{p}\}$, obtained from solution of the following:

$\label{eq:NSbulk} (3) \qquad \left\{\begin{array}{rcl} \mat{X}\, \bar{\vec{u}} & = & \mat{Y}\, \vec{u}^q + \vec{N}^{q+\frac{1}{2}} - \bnabla \bar{p} \, , \\ \nabla^2 \bar{p} & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q + \vec{N}^{q+\frac{1}{2}}) , \end{array}\right.$

with boundary conditions $\bar{\vec{u}}=\vec{0}$ and $\partial_{r}\bar{p}=0$. Introduce the following systems:

$\label{eq:NSp0} (4) \qquad \left\{\begin{array}{rcl} \mat{X}\,\vec{u}' & = & -\bnabla p' \, , \\ \nabla^2 p' & = & 0, \end{array}\right.$

with boundary conditions $\vec{u}'=\vec{0}$ and $\partial_{r}p'=1$ on $r=R$, and

$\label{eq:NSu0} (5) \qquad \left\{\begin{array}{rcl} \mat{X}\, \vec{u}' & = & \vec{0}, \end{array}\right.$

with boundary conditions $u'_+=1$, $u'_-=1$, $u'_z=\mathrm{i}$ on $r=R$ (see Decoupling_the_equations). The system $(4)$ provides a linearly independent function $\vec{u}'_1$ that may be added to $\bar{\vec{u}}$ without affecting the right-hand side in $(3)$, but altering (to correct) the boundary condition applied. Similarly the system $(5)$ provides a further three functions $\vec{u}'_j$ (where $j=2,3,4$; the +,- and $z$ components may be considered separately). The superposition

$\label{eq:usuperpos} (6) \qquad \vec{u}^{q+1} = \bar{\vec{u}} + \sum_{j=1}^4 a_j\, \vec{u}'_j \,$

may be formed in order to satisfy the four boundary conditions, $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$. Let $\vec{g}(\vec{u})$ be a 4-vector composed of these boundary conditions, such that they are satisfied when $\vec{g}(\vec{u})=\vec{0}$. Substituting the superposition $(6)$ into the boundary conditions, they may be written

$\label{eq:velBCs} (7) \qquad \mat{A}\,\vec{a} = -\vec{g}(\bar{\vec{u}}) ,$

where $\mat{A}=\mat{A}(\vec{g}(\vec{u}'))$ is a 4$\times$4 matrix and the 4-vector $\vec{a}$ is composed of the $a_j$. Thus, the appropriate coefficients required to satisfy the boundary conditions are recovered from solution of this small system for $\vec{a}$.

The error in the boundary conditions $g_j(\vec{u}^{q+1})$ using the influence-matrix technique is at the level of the machine epsilon, typically order 1e-16. The functions $u'_j(r)$, the matrix $\mat{A}$ and its inverse may all be precomputed. The boundary conditions for $\vec{u}'$ have been chosen so that that $u'_\pm$ are pure real, $u'_z$ is pure imaginary, and $\mat{A}$ is real. For each timestep, this application of the influence matrix technique requires only evaluation of the deviation from the boundary condition, multiplication by a 4$\times$4 real matrix, and the addition of only two functions to each component of $\vec{u}$, each either pure real or pure imaginary. Compared to the evaluation of nonlinear terms, the computational overhead is negligible.

If needed, the pressure may be calculated using the same $a_1$ as in $(6)$ using the $p'$ from $(4)$:

$\label{eq:pressure} (8a) \qquad \tilde{p} = \bar{p} + a_1\, p' \, \qquad$ (for each Fourier mode, then)

$ (8b) \qquad p = \tilde{p} - \frac{1}{2}\vec{u}\cdot\vec{u}\, \qquad$ (in physical space).

The adjustment in $(8b)$ arises when the Navier-Stokes equations are in rotational form (see Equations_and_parameters).