Hi Maneshi,
the only thing i can say is that i found it unphysical for a flow not to converge to a steady state condition. This is because each system has a time constant Tau, and a solution can be considered as converged after around 3 *Tau or 4*Tau approximately.
I think the problem could be connected to the impossibility to simulate a real discharge condition at the moment with this solver, as you probably have seen in the previous posts by Archy and me.
It could be also something due to some errors in setting up the boundary conditions for your case.

I once did a simulation with a tube having a bubble inside, where i set a cyclic condition at the inlet and outlet of my domain. I had the same result of yours, meaning that i did not reach a steady state condition, but this was because the patches in which i had the cyclic boundary were not isobaric (i.e they were normal to the gravity vector). This way, when the bubble exited the outlet patch and came back to the inlet one, it was subjected to a jump in pressure, which prevented the solution from reaching a real steady state.

Could this be interesting for you?

I'm sorry, but I do not have any other ideas at the moment...

Marta

The case you mentioned is quite similar to my problem, but here for a nano channel with infinite length the patches of inlet and outlet are isobaric. aren't they ?
And the other interesting thing is that velocity distributions which are supposed to be of a quadratic form are not like that, they are not even symmetric! the velocity jumps and density distribution are not even similar as those anticipated! where is the problem?
I thought maybe the problem is about not reaching a steady state within 5000 time steps with 1e-15 second time step, but no improvements can be seen even after 50000 time steps elapsed!!!
Any ideas what point I'm missing?

Thanks

Mehdi

Dear Maneshi, i would like to help you more, but i think i will need to see your test case to make a couple of runs and see what's up...
I will be happy to help you analysing your test case if you can send me something. You can also use my email account, if you want.

Best Regards

Marta

pressure in dsmcFOAM
 

Hi Chris - I am trying to run a dsmcFOAM and I noticed in your post that you were at one point using pressure based boundary conditions. I was wonder how you implemented these types of conditions. Did you simply use the native fD boundary condition or is there some other way to specify a pressure. Thanks for any and all help.

Dear Hyperion, as far as I know there are no explicit pressure conditions to be applied in dsmcFoam. You can extrapolate pressure by calculating k*rhoN*T, where k is the Boltzmann constant.

Marta

Hi Marta - Thanks for the response. I am trying to produce a low density jet into a vacuum, and when I try specify temperature and number density at the boundaries I get some odd results. Have you seen any simulations of jets into a low pressure background that I can look at?

Also, do you know a way to calculate the number density on a per cell basis (#/cell) as a quality metric.

Also, do you happen to know how the mean fields are calculated? Are they just time averaged over the entire simulation?

Also, I just assumed that the wall interaction models automatically apply to patches defined as walls. Is that the case.

Thanks for all the help.

Hi Hyperion! I will start answering a couple of your questions, for the others i need to have a quick look at the source files, so i will answer after having checked.

Let's start with the test case you are trying to run. I have something similar to simulate: a flow inside a given geometry with an outlet at a very low pressure, almost in vacuum.
At the moment i don't know how to impose boundary conditions in order to simulate such a flow correctly. I'm studying the inflowBoundaryModel, but it will take me another week, i think (more or less).

Concerning how the averages are calculated, refer to the controlDict file. There you see the 'functions' keyword. The average fields are calculated in the fieldAverage.C file (look at the private member functions).
They are time averages between the start of the simulation and the time when you save your data, as far as i understood.

Regarding the average value of rhoN in each cell of the domain, i'm checking a couple of possibilities.

I will write back soon!

Marta

Hi Hyperion, this is to answer about the wallInteractionModel.
To find the source code you can have a look in:

src/lagrangian/dsmc/submodels/wallInteractionModel/wallInteractionModel

In the WallInteractionModel.H file contained in that folder, you will find this virtual function:

virtual void correct
(
const wallPolyPatch& wpp,
const label faceId,
vector& U,
scalar& Ei,
label typeId
) = 0;

I think that since a pointer is passed from wallPolyPatch, this particular models for particles behaviour near walls can be applied only to patches defined as 'wall'.

Moreover, if you do not specify which particular wall model you want to apply to your system in the 'dsmcProperties' file (contained in the 'constant' folder), then no wall model is automatically applied and the solver will complain that no wall model is specified or, if you eliminate the keyword completely, it will say this:

keyword WallInteractionModel is undefined in dictionary "/home/marta/OpenFOAMRelease/marta-1.6/run/tutorials/discreteMethods/dsmcFoam/Portata8Fori/constant/dsmcProperties"

Marta

Hi Hyperion. Concerning the rhoN calculated based on cells, i think this is already done in the dsmc code.
I made this check for you.

1) i looked at my geometry and saw the number of cells in the domain (internalField). The resulting parameters are:

