|
|
| (3 intermediate revisions by the same user not shown) |
| Line 1: |
Line 1: |
| − | <math>
| |
| | | | |
| − |
| |
| − | \documentclass[12pt]{book}
| |
| − | \usepackage{amsmath}
| |
| − | \usepackage{epsfig}
| |
| − | %\usepackage{graphicx}
| |
| − | \topmargin0truemm
| |
| − | \oddsidemargin0truemm
| |
| − | \textheight=210truemm
| |
| − | \flushbottom
| |
| − | \footskip15truemm
| |
| − | \textwidth=150truemm
| |
| − | \renewcommand{\baselinestretch}{1.5}
| |
| − | \newcommand{\pd}{\partial}
| |
| − | \newcommand{\beq}{\begin{equation}}
| |
| − | \newcommand{\beqa}{\begin{eqnarray}}
| |
| − | \newcommand{\eeq}{\end{equation}}
| |
| − | \newcommand{\eeqa}{\end{eqnarray}}
| |
| − | \newcommand{\mi}{\mbox{i}}
| |
| − | \newcommand{\kp}{\kappa}
| |
| − | \newcommand{\md}{\mbox{d}}
| |
| − | \newcommand{\bs}{\boldsymbol}
| |
| − | \newcommand{\ba}{\begin{array}}
| |
| − | \newcommand{\ea}{\end{array}}
| |
| − | \providecommand{\D}{\Delta}
| |
| − | \providecommand{\nn}{\nonumber}
| |
| − | \title{\sc SELFE Theory Manual}
| |
| − | \author{Yinglong Joseph Zhang}
| |
| − | % \\Center for Coastal Margin Observation \and Prediction (CMOP)\\
| |
| − | % \\Oregon Health \and Science University}
| |
| − | \date{\today}
| |
| − |
| |
| − | \begin{document}
| |
| − | \maketitle
| |
| − |
| |
| − | \chapter{Transport}
| |
| − | \section{Notation on symbols used in this manual}
| |
| − |
| |
| − | We typically use bold characters to denote vectors and matrices, and unbold
| |
| − | characters to denote scalars.
| |
| − |
| |
| − | \section{Transport equation}
| |
| − | \beqa
| |
| − | \frac{D T}{D t}&=&\frac{\pd }{\pd z} (\kappa \frac{\pd T}{\pd z}) +\hat{Q},\label{tr1}\\
| |
| − | \kappa \frac{\pd T}{\pd z} &=& \hat{T}, \mbox{ at } z=\eta, \label{tr2} \\
| |
| − | \kappa \frac{\pd T}{\pd z} &=& \hat{T_b}, \mbox{ at } z=-h, \label{tr3}
| |
| − | \eeqa
| |
| − |
| |
| − | where $T$ is a generic scalar variable, and the total derivative is written in a conservative form:
| |
| − | \beq
| |
| − | \frac{D }{D t}=\frac{\pd }{\pd t}+\nabla \cdot (\bs{u}T) \label{tr4}
| |
| − | \eeq
| |
| − |
| |
| − | and the 3D velocity $\bs{u}$ is divergence free:
| |
| − | \beq
| |
| − | \nabla \cdot (\bs{u})=0 \label{tr5}
| |
| − | \eeq
| |
| − |
| |
| − | Eq. (\ref{tr5}) guarentees the constancy condition for the transport equation i.e. $T=const$ initially
| |
| − | will remain so in the absence of sinks/sources.
| |
| − |
| |
| − |
| |
| − | \section{Finite volume for 3D continuity equation}
| |
| − | We use F.V. method to solve the 3D continuity eq. (\ref{tr5}), not only to enforce volume conservation,
| |
| − | but also serves as a foundation for mass conservation due to the reason stated above.
| |
| − |
| |
| − | The continuity eq. is solved prism-wise from bottom to surface; note that the top and bottom faces
| |
| − | of the prism are not necessarily horizontal, and so the normal velocity there has horizontal
| |
| − | components in it (cf. Fig. \ref{tr6}).
| |
| − | \beqa
| |
| − | & &\sum_{m=1}^3 \hat{P}_{js(i,m),l} s_{i,m} \frac{q_{js(i,m),l}+q_{js(i,m),l+1}}{2} +
| |
| − | \hat{S}_{i,l+1} (\bar{u}_{i,l+1} n_{l+1}^x +\bar{v}_{i,l+1} n_{l+1}^y+\bar{w}_{i,l+1} n_{l+1}^z)- \nn \\
| |
| − | & & \hat{S}_{i,l}(\bar{u}_{i,l} n_{l}^x +\bar{v}_{i,l} n_{l}^y+\bar{w}_{i,l} n_{l}^z)=0 \nn \\
| |
| − | & & (l=k_b^e,...,N_z-1) \label{tr7}
| |
| − | \eeqa
| |
| − | where $\hat{P}$ and $\hat{S}$ are areas of horizontal and vertical faces of the prism, $i$ is the element number,
| |
| − | $js(i,m)$ is the side number, $s_{i,m}$ is a sign
| |
| − | related to the outer normal, $k_b^e$ is the bottom index of the element,
| |
| − | $(n_l^x,n_l^y,n_l^z)$ is the unit upward normal on the 2 vertical faces,
| |
| − | $q$ is the normal velocity at side (with positive direction following internal convention of SELFE, and thus
| |
| − | the sign $s_{i,m}$), and $(\bar{u},\bar{v})$ is the horizontal velocity, averaged from 3 sidecenter velocities
| |
| − | (which are the 'original' velocities defined in SELFE), on 2 vertical faces. The vertical velocity is then
| |
| − | solved from bottom to surface starting from the bottom no-flux b.c.
| |
| − |
| |
| − | The eq. (\ref{tr7}) can be written in the following compact statement for volume conservation:
| |
| − | \beq
| |
| − | \sum_{j\in S^+} |Q_j|=\sum_{j\in S^-} |Q_j| \label{tr8}
| |
| − | \eeq
| |
| − | where $Q=u_n A$ is the flux, $A$ is the area of the face, and $S^\pm$ are the out/inflow faces.
| |
| − | This eq. simply states that the incoming and outging fluxes must be balanced out against each
| |
| − | other; this will be used below.
| |
| − |
| |
| − | \section{Upwind scheme}
| |
| − |
| |
| − | A F.V. discretization of eq. (\ref{tr1}) for a prism shown in Fig. \ref{tr6} is:
| |
| − | \beqa
| |
| − | T_i^{m+1} &=& T_i^{m}-\frac{\D t'}{V_i}\sum_{j\in S}Q_j T_{j*} +\hat{Q}_{i,k}^{n+1}\D t'+ \nn \\
| |
| − | && \frac{A_i \D t'}{V_i}[\kappa_{i,k}\frac{T_{i,k+1}^{m+1}-T_{i,k}^{m+1}}{\D z_{i,k+1/2}}-
| |
| − | \kappa_{i,k-1}\frac{T_{i,k}^{m+1}-T_{i,k-1}^{m+1}}{\D z_{i,k-1/2}}] \nn \\
| |
| − | && (k=k_b^e+1,...,N_z)\label{tr9}
| |
| − | \eeqa
| |
| − | where $\D t' \neq \D t$ is the trasnport step (subject to Courant condition below), $V_i$ is the volume of the prism,
| |
| − | $T_i$ is a shorthand for $T_{i,k}$ (i.e., concentration at prism $(i,k)$).
| |
| − | Note that we have treated the diffusion term implicitly. For the sake of brevity
| |
| − | we'll drop the source and diffusion terms from now on and focus on the advection term, where different methods diverge.
| |
| − |
| |
| − | The upwind concentration is defined as:
| |
| − | \beq
| |
| − | T_{j*}\equiv \tilde{T}_{j*}=
| |
| − | \left\{
| |
| − | \ba{ll}
| |
| − | T_i, & j\in S^+ \\
| |
| − | T_j, & j\in S^-
| |
| − | \ea
| |
| − | \right.
| |
| − | \eeq
| |
| − | where again we have used shorthand $T_j$ for concentration at prism $j$ (i.e. the prism adjacent to $(i,k)$
| |
| − | from face $j$).
| |
| − |
| |
| − | Eq. (\ref{tr9}) then becomes (only retaining the advective term):
| |
| − | \beqa
| |
| − | T_i^{m+1} &=& T_i^{m}-\frac{\D t'}{V_i}[\sum_{j\in S^+}|Q_j|T_i-\sum_{j\in S^-}|Q_j|T_j] \nn \\
| |
| − | &=& \left( 1-\frac{\D t'}{V_i}\sum_{j\in S^+}|Q_j|\right) T_i + \frac{\D t'}{V_i} \sum_{j\in S^+}|Q_j|T_j \label{tr10}
| |
| − | \eeqa
| |
| − | ($T$ without superscript are evaluated at step $m$, i.e. explicitly). We have used the volume conservation
| |
| − | statement (\ref{tr8}) here. Therefore the Courant condition is:
| |
| − | \beq
| |
| − | 1-\frac{\D t'}{V_i}\sum_{j\in S^+}|Q_j| \geq 0 \label{tr11}
| |
| − | \eeq
| |
| − | SELFE uses this eq. to dynamically adjust time step for transport for each step. Moreover, to improve
| |
| − | efficiency, the vertical flux terms in (\ref{tr10}) are treated implicitly, and the corresponding
| |
| − | terms are then removed in the Courant condition (\ref{tr11}). This is allowable because upwind
| |
| − | is a linear method.
| |
| − |
| |
| − | \begin{figure} [100]
| |
| − | \centering
| |
| − | \includegraphics{transport-prism.eps}
| |
| − | \caption{Prism and variable locations in SELFE}
| |
| − | \label{tr6}
| |
| − | \end{figure}
| |
| − |
| |
| − | \section{TVD scheme}
| |
| − |
| |
| − | The only difference between TVD and upwind schemes lies in the evaluation of the interfacial concentation:
| |
| − | \beq
| |
| − | T_{j*}=\tilde{T}_{j*}+\frac{\phi_j}{2}(T_{jD}-\tilde{T}_{j*}) \label{tr12}
| |
| − | \eeq
| |
| − | where $T_{jD}$ is the downstream concentration, and $0\leq \phi_j \leq 2$ is a limiter function. TVD
| |
| − | scheme nominally has 2nd order accuracy due to the anti-diffusion term in (\ref{tr12}).
| |
| − |
| |
| − | The final eq. for TVD is:
| |
| − | \beqa
| |
| − | T_i^{m+1} &=& T_i^{m}+\frac{\D t'}{V_i}\sum_{j\in S^-}|Q_j|(T_j-T_i) + \nn \\
| |
| − | &&\frac{\D t'}{V_i}\sum_j |Q_j|\frac{\phi_j}{2}(T_i-T_j)+\mbox{souce}+\mbox{diffusion} \label{tr13}
| |
| − | \eeqa
| |
| − | and the Courant condition is (option 'AA' in param.in):
| |
| − | \beq
| |
| − | \D t' \leq \frac{V_i}{\sum_{j\in S^-} |Q_j| (1-\frac{\phi_j}{2}+\delta_i)} \label{tr14}
| |
| − | \eeq
| |
| − | where
| |
| − | \beqa
| |
| − | \delta_i &=& \sum_{p\in S^+} \frac{\phi (r_p)}{2r_p} \nn \\
| |
| − | r_p &=& \frac{\sum_{m\in S^-} |Q_m| (T_m-T_i)}{|Q_p|(T_i-T_p)}, p\in S_i^+ \label{tr15}
| |
| − | \eeqa
| |
| − | The 2nd eq. involves upwind of upwind neighboring prism.
| |
| − |
| |
| − | Again SELFE automatically calculates the time step accroding to the Courant condition above.
| |
| − | The choices for the limiter function include: MINMOD, OSHER, van Leer, Super Bee etc.
| |
| − |
| |
| − | Since TVD is a nonlinear method, we cannot treat the vertical fluxes implicitly, and so all fluxes
| |
| − | have to be treated explicitly. TVD method is therefore more expensive than upwind. A hybrid upwind/TVD,
| |
| − | with TVD being used in the deeper depths and upwind in the shallow depths, has been implemented in SELFE
| |
| − | to improve efficiency. The user can also manually specify upwind/TVD zones in the domain.
| |
| − |
| |
| − | \section{Boundary conditions}
| |
| − |
| |
| − | In either upwiund or TVD schemes, the concentration at the neighboring prism $T_j$ at the open boundary
| |
| − | is known. For outflow, $T_j=T_i$ and the signal is advected out of the domain without hindrance.
| |
| − | For incoming flow, $T_j$ is specified by the b.c., and SELFE nudges to this value with a relaxation constant,
| |
| − | in order to prevent sharp gradient there.
| |
| − |
| |
| − | \section{Interface for the generic tracer transport module}
| |
| − |
| |
| − | In SELFE v3.2.x, it's easy to use SELFE's generic tracer transport module by using flag\_model=-1 in param.in. The
| |
| − | 3 arrays in selfe\_step.F90, $bdy\_frc,flx\_sf,flx\_bt$ correspond to the source, surface/bottom b.c. in Eqs.
| |
| − | (\ref{tr1}-\ref{tr3}). You may use SELFE's conventional way of specifying i.c. (flag\_ic in param.in) or
| |
| − | choose to specify your own by modifying sections in selfe\_init.F90.
| |
| − |
| |
| − | \end{document}
| |
| − |
| |
| − | </math>
| |
