0 Folder for a dsmcFoam simulation
 

Hi, I'm new in using the dsmc solver and I#m having the problem with the definition of the 0 folder because I don't know the exact meaning of the properties that you are defininig:

boundaryT (temperature)
boundaryU (velocity)
dsmcRhoN (numberdensity??)
dsmcRhoM (mass density??)
fD (Pressure??)
iDof (???)
internalE
internalKE
momentum
q
rhoM (Why these rhoM and rhoN when you have already defined dsmcRhoN and RhonM?)
rhoN

Is it really necessary to define all of them?

I would appreciate If someone could clarify my doubts and recommend me a tutorial, paper, link, book or whatever that could help me to understand this solver.

I am still reading Bird's book, so I cannot help you a lot.

You find a description of the fields in the documentation of the DsmcCloud class:

http://foam.sourceforge.net/docs/cpp/a00438.html

To clarify your confusion about rhoN and dsmcRhoN, notice that rhoN is the number density (the actual one), while dsmcRhoN is the number density of the DSMC particles.

Best,

You have to define boundaryU and boundaryT according to the corresponding patch Maxwellian distribution properties.

The rest of them may be set to zeroGradient for all patches except symmetryPlanes and empty (it does not make a difference in the simulation).

Thank you all for your answers, definitely it is totally recommendable and form my point of view necessary to read Bird's book and the whole C++ code to understand the operation of this solver

Hi all,

I'm also starting to use the dsmcFoam solver for my OpenFOAM simulations and to get more familiar with the solver, I set up an easy case:

2-D rectangle (1x10 mm) with simple Grading 10x100 cells; inlet and outlet on the short sides and all other sides set as walls. The number density is set to 2.6867774x10^25 for ambient conditions. For an average of 6 dsmcParticles in each cell, the nEquivalentParticles in the dsmcProperties was set to 4.47796x10^12. deltat was set to 10^-7 s to ensure that the gas molecules stayed approximately 3 iteration steps in one cell (assuming a molecule velocity of the speed of sound). boundaryT was set to 300 K and the velocity at the inlet to 50 m/s. All other files in the 0-folder were set to the appropriate values (zeroGradient for patches and calculated with value of 0 for the walls).

The calculated dsmcParticles and their movement as well as the pressure change with time look fine to me but I have had some difficulties understanding the other calculated values.

1. What is the difference between overallT, translationalT and internalT? Are they different temperature values or can they be used to calculate the temperature (units being kg*m^2/s^2 instead of K)?

2. I noticed that the boundary velocity at the inlet seems to be only 25 m/s even though it should be 50 m/s (same thing happens to boundaryT, e.g. when I set it to 700 K, paraView shows only 550 K). Additionally, the calculated velocity is unreasonable where it is 8 m/s within the first cell, then a step to 0 m/s which stays constant within the rectangle until the outlet where it goes back to 8 m/s.

Does anybody see where the problem(s) may be?

Thank you very much for the help

A few things to get you started.

The units you wrote (kg*m^2/s^2) correspond to energy. Maybe you mixed internalT with internalE?

