Strange problem with Direct Simulation Monte Carlo

[Declan]

Hi everyone!

I'm running and modifying some code from the Alejandro Garcia book, "Numerical Methods for Physics". I generalized it to three dimensions, but noticed a problem when calculating the average temperature. Then I went back to the original, untouched program, and found that the problem was already there.

The program essentially uses "representative particles", meaning each particle represents a certain number of particles in roughly the same position and velocity in a gas. It then divides up the volume it's in into cells, for processing collisions and a few other things. The cell widths should be a fraction of the mean free path to "ensure that a cell is locally homogeneous". It goes through a user submitted number of time steps, in which it moves the particles, processes collisions among themselves and the walls, sorts them into their new cells, and takes statistical averages.

Basically, my problem is that the average temperature calculated by the program seems to get higher as there are less cells in the system, and lower as there are more cells. I think this happens because of how he's calculating the average temperature. Every time step a few averages are calculated for each cell. For a cell, it calculates ave_vx, ave_vy, ave_vz, and ave_v2, where the first three are averages of the x, y, and z velocities of the particles in that cell, and ave_v2 is the average of vx^2+vy^2+vz^2 for the particles in that cell. Now, the way it calculates temperature for each cell is by doing ave_v2-((ave_vx)^2+(ave_vy)^2+(ave_vz)^2).

The book mentions that this isn't a proper measurement of temperature: "If you've ever had a rigorous training in equilibrium statistical mechanics, you should instinctively cringe at [this equation]. The proper thermodynamic definition of temperature is based on entropy and involves an average over an ensemble of states. All the same, it is useful to extend definitions of thermodynamic quantities, such as temperature, to nonequilibrium systems by using equilibrium identities, such as the equipartition theorem."

Is this the reason why? Can anyone help me?

Thanks!

 

[Alej Garcia]

Hello,

It's a tricky problem and, unfortunately, the answer is a little complicated. It's all spelled out in this paper:

"Estimating Hydrodynamic Quantities in the Presence of Microscopic Fluctuations", A. Garcia, Communications in Applied Mathematics and Computational Science 1 53-78 (2006).

You can find a copy on my website:
http://www.algarcia.org/Pubs/index.html

Cheers,
Alejandro Garcia

 

[Declan]

Wow, are you really the author?? If so, awesome!

Anyway, I'll read that paper right now. The conclusion we came to (before I read what your paper says) was that with too few particles, the center of mass velocity is too powerful. The way I have it set up, I have the particles simply in a 3D cube, so the CoM velocity should be zero. I tried with varying numbers of particles and it seems that the CoM velocity rounds off to zero as it should be when there are a lot of particles per cell.

This also seemed to be reflected in G.A. Bird's Fortran code we have, in which when finding the temperature the way you did in the dsmcne.cpp program, he takes the average of the particles in groups of cells, rather than one cell.

But that was just our guess, and I'll actually read the right answer in your paper now. Thanks!

 

[Alej Garcia]

Yeah, it's me. The brief answer is that you can't define an instantaneous temperature that has the right mean value. But if you measure the mean values of the conserved quantities (mass, momentum, energy) then from those mean values you get the right temperature. The latter is what Bird does in his Fortran code. All that said, with enough particles you can get an instantaneous temperature that's not too bad. The funny thing is that the same phenomenon occurs for the fluid velocity. Anyway, let me know if you have questions after reading the paper.

 

[Declan]

Ok, so I read the paper. I understood the first three pages pretty well, but started getting lost after that. It's very clearly written; there are just some big gaps in my knowledge of this stuff.

For example, I don't understand equation (8):

Is that a Taylor expansion? I understand what <H(M)>_S is, I just don't understand how it's equal to that.

I also don't understand the modification you do here to improve the definition of temperature based on the instantaneous internal energy per particle:

I see that you factored the rho out, but why do you subtract m/V from it? Is that to correct for the samples with N=0 or N=1 that we don't use?

Thanks for your time!

 

[Declan]

Hmmm. I went through it more thoroughly and I think I understood most of the math this time, but I tried implementing what I got from the paper. So I calculated temperature at the end based one the time averages of the mechanical variables I got. But even with what seems to be a decent number of particles (100,000 particles for 1,000 cells, doing 10,000 timesteps), the average temperature is still falling low.

For large numbers of samples, this method of finding the average temperature in a cell should be unbiased, even out of equilibrium, right?

I'm also wondering something else now. In the original code of dsmcne.cpp, the temperature seems to be calculated as the time average of the instantaneous temperature of a cell, which is biased out of equilibrium. Why did you do it this way?

 

[Declan]

Durr, I made a very stupid mistake with the velocity distribution for particles colliding off walls, when I was generalizing to 3d...

Everything is dandy at equilibrium. But now I set up a gas flow in one direction, and the temperature is dropping. I'm using the method from the paper that is good for nonequilibrium.

 

[Alej Garcia]

Hello again,
Sorry for not getting back to you sooner but I've been traveling most of the summer. Any luck getting things working?
Alej Garcia

 

[Declan]

Hello again,
Sorry for not getting back to you sooner but I've been traveling most of the summer. Any luck getting things working?
Alej Garcia

Haha, no problem. I wish I could travel too.

Lots of luck! This particular problem was due to a few things, all stupid errors on my part. Since then, I have made it so the gas can flow over any surface that is formed from arbitrary rectangles (which, in the limit as the rectangles get very small, can be most continuous surfaces).

I'm actually curious as to how other people have done this. Most papers I've tried to find on the subject of 3d DSMC have talked about more physics related (and I guess more important) subjects. But how have others had the gas deal with surfaces which are in 3 dimensional space?

I'm also working on getting it to work with multiple processors using OpenMP right now.