Changeset 4265
- Timestamp:
- Feb 16, 2007, 5:26:04 PM (18 years ago)
- Location:
- anuga_core/documentation/user_manual
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/documentation/user_manual/anuga_user_manual.tex
r4258 r4265 2660 2660 Elevation Model) and converts it to PTS format. 2661 2661 \end{funcdesc} 2662 2663 2664 2665 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2666 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2667 2668 \chapter{\anuga mathematical background} 2669 \label{cd:mathematical background} 2670 2671 \section{Introduction} 2672 2673 This chapter outlines the mathematics underpinning \anuga. 2674 2675 2676 2677 \section{Model} 2678 \label{sec:model} 2679 2680 The shallow water wave equations are a system of differential 2681 conservation equations which describe the flow of a thin layer of 2682 fluid over terrain. The form of the equations are: 2683 \[ 2684 \frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial 2685 x}+\frac{\partial \GG}{\partial y}=\SSS 2686 \] 2687 where $\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 2691 entering the system are bed elevation $z$ and stage (absolute water 2692 level) $w$, where the relation $w = z + h$ holds true at all times. 2693 The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given 2694 by 2695 \[ 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] 2705 \] 2706 and the source term (which includes gravity and friction) is given 2707 by 2708 \[ 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] 2714 \] 2715 where $S_f$ is the bed friction. The friction term is modelled using 2716 Manning's resistance law 2717 \[ 2718 S_{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}} 2720 \] 2721 in which $\eta$ is the Manning resistance coefficient. 2722 2723 As demonstrated in our papers, \cite{modsim2005,Rob99l} these 2724 equations provide an excellent model of flows associated with 2725 inundation such as dam breaks and tsunamis. 2726 2727 \section{Finite Volume Method} 2728 \label{sec:fvm} 2729 2730 We use a finite-volume method for solving the shallow water wave 2731 equations \cite{Rob99l}. The study area is represented by a mesh of 2732 triangular cells as in Figure~\ref{fig:mesh} in which the conserved 2733 quantities of water depth $h$, and horizontal momentum $(uh, vh)$, 2734 in each volume are to be determined. The size of the triangles may 2735 be varied within the mesh to allow greater resolution in regions of 2736 particular interest. 2737 2738 \begin{figure} 2739 \begin{center} 2740 \includegraphics[width=5.0cm,keepaspectratio=true]{step-five} 2741 \caption{Triangular mesh used in our finite volume method. Conserved 2742 quantities $h$, $uh$ and $vh$ are associated with the centroid of 2743 each triangular cell.} \label{fig:mesh} 2744 \end{center} 2745 \end{figure} 2746 2747 The equations constituting the finite-volume method are obtained by 2748 integrating the differential conservation equations over each 2749 triangular cell of the mesh. Introducing some notation we use $i$ to 2750 refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the 2751 set of indices referring to the cells neighbouring the $i$th cell. 2752 Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is 2753 the length of the edge between the $i$th and $j$th cells. 2754 2755 By applying the divergence theorem we obtain for each volume an 2756 equation which describes the rate of change of the average of the 2757 conserved quantities within each cell, in terms of the fluxes across 2758 the edges of the cells and the effect of the source terms. In 2759 particular, rate equations associated with each cell have the form 2760 $$ 2761 \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i 2762 $$ 2763 where 2764 \begin{itemize} 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, 2767 and 2768 \item $\HH_{ij}$ is the outward normal flux of 2769 material across the \textit{ij}th edge. 2770 \end{itemize} 2771 2772 2773 %\item $l_{ij}$ is the length of the edge between the $i$th and $j$th 2774 %cells 2775 %\item $m_{ij}$ is the midpoint of 2776 %the \textit{ij}th edge, 2777 %\item 2778 %$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing 2779 %normal along the \textit{ij}th edge, and The 2780 2781 The flux $\HH_{ij}$ is evaluated using a numerical flux function 2782 $\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow 2783 water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$ 2784 $$ 2785 H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 . 2786 $$ 2787 2788 Then 2789 $$ 2790 \HH_{ij} = \HH(\UU_i(m_{ij}), 2791 \UU_j(m_{ij}); \mathbf{n}_{ij}) 2792 $$ 2793 where $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 2796 T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and 2797 neighbouring cells. 2798 2799 We use a second order reconstruction to produce a piece-wise linear 2800 function construction of the conserved quantities for all $x \in 2801 T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This 2802 function is allowed to be discontinuous across the edges of the 2803 cells, but the slope of this function is limited to avoid 2804 artificially introduced oscillations. 2805 2806 Godunov's method (see \cite{Toro-92}) involves calculating the 2807 numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly 2808 solving the corresponding one dimensional Riemann problem normal to 2809 the edge. We use the central-upwind scheme of \cite{KurNP2001} to 2810 calculate an approximation of the flux across each edge. 2811 2812 \begin{figure} 2813 \begin{center} 2814 \includegraphics[width=5.0cm,keepaspectratio=true]{step-reconstruct} 2815 \caption{From the values of the conserved quantities at the centroid 2816 of the cell and its neighbouring cells, a discontinuous piecewise 2817 linear reconstruction of the conserved quantities is obtained.} 2818 \label{fig:mesh:reconstruct} 2819 \end{center} 2820 \end{figure} 2821 2822 In the computations presented in this paper we use an explicit Euler 2823 time stepping method with variable timestepping adapted to the 2824 observed CFL condition. 2825 2826 2827 \section{Flux limiting} 2828 2829 The shallow water equations are solved numerically using a 2830 finite volume method on unstructured triangular grid. 2831 The upwind central scheme due to Kurganov and Petrova is used as an 2832 approximate Riemann solver for the computation of inviscid flux functions. 2833 This makes it possible to handle discontinuous solutions. 2834 2835 To alleviate the problems associated with numerical instabilities due to 2836 small water depths near a wet/dry boundary we employ a new flux limiter that 2837 ensures that unphysical fluxes are never encounted. 2838 2839 2840 Let $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 2843 same height datum. 2844 The conserved quantities tracked by ANUGA are momentum in the 2845 $x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$) 2846 and depth ($h = w-z$). 2847 2848 The flux calculation requires access to the velocity vector $(u, v)$ 2849 where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively. 2850 In the presence of very small water depths, these calculations become 2851 numerically unreliable and will typically cause unphysical speeds. 2852 2853 We have employed a flux limiter which replaces the calculations above with 2854 the limited approximations. 2855 \begin{equation} 2856 \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h}, 2857 \end{equation} 2858 where $h_0$ is a regularisation parameter that controls the minimal 2859 magnitude of the denominator. Taking the limits we have for $\hat{u}$ 2860 \[ 2861 \lim_{h \rightarrow 0} \hat{u} = 2862 \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0 2863 \] 2864 and 2865 \[ 2866 \lim_{h \rightarrow \infty} \hat{u} = 2867 \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u 2868 \] 2869 with similar results for $\hat{v}$. 2870 2871 The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator) 2872 \[ 2873 1 - h_0/h^2 = 0 2874 \] 2875 or 2876 \[ 2877 h_0 = h^2 2878 \] 2879 2880 2881 ANUGA has a global parameter $H_0$ that controls the minimal depth which 2882 is considered in the various equations. This parameter is typically set to 2883 $10^{-3}$. Setting 2884 \[ 2885 h_0 = H_0^2 2886 \] 2887 provides a reasonable balance between accurracy and stability. In fact, 2888 setting $h=H_0$ will scale the predicted speed by a factor of $0.5$: 2889 \[ 2890 \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = \frac{\mu}{2 H_0} 2891 \] 2892 In general, for multiples of the minimal depth $N H_0$ one obtains 2893 \[ 2894 \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} = 2895 \frac{\mu}{H_0 (1 + 1/N^2)} 2896 \] 2897 which converges quadratically to the true value with the multiple N. 2898 2899 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. 2901 2902 2903 2904 2905 2906 \section{Slope limiting} 2907 A 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. 2908 2909 However 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. 2910 2911 2912 Let $w, z, h$ be the stage, bed elevation and depth at the centroid and 2913 let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$. 2914 Define the minimal depth across all vertices as $\hmin$ as 2915 \[ 2916 \hmin = \min_i h_i 2917 \] 2918 2919 Let $\tilde{w_i}$ be the stage obtained from a gradient limiter 2920 limiting on stage only. The corresponding depth is then defined as 2921 \[ 2922 \tilde{h_i} = \tilde{w_i} - z_i 2923 \] 2924 We would use this limiter in deep water which we will define (somewhat boldly) 2925 as 2926 \[ 2927 \hmin \ge \epsilon 2928 \] 2929 2930 2931 Similarly, let $\bar{w_i}$ be the stage obtained from a gradient 2932 limiter limiting on depth respecting the bed slope. 2933 The corresponding depth is defined as 2934 \[ 2935 \bar{h_i} = \bar{w_i} - z_i 2936 \] 2937 2938 2939 We introduce the concept of a balanced stage $w_i$ which is obtained as 2940 the linear combination 2941 2942 \[ 2943 w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i} 2944 \] 2945 or 2946 \[ 2947 w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} 2948 \] 2949 where $\alpha \in [0, 1]$. 2950 2951 Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope 2952 is ignored we have immediately that 2953 \[ 2954 \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0 2955 \] 2956 %where the maximal bed elevation range $dz$ is defined as 2957 %\[ 2958 % dz = \max_i |z_i - z| 2959 %\] 2960 2961 If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that 2962 no negative depths occur. Formally, we will require that 2963 \[ 2964 \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i 2965 \] 2966 or 2967 \begin{equation} 2968 \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i 2969 \label{eq:limiter bound} 2970 \end{equation} 2971 2972 There are two cases: 2973 \begin{enumerate} 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$. 2978 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 \] 2996 \end{enumerate} 2997 2998 Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one 2999 arrives at the definition 3000 \[ 3001 \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}} 3002 \] 3003 which 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. 3006 3007 %Furthermore, 3008 %dropping the $\epsilon$ ensures that alpha is always positive and also 3009 %provides a numerical safety {??) 3010 3011 2662 3012 2663 3013 -
anuga_core/documentation/user_manual/definitions.tex
r2895 r4265 10 10 %\newcommand{\anugav}[1]{\textbf{ANUGA}$^\copyright$ V#1\ } 11 11 12 13 \newcommand{\Python}{\textsc{Python}} 14 \newcommand{\VPython}{\textsc{VPython}} 15 \newcommand{\pypar}{\textsc{mpi}} 16 \newcommand{\Metis}{\textsc{Metis}} 17 \newcommand{\mpi}{\textsc{mpi}} 18 19 \newcommand{\UU}{\mathbf{U}} 20 \newcommand{\VV}{\mathbf{V}} 21 \newcommand{\EE}{\mathbf{E}} 22 \newcommand{\GG}{\mathbf{G}} 23 \newcommand{\FF}{\mathbf{F}} 24 \newcommand{\HH}{\mathbf{H}} 25 \newcommand{\SSS}{\mathbf{S}} 26 \newcommand{\nn}{\mathbf{n}} 27 28 %\newcommand{\code}[1]{\texttt{#1}} 29 30 31 32 33 34 12 35 %Used with \code. 13 36 % -
anuga_core/documentation/user_manual/old_pyvolution_documentation/limiting.tex
r4236 r4265 1 1 2 \documentclass[12pt]{article} 2 3 \usepackage{graphicx} 3 4 5 6 % THIS HAS BEEN MOVED INTO THE USER MANUAL (16 Feb 2007) 4 7 5 8 \newcommand{\hmin}{h_{\mbox{\tiny min}}}
Note: See TracChangeset
for help on using the changeset viewer.