These temperatures are calculated in dsmcFields.C
(http://foam.sourceforge.net/docs/cpp/a07456_source.html)
in lines 168-208. With a little patience, you can see their definition and verify that their units are indeed Kelvin. Keep in mind that many quantities (such as iDofMean, linearKEMean, etc) are naturally defined per unit of volume.

How did you define your inlet and outlet conditions?
In the default version of dsmcFoam you can only give a number density for all open boundary surfaces (unless they worked on that in 2.1, I didn't check that but I doubt it).

See the following:
http://www.cfd-online.com/Forums/ope...am-solver.html
http://www.cfd-online.com/Forums/ope...am-vacuum.html

When you say that you have incorrect velocities and temperatures in ParaFoam, are you visualising boundaryT and boundaryU or overallT and overallU? Because the first ones are not very meaningful, they only serve as boundary conditions and their inner field is not calculated.

The unit of overallT, translationalT and internalT is K, as anon_a said, you can check it in the C++ code.

You can see the definition of these temperatures, for example, in the first chapter (page 14) of Bird's book: Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Oxford Engineering Science Series.

The translationalT is calculated from the velocities of a set of molecules. The internalT is calculated taking in account the vibrational and rotational energy of the molecules (so its value is zero for a monoatomic gas).

The overallT is defined as a weighted mean (taking in account the internal degrees of freedom) of the other two.

Thank you very much for your answers! I’m not so familiar with C++ so I only looked in the calculated files for a specific time step (e.g. internalT). Here, the dimension is [1 2 -2 0 0 0 0], which confused me a bit (expected [0 0 0 1 0 0 0]). Is that always the case with temperatures (dimensions and actually calculated units are not the same) or is there something wrong with my case? Due to your comment I checked the dimensions of internalE, which would be [1 -1 -2 0 0 0 0], thus energy per volume? Seems like I have to learn a bit more about C++…


Do you mean the inlet and outlet conditions in the 0-folder? They are defined as patches in the blockMesh-file with these parameters in the 0-folder: boundaryT=300 K and boundaryU=(0 0 50) m/s and boundaryU=(0 0 0) for the inlet and outlet, respectively. Everything else is set to zeroGradient.


The overallU and overallT seem unreasonable. At the inlet, where there should be a velocity close to the boundary conditions (+50 m/s; shown boundaryU is +25 m/s) a value for UMean_z of -8 m/s is calculated (the opposite for the outlet: a velocity of +8). Within the geometry, the velocity is 0 after a while.


The calculated parameters of the tutorial case “freeSpaceStream” look strange to me, too. boundaryT is set to 700 K but overallT does only reach about 500. And shouldn’t the whole cube reach a constant temperature around 700 K after a while (all other boundaryT conditions are set to walls with zeroGradient).


Thanks a lot and regards

First off, sorry for not checking out the temperature files first but I was not in my PC. Second, to be honest, I am not exactly a C++ master myself too but it is worth investing time on the code.

You are right about the temperature units. I would guess that the units of the Boltzmann constant are not passed through the .value() function (even though it is correctly set up in etc/controlDict). If k was dimensionless, the units would match those given in the temperature files. The same thing would happen with pressure p but in this case k is both in the numerator and denominator and therefore it's ok. I wouldn't worry too much about that anyway.

That is correct.

In order to fully determine the distribution function of particles in a boundary patch (and assuming it is Maxwellian), you need the temperature, velocity and number density. The first two are imposed through the "boundaryT/U" files. The third one is imposed through the constant/dsmcProperties file, in the Inflow Boundary Model section (if you were changing only system/dsmcInitialiseDict, this refers to the initial conditions). However, in the default inflow models you can only impose a specific number density to all open boundaries simultaneously. See the links I sent for some hints on how to change that.

Even if you change the patch inflow properties, you will not get exactly what you prescribe at these boundaries. On a boundary patch for example, the macroscopic quantities are calculated using the molecular properties of particles crossing the patch both ways and you can not be sure (in fact it rarely happens) that molecules moving in opposite directions have the properties you want. Also, their distribution is probably not Maxwellian. It is very common to include a part of the upstream and downstream reservoirs in the computational domain (see e.g. http://www.sciencedirect.com/science...2037961000253X) to solve this problem (but be warned that the computational effort is increased).

Also, the fields may be lost in noise if you have very slow flows and very low Kn.

Sorry but I don't know what you are talking about: In my version of the freeSpaceStream tutorial, the boundaryT field is set to 300K and this is achieved in the overallT temperature. Also, I would highly recommend a small time step to eliminate non-uniformities at the boundaries (if you set a "large" time step, many molecules may enter in "pulses" and then exit with high velocity, causing some unnatural distributions of density and velocity at the boundaries. This is frequently mistaken as a coding error).

I hope these help.

Hey everybody,


First of all, thanks a lot for the explanations, anon_a and Jesgonse. I’ve done some simulations over the weekend and the results look much better now (used a reservoir at the outlet).


Sorry about the freeSpaceStream tutorial. I mixed something up there with the temperatures.


Another case that I am running right now is about heat transfer. The geometry is a simple 2-D rectangle with walls on all sides. boundaryT is 700 K on one side, internalT is 300 K and all other sides have “zeroGradient” as boundaryT. All other parameters in the 0-folder are set to “zeroGradient” and InflowBoundaryModel is “none” due to the absence of patches. This way the whole rectangle should heat up to a constant temperature after some time, right? Unfortunately, the rectangle does not heat up. Only in a very small region close to the heated wall does the overallT reach an increased temperature (about 500 K) but stays at about 300 K everywhere else. For me, it seems like the heat transfer between the wall and the molecules or between the molecules themselves somehow do not work. BinaryCollisionModel is set to “LarsenBorgnakkeVariableHardSphere” and WallInteractionModel to “MaxwellianThermal”. Does anybody have an idea what’s wrong with my case?


Thanks in advance.

You are confusing boundaryT, internalT, internalField and initial temperature.

internalField: This is an entry existing in all quantities. It is the non-boundary value of each quantity, the internal field as the name says. Normally, if you wanted to impose a temperature in the internal field you would change this value in the file T. But dsmcFoam is an exception.

boundaryT: The temperature of the boundary, according to which the boundary Maxwellian is produced. The value given in the "internalField" of this file has no meaning, it is not used anywhere, only the boundary terms make sense.

internalT: This is the part of the temperature associated with internal energy. It is not used for the initial generation of the field.

initial conditions: These are imposed through the system/dsmcInitialiseDict. This is your initial distribution of particles.

Also, there might be a problem in using zeroGradient for boundaryT in DSMC.
Please check that everything is in order and we'll talk again.

Hey anon_a,


I just checked all of the temperature settings and they are as follows:


In the 0-folder, for boundaryT: one wall is set at 700 K, all other wall are at 300 K. The internalField in this file is set to 300 K but as you said, it has no meaning/is not used.


In the constant-folder, for dsmcProperties: Tref of the BinaryCollisionModel is 273 (although I don’t know the meaning of this value)


In the system-folder, for dsmcInitialiseDict: temperature is set as 300K.


Other than these, I didn't find anything else with regard to temperature settings. Did I perhaps miss any, for example where can I find internalT? internalField is also not present in the dsmc solver (except of course in the boundaryT file)?


Is there a way to check if zeroGradient is working properly for boundaryT or can I simply use another condition for my case?


Thank you very much

How small is this affected area? can you post some pictures?
If it is extremely small, I would bet money that you have used specular boundary conditions in constant/dsmcProperties, especially since you started with the freeSpacePeriodic tutorial. Check for this line:
WallInteractionModel SpecularReflection;
and change it to
WallInteractionModel MaxwellianThermal;
Also study a little on the literature for the kinetic boundary conditions.

Regarding Tref, take a look at Appendix A in Bird's book.

Just forget about setting InternalT, you don't need to do that.
InternalT comes as a result, it is not a prerequisite.

I have only tried to use zeroGradient as a BC for boundaryT a long time ago and I remember it gave me an error (didn't it give you an error?). But it also makes sense: when zeroGradient is used, the boundary value is obtained using values from the internal field (which is meaningless for this file). Therefore, it would produce erroneous results.

I just noticed that you had mentioned using Maxwellian BCs (yes, 3 days after my post).
Sorry about that! It's a good thing noone called my bet :-)

Did you have any progress after all?

Hey anon_a,
I did some more simulations with your comments in mind. The cell size is 0.01 mm^2 (see picture) where the inlet, outlet, and sides are defined as walls and boundaryT at the inlet is set to 700 K and dsmcInitialiseT is 300 K.
http://www.cfd-online.com/Forums/mem...7-overallt.png

After applying “MaxwellianThermal”(see 1st picture), the temperature near the inlet gets close to the “boundaryT”, but the rest of the rectangle does not heat up. For me, it looks like the molecules at the inlet receive thermal energy from the boundary (as they were initialized at only 300 K) but somehow do not transfer the thermal energy to the rest of the geometry. This might make sense because as I understand, the molecules do not get reflected from the wall when “MaxwellianThermal” is used since this condition is usually for very small geometries/low pressures. Or the “boundaryU” condition (set to “zeroGradient”) could be wrong? Or even the particle-particle interaction (“LarsenBorgnakkeVariableHardSphere” with “Tref” set at 273 K and “relaxationCollisionNumber” at 5.0) could be wrong (see 2nd picture)? In these centerline “overallT” profiles along the length of the rectangle, one can see that initially there is quite a bit of variation in the temperature within the geometry but after some time it reaches equilibrium around 300 K.




http://www.cfd-online.com/Forums/mem...rallt-plot.png

For comparison, I applied “SpecularReflection”, but there is no interaction with the heated wall at all (remains at 300 K near the inlet). It looks like the molecules are reflected without absorbing any thermal energy.
I would also like to try the “MixedDiffuseSpecular” condition for my case but I don’t know the syntax. It requires a fraction of Maxwell and reflection, otherwise I get an error. Does anybody know how to use “MixedDiffuseSpecular”?
Also trying to get a copy of Bird’s book, but have not been successful as of yet.
Thank you and have a nice weekend

Hmmm... it seems that your field is almost not influenced at all by your wall (observe that only the first cell is practically warm). What is your Kn number? If it is too high, then a large temperature jump is justified on the wall. I would run the same simulation for Kn = 0.1, 1, 10 and observe the differences. You can control Kn by adjusting the number density, but make sure you also keep the same number of particles (by adjusting nEquivalentParticles).

Also, I would use fixedValue BCs in the boundaryT and boundaryU files and then, after being sure about my results, I would try zeroGradient (but, again, I seriously doubt that zeroGradient works the correct way).

The "variations" (a.k.a. fluctuations) you observe are due to the statistical noise and are to be expected in DSMC. There is nothing wrong with the models. In order to correct this, you usually average over a lot of time steps in steady state problems (this is what you currently do in overallT, even without knowing it :-) that's why the field is smoother at step 125 than at step 1).

The results regarding specular reflection are, off course, natural.
MixedDiffuseSpecular is practically a combination of MaxwellianThermal and SpecularReflection and its syntax should be something like:

type MixedDiffuseSpecular;
value uniform 300;
diffuseFraction 0.9;

(I did not confirm this before writing). diffuseFraction represents the percentage of molecules interacting diffusely with the wall, while the rest behave specularly.

By the way, what is the purpose of the study?
Have a nice weekend

Hey anon_a,


Thank you very much for your help. I will run a few simulations with different Kn numbers and boundary conditions. As soon as I have any results, I will get back to you.


The future goal will be to simulate gas flows thru nano-porous layers. I have not worked with the dsmc solver before, so I wanted to start with a quite easy case to get familiar with the parameters.


Regards

Hey,


I have run one case with Kn 0.1 and one with Kn 1, where all other conditions were kept the same, such as wall interaction model (Maxwell), collision model (VariableHardSphere) and number of dsmcParticles per cell (4 in average). boundaryT was set to fixedValue as was suggested (or did you mean all boundary conditions should not be zeroGradient?) but there still seems to be no heat transfer from the wall to the rest of the geometry.


Right now, I’m running two more simulations: one with Kn 10 and one with Kn 0.1 and the MixedDiffuseSpecular wall interaction model. Hopefully, the results will be better this time with the syntax (thanks for your help, anon_a):


WallInteractionModel MixedDiffuseSpecular;
MixedDiffuseSpecularCoeffs
{
value uniform 300;
diffuseFraction 0.9;
}


In the meantime, a new question came up. Is there an easy way to calculate partial pressures, concentrations of a specific gas, or plot moving molecules? A similar question was already posted some time ago (http://www.cfd-online.com/Forums/ope...-dsmcfoam.html) and partial pressures have also been discussed for fluent (http://www.cfd-online.com/Forums/flu...re-oxygen.html). Does anybody know how this works with dsmcFoam? I would like to set a different gas composition at the boundary as compare to within the geometry in order to simulate the diffusion of different gas species within the structure. For me, it sounds like the positions of molecules are stored in each time step at lagrangian/dsmc/positions (fits to number of dsmcParticles) and their type in the same folder in typeld (e.g. 0 for N2 and 1 for O2).


Thank you very much

Could you please post your case here? I am getting very curious about what's going on and can not help you any further otherwise. The mixed boundary condition will not provide anything different since you add specular interaction and this should make things worse.

I have not been really involved with mixtures, so I can not help you with the partial pressures. The particles are indeed stored where you mention. In order to visualize them in with the paraFoam script, activate dsmc in the Mesh parts and your quantity of interest in the Lagrangian fields.

Hey anon_a,


Sorry for my late reply, I was on vacation last week.


I’m trying to simulate a closed box, where one of the walls is heated. After some time, the entire box should reach a constant temperature (the same as the heated wall) due to heat accumulation. The case I used so far was a box (all sides set to “wall”) of 0.5mm x 0.1mm x 1 mm with ambient conditions (300 K and 1 atm) inside while the bottom of the box was set to 700 K. Similar to the images previously posted, the temperature inside the box does not seem to change at all and I don’t know why.


By the way, thanks a lot for your suggestion with paraView. It works fine: set Mesh part to “dsmc” to visualize the dsmc particle and colored them by the Lagrangian Field “typeld” to see their type (e.g. N2 or O2).


I am still working on the partial pressures/number densities. Does anybody know what “iDof” is? In the readme it is mentioned that” iDof” is the number density of extensive quantities, this could be used to calculate the concentrations of a species, right? Unfortunately, “iDof” shows higher numbers than “rhoN”, which makes no sense to me.


Thanks a lot and cheers,

Hello again.

What I meant by "please post your case" is "could you please make a zip with your case directory and attach it here?". I will then take a look and let you know if I find something unusual.

Dear anon_a,


You can download the case I’m working on from the following link:


http://dl.dropbox.com/u/14930252/HeatKn0.1.zip



It is just a simple closed case (all sides are walls) in 2-D geometry with one heated side. As described earlier, the temperature does not increase over time. Hope you can find the problem.


Cheers

I would do the following:
(1) You need more particles per cell for a better accuracy, at least 10. Thus, I would increase the particles to 5000.
(2) Your Knudsen is definitely incorrect. The number of collisions is HUGE for this amount of particles (and that is the reason it is so sloooow). Therefore, I would go back to something simpler, i.e. single gas, before doing this. Which reference length do you use for your Knudsen? I would use the smallest dimension. If the collisions are not so many, you wouldn't even need to run it in parallel.
(3) Your geometry is not in the order of mm. Check the convertToMeters parameter
(4) Is your box really like that or are you trying to make something 2D? Just asking, because the long aspect ratio, checking if you know about "empty" boundary conditions.

You really need to check (2) and (3) with a single gas.

Dear anon_a,


Thank you very much for your suggestions.


The number density I used in dsmcInitialiseDict is 1.836E+25, so the mean free path is 100nm. The smallest distance is 1000nm (thickness of the block), so the Knudsen number should be 0.1. Or is there an error in my reasoning?


I actually want to simulate a 3D case but used a 2D case to start with. Actually, I don’t know what an “empty” boundary condition is for or what it does?


Just started a new simulation with the same geometry but higher number of dsmc particles (45899 overall, thus in average 11.5 per cell) as you suggested, 3D (grading of 2 for the thickness), and only a single gas specimen but still, the temperature seems to stay constant.


Really no idea where the problem is…


Do you have a similar case where the heat transfer works?


Happy Easter

I would recommend to read the manual and apply the tutorials, even if they do not seem directly related with DSMC. The empty BC makes your simulation 2D. You will find this in the early chapters.


Using the hard sphere mean free path definition
lambda = 1 / (sqrt(2) * n * (pi*d^2)) with n=1.836E+25
I get lambda = 7.051E-08 m = 70.5nm, which is close to what you say. Your geometry has indeed 10^(-6)m as the smallest dimension. I just got confused because in
http://www.cfd-online.com/Forums/ope...tml#post352121
you mention that it is in the millimeter (instead of micrometer) range. Furthermore, you have different numbers in that post than in the blockMesh
(0.5x0.1x1 vs. 5x0.1x1). Finally, since you are trying a 2D problem, it does not make sense to use the depth dimension as a reference length since you can adjust it to whatever you want. So, if you chose the 1 micrometer in this case, your Knudsen is around 0.007.


I realized that this number is the one appearing in your dsmcProperties dictionary. You got it wrong: nEquivalentParticles is the number of real molecules that a dsmc particle represents. It is not the expected number of dsmc particles in the simulation. Please make sure you understand the difference and you include this change to your particle number calculations.



I ran a simulation taking into account these remarks for nitrogen at 293.15 and Kn around 0.7 by changing the number density. It was done rather fast, so I apologize if I made a mistake in the process. You can see the results in the attached zip file. I also removed some files to conform to the site attachment size rules but you will see that a temperature distribution is indeed developed near the hot wall. Let me know if there is something you don't understand.

Happy Easter! :-)

moleculeProperties in dsmcFoam
 

Hello everybody,

My question regards the dsmc properties rather than the 0 folder. However, I hope somebody can help me.

In the file dsmcProperties located in the constant folder “moleculeProperties” are defined as:

N2
{
mass 46.5e-27;
diameter 4.17e-10;
internalDegreesOfFreedom 2;
omega 0.74;
}

O2
{
mass 53.12e-27;
diameter 4.07e-10;
internalDegreesOfFreedom 2;
omega 0.77;
}

How does one get to these values? E.g. in the case of N2 the covalent radius is 71 pm, which would lead to a max. diameter of 284 pm=2.84e-10
For a triple bond N atom, the covalent radius is 54 pm resulting in a max. diameter of 2.16e-10 m
As the molecule is not spherical, the other diameter would be even smaller (about half perpendicular to the max. diameter). In any case, I get diameters much smaller than defined in dsmcProperties (similar for O2).

And what is the unit of the mass? It seems like it is not the atomic mass, right (would be something with 28 and 32 for N2 and O2, respectively)?

Cheers and thank you very much,
j3r

Hi there

Check Bird's book, Appendix A, Tables A1 and A2.
Quoting Bird, they are "reference diameters for the VHS molecular model".

The units are simply kg.

Hi,

The reference diameter is calculated using the value of a reference viscosity at a reference temperature. For a variable soft sphere model see the equation 4.62 in Bird's book.
The mass represents the mass of a molecule in kg.

Exactly what I meant but more clearly stated :-)

Thank you very much for your fast answers. That makes sense…In that case, there must be some other problem.

Actually, I’m trying to simulate diffusion on the nano-scale. I started with a very simple case (all at ambient conditions):
Cube with 200 nm side length; one face is defined as wall, the opposite side is the inlet and all other ones are cyclic (this condition is defined for dsmc, right?). Initially, there is 100% N2 in the cube and a 20% O2 80% N2 mixture at the inlet. nEquivalentParticles is set to 1 and I’m using the VHS right now.

After the simulation I compered the concentrations at different time steps with the analytical solution. No matter which parameters I chose (different cell size, other collision model, etc.) the analytical solution is always much faster than the dsmc simulation results. That’s why I thought that the molecule diameter might be too big and thus, too many collisions are calculated which seems to be not true.

Any idea what the problem could be?

Thank’s a lot,
j3r

Some comments and questions:
- If you set up your DSMC case with nEquivalentParticles = 1, then what about using molecular dynamics?
- Time dependent DSMC calculations are supposed to be performed taking a number of realizations with different random seeds and averaging. One realization is usually not accurate enough, unless you take a lot of particles into account. If you don't do that, then you may get the order of magnitude right but some noise will be present. But if your flow is not very violent (as far as I understand, the macroscopic velocity at the inlet is zero and the temperature is also the same as the wall's) then the noise becomes comparable to the solution.
- Have you checked your wall interaction model? Do you use diffuse or specular BCs? (Check WallInteractionModel in dsmcProperties file)
- Could you state your pressure and Knudsen number?

As well I wouldn't simulate diffusion of multispecies gases using the VHS model because several reviews of DSMC state that the model cannot represent properly the diffusion coefficient.

See:
Recent Advances and Current Challenges for DSMC. G.A.Bird
http://www.sciencedirect.com/science...9812219700254X

"The VHS model that reproduces the viscosity coefficient does not, however, lead to the correct
diffusion coefficient in a gas mixture"

anon_a
-So far I only worked with dsmcFoam and never tried another solver. Which one would you suggest for my case, i.e. molecular dynamics?
-I’m using the MaxwellianThermal model right now. However, the wall interaction model should be no problem as the concentration profile looks different from the analytical solution before there are any O2 molecules close to the wall.
-The number density is 2.7e25, thus about 1 atm. The Knudsen number should be about 0.34 applying the distance between the wall and the inlet as representative length.
-Yes, pressure and temperature at the inlet are equal to the initial/internal values and the velocity is set to zero. Right now I simulate a case with bigger dimensions. This way the Knudsen number is decreased, number of molecules is increased and wall interaction model shouldn’t matter. Maybe the simulation results are only different from the analytical solution as they are in the transition regime (will see in a few hours).

Jesgonse:
Do you know if one of the other implemented dsmc collision models is better suited? Or would you also suggest a different OpenFOAM solver for my case?

Thank’s a lot for your help,
j3r

-I have not used molecular dynamics in the past (mdFoam is the name of the solver) but as far as I know, the main advantage of DSMC is to reduce the number of particles. If each particle represents a single molecule (and if MD is not much more computationally intensive for the same number of particles) then I see no point in using DSMC.

-I would try to switch to LarsenBorgnakkeVariableHardSphere to take into account other degrees of freedom. But I hardly believe that this will change your solution significantly, Jesgonse's answer about the VHS inaccuracy is probably much more relevant here than setup.

-Would you also mind posting your case here? Just to make sure that nothing else is wrong.

I agree with anon_a, as long as you are simulating diatomic gases, internal degrees of freedom must be taken into account, and right now only the LarsenBorgnakkeVariableHardSphere model is avaliable. I have never worked with internal degrees of freedom, so I don't really know that model very well

As an idea I would suggest you to try to program the variable soft sphere model and run the simulation with a mixture of monoatomic gases, the diffusion should be properly represented

Hey everybody,

Sorry for my late reply. I tried to verify your suggestions and there were some quite time intensive simulations required.

In the attached file you can find an example case as well as an image with my latest results (unfortunately, the upload limit was too low for e.g. pdf).

Summarized, the diffusion simulated with dsmcFoam is always slower than the analytic solution. I tried varying the length and thus Kn-number as well as wall influence, simulated with VHS and LBVHS model, and increased the diameter to avoid difficulties with the cyclic boundary condition (LBVHS_big in the graphs). When I used Ar and He (applied values from Bird’s book for moleculeProperties as well as an article from Taylor and Cain; J. Chern. Phys. 78(10) 1983 for diffusion coefficient for analytic calculation) instead of O2 and N2 to decrease the degrees of freedom, the results looked similar.

Do you see any mistake in the example case or is there an error in reasoning? I could try to vary the relaxationCollisionNumber for the LBVHS model but that seems to be more like fitting until the results look good and I am not sure which numbers are physically reasonable for my case (not sure if it makes much difference either). Do you have any other suggestions? Do you think the earlier suggested soft sphere model will improve the results or is another solver necessary (e.g. the suggested molecular one)?

Thank you very much for your help,
j3r

Attachment 14387

Uploaded the file once again as some of you had difficulties to open the previous zip-file.

Cheers,
j3r

Implementing the VSS model
 

Hello everybody,

I tried to implement the VSS model as you suggested but have same difficulties. Hopefully one of you can help me. First I changed everything at once, the compilation worked fine but the simulations looked strange (number of dsmcParticles changed and only one gas type was generated at the boundary). Thus, I started all over again and checked the simulation after each step:

1) Copy and paste the VHS folder in openfoam/src/lagrangian/dsmc/submodels and rename the copied folder and containing .H/.C files in VSS. Replace “Hard” in the .H/.C files by “Soft”
2) Add the VSS model in the makeDsmcParcelBinaryCollisionModels.C file in the openfoam/src/lagrangian/dsmc/parcels/derived/dsmcParcel/ folder
3) wmake libso in the openfoam/src/lagrangian/basic/ and openfoam/src/lagrangian/dsmc/ folder as well as wmake all in the openfoam/ folder


