Changeset 4265

Feb 16, 2007, 5:26:04 PM (17 years ago)

Move mathematical destcription from linuxconf paper and limiting.tex into user manual.
Have not tried to compile it though.

3 edited


  • anuga_core/documentation/user_manual/anuga_user_manual.tex

    r4258 r4265  
    26602660  Elevation Model) and converts it to PTS format.
    26612661  \end{funcdesc}
     2668\chapter{\anuga mathematical background}
     2669\label{cd:mathematical background}
     2673This chapter outlines the mathematics underpinning \anuga.
     2680The shallow water wave equations are a system of differential
     2681conservation equations which describe the flow of a thin layer of
     2682fluid over terrain. The form of the equations are:
     2684\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
     2685x}+\frac{\partial \GG}{\partial y}=\SSS
     2687where $\UU=\left[ {{\begin{array}{*{20}c}
     2688 h & {uh} & {vh} \\
     2689\end{array} }} \right]^T$ is the vector of conserved quantities; water depth
     2690$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
     2691entering the system are bed elevation $z$ and stage (absolute water
     2692level) $w$, where the relation $w = z + h$ holds true at all times.
     2693The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
     2696\EE=\left[ {{\begin{array}{*{20}c}
     2697 {uh} \hfill \\
     2698 {u^2h+gh^2/2} \hfill \\
     2699 {uvh} \hfill \\
     2700\end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c}
     2701 {vh} \hfill \\
     2702 {vuh} \hfill \\
     2703 {v^2h+gh^2/2} \hfill \\
     2704\end{array} }} \right]
     2706and the source term (which includes gravity and friction) is given
     2709\SSS=\left[ {{\begin{array}{*{20}c}
     2710 0 \hfill \\
     2711 -{gh(z_{x} + S_{fx} )} \hfill \\
     2712 -{gh(z_{y} + S_{fy} )} \hfill \\
     2713\end{array} }} \right]
     2715where $S_f$ is the bed friction. The friction term is modelled using
     2716Manning's resistance law
     2718S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy}
     2719=\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}
     2721in which $\eta$ is the Manning resistance coefficient.
     2723As demonstrated in our papers, \cite{modsim2005,Rob99l} these
     2724equations provide an excellent model of flows associated with
     2725inundation such as dam breaks and tsunamis.
     2727\section{Finite Volume Method}
     2730We use a finite-volume method for solving the shallow water wave
     2731equations \cite{Rob99l}. The study area is represented by a mesh of
     2732triangular cells as in Figure~\ref{fig:mesh} in which the conserved
     2733quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
     2734in each volume are to be determined. The size of the triangles may
     2735be varied within the mesh to allow greater resolution in regions of
     2736particular interest.
     2741\caption{Triangular mesh used in our finite volume method. Conserved
     2742quantities $h$, $uh$ and $vh$ are associated with the centroid of
     2743each triangular cell.} \label{fig:mesh}
     2747The equations constituting the finite-volume method are obtained by
     2748integrating the differential conservation equations over each
     2749triangular cell of the mesh. Introducing some notation we use $i$ to
     2750refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
     2751set of indices referring to the cells neighbouring the $i$th cell.
     2752Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
     2753the length of the edge between the $i$th and $j$th cells.
     2755By applying the divergence theorem we obtain for each volume an
     2756equation which describes the rate of change of the average of the
     2757conserved quantities within each cell, in terms of the fluxes across
     2758the edges of the cells and the effect of the source terms. In
     2759particular, rate equations associated with each cell have the form
     2761 \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i
     2765\item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell,
     2766\item $\SSS_i$ is the source term associated with the $i$th cell,
     2768\item $\HH_{ij}$ is the outward normal flux of
     2769material across the \textit{ij}th edge.
     2773%\item $l_{ij}$ is the length of the edge between the $i$th and $j$th
     2775%\item $m_{ij}$ is the midpoint of
     2776%the \textit{ij}th edge,
     2778%$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing
     2779%normal along the \textit{ij}th edge, and The
     2781The flux $\HH_{ij}$ is evaluated using a numerical flux function
     2782$\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow
     2783water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$
     2785H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
     2790\HH_{ij}  = \HH(\UU_i(m_{ij}),
     2791\UU_j(m_{ij}); \mathbf{n}_{ij})
     2793where $m_{ij}$ is the midpoint of the \textit{ij}th edge and
     2794$\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the
     2795\textit{ij}th edge. The function $\UU_i(x)$ for $x \in
     2796T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and
     2797neighbouring  cells.
     2799We use a second order reconstruction to produce a piece-wise linear
     2800function construction of the conserved quantities for  all $x \in
     2801T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This
     2802function is allowed to be discontinuous across the edges of the
     2803cells, but the slope of this function is limited to avoid
     2804artificially introduced oscillations.
     2806Godunov's method (see \cite{Toro-92}) involves calculating the
     2807numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly
     2808solving the corresponding one dimensional Riemann problem normal to
     2809the edge. We use the central-upwind scheme of \cite{KurNP2001} to
     2810calculate an approximation of the flux across each edge.
     2815\caption{From the values of the conserved quantities at the centroid
     2816of the cell and its neighbouring cells, a discontinuous piecewise
     2817linear reconstruction of the conserved quantities is obtained.}
     2822In the computations presented in this paper we use an explicit Euler
     2823time stepping method with variable timestepping adapted to the
     2824observed CFL condition.
     2827\section{Flux limiting}
     2829The shallow water equations are solved numerically using a
     2830finite volume method on unstructured triangular grid.
     2831The upwind central scheme due to Kurganov and Petrova is used as an
     2832approximate Riemann solver for the computation of inviscid flux functions.
     2833This makes it possible to handle discontinuous solutions.
     2835To alleviate the problems associated with numerical instabilities due to
     2836small water depths near a wet/dry boundary we employ a new flux limiter that
     2837ensures that unphysical fluxes are never encounted.
     2840Let $u$ and $v$ be the velocity components in the $x$ and $y$ direction,
     2841$w$ the absolute water level (stage) and
     2842$z$ the bed elevation. The latter are assumed to be relative to the
     2843same height datum.
     2844The conserved quantities tracked by ANUGA are momentum in the
     2845$x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$)
     2846and depth ($h = w-z$).
     2848The flux calculation requires access to the velocity vector $(u, v)$
     2849where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively.
     2850In the presence of very small water depths, these calculations become
     2851numerically unreliable and will typically cause unphysical speeds.
     2853We have employed a flux limiter which replaces the calculations above with
     2854the limited approximations.
     2856  \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h},
     2858where $h_0$ is a regularisation parameter that controls the minimal
     2859magnitude of the denominator. Taking the limits we have for $\hat{u}$
     2861  \lim_{h \rightarrow 0} \hat{u} =
     2862  \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0
     2866  \lim_{h \rightarrow \infty} \hat{u} =
     2867  \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u
     2869with similar results for $\hat{v}$.
     2871The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator)
     2873  1 - h_0/h^2 = 0
     2877  h_0 = h^2
     2881ANUGA has a global parameter $H_0$ that controls the minimal depth which
     2882is considered in the various equations. This parameter is typically set to
     2883$10^{-3}$. Setting
     2885  h_0 = H_0^2
     2887provides a reasonable balance between accurracy and stability. In fact,
     2888setting $h=H_0$ will scale the predicted speed by a factor of $0.5$:
     2890  \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = \frac{\mu}{2 H_0}
     2892In general, for multiples of the minimal depth $N H_0$ one obtains
     2894  \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} =
     2895  \frac{\mu}{H_0 (1 + 1/N^2)}
     2897which converges quadratically to the true value with the multiple N.
     2900%The developed numerical model has been applied to several test cases as well as to real flows. Numerical tests prove the robustness and accuracy of the model.
     2906\section{Slope limiting}
     2907A multidimensional slope-limiting technique is employed to achieve second-order spatial accuracy and to prevent spurious oscillations. This is using the MinMod limiter and is documented elsewhere.
     2909However close to the bed, the limiter must ensure that no negative depths occur. On the other hand, in deep water, the bed topography should be ignored for the purpose of the limiter.
     2912Let $w, z, h$  be the stage, bed elevation and depth at the centroid and
     2913let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$.
     2914Define the minimal depth across all vertices as $\hmin$ as
     2916  \hmin = \min_i h_i
     2919Let $\tilde{w_i}$ be the stage obtained from a gradient limiter
     2920limiting on stage only. The corresponding depth is then defined as
     2922  \tilde{h_i} = \tilde{w_i} - z_i
     2924We would use this limiter in deep water which we will define (somewhat boldly)
     2927  \hmin \ge \epsilon
     2931Similarly, let $\bar{w_i}$ be the stage obtained from a gradient
     2932limiter limiting on depth respecting the bed slope.
     2933The corresponding depth is defined as
     2935  \bar{h_i} = \bar{w_i} - z_i
     2939We introduce the concept of a balanced stage $w_i$ which is obtained as
     2940the linear combination
     2943  w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i}
     2947  w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i}
     2949where $\alpha \in [0, 1]$.
     2951Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope
     2952is ignored we have immediately that
     2954  \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0
     2956%where the maximal bed elevation range $dz$ is defined as
     2958%  dz = \max_i |z_i - z|
     2961If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that
     2962no negative depths occur. Formally, we will require that
     2964  \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i
     2968  \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i
     2969  \label{eq:limiter bound}
     2972There are two cases:
     2974  \item $\bar{h_i} \le \tilde{h_i}$: The deep water (limited using stage)
     2975  vertex is at least as far away from the bed than the shallow water
     2976  (limited using depth). In this case we won't need any contribution from
     2977  $\bar{h_i}$ and can accept any $alpha$.
     2979  E.g.\ $\alpha=1$ reduces Equation \ref{eq:limiter bound} to
     2980  \[
     2981    \tilde{h_i} > \epsilon 
     2982  \]
     2983  whereas $\alpha=0$ yields
     2984  \[
     2985    \bar{h_i} > \epsilon       
     2986  \]
     2987  all well and good.
     2988  \item $\bar{h_i} > \tilde{h_i}$: In this case the the deep water vertex is
     2989  closer to the bed than the shallow water vertex or even below the bed.
     2990  In this case we need to find an $alpha$ that will ensure a positive depth.   
     2991  Rearranging Equation \ref{eq:limiter bound} and solving for $\alpha$ one
     2992  obtains the bound
     2993  \[
     2994    \alpha < \frac{\epsilon - \bar{h_i}}{\tilde{h_i} - \bar{h_i}}, \forall i
     2995  \]
     2998Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one
     2999arrives at the definition
     3001  \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}}
     3003which will guarantee that no vertex 'cuts' through the bed. Finally, should
     3004$\bar{h_i} < \epsilon$ and therefore $\alpha < 0$, we suggest setting
     3005$alpha=0$ and similarly capping $\alpha$ at 1 just in case.
     3008%dropping the $\epsilon$ ensures that alpha is always positive and also 
     3009%provides a numerical safety {??)
  • anuga_core/documentation/user_manual/definitions.tex

    r2895 r4265  
    1010%\newcommand{\anugav}[1]{\textbf{ANUGA}$^\copyright$ V#1\ }
    1235%Used with \code.
  • anuga_core/documentation/user_manual/old_pyvolution_documentation/limiting.tex

    r4236 r4265  
    58\newcommand{\hmin}{h_{\mbox{\tiny min}}}
Note: See TracChangeset for help on using the changeset viewer.