maelstrom.navier_stokes¶
Numerical solution schemes for the Navier–Stokes equation in cylindrical coordinates,
cf. https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Representations_in_3D. In the weak formulation, we consider integrals in pseudo 3D, resulting in a weighting with \(2\pi r\) of the equations. (The volume element is \(2\pi r\,\text{d}x\).)
The order of the variables is taken to be \((r, z, \theta)\). This makes sure that for planar domains, the \(x\)- and \(y\)-coordinates are interpreted as \(r\), \(z\).
Note that [LMW+12] contains a chapter with an extensive comparison of solution methods for the transient Navier–Stokes equations for many different problems. It is found that IPCS is the best (fastest, most accurate) solution method.
An overview of projection methods for incompressible flow can be found in [GMS06] and [Joh16].
-
class
maelstrom.navier_stokes.IPCS(time_step_method='backward euler')¶ Bases:
objectIncremental pressure correction scheme; for details see [GMS06].
-
order= {'pressure': 1, 'velocity': 1}¶
-
step(dt, u, p0, W, P, u_bcs, p_bcs, rho, mu, f, verbose=True, tol=1e-10, my_dx=<Mock name='mock.dx' id='139755140315512'>)¶
-
-
maelstrom.navier_stokes.compute_pressure(P, p0, mu, ui, u, my_dx, p_bcs=None, rotational_form=False, tol=1e-10, verbose=True)¶ Solve the pressure Poisson equation
\[\begin{split}\begin{align} -\frac{1}{r} \div(r \nabla (p_1-p_0)) = -\frac{1}{r} \div(r u),\\ \text{(with boundary conditions)}, \end{align}\end{split}\]for \(\nabla p = u\).
The pressure correction is based on the update formula
\[\begin{split}\frac{\rho}{dt} (u_{n+1}-u^*) + \begin{pmatrix} \text{d}\phi/\text{d}r\\ \text{d}\phi/\text{d}z\\ \frac{1}{r} \text{d}\phi/\text{d}\theta \end{pmatrix} = 0\end{split}\]with \(\phi = p_{n+1} - p^*\) and
\[ \frac{1}{r} \frac{\text{d}}{\text{d}r} (r u_r^{(n+1)}) + \frac{\text{d}}{\text{d}z} (u_z^{(n+1)}) + \frac{1}{r} \frac{\text{d}}{\text{d}\theta} (u_{\theta}^{(n+1)}) = 0\]With the assumption that u does not change in the direction \(\theta\), one derives
\[\begin{split}- \frac{1}{r} \div(r \nabla \phi) = \frac{1}{r} \frac{\rho}{dt} \div(r (u_{n+1} - u^*))\\ - \frac{1}{r} \langle n, r \nabla \phi\rangle = \frac{1}{r} \frac{\rho}{dt} \langle n, r (u_{n+1} - u^*)\rangle\end{split}\]In its weak form, this is
\[\int r \langle\nabla\phi, \nabla q\rangle \,2 \pi = - \frac{\rho}{dt} \int \div(r u^*) q \, 2 \pi - \frac{\rho}{dt} \int_{\Gamma} \langle n, r (u_{n+1}-u^*)\rangle q \, 2\pi.\](The terms \(1/r\) cancel with the volume elements \(2\pi r\).) If the Dirichlet boundary conditions are applied to both \(u^*\) and \(u_n\) (the latter in the velocity correction step), the boundary integral vanishes.
If no Dirichlet conditions are given (which is the default case), the system has no unique solution; one eigenvalue is 0. This however, does not hurt CG convergence if the system is consistent, cf. [vdV03]. And indeed it is consistent if and only if
\[\int_\Gamma r \langle n, u\rangle = 0.\]This condition makes clear that for incompressible Navier-Stokes, one either needs to make sure that inflow and outflow always add up to 0, or one has to specify pressure boundary conditions.
Note that, when using a multigrid preconditioner as is done here, the coarse solver must be chosen such that it preserves the nullspace of the problem.
-
maelstrom.navier_stokes.compute_tentative_velocity(time_step_method, rho, mu, u, p0, dt, u_bcs, f, W, my_dx, tol)¶ Compute the tentative velocity via
\[\rho (u_0 + (u\cdot\nabla)u) = \mu \frac{1}{r} \div(r \nabla u) + \rho g.\]
-
maelstrom.navier_stokes.compute_velocity_correction(ui, p0, p1, u_bcs, rho, mu, dt, rotational_form, my_dx, tol, verbose)¶ Compute the velocity correction according to
\[U = u_0 - \frac{dt}{\rho} \nabla (p_1-p_0).\]