Mesh stats
points: 166789
internal points: 137980
faces: 1845299
internal faces: 1787685
cells: 908246
boundary patches: 11
point zones: 0
face zones: 1
cell zones: 1

Overall number of cells of each type:
hexahedra: 0
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 908246
polyhedra: 0

So, as you can see, the cells are 908246.

2) At this point i controlled a rhoN file of a generic time step different from 0. What it is possible to see, is:

FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0.0014";
object rhoN;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 -3 0 0 0 0 0];

internalField nonuniform List<scalar>
908246
(
....
....
...
)

So, apparently from the dimensions of this parameter and considering the fact that you have as many resulting rhoN values as the cells in your domain, i think it already does what you expect.

Hope this helps!

Marta

Hi Marta - Thanks for all of your excellent answers. I too will continue to try to get a better understanding of how to simulate internal flows more readily. I have has experience with another 2D dsmc program that used chemistry models at the boundaries that enabled me to specify what fraction of each species would be absorbed or reflected and other reactions at the wall. I used this to force every particle of a specific species that encountered the wall to exit the computational domain. This worked rather well, but now I need to simulate 3D flows. Hence, I am now trying to use dsmcFOAM. I noticed that another openFOAM module uses similar (maybe the exact same) chemistry files, but I don't know how to incorporate them into dsmcFOAM.

Once again, thanks for all your help.

Hi Marta - I just saw your answer regarding the rhoN internal field. What I'm trying to calculated is the number of particles per cell, so if I had the cell a similar list as the one you directed me to but with corresponding cell volumes I could determine how many particles were in each cell. I want to do this because some of the cells in my mesh are extremely small (~10E-12 m^3). Thus, it is possible that in these cells there might only be 1 or 2 dsmc particles. I want to make sure all cells have enough particles to produce statistically accurate results.

Once again, thanks.

Ok...now i understand your point, i'll do a verification and let you know. =)

Marta

I think you can create your own functionObject. To do this, you can copy for example forceCoeffs.C and forceCoeffs.H and see how they work, or you can use any of the files used to create the dynamic libraries, that you can find in the controlDict file of some of the tutorials in OpenFOAM.

The piece of code doing this should look like:


functions
{
forces
{
type forceCoeffs;
functionObjectLibs ( "libforces.so" );
outputControl timeStep;
outputInterval 1;
patches
(
intra
extra
);
pName p;
UName U;
log true;
rhoInf 1;
CofR ( 0 0 0 );
liftDir ( -0.239733 0.970839 0 );
dragDir ( 0.970839 0.239733 0 );
pitchAxis ( 0 0 1 );
magUInf 618.022;
lRef 1;
Aref 1;
}
}

You can create a folder where you should put a Make directory, a lnInclude directory and the directory containing your .C and .H files. Since the linked libraries are part of a list, you should link them this way, if you have more than one library:

functionObjectLibs ( "libforces.so lib2.so lib3.so" );

You should create your own private member function, calculating the particles in each cell and then modify the two original files accordingly.

There are more than one way to access cell volumes, i know the following ones:

1) const dimensionedField<scalar,volMesh> & V() const;
2) const scalarField & cellVolumes() const;
3) mesh.V() [cellID]

but you rhoN is a volScalarField, so i don't know if they work properly, anyway you can try doing a cast, if this is possible.

Marta

Dear All, after trying different boundary conditions, i have improved results a little bit, but still they are not acceptable.
For example, here you see a sort of mass accumulation at the end of the domain...have you ever observed something similar when trying to simulate a flow inside a geometry with inlet - outlet conditions?

Thanks

Marta

Hi Marta - I have obtained similar results. I modified the supersonicCorner tutorial by making a domain of the same size and including a small( r= 5mm) cylindrical nozzle at one end. The outlet side consists of all walls except the surface opposite to the nozzle. Check out the attached pdf for a more compete synopsis. I also ran at higher inlet and outlet velocity (1936 m/s) and obtained no accumulation at the outlet.

Hi Hyperion, thank you for your PDF. I tried your BC combination, and it worked because no mass accumulation was visible. But i have to simulate something different, with no free stream conditions.

