Algorithm for modeling expanding gases on a 2D grid

Question

I have a simple program, at it's heart is a two dimensional array of floats, supposedly representing gas concentrations, I have been trying to come up with a simple algorithm that will model the gas expanding outwards, like a cloud, eventually ending up with the same concentration of the gas everywhere across the grid.

For example a given state progression could be: (using ints for simplicity)

starting state

00000
00000
00900
00000
00000

state after 1 pass of algorithm

00000
01110
01110
01110
00000

one more pas should give a 5x5 grid all containing the value 0.36 (9/25).
I've tried it out on paper but no matter how I try, I cant get my head around the algorithm to do this.

So my question is, how should I set about trying to code this algorithm? I've tried a few things, applying a convolution, trying to take each grid cell in turn and distributing it to its neighbours, but they all end up having undesirable effects, such as ending up eventually with less gas than I originally started with, or all of gas movement being in one direction instead of expanding outwards from the centre. I really can't get my head around it at all and would appreciate any help at all.

Answer 1

In the example you give, your second stage has a core of 1's. Usually diffusion requires a concentration gradient, so most diffusion related techniques won't change the 1 in the middle on the next iteration (nor would they have got to that state after the first one, but it's a bit easier to see once you've got a block of equal values). But as the commenters on your post say, that's not likely to be the cause of a net movement. Reducing the gas may be edge effects, but can also be a question of rounding errors - set the cpu to round half even, and total the gas and apply a correction now and again.

Answer 2

It's either a diffusion problem if you ignore convection or a fluid dynamics/mass transfer problem if you don't. You would start with equations for conservation of mass and momentum for an Eulerian (fixed control volume) viewpoint if you were solving from scratch.

It's a transient problem, so you need to perform an integration to advance the state from time t(n) to t(n+1). You show a grid, but nothing about how you're solving in time. What integration scheme have you tried? Explicit? Implicit? Crank-Nicholson? If you don't know, you're not approaching the problem correctly.

One book that I really liked on this subject was S.W. Patankar's "Numerical Heat Transfer and Fluid Flow". It's a little dated now, but I liked the treatment. It's still good after 29 years, but there might be better texts since I was reading on the subject. I think it's approachable for somebody looking into it for the first time.

Answer 3

It looks like you're trying to implement a finite difference solver for the heat equation with Neumann boundary conditions (insulation at the edges). There's a lot of literature on this kind of thing. The Wikipedia page on finite difference method describes a simple but stable method, but for Dirichlet boundary conditions (constant density at edges). Modifying the handling of the boundary conditions shouldn't be too difficult.

Answer 4

It looks like what you want is something like a smoothing algorithm, often used in programs like Photoshop, or old school demo effects, like this simple Flame Effect.

Whatever algorithm you use, it will probably help you to double buffer your array.

A typical smoothing effect will be something like:

begin loop forever
    For every x and y
    {
        b2[x,y]  = (b1[x,y] + (b1[x+1,y]+b1[x-1,y]+b1[x,y+1]+b1[x,y-1])/8) / 2
    }
    swap b1 and b2
end loop forever

Answer 5

See Tom Forsyth's Game Programming Gems article. Looks like it fulfils your requirements, but if not then it should at least give you some ideas.

Answer 6

The concentration for each iteration given a starting concentration can be obtained by the equation:

concentration = startingConcentration/(2*iter + 1)**2

iter is the time iteration. So for your example.

startingConcentration = 9
iter = 0
concentration = 9/(2*0 + 1)**2 = 9

iter = 1
concentration = 9/(2*1 + 1)**2 = 1

iter = 2
concentration = 9/(2*2 + 1)**2 = 9/25 = .35

you can set the value of the array after each "time step"

Answer 7

Here's a solution in 1D for simplicity:

The initial setup is with a concentration of 9 at the origin (), and 0 at all other positive and negative coordinates.

initial state: 0 0 0 0 (9) 0 0 0 0

The algorithm to find next iteration values is to start at the origin and average current concentrations with adjacent neighbors. The origin value is a boundary case and the average is done considering the origin value, and its two neighbors simultaneously, i.e. average among 3 values. All other values are effectively averaged among 2 values.

after iteration 1: 0 0 0 3 (3) 3 0 0 0

after iteration 2: 0 0 1.5 1.5 (3) 1.5 1.5 0 0

after iteration 3: 0 .75 .75 2 (2) 2 .75 .75 0

after iteration 4: .375 .375 1.375 1.375 (2) 1.375 1.375 .375 .375

You do these iterations in a loop. Outputting the state every n number of iterations. You may introduce a time constant to control how many iterations represent one second of clock-on-the-wall time. This is also a function of what length units the integer coordinates represent. For a given H/W system, you can tune this value empirically. You may also introduce a steady state tolerance value to control when the program says " all neighbor values are within this tolerance" or "no value changed between iterations by more than this tolerance" and so the algorithm has reached a steady state solution.

Answer 8

You could try a search on "cellular automata" ..