sandbox/HCH/document/bcg-doc

Introduction and Background

bcg.h file contains two functions: tracer_fluxes and advection whose major purpose is to construct a solver for advective equation: \displaystyle \frac{\partial \Phi}{\partial t} + ( \mathbf{u} \cdot \nabla)\Phi = 0 \quad (1) where \Phi is the scalar and \mathbf{u} is the velocity. The discrete time formulation of equation (1) reads: \displaystyle \frac{\Phi^{n+1}-\Phi^{n}}{\Delta t} + \mathbf{A}^{n+\frac{1}{2}} = 0 \quad (2) Where \mathbf{A}^{n+ \frac{1}{2}} is the abbreviation of advection term in BCG (see Bell et al, 1989 and Popinet, 2003) scheme (and is the reason why the name of this file called ‘bcg.h’), and is intergrated within the controlled volume \displaystyle \int_{\Gamma} A^{n+ \frac{1}{2}} = \int_{\Gamma} [( \mathbf{u}\cdot \nabla)\Phi]^{n+ \frac{1}{2}} = \int_{\Gamma} [\nabla\cdot( \mathbf{u} \Phi)]^{n+ \frac{1}{2}} = \int_{\partial \Gamma} ( \mathbf{u}^{n + \frac{1}{2}} \cdot \mathbf{n}) \Phi^{n+ \frac{1}{2}}\quad (3)

Equation test1 Equation test2 The second step above is achieved by combining divergence free constraints \nabla\cdot \mathbf{u}=0 and the \mathbf{n} represents the normal direction of cell interface. For quadtree or octree cells, the overall calculation turns out to be \displaystyle \Delta A^{n+ \frac{1}{2}} = \sum_{i=d} u_d^{n+ \frac{1}{2}}\Phi_d^{n + \frac{1}{2}} \quad (4) where u_d is the component of \mathbf{u} on normal direction of interface, Delta is the length of the cell and \Phi_d is the corresponding value at the interface. Then the critical is how to obtain \Phi_d^{n + \frac{1}{2}}. According to Popinet, 2003, using Taylor expansion we have \displaystyle \Phi^{n+ \frac{1}{2}}_d = \Phi^n+ \frac{\Delta}{2} \frac{\partial \Phi^n}{\partial x_d} + \frac{\Delta t}{2} \frac{\partial \Phi^n}{\partial t} + O(\Delta^2,\Delta t^2)\quad (5) Replacing \frac{\partial \Phi^n}{\partial x_d} with equation (1) yielding \displaystyle \Phi^{n+ \frac{1}{2}}_d = \Phi^n + [\frac{\Delta}{2}- \frac{\Delta t}{2} \mathbf{u}^n \cdot \mathbf{n}_d] \frac{\partial \Phi^n}{\partial x_d} - \frac{\Delta t}{2} \mathbf{u}^n\cdot \mathbf{n}_e \frac{\partial \Phi^n}{\partial x_e} - \frac{\Delta t}{2} \mathbf{g}^n \quad (6)
Subscript e represents direction other than current direction d. Compromise is taken for sake of convenience (maybe) that \mathbf{u}^n in equation (6) is replaced by \mathbf{u}^{n + \frac{1}{2}} so that \mathbf{u}^n\cdot \mathbf{n}_d can be computed by u_d^{n + \frac{1}{2}}.

tracer_fluxes

Data of face vector type is stored on staggered mesh Harlow and Welch 1965. Take as an example, the storage on a 2D cell displays in figure 1.

Figure 1: Staggered mesh.

Components are saved in corresponding face with single value. The direction of the component express itself in the feature of the cell face. Therefore, \mathbf{u}^{n + \frac{1}{2}} \cdot \mathbf{n}_x can simply be \frac{1}{2}(uf_x[i]+uf_x[i+1]).