I think that in the case of your simulation, the density gets lower because the velocity is maintained constant between the inlet and outlet (free stream), so it's just a matter of continuity.

In my case, i do not know which is the discharge velocity, it should be calculated from the code. Maybe i shall just modify the inflowBoundaryModel...

Thank you very much again

Marta

Hi Marta - I have had some success with the vacuum boundary by setting the outlet normal boundaryU to extremely high values (~2000 m/s) while setting the inlet boundary to something more realistic (400 m/s). The molecules end up leaving the computational domain at around 700 m/s, which is pretty close to the adiabatic speed limit. I still haven't found a good way to enforce a specific pressure at the outlet, but for my application that's OK.

I did have a question about an earlier post, however. You mentioned that you can access cell volume using mesh.V, and cellVolumes(). How do access these? It would be ideal to be to be able to get something that writes the cell volumes to a file.

Thanks again for all the help.

Hi Hyperion,

Marta and I have talked about this lately and the situation is the following:

For dsmcFoam only boundaryU and boundaryT are real input files, all the others are merely output fields ( yes also rhoN and rhoM ).

how are the number densities defined then ?
they are defined by the inflowBoundaryModel. for now there are two models NoInflow and FreeStream, where freestream sees all non-wall patches as openings with the same numberdensity defined in dsmcProperties under FreeStreamCoeffs.

So an outlet as we wish to use is not implemented yet ( but i think you could define a negative speed at your outlet if you know this parameter ).

I wanted to use a pressure condition so I tried to create a new inflow model "inOutflow" which is basically a copy of FreeStream but gives you the possibility to define a numberDensity for each opening patch.
Therefore you will have to put the patch name ( look it up in the boundary file ) in the inOutflowCoeffs dict like

inOutflowCoeffs
{
inlet_patch
{
numberDensities
{
Ar 1e20;
O2 1.21e20;
}
}

outlet_patch
{
numberDensities
{
Ar 1e5;
O2 1.21e5;
}
}
}


I am not completely sure if its doing the job but you can have a try if you want.

The compressed folder is actually tar.bz2 so rename the extension first.
wclean and "wmake libso" it
which will create a new libdsmc.so in your user_libbin folder.
this lib has the same name as the orginal library in the openfoam directory but user libraries are ranked first, so the one in your user_libbin will be used.

Let me know if it works with your case.

Cheers,
Archy

Hi Archy - First, I want to thank you and Marta for all the good work in trying to implement pressure boundaries. When I specify a large velocity at the outlet it seems to have the effect of a vacuum pump. I think this is because the freestream model uses boundaryU to calculate number flux on patches of type "patch." According to the code it does this according to equation 4.22 of Bird's book. Like you said the number flux only depends on boundaryT and boundaryU, so it seems that people have been varying one or the other to obtain maximum outflow. I know that in some simulations that I've run this resulted in an accumulation of particles at the outflow patches, but only when I specified a number density at those patches. When I leave the freestream outflow patches with the rhoN as calculated or zeroGradient then the number density decays radially and axially from my nozzle with no accumulation at the boundaries.

One annoyance that I have encountered as a result of the high velocity outlet is that the code won't output dsmc fields (e.g. Umean, overallT) due to rhoN being too small somewhere in the domain. Have you experienced this?

I am thinking I might modify the dsmcFields.C to remove the requirement of large rhoN. I also want to add some additional fields such as mean free path and number of particles per each cell to ensure quality of results.

Thanks again for providing the changes to get a good vacuum boundary. I'll try to implement it soon.

Hi Hyperion, i have experienced this too. In my case the reason was not due to the high boundary velocity itself, but to a wrong initialisation of the velocity inside the domain (in other words it was too high with respect to the boundaries i specified).

In general this happens when a cell is detected with a very small density in the meandensity field during the evaluation of the intensive variable fields (U, T etc)by the dsmcFields function object. Are you sure that there are no divisions by zero?
You can increase the averaging period of the case (using longer write interval) and/or make
sure that you have the fieldAverage functionObject set to:

type fieldAverage;
...
resetOnOutput off;

so that the averages accumulate rather than being reset at each write to disk.
Once all cells have a suitably high density, this warning will disappear and the
intensive fields will be calculated and output.

Anyway, by using a more appropriate internal initial velocity, the warning disappeared in my case.

PS: I will post here my code with the possibility to initialise particles with different densities in different areas as sonn as i will have it finished!