How could I go about solving the continuity equation below in 2D?

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$$

I saw that a similar question was posted here: A good finite difference for the continuity equation, but that only discussed the 1-D case.

In it, one of the answers defines incoming, and outgoing fluxes, as $$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i} + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1}\\ \Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1} + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i}$$

With these, one could use the finite difference scheme

$$\dfrac{\partial(\rho)}{\partial t}=-\dfrac{\partial(\rho u)}{\partial x} = -\dfrac{\partial \Phi}{\partial x} \approx -\dfrac{\Delta \Phi}{\Delta x} \approx -\dfrac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x}$$

The problem is that when I use a naive generalization of this, namely defining

$$\Phi_{i+1/2,j;x} = \dfrac{u_{i+1/2,j;x}+|u_{i+1/2,j;x}|}{2}\rho_{i,j} + \dfrac{u_{i+1/2,j;x}-|u_{i+1/2,j;x}|}{2}\rho_{i+1,j}\\ \Phi_{i,j+1/2;y} = \dfrac{u_{i,j+1/2;y}+|u_{i,j+1/2;y}|}{2}\rho_{i,j} + \dfrac{u_{i,j+1/2;y}-|u_{i,j+1/2;y}|}{2}\rho_{i,j+1}$$

and taking

$$\dfrac{\partial (\rho u_x)}{\partial x} \approx \frac{\Phi_{i+1/2,j;x}-\Phi_{i-1/2,j;x}}{\Delta x}$$ and $$\dfrac{\partial (\rho u_y)}{\partial y} \approx \frac{\Phi_{i,j+1/2;y}-\Phi_{i,j-1/2;y}}{\Delta y}$$ my density diverges. What scheme could I use to solve this equation properly in 2D?

  • $\begingroup$ What you have written is an upwind scheme which is very robust. But are you solving this coupled to some other equations ? What is the time scheme, time step ? What does the velocity field look like ? $\endgroup$ – cfdlab Jun 20 '18 at 12:24
  • $\begingroup$ Hey. So I'm using a Runge-Kutta 4th order time integrator, with a time step that's $\Delta t = \Delta x/4$. The density is indeed coupled to another set of equations (Maxwell's equations) that I'm solving using a 4th order in space finite difference scheme (I'm comparing with an steady-state analytical solution, and if I plug in the correct steady-state density, the Maxwell solver works well). The velocity is a function of the electric and magnetic field, $\propto (EB-E-B)/(E^2)$. $\endgroup$ – grizzlyjoker Jun 20 '18 at 15:03
  • $\begingroup$ where $E > 0$ always $\endgroup$ – grizzlyjoker Jun 20 '18 at 19:06
  • $\begingroup$ Does your time step satisfy CFL condition ? Make sure that $\Delta t \le \frac{1}{\max(|u_x|/\Delta x + |u_y|/\Delta y)}$ where the max is taken over the whole grid. Also your spatial scheme is first order accurate, so you dont need RK4 in time. Moreover, the first order upwind scheme is too dissipative. $\endgroup$ – cfdlab Jun 21 '18 at 4:11
  • $\begingroup$ My time step does indeed satisfy the CFL condition. Ok, I'll try getting a higher order upwind scheme to work. $\endgroup$ – grizzlyjoker Jun 21 '18 at 7:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.