Algorithm for modeling expanding gases on a two-dimensional grid

Question

I have a simple program, her heart is a two-dimensional array of floats, presumably representing gas concentrations. I tried to come up with a simple algorithm that would simulate a gas expanding outward like a cloud, ending up with the same gas concentration throughout the grid.

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

initial state

00000
00000
00900
00000
00000

after 1 pass of the algorithm

00000
01110
01110
01110
00000

another pas should give a 5x5 grid containing a value of 0.36 (9/25).
I tried this on paper, but no matter how hard I try, I can't get around the algorithm to do this.

So my question is: how should I try to code this algorithm? I tried several things, using convolution, trying to take each grid cell in turn and distribute its neighbors, but they all turn out to have undesirable effects, such as, in the end, with less gas than I originally started, or all the gas movements are in one direction rather than expanding outward from the center. I really can't get around it at all and would appreciate any help.

Answer 1

This is either a diffusion problem if you ignore convection or a fluid dynamics / mass transfer problem if you do not. You would start with equations for the conservation of mass and momentum for the Euler (fixed control value) point of view, if you were to solve from scratch.

This is a temporary problem, so you need to perform the integration to advance the state from time t (n) to t (n + 1). You show the grid, but nothing about how you decide on time. What integration scheme have you tried? Explicit? Implicit? Crank-Nicholson? If you do not know, you do not approach the problem very correctly.

One book that I really liked on this subject was SW Patankar Numerical Heat Transfer and Fluid Flow . This is a little from now on, but I liked the treatment. He is still good after 29 years, but since I read this topic, there may be better texts. I think this is available to someone who is looking at him for the first time.

Answer 2

In the example you give, your second stage has the core of the 1st. Typically, diffusion requires a concentration gradient, so most diffusion-related methods will not change 1 in the middle at the next iteration (and they will not fall into this state after the first, but it is a little easier to see once you get a block of equal values). But, as the commentators at your post say, this is unlikely to cause pure movement. Reducing gas can be an edge effect, but it can also be associated with rounding errors - install the processor by half to half, and also sum the gas and apply the correction again and again.

Answer 3

It seems that you are trying to implement a finite difference solver for the heat equation with Neumann boundary conditions (isolation along the edges). There is a lot of literature on this subject. The Wikipedia page on the finite difference method describes a simple but stable method, but for the Dirichlet boundary conditions (constant density at the edges). Changing the processing of boundary conditions should not be too complicated.

Answer 4

It looks like you want something like a smoothing algorithm, often used in programs such as Photoshop, or old-school demos like this simple Flame Effect .

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

A typical smoothing effect would 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 the Tom Forsyth Game Programming Gems release article . It seems to meet your requirements, but if not, it should at least give you some ideas.

Answer 6

Here's a solution in 1D for simplicity:

Initial setting with a concentration of 9 at the origin () and 0 for all other positive and negative coordinates.

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

The algorithm for searching for the following iterative values ​​should start from the beginning and average current concentrations with adjacent neighbors. The initial value is a boundary case, and the average value is performed taking into account the value of the beginning and its two neighbors at the same time, i.e. The average of 3 values. All other values ​​are effectively averaged between 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 0,375 0,375

You do these iterations in a loop. Status output every n iterations. You can enter a time constant to control how many iterations are one second of the clock time on the wall. It is also a function of which length units are integer coordinates. For a given H / W system, you can set this value empirically. You can also enter the value of the tolerances in the steady state for control, when the program says that “all values ​​of the neighbor are within this tolerance” or “there is no value changed between iterations by more than this tolerance”, and therefore the algorithm has reached a steady state solution.

Answer 7

The concentration for each iteration, taking into account the initial concentration, can be obtained by the formula:

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

iter is an iteration of time. 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"