Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerical solution of continuity equation, implicit scheme, staggered grid

  1. Nov 29, 2011 #1

    I'm trying to implement an implicit scheme for the continuity equation.
    The scheme is the following:
    http://img28.imageshack.us/img28/3196/screenshot20111130at003.png [Broken]
    With [TEX]\rho[/TEX] being the density, [TEX]\alpha[/TEX] is a weighing constant. d is a parameter that relates the grid spacing to the velocity of flow. j is the spatial grid and n is the time grid.

    The problem with this scheme is the fact that it is implicit and effectively I have no idea how to successfully implement it. I tried by doing the following:
    Assume have initial time n=0 spatial (j) grid full. Also assume I know spatial boundaries (j=-1/2 and j=max) at time n+1.

    Then I set j=1/2 equal to a variable x.
    Next rearrange so that have j=3/2 in terms of variable x and boundary.
    Do likewise for all the (half)grid points so that they can be written in terms of x. then when I reach the opposite boundary of grid (j=max) solve for x. So then I have j=1/2 and hence can substitute x into all the other equations to fill the grid points with data.

    This can't work since:
    When writing expression for any grid point in terms of x, i divide by [tex]0.5d(1-\alpha)[/tex]. Which is roughly 0.5. So going through all the grid points I end up with 0.5^400 (400 grid points) in the denominator which is bound to kill the calculation.

    Can anyone tell me how to deal with this implicit scheme?
    Last edited by a moderator: May 5, 2017
  2. jcsd
  3. Dec 5, 2011 #2
    Are you sure you want to work at the 1/2 nodes? Shifting everything to the left by 1/2 and you will get a scheme with unknowns at j-1,j,j+1, which will be MUCH more convenient.

    First of all, assume everything at timestep n known. When n=0, you can substitute the initial conditions. Put all unknowns (the solution at n+1) to the left and all known terms to the right. This is generally a good idea for any numerical scheme.

    I will now for simplicity move your scheme 1/2 node to the left. The result is:
    [itex]\rho_{j} + 0.5d(1-\alpha)(\rho_{j+1}-\rho_{j-1}) = -0.5d(...)[/itex]

    When you have N nodes to discretize space, from j=0..N, then j=0 and j=N are on the boundaries.
    Suppose j=1. We then get:
    [itex]\rho_{1} + 0.5d(1-\alpha)(\rho_{2}-\rho_{0}) = -0.5d(...)[/itex]

    [itex]\rho_{0}[/itex] is on a boundary, so if a boundary condition was imposed on x=0, then you can substitute its value (and move it to the right of the equal sign).

    As long as the boundary conditions have not been substituted yet, you will have N-2 equations (j=1..N-1) for N unknowns (j=0..N).
    So you have a matrix-vector system of the form [itex][A][\rho^{n+1}]=[/itex]
    After substituting the boundary conditions, you should have a system of N-2 equations with N-2 unknowns.
    This matrix-vector system is tri-diagonal for which very efficient algorithms exist.
    If you solve this, you will get the solution at time n+1.

    Hope this helps
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook