vp-

We describe the numerical scheme for incompressible / weakly compressible. A document containing the scheme is found in DOCUMATH_GeneralNumericalScheme.pdf. The timestep starts with an (explicit!) movement of the MESHFREE points. \begin{align} {\bf x}^{n+1}_i = {\bf x}^{n}_i + \Delta t \cdot {\bf v}^{n}_i \end{align} The new positions of time level \( n+1\) are found in %ind_x(1)% , %ind_x(2)% , %ind_x(3)% . The old positions of time level \( n\) are kept in %ind_x0(1)% , %ind_x0(2)% , %ind_x0(3)% . Compute all necessary material data. Especially see %ind_r%, %ind_ETA%, %ind_LAM%, %ind_MUE%, %ind_betaDarcy%, %ind_v0Darcy(1)% ... %ind_v0Darcy(3)%, %ind_SIG%, ... \begin{align} \begin{array}{*{35}{l}} \rho & =\rho \left( t^{n+1},p_{hyd}^{n}+p_{dyn}^{n},T^{n},A_{v}^{n} \right) \\ \eta & =\eta \left( t^{n+1},p_{hyd}^{n}+p_{dyn}^{n},T^{n},A_{v}^{n} \right) \\ \lambda & =\lambda \left( t^{n+1},p_{hyd}^{n}+p_{dyn}^{n},T^{n},A_{v}^{n} \right) \\ \mu & =\mu \left( t^{n+1},p_{hyd}^{n}+p_{dyn}^{n},T^{n},A_{v}^{n} \right) \\ k_{D} & =k_{D}\left( t^{n+1},p_{hyd}^{n}+p_{dyn}^{n},T^{n},A_{v}^{n} \right) \\ ... & {} \\ \end{array}\end{align} Also, compute derived data, for example: the compressibility, see %ind_R_P%, also %ind_DiagPcorr% . \begin{align} \frac{\partial \rho}{\partial p}\end{align} Compute the effective viscosity, see %ind_ETA_sm% and VelocityAlgorithm. \begin{align} \hat{\eta } = \eta + C\cdot\Delta t \cdot \mu + c_{\mu }\cdot \frac{k^{2}}{\varepsilon }\end{align} Compute the effective body forces, see %ind_g(1)% ... %ind_g(3)% \begin{align} \mathbf{\hat{g}}^n=\mathbf{g}^{n+1}+\frac{1}{\rho }\nabla \mathbf{S}_{s}^{n}\end{align} Solve the hydrostatic pressure. See HydrostaticPressureAlgorithm. See also LIQUID.%ind_p% . \begin{align} -\frac{1}{\Delta t^{2}}\left( \frac{1}{\rho }\frac{\partial \rho }{\partial p}\tilde{p}_{hyd}^{n+1} \right)+\nabla ^{T}\left( \frac{1}{\rho }\nabla \tilde{p}_{hyd}^{n+1} \right)=-\frac{1}{\Delta t^{2}}\left( \frac{1}{\rho }\frac{\partial \rho }{\partial p}p_{hyd}^{n} \right)+\nabla ^{T}\mathbf{\hat{g}}^n\end{align} Solve the temperature. See TemperatureAlgorithm. See %ind_T% . \begin{align} \rho \cdot c_{v}\frac{T^{n+1}-T^{n}}{\Delta t}=q+\nabla ^{T}\left( \lambda \cdot \nabla T^{n+1} \right)\end{align} Set up the preliminary dynamic pressure for the momentum equation. \begin{align} \tilde{p}^n_{dyn} = \mathcal{C} \cdot p^n_{dyn} \end{align} The preliminary value is stored in %ind_p_dyn% . The original value of the dynamic pressure at time level \( n\) is stored in %ind_p_dyn_0% . Remember, that the parameter \( \mathcal{C}\) is set by the input-file-parameter damping_p_corr . Compute the nominal divergence of velocity, needed for the desired divergence of velocity in CorrectionPressureAlgorithm, see especially DesiredAndNominalDivergenceOfVelocity. Temporarily saved in %ind_div_bar% . Accessible for post-processing via %ind_div_bar_c% . \begin{align} \left( \overline{\nabla ^{T}\mathbf{v}} \right)_{c}^{n+1}\equiv -\frac{1}{\Delta t}\left( \log \left( \rho \left( t^{n+1},p_{hyd}^{n+1}+\tilde{p}_{dyn}^{n},T^{n+1},A_{v}^{n} \right) \right)-\log \left( \rho_{p_{dyn}}^{n} \right) \right)\end{align} Solve the velocity and the correction pressure in one big system of equations. See %ind_v(1)% ... %ind_v(3)% as well as %ind_v_tild(1)% ... %ind_v_tild(3)% . See %ind_c% (correction pressure). \begin{align} \left\{ \begin{matrix} \frac{\mathbf{\tilde{v}}^{n+1}-\mathbf{v}^{n}}{\Delta t}+\frac{1}{\rho }\nabla \tilde{p}_{hyd}^{n+1}+\frac{1}{\rho }\nabla \tilde{p}_{dyn}^{n}+\frac{1}{\rho }\nabla c=\frac{1}{\rho }\nabla \mathbf{S}_{v}\left( \mathbf{\tilde{v}}^{n+1} \right)+\mathbf{\hat{g}}^n-\beta \cdot ( \mathbf{\tilde{v}}^{n+1} - \mathbf{v}_{\beta}) \\ -\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{matrix} \right\} \end{align} Update the dynamic pressure. See %ind_p_dyn% \begin{align} \tilde{p}_{dyn}^{n+1}=\tilde{p}_{dyn}^{n}+c\end{align} Correct the velocity with the help of the correction pressure if VP0_VelocityCorrection is switched on. Result in %ind_v(1)% ... %ind_v(3)% . \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} Compute the new density as a backup for the next time step. See %ind_r_c% . \begin{align}\rho _{c}^{n+1}=\rho \left( t^{n+1},p_{hyd}^{n+1}+\tilde{p}_{dyn}^{n+1},T^{n+1},A_{v}^{n} \right)\end{align} Compute the stress tensor at time level \( n+1\) by the stress tensor algorithm, i.e. \begin{align} {\bf S}^{n+1} = f \left( \Delta t, {\bf S}^{n}, {\bf v}^{n} \right)\end{align} see the StressTensorAlgorithm . Update turbulence values for k-epsilon. See %ind_k% and %ind_eps% . \begin{align} k^{n+1}=k^{n}+...\end{align} \begin{align} \varepsilon ^{n+1}=\varepsilon ^{n}+...\end{align} Recompute the resulting body forces. See %ind_g(1)% ... %ind_g(3)% . \begin{align} \mathbf{\hat{g}}^{n+1}=\frac{1}{\rho }\nabla \mathbf{S}_{s}^{n+1}+\mathbf{g}\end{align} Recompute, if needed, the hydrostatic pressure. See HydrostaticPressureAlgorithm . \begin{align} -\frac{1}{\Delta t^{2}}\left( \frac{1}{\rho }\frac{\partial \rho }{\partial p}p_{hyd}^{n+1} \right)+\nabla ^{T}\left( \frac{1}{\rho }\nabla p_{hyd}^{n+1} \right)=-\frac{1}{\Delta t^{2}}\left( \frac{1}{\rho }\frac{\partial \rho }{\partial p}p_{hyd}^{n} \right)+\nabla ^{T}\left( {\mathbf{\hat{g}}}^{n+1} \right)\end{align} Nominal divergence of velocity, motivated by dynamic pressure. Temporarily resulting in %ind_div_bar% . Accessible for post-processing via %ind_div_bar_pDyn% . \begin{align} \left( \overline{\nabla ^{T}\mathbf{v}} \right)_{dyn}^{n+1}\equiv -\frac{1}{\Delta t}\left( \log \left( \rho \left( ...,p_{hyd}^{n}+p_{dyn}^{n},... \right) \right)-\log \left( \rho _{{{p}_{dyn}}}^{n} \right) \right)\end{align} Out of the velocity field, compute the consistent dynamic pressure. See %ind_p_dyn% . \begin{align} \begin{array}{*{35}{l}} -\frac{1}{\Delta t^{2}}\left( \frac{1}{\rho }\frac{\partial \rho }{\partial p}p_{dyn}^{n+1} \right)+\nabla ^{T}\left( \frac{1}{\rho }\nabla \left( p_{dyn}^{n+1} \right) \right)= & -\frac{1}{\Delta t}\left( \left( \overline{\nabla ^{T}\mathbf{v}} \right)_{dyn}^{n+1}-\left( \nabla ^{T}\mathbf{v} \right)^{n} \right) \\ {} & -\frac{1}{\Delta t^{2}}\left( \frac{1}{\rho }\frac{\partial \rho }{\partial p}p_{dyn}^{n} \right) \\ {} & + \Psi \left( \mathbf{v}^{n+1} \right)-\Theta \left( \mathbf{v}^{n+1} \right)-\Phi \left( \mathbf{v}^{n+1} \right) \\ \end{array}\end{align} See especially the particular philosophy of providing boundary conditions to the solution of \( p_{dyn}^{n+1}\) ! Refer to %ind_p_dyn% . Switch off this algorithm by placing "PDYN:NONE" in the KindOfProblem definition ( KOP ). In this case, the dynamic pressure is taken from the CorrectionPressureAlgorithm in the sense \( p_{dyn}^{n+1}= \mathcal{C} \cdot p^n_{dyn} + c\) backup the density after computing the dynamic pressure \begin{align} \rho_{p_{dyn}}^{n+1}=\rho \left( t^{n+1},p_{hyd}^{n+1}+p_{dyn}^{n+1},T^{n+1},A_{v}^{n} \right)\end{align} compute the target divergence of velocity as a backup for the next time cycle. Resulting both in %ind_div_bar_0% and %ind_div_bar%. \begin{align} \left( \nabla ^{T}\mathbf{v} \right)_{target}^{n+1}=-\frac{1}{\Delta t}\left( \log \left( \rho_{p_{dyn}}^{n+1} \right)-\log \left( \rho_{p_{dyn}}^{n} \right) \right)\end{align} Compute the residuals for the velocity. See %ind_v_residual(1)% ... %ind_v_residual(3)% . Compute the residuals for the density. See %ind_r_residual% integrate all additional variables defined by the CODI commands. See also %ind_addvar% ... \begin{align} A_{v}^{n+1}=F\left( A_{v}^{n},t^{n+1},\text{v}^{n+1},T^{n+1},p_{hyd}^{n+1},p_{dyn}^{n+1},k^{n+1},\varepsilon ^{n+1} \right)\end{align}