~fluidity-core/fluidity/shallow-water-dev

« back to all changes in this revision

Viewing changes to manual/numerical_discretisation.tex

  • Committer: colin.cotter at ac
  • Date: 2012-02-16 17:35:27 UTC
  • mfrom: (3565.1.365 fluidity)
  • Revision ID: colin.cotter@imperial.ac.uk-20120216173527-wpnlu4v9v1kmf4xz
MergeĀ fromĀ trunk.

Show diffs side-by-side

added added

removed removed

Lines of Context:
1451
1451
 
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.
1453
1453
 
 
1454
\subsection{Porous Media}
 
1455
\label{sec:time_disc_scalar_field_porous}
 
1456
\index{porous media darcy flow}
 
1457
 
 
1458
In this section the modification to the temporal discretisation of the general scalar transport equation to include porosity is described. 
 
1459
 
 
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:
 
1461
\begin{equation}
 
1462
   \phi \ppt{c} + \ppt{\phi} c + \nabla\cdot(\vec{u}_{\phi} c) = 0.
 
1463
\end{equation}
 
1464
This is discretised using the classical $\theta$ scheme to give:
 
1465
\begin{equation}
 
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}.
 
1467
\end{equation}
 
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:
 
1469
\begin{equation}
 
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}).
 
1473
\end{equation}
 
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:
 
1475
\begin{equation}
 
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},
 
1478
\end{equation}
 
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
 
1480
 \begin{equation}
 
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.
 
1483
\end{equation}
 
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.
 
1485
 
1454
1486
\section{Momentum equation}
1455
1487
\index{momentum equation!discretised}
1456
1488
 
1532
1564
\end{split}
1533
1565
\end{equation*}
1534
1566
 
 
1567
\subsection{Porous Media Darcy Flow}
 
1568
\label{sec:porous_media_darcy_flow_discretisation}
 
1569
\index{porous media darcy flow}
 
1570
 
 
1571
\subsubsection{Single Phase}
 
1572
 
 
1573
The discretisation of the Darcy velocity equation \eqref{eq:mom_press_single_phase_darcy} is given by:
 
1574
\begin{equation}
 
1575
  \mat{M}_{\sigma} \dvec{u}_{\phi}
 
1576
    = - \matC \dvec{p} + \dvec{b}(\rho),
 
1577
\end{equation}
 
1578
where (ignoring boundary conditions):
 
1579
\begin{equation*}
 
1580
\begin{split}
 
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.
 
1584
\end{split}
 
1585
\end{equation*}
 
1586
 
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.
1749
1801
 
 
1802
\subsection{Porous Media Darcy Flow}
 
1803
\label{sec:porous_media_darcy_flow_pressure_correction}
 
1804
\index{porous media darcy flow}
 
1805
 
 
1806
\subsubsection{Single Phase}
 
1807
 
 
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:
 
1815
\begin{equation}
 
1816
  \mat{M}_{\sigma} \dvec{u}_{\phi_{*}}
 
1817
    = - \matC \dvec{p}_{*} + \dvec{b}(\rho).
 
1818
\end{equation}
 
1819
A pressure correction is then solved for via the derived equation
 
1820
\begin{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}.
 
1823
\end{equation}
 
1824
After obtaining the pressure correction the pressure and velocity can
 
1825
be updated by:
 
1826
\begin{equation*}
 
1827
\begin{split}
 
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.
 
1831
\end{split}
 
1832
\end{equation*}
 
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.
 
1837
 
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