1452
1452
In \fluidity\ the low order pivot value uses first order upwinding (see section \ref{sec:cv_fou}) and the pivot implicitness factor, $\theta_p$, defaults to $1$. This implicit pivot is then used to overcome the timestep restriction and an extra advection iteration loop is introduced to update the values of $\tilde{c}_{f_k}^{HO}$ and $\tilde{c}_{f_k}^{LO}$. If this converges, then after a number of iterations $c_{f_k}^{LO n+1} \approx \tilde{c}_{f_k}^{LO}$ and $c_{f_k} \approx \theta c_{f_k}^{HO n+1} + \theta c_{f_k}^{HO n}$. Hence an implicit high order face value is achieved. If the iteration does not converge, either due to too few advection iterations being selected or as a result of non-convergent behaviour in the face value limiter, then $c_{f_k}$ is just a linear combination of low order and high order face values. This still results in a valid face value, although as is generally the case for implicit advection methods, it is not possible to guarantee the boundedness of this scheme.
1454
\subsection{Porous Media}
1455
\label{sec:time_disc_scalar_field_porous}
1456
\index{porous media darcy flow}
1458
In this section the modification to the temporal discretisation of the general scalar transport equation to include porosity is described.
1460
Using the product rule for the time partial derivative of the modified scalar transport equation \eqref{eq:general_scalar_eqn_single_phase_darcy} gives:
1462
\phi \ppt{c} + \ppt{\phi} c + \nabla\cdot(\vec{u}_{\phi} c) = 0.
1464
This is discretised using the classical $\theta$ scheme to give:
1466
\phi^{n} \frac{c^{n+1}-c^{n}}{\Delta t} + \nabla\cdot(\urelax_{\phi} c^{n+\theta}) = - c^{n} \frac{\phi^{n+1}-\phi^{n}}{\Delta t}.
1468
This discretisation has chosen to explicitly represent $\phi$ when associated with the discretised time derivative of $c$ and vice versa. In a similar manner to that derived in section \ref{sec:ND_time_disct_adv_diff} this can be expressed in semi discrete form as:
1470
\left ( \mat{M}_{\phi^{n}} + \thetac \Delta t \mat{A}(\urelax_{\phi}) \right ) \cnew =
1471
\left ( \mat{M}_{\phi^{n}} - ( 1 - \thetac ) \Delta t \mat{A}(\urelax_{\phi}) \right ) \cold -
1472
\mat{M}_{c^{n}} (\dvec{\phi}^{n+1}-\dvec{\phi}^{n}).
1474
Taking $\psi$ to represent the test and basis function space associated with $c$ and $\zeta$ to represent the basis function space associated with the porosity $\phi$, the modified mass matrices take the form:
1476
\mat{M}_{\phi^{n}_{ij}} = \int_\Omega \psi_i \psi_j \phi^{n}, \quad
1477
\mat{M}_{c^{n}_{ij}} = \int_\Omega \psi_i \zeta_j c^{n},
1479
where $\phi^{n}$ and $c^{n}$ are known functions (which are thus numerically evaluated at the quadrature points during the integration sum). If the rate of change of porosity with time is negligible then the discretised equation for $c$ reduces to
1481
\left ( \mat{M}_{\phi^{n}} + \thetac \Delta t \mat{A}(\urelax_{\phi}) \right ) \cnew =
1482
\left ( \mat{M}_{\phi^{n}} - ( 1 - \thetac ) \Delta t \mat{A}(\urelax_{\phi}) \right ) \cold.
1484
This form of the porous media scalar transport equation is the one permitted in \fluidity for control volume discretisations with $\zeta$ equal to $\psi$ or piece wise constant over an element. A similar derivation is used for the transport of the individual components of the metric tensor associated with mesh adaptivity.
1454
1486
\section{Momentum equation}
1455
1487
\index{momentum equation!discretised}
1533
1565
\end{equation*}
1567
\subsection{Porous Media Darcy Flow}
1568
\label{sec:porous_media_darcy_flow_discretisation}
1569
\index{porous media darcy flow}
1571
\subsubsection{Single Phase}
1573
The discretisation of the Darcy velocity equation \eqref{eq:mom_press_single_phase_darcy} is given by:
1575
\mat{M}_{\sigma} \dvec{u}_{\phi}
1576
= - \matC \dvec{p} + \dvec{b}(\rho),
1578
where (ignoring boundary conditions):
1581
\mat{M}_{\sigma_{ij}}=\int_\Omega \bmphi_i\cdot\bmphi_j \tensor{\sigma}, \quad
1582
\matC_{ij}=\int_\Omega \bmphi_i\cdot\grad\psi_j, \quad
1583
\dvec{b}_i(\rho)=\int_\Omega \bmphi_i\cdot\vec g\rho.
1535
1587
\section{Pressure equation for incompressible flow}
1536
1588
\index{pressue!discretisation}
1537
1589
\label{sec:ND_pressure_equation}
1747
1799
is the same as that would be obtained from a direct $P_{n+1}$
1748
1800
discretisation of the continuous pressure poisson equation.
1802
\subsection{Porous Media Darcy Flow}
1803
\label{sec:porous_media_darcy_flow_pressure_correction}
1804
\index{porous media darcy flow}
1806
\subsubsection{Single Phase}
1808
To solve the single phase porous media Darcy force balance \eqref{eq:mom_press_single_phase_darcy}
1809
and incompressible mass conservation \eqref{eq:cty_single_phase_darcy} the pressure correction
1810
algorithm described in section \ref{sec:pressure_correction} requires modification. As there is
1811
no momentum time term the absorption term is instead included in the correction step. The pressure
1812
correction algorithm then takes the same steps as before but solving for the Darcy velocity rather
1813
than the interstitial velocity. The first step is to solve for a preliminary $\dvec{u}_{\phi}^{*}$
1814
from an initial pressure guess $\dvec{p}^{*}$ via:
1816
\mat{M}_{\sigma} \dvec{u}_{\phi_{*}}
1817
= - \matC \dvec{p}_{*} + \dvec{b}(\rho).
1819
A pressure correction is then solved for via the derived equation
1821
\matB^T \mat{M_{\sigma}^{-1}}\matC~\Delta\dvec p
1822
= \matB^T \dvec{u}_{\phi_{*}}-\vec{\mat{M}}_{D} ~\dvec{g}_{D}.
1824
After obtaining the pressure correction the pressure and velocity can
1828
\dvec{\tilde p} &= \dvec p_{*} + \Delta\dvec p, \\
1829
\dvec{\tilde u}_{\phi} &= \dvec{u}_{\phi_{*}}
1830
- M_{\sigma}^{-1} \matC \Delta\dvec p.
1833
For incompressible flow as there is no Darcy velocity time dependence and because all the
1834
Darcy velocity terms (being only the absorption) are included in the pressure correction step
1835
there is no requirement for multiple non linear iterations, unless there is a non linear absorption
1836
coefficient or source term.
1750
1838
\section{Velocity and pressure element pairs}
1751
1839
\label{sec:velocity_pressure_element_pairs}
1752
1840
As we have seen, various methods are available in \fluidity for the discretisation of