So far everything works fine and I get no error messages during compilation or running the simulations. Now I tried to implement alpha, which is required for the VSS model (see e.g. Eq. 2.36 in Bird’s book). This has to be defined for each gas type in dsmcProperties in moleculeProperties (see A3 in Bird’s book).

4) Add alpha in DsmcParcel.H and DsmcParcell.H in the openfoam/src/lagrangian/dsmc/parcels/Templates/DsmcParcel folder (just copy/paste all sections with e.g. omega and rename to alpha).
5) wmake libso in the openfoam/src/lagrangian/basic/ and openfoam/src/lagrangian/dsmc/ folder as well as wmake all openfoam/

Again, I don’t see any error messages during compilation. dsmcInitilizes gives an error massage if alpha is not defined under moleculeProperties, which makes sense. However, the simulation doesn’t work correctly any more even though I still apply the VHS model and alpha is not used anywhere. This is what happens:

When I only initialize one gas type and have two different ones at the boundary, only the one that was initialized is generated at the boundary. The number of molecules decreases during the simulation.

When I initialize two different gas types dsmcFoam gives the following floating point exception error message:

Create mesh for time = 0


Constructing dsmcCloud
Selecting BinaryCollisionModel VariableHardSphere
Selecting WallInteractionModel MaxwellianThermal
Selecting InflowBoundaryModel FreeStream