Advection term employs BCG schemeMartin et al, 2000, a second order upwinded scheme. It first requires face-centered approximation of velocity u_d^n,v_e^n(w_f^n). Then equation (6) is adapted to a upwinded scheme \displaystyle \Phi_d^{n+ \frac{1}{2}}= \left \{ \begin{array}{cc} \bar{\Phi}^{L,n+ \frac{1}{2}}_d\quad &if\, u_d^n>0,\\ \bar{\Phi}^{R,n+ \frac{1}{2}}_d\quad &if\, u_d^n<0,\\ \frac{1}{2}(\bar{\Phi}^{L,n+ \frac{1}{2}}_d+\bar{\Phi}^{R,n+ \frac{1}{2}}_d)\quad &if \,u_d^n=0. \end{array} \right. (7) where \displaystyle \begin{aligned} \bar{\Phi}^{L,n+ \frac{1}{2}}_d &= \Phi^n[i] + \frac{\Delta}{2}min[1- \frac{\Delta t}{2\Delta}(u_d^n[i]+u_d^n[i+1]),1] \frac{\partial \Phi}{\partial x_d}[i]- \frac{\Delta t}{2} \mathbf{g}^n-flux_e[i]\quad (8)\\ \bar{\Phi}^{R,n+ \frac{1}{2}}_d &= \Phi^n[i+1] - \frac{\Delta}{2}min[1- \frac{\Delta t}{2\Delta}(u_d^n[i]+u_d^n[i+1]),1]\frac{\partial \Phi}{\partial x_d}[i+1]- \frac{\Delta t}{2} \mathbf{g}^n-flux_e[i+1]\quad (9) \end{aligned} flux_e represents contribution from direction other than d, which also takes upwinded scheme \displaystyle flux_e = \frac{\Delta t}{2} \mathbf{u}^n\cdot\mathbf{n}_e \frac{\partial \Phi^n}{\partial x_e}=\left \{ \begin{array}{cc} \frac{\Delta t}{2\Delta} v_{e}^{trans}(u^n[i,j]-u^n[i,j-1])\quad if\,v_{e}^{trans}>0,\\ -\frac{\Delta t}{2\Delta} v_e^{trans}(u^n[i,j]-u^n[i,j+1])\quad if\,v_e^{trans}<0. \end{array} \right. (10) where v_e^{trans} = \frac{1}{2}(v_e^{trans}[i,j]+v_e^{trans}[i,j+1])

advection

The input of this function (tracers,scr) can be vector type, e.g. \mathbf{u}, where components of each direction is deemed as scalar type data then is finally assembled as a vector.
Another thing needs to be concerned is the direction of the normal unity associates with coordinate axis. Positive direction of coordinate remains unchanged will face normal direction varies with respect to current cell. Take figure 2 as an example.

Figure 2: Facenormal.

For cell A, direction of highlighted face is opposite to coordinate while positive for that of cell B. For the sake of clarity and convenience, all value involved compute in function tracer_fluxes take default positive direction. Which means for further computation for \sum_{i=d}u_d^{n+ \frac{1}{2}}\Phi_d^{n+ \frac{1}{2}} of cell A, the value added at highlighted face should take negative value reads \displaystyle u_d^{n+ \frac{1}{2}}\Phi_d^{n + \frac{1}{2}}(i)+u_d^{n+ \frac{1}{2}}\Phi_d^{n + \frac{1}{2}}(i+1) = -flux.x[i]+flux.x[i+1]\quad (11)

Appendix: Calculation of Face Centered Normal Velocity

It can be seen from previous discussion, the computation of face centered normal velocity u_d^{n+ \frac{1}{2}} is of great importance to construct advection term. The detailed procedure of extrapolation is shown in documentation of centered,h. The output of such variable follows a similar step as equation (6) but without pressure term since the face centered velocity will be projected with an edge-centered projection. This additional projection ensures the feature of conservative of corresponding method.

When it comes to the resolution of NS equation, \Phi in governing equation (1) is the cell-centered vector \mathbf{u}^n. According to Martin et al, 2000, normal components calculated by equation (6) can simply be replaced by u_d^{n+ \frac{1}{2}}, which means we only need to compute tangential fluxes. However, after series of tests, reuse of normal fluxes will lead to unstable at sharp angles. Thus all components is recomputed in Basilisk (see Popinet, 2003).

References