I'm trying to improve the advection method in a 2D-windfield. The Navier-Stokes Equations (NSE) are currently used for the influence of pressure, viscosity,… I am just focusing on the convective acceleration term $(u\cdot\nabla)u$. So after rearranging the NSE I get $$\left(\frac{\partial \vec v}{\partial t}\right)_{\text{Adv}}=-(u\cdot\nabla)u.$$ For my case I can write this as: $$\begin{pmatrix} \frac{\partial u}{\partial t} \\ \frac{\partial v}{\partial t} \end{pmatrix}=
\begin{pmatrix} -u(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}) \\ -v(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}) \end{pmatrix}.$$
My question is, can I now apply backward differencing $$\left(\frac{\partial \Phi}{\partial x}\right)_{i}=\frac{\Phi_i-\Phi_{i-1}}{\Delta x},\quad\left(\frac{\partial u}{\partial t}\right)_{i}=\frac{u_i^{t+1}-u_i^t}{\Delta t}~?$$ My final two equation then:
$$u_{x,y}^{t+1}=u_{x,y}^{t}-\frac{u_{x,y}^{t}\Delta t}{\Delta x}(u_{x,y}^{t}-u_{x-1,y}^{t}+v_{x,y}^{t}-v_{x,y-1}^{t}) $$
$$v_{x,y}^{t+1}=v_{x,y}^{t}-\frac{v_{x,y}^{t}\Delta t}{\Delta x}(u_{x,y}^{t}-u_{x-1,y}^{t}+v_{x,y}^{t}-v_{x,y-1}^{t}). $$
Now I have an equation for the velocity $u$ in $x$-direction and $v$ in $y$-direction.
It's not working properly in the simulation. For example, if the direction of the wind should change over time, the direction stays the same but the velocity on the edges of the field (where the initial conditions of $u$ and $v$ are set) are getting bigger or smaller, depending on the "direction of the change". It seems like the information is not spreading properly over the whole wind field.
Best Answer
I advocate the view that the advection component is an operator, $$\left(\mathbf{u}\cdot\nabla\right)\mathbf{u}=\operatorname{adv}\mathbf{u};$$ hence you have, depending on what is most natural to you, $$\operatorname{adv}\mathbf{u}\equiv\operatorname{adv}\left(\begin{array}{c}u \\ v \end{array}\right)\equiv\left(\begin{array}{c}\operatorname{adv}u \\ \operatorname{adv}v \end{array}\right)\equiv \operatorname{adv}u_i\mathrm{e}_i.\tag{1}$$ Since $\operatorname{adv}=\mathbf{u}\cdot\nabla$, then Eq (1) is, $$\operatorname{adv}\mathbf{u}\equiv\left(u\partial_x+v\partial_y\right)\left(\begin{array}{c}u \\ v \end{array}\right)\equiv\left(\begin{array}{c}\left(u\partial_x+v\partial_y\right)u \\ \left(u\partial_x+v\partial_y\right)v \end{array}\right).\tag{2}$$ Then applying the backwards finite difference to the derivatives, $$\operatorname{adv}\mathbf{u}=\left(\begin{array}{c} u_{i,j}\left(u_{i,j}-u_{i-1,j}\right)/\Delta x +v_{i,j}\left(u_{i,j}-u_{i,j-1}\right)/\Delta y\\ u_{i,j}\left(v_{i,j}-v_{i-1,j}\right)/\Delta x +v_{i,j}\left(v_{i,j}-v_{i,j-1}\right)/\Delta y \end{array}\right)$$ So along with your time derivative component, you'd see something like, $$u_{i,j}^{n+1}=u_{i,j}^n-\frac{\Delta t}{\Delta x}u_{i,j}^n\left(u_{i,j}^n-u_{i-1,j}^n\right)-\frac{\Delta t}{\Delta y}v_{i,j}^n\left(u_{i,j}^n-u_{i,j-1}^n\right)\tag{3}$$ and similarly for the $v$ component of the velocity.
If you wanted to increase the accuracy by using a 2nd order scheme (still with backwards differencing), you'd have something like, $$u_{i,j}^{n+1}=u_{i,j}^n-\frac{\Delta t}{\Delta x}u_{i,j}^n\left(\frac{3}{2}u_{i,j}^n-2u_{i-1,j}^n+\frac{1}{2}u_{i-2,j}^n\right)-\frac{\Delta t}{\Delta y}v_{i,j}^n\left(\frac{3}{2}u_{i,j}^n-2u_{i,j-1}^n+\frac{1}{2}u_{i,j-2}^n\right)$$ It might also be useful to look over some of the predictor-corrector methods, such as the MacCormack method, wherein you compute an intermediate state, $u^p_{i,j}$, using Equation (3) and then use that state to evaluate the next step using, $$u^{n+1}_{i,j}=\frac{1}{2}\left(u^n_{i,j}+u^p_{i,j}\right)-\frac{\Delta t}{\Delta x}u_{i,j}^p\left(u_{i,j}^p-u_{i-1,j}^p\right)-\frac{\Delta t}{\Delta y}v_{i,j}^p\left(u_{i,j}^p-u_{i,j-1}^p\right)$$