Constructing constant properties for
N2
O2

Starting time loop

Reading/calculating field rhoNMean

Reading/calculating field rhoMMean

Reading/calculating field dsmcRhoNMean

Reading/calculating field momentumMean

Reading/calculating field linearKEMean

Reading/calculating field internalEMean

Reading/calculating field iDofMean

Reading/calculating field qMean

Reading/calculating field fDMean

fieldAverage: starting averaging at time 0

--> FOAM Warning :
From function polyBoundaryMesh::patchSet(const wordReList&, const bool) const
in file meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C at line 572
Cannot find any patch names matching obstacle
Time = 1e-12

Particles inserted = 27
#0 Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 in "/lib/libc.so.6"
#3 in "/lib/libm.so.6"
#4 Foam::DsmcCloud<Foam::DsmcParcel<Foam::particle> >::equipartitionInternalEnergy(double, double) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libdsmc.so"
#5 Foam::MaxwellianThermal<Foam::DsmcCloud<Foam::Dsmc Parcel<Foam::particle> > >::correct(Foam::DsmcParcel<Foam::particle>&, Foam::wallPolyPatch const&) in "/opt/openfoam201/platforms/linux64GccDPOpt/lib/libdsmc.so"
#6
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
#7
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
#8
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
#9
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
#10
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
#11
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
#12 __libc_start_main in "/lib/libc.so.6"
#13
in "/home/dreyer/OpenFOAM/dreyer-2.0.1/platforms/linux64GccDPOpt/bin/dsmcFoam"
Floating point exception

Does anybody have an idea where the problem is? For me it looks like something is wrong with gas type 1 (whichever one is defined as second type) while type 0 works fine. Maybe I overlooked a file in which adjustments are required when the number of moleculeProperties is increased?

Thank you very much for your help,
j3r

Hello everybody,

I’m still trying to implement the variable soft sphere (VSS) model in dsmcFoam but am stuck with the same problem as described before. Just can’t figure out where the problem is. For me it seems like only changes in DsmcParcel.H and DsmcParcell.H are required or is there any other file I have to modify when I want to add a parameter to molculeProperties? I could imagine that DsmcCloud.C also needs adjustments but my C++ knowledge is too limited.

Any hint would be highly appreciated, I’m really lost here. If you need additional information please let me know.

Thank you so much,
j3r

Hi J3r,

I am right now also thinking of implementing VSS in dsmcFoam. Did you finally managed to get it work?

Cheers,
Wenjie