|
version 0.6.0
|
|
|
version 0.6.0
|
Prepare and solve the Navier-Stokes equations. More...
Namespaces | |
| module | mod_compressible_predict_density |
| Predict the density at the next time step. | |
| module | mod_compute_pressure_gradient |
| Compute pressure gradient on the face grid. | |
| module | enum_navier_velocity_pressure_method |
| Enumerations associated to the Navier-Stokes equations. | |
| module | fields_navier |
| Declaration of the field arrays associated to the Navier-Stokes equations. | |
| module | mod_navier_predict_density |
| Predict the density at the next time step. | |
| module | mod_solve_pressure |
| This module solves the correction step of the pressure correction method. | |
| module | variables_navier |
| Define variables (time step, schemes, solvers, etc.) associated to the Navier-Stokes equations. | |
Functions | |
| subroutine | mod_apply_bc_on_pressure_increment::apply_bc_on_pressure_increment () |
| Set boundary conditions on pressure increment as a function of the ones defined for velocity field. | |
| subroutine | mod_initialize_hydrostatic_pressure::initialize_hydrostatic_pressure () |
| Initialize hydrostatic pressure. | |
| subroutine | mod_project_velocity::project_velocity () |
| Project the velocity on incompressible velocities This is useful when one initialize the velocity field with a non divergence free field. | |
| subroutine | mod_set_brinkman_term_coefficient::set_brinkman_term_coefficient (permeability, viscosity, brinkman_coef) |
| Set Brinkman term coefficient \( \mu/K \). | |
| subroutine | solve_navier () |
| Solve the Navier-Stokes equations. | |
| subroutine | mod_solve_pressure::solve_pressure (density_star_face, brinkman_coef, carman_kozeny_coef) |
| This routine solves the correction step of the pressure correction method. | |
This directory contains the routines necessary to prepare the code to solve the Navier-Stokes equations. For incompressible flows, the Navier-Stokes has the following form: \( \operatorname{div} \mathbf{u} = 0 \tag{Continuity equation} \\ \)
\begin{align} \rho \frac{\partial \mathbf{u}}{\partial t} =& \\ &- \nabla p \tag{Pressure gradient term} \\ &+ \operatorname{div} \Big( \mu \big( \nabla\mathbf{u} + \nabla\mathbf{u}^{\intercal} \big) \Big) \tag{Stress term} \\ &- \rho \operatorname{div} \big( \mathbf{u} \otimes \mathbf{u} \big) \tag{Inertial term} \\ &+ \nabla r \big( \operatorname{div} \mathbf{u} \big) \tag{Grad-div term - optional} \\ &- \frac{\mu}{K} \mathbf{u} \tag{Brinkman term - optional} \\ &+ \mathbf{f} \tag{Source term - optional} \\ &+ \mathbf{g} \tag{Gravity term - optional} \\ &- \mathbf{L} \mathbf{u} \tag{Linear term - optional} \\ &+ \mathbf{a} \tag{Capillarity term - optional} \end{align}
To solve the velocity/pressure coupling, the pressure-correction methods are used. If viscosity is constant in the domain, which corresponds to monophasic flows in Notus (since viscosity is constant for a given fluid so far), the rotational pressure-correction method of Timmermans [2] is used. For mulptiphase flows, the Goda method [1] is used. These method are based on two steps:
This step computes a non-solenoidal predicted velocity \( \mathbf{u}^{*,n+1} \):
\[ \rho \left( \frac {\alpha \mathbf{u}^{*,n+1} + \beta \mathbf{u}^n + \gamma \mathbf{u}^{n-1}} {\Delta t} + \operatorname{div} \Big( (\zeta \mathbf{u}^n + \eta \mathbf{u}^{n-1}) \otimes \tilde{\mathbf{u}^*} \big) \right) = - \nabla p^n + \nabla \cdot \left[ \mu \left( {\nabla \mathbf{u}^{*,n+1} + \nabla^\top \mathbf{u}^{*,n+1}} \right) \right] \]
\[ B.C.(\mathbf{u}^{*,n+1}) \]
where values of \( \alpha, \beta, \gamma \) helps to switch from the 1st order Euler time discretization scheme to the 2nd order backward differential one:
\( \alpha = 1, \beta = -1, \gamma = 0 \) for the Euler scheme
\( \alpha = \frac {3} {2}, \beta = -2, \gamma = \frac {1} {2} \) for 2nd order BDF
The non linear term of the momentum equation is linearized in time at first or second order depending of the previous time order discretization scheme:
\( \zeta = 1, \eta = 0 \) if the Euler scheme is used
\( \zeta = 2, \eta = -1 \) (second time order extrapolation) if the 2nd order BDF is used
This term can be discretized explicitly \( \tilde{\mathbf{u}^*} = \mathbf{u}^{*,n} \) or implicitly \( \tilde{\mathbf{u}^*} = \mathbf{u}^{*,n+1} \)
This step computes the pressure increment \( \varphi^{n+1} = p^{n+1} - p^n \).
\[ \nabla \cdot \left( \frac {\Delta t} {\alpha \rho} \nabla \varphi^{n+1} \right) = \nabla \cdot \mathbf{u}^{*,n+1} \]
\[ B.C.(\varphi^{n+1}) \]
Pressure and solenoidal velocity are then updated according to:
\[ p^{n+1} = p^n + \varphi^{n+1} - \chi \mu \nabla \cdot \mathbf{u}^{*,n+1} \]
\[ \mathbf{u}^{n+1} = \mathbf{u}^{*,n+1} - \frac {\Delta t} {\rho \alpha} \nabla \varphi^{n+1} \]
with \( \chi = 0\) (Goda method) or \(\chi = 1\) (Timmermans method).
The corresponding boundary conditions between velocity and pressure increment is the following: where a Dirichlet boundary condition is set to velocity, a Neumann one is set to pressure increment; where a Neumann boundary condition is set to velocity, a Dirichlet one of zero is set to pressure increment.
Regardless of the spatial discretization schemes used for the implicit or explicit predictor step, the Poisson equation of the pressure-correction methods must be solved with a linear solver. The iterative methods of BiCGStab(2) (or GMRES) with the geometric multigrid preconditioners PFMG or SMG of the high performance preconditioners and solvers library Hypre [3] is proposed in Notus. This computationally intensive calculation may be further complicated by a non-constant coefficient matrix, as is obtained with the non-uniform density flows. To reduce the impact of this computationally intensive calculation, the approach proposed by Dodd and Ferrante [4] for the Chorin pressure-correction method is extended to the pressure increment Poisson equation [6]. As a result, the variable coefficient matrix is reformulated into a constant coefficient matrix following
\[ \frac {\nabla \varphi^{n+1}} {\rho^n} \approx \frac {\nabla \varphi^{n+1}} {\rho^0} + \left( \frac {1}{\rho^n} - \frac {1}{\rho^0} \right) \nabla \varphi^n \]
where \( \rho^0 = \min(\rho^n) \). The resulting Poisson pressure increment equation follows,
\[ \Delta \varphi^{n+1} = \nabla \cdot \left[ \left( 1 - \frac {\rho^0} {\rho^n} \right) \nabla \varphi^n \right] + \frac {\rho^0}{\Delta t} \nabla \cdot \mathbf{u}^{*,n+1} \]
This method may yield significant performance increases. However, if density ratio is high, the method yields to a decrease of the time step to ensure stability.
Another (experimental) method exists in Notus to reduce computational time of the correction step in case of immersed domains into the Cartesian mesh. Frantzis and Grigoriadis [5] proposed to avoid applying the Neumann boundary condition of pressure increment on the immersed interface. Even if it is not mathematically consistent, it is shown in [5] that, with minor adaptation, the solution remains very close to the unmodified system ones. A modification to this method is implemented [6], wherein the correction step is solved on the whole domain including the inside of the immersed subdomain.
The 2D/3D fields associated to this equation are defined in fields.f90 file.
The scalar variables such as the time step, the schemes, the residuals, etc. are defined in variables.f90 file.
Enumarations that help to do some choices (schemes, terms, etc.) are declared in enum.f90.
The time splitting method is done in solve_navier.f90 file that is separated into 2 parts: the first one solves the momentum equation, the second one solves the pressure increment thanks to the solve_pressure.f90 routine.
A set of upper level routines helps to discretize the Navier-Stokes equations: compute_pressure_gradient.f90, compute_surface_tension_term.f90, set_brinkman_term_coefficient.f90,
The directory contains some utilities: initialize_hydrostatic_pressure.f90, project_velocity.f90
[1] K. Goda, A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows, Journal of Computational Physics 30 (1979) 76–95.
[2] L. J. P. Timmermans, P. D. Minev, F. N. Van De Voss, An approximate projection scheme for incompressible flow using spectral elements, International Journal for Numerical Methods in Fluids 22 (1996) 673–688.
[3] R.D. Falgout, J.E. Jones, U.M. Yang, The Design and Implementation of hypre, a Library of Parallel High Performance Preconditioners, in: A.M. Bruaset, A. Tveito (Eds.), Numer. Solut. Partial Differ. Equ. Parallel Comput., Springer, Berlin, Heidelberg, 2006: pp. 267–294.
[4] M.S. Dodd, A. Ferrante, A fast pressure-correction method for incompressible two-fluid flows, J. Comput. Phys. 273 (2014) 416–434.
[5] C. Frantzis, D.G.E. Grigoriadis, An efficient method for two-fluid incompressible flows appropriate for the immersed boundary method, J. Comput. Phys. 376 (2019) 28–53.
[6] A. M. D. Jost, S. Glockner, A. Erriguible, Direct numerical simulations of fluids mixing above mixture critical point, Submitted to Journal of Supercritical Fluids, 2020
| subroutine mod_initialize_hydrostatic_pressure::initialize_hydrostatic_pressure |
This routine initialzes hydrostatic pressure inside the domain by solving hydrostatic equilibrium equation.
| subroutine mod_set_brinkman_term_coefficient::set_brinkman_term_coefficient | ( | double precision, dimension(nx,ny,nz), intent(in) | permeability, |
| double precision, dimension(nx,ny,nz), intent(in) | viscosity, | ||
| type(t_face_field), intent(inout) | brinkman_coef ) |
This term is defined on velocity nodes. Permeabilty and viscosity are interpolated from cell center to faces thanks to the harmonic interpolation.
| subroutine solve_navier |
First part of the routine builds rhs and matrix of the prediction step, term by term. Then, the linear system is solved thanks to Hypre or MUMPS libraries.
The second part solves the correction step. This is done in solve_pressure.f90 that is called at the end of solve_navier.
| subroutine mod_solve_pressure::solve_pressure | ( | type(t_face_field), intent(in) | density_star_face, |
| type(t_face_field), intent(in) | brinkman_coef, | ||
| type(t_face_field), intent(in) | carman_kozeny_coef ) |
This routine fills the linear system relative to the correction step on pressure increment:
velocity_correction.f90) 
1.11.0