CorrectionPressureAlgorithm

Let us suppose we have a velocity field \( \mathbf{\tilde{v}}^{n+1}\) that stems from the numerical integration of the momentum equation (see EquationsToSolve), i.e. we have computed \begin{align} \frac{\mathbf{\tilde{v}}^{n+1}-\mathbf{v}^{n}}{\Delta t}+\frac{1}{\rho }\nabla p^{n}=\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{s}^{n} \right)^{T}+\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{v}\left( \mathbf{\tilde{v}}^{n+1} \right) \right)^{T}+\mathbf{g}-\beta \cdot \mathbf{\tilde{v}}^{n+1}\end{align} The resulting velocity field \( \mathbf{\tilde{v}}^{n+1}\) does most probably not provide the correct value of divergence of velocity \( \nabla ^{T}\mathbf{\tilde{v}}^{n+1}\) . Let us suppose there is a correction to the pressure \( c\) that exactly leads to the correct divergence, that is \begin{align} \frac{\mathbf{v}^{n+1}-\mathbf{v}^{n}}{\Delta t}+\frac{1}{\rho }\nabla p^{n}+\frac{1}{\rho }\nabla c=\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{s}^{n} \right)^{T}+\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{v}\left( \mathbf{v}^{n+1} \right) \right)^{T}+\mathbf{g}-\beta \cdot \mathbf{v}^{n+1}\end{align} With the presumption that \( \nabla ^{T}\mathbf{v}^{n+1}\) has the correct value. In order to find the correction pressure \( c\) we subtract the two equations from one another, that is \begin{align} \left( \frac{1+\Delta t\beta }{\Delta t} \right)\cdot \left( \mathbf{v}^{n+1}-\mathbf{\tilde{v}}^{n+1} \right)+\frac{1}{\rho }\nabla c=+\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{v}\left( \mathbf{v}^{n+1} \right) \right)^{T}-\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{v}\left( \mathbf{\tilde{v}}^{n+1} \right) \right)^{T}\end{align} Written in another way, we have \begin{align} \left( \mathbf{v}^{n+1}-\mathbf{\tilde{v}}^{n+1} \right)+\left( \frac{\Delta t}{1+\Delta t\beta } \right)\frac{1}{\rho }\nabla c=\left( \frac{\Delta t}{1+\Delta t\beta } \right)\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{v}\left( \mathbf{v}^{n+1} \right) \right)^{T}-\left( \frac{\Delta t}{1+\Delta t\beta } \right)\frac{1}{\rho }\left( \nabla ^{T}\mathbf{S}_{v}\left( \mathbf{\tilde{v}}^{n+1} \right) \right)^{T}\end{align} For incompressible problems with constant viscosity, we can simplify \begin{align} \left( \mathbf{v}^{n+1}-\mathbf{\tilde{v}}^{n+1} \right)+\left( \frac{\Delta t}{1+\Delta t\beta } \right)\frac{1}{\rho }\nabla c=\left( \frac{\Delta t}{1+\Delta t\beta } \right)\frac{\eta }{\rho }\Delta \left( \mathbf{v}^{n+1}-\mathbf{\tilde{v}}^{n+1} \right)\end{align} The correction pressure stems from the simplified correction ansatz of a given velocity field (marked by tilde) towards a velocity field with a desired divergence of velocity. \begin{align} \mathbf{v}^{n+1}=\mathbf{\tilde{v}}^{n+1}-\frac{\Delta t_{virt}}{\left( 1+\Delta t_{\beta }\cdot \beta \right)}\frac{1}{\rho }\nabla c\end{align} For the term \( \Delta t_{virt}\) see VirtualTimeStepSize. By application of the divergence operator from left, we obtain \begin{align} \left( \nabla ^{T}\mathbf{v}^{n+1} \right)_{target}-\nabla ^{T}\mathbf{\tilde{v}}^{n+1}+\nabla ^{T}\left( \frac{\Delta t_{virt}}{\left( 1+\Delta t_{\beta }\cdot \beta \right)}\frac{1}{\rho }\nabla c \right)=0\end{align} The desired divergence of the velocity is depending on the compressibility of the fluid as well as on temporal changes of the density due to other effects such as chemical reaction, expansion due to heating, etc. Derivation of this term is found in DesiredAndNominalDivergenceOfVelocity. Having a formulation for the divergence of velocity, the equation to be solved for the correction pressure is \begin{align} -\frac{1}{\rho }\frac{\partial \rho }{\partial p}\frac{1}{\Delta t}c+\nabla ^{T}\left( \frac{\Delta t_{virt}}{\left( 1+\Delta t_{\beta }\cdot \beta \right)}\frac{1}{\rho }\nabla c \right)-\nabla ^{T}\mathbf{\tilde{v}}^{n+1}=-\left( \overline{\nabla ^{T}\mathbf{v}} \right)_{c}^{n+1}\end{align} which numerically leads to the (linear) system to be solved \begin{align} -\frac{1}{\rho }\frac{\partial \rho }{\partial p}\frac{1}{\Delta t}c_{i}+\sum\limits_{j}{d_{ij}^{\nabla ^{T}\left( \frac{\Delta t_{virt}}{\left( 1+\Delta t_{\beta }\cdot \beta \right)}\frac{1}{\rho }\nabla \right)}}\cdot c_{ij}-\nabla ^{T}\mathbf{\tilde{v}}^{n+1}=-\left( \overline{\nabla ^{T}\mathbf{v}} \right)_{c}^{n+1}\end{align} The result of this equation is stored in %ind_c% .