[abminternet]

Hi everybody,

I am trying to set up a case in OpenFOAM 1.6.x and run it using sonicDyMFoam, but I can't get past the thermophysicalProperties. when I run the solver it gives me the following error:

Selecting thermodynamics package hPsiThermo<pureMixture<constTransport<specieThermo <hConstThermo<perfectGas>>>>>


--> FOAM FATAL ERROR:
Not implemented

From function basicThermo::e()
in file basicThermo/basicThermo.C at line 355.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/software/openfoam/public/OpenFOAM/Git/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/opt/software/openfoam/public/OpenFOAM/Git/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 Foam::basicThermo::e() in "/opt/software/openfoam/public/OpenFOAM/Git/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libbasicThermophysicalModels.so"
#3 main in "/opt/software/openfoam/public/OpenFOAM/Git/OpenFOAM/OpenFOAM-1.6.x/applications/bin/linux64GccDPOpt/sonicDyMFoam"
#4 __libc_start_main in "/lib64/libc.so.6"
#5 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/opt/software/openfoam/public/OpenFOAM/Git/OpenFOAM/OpenFOAM-1.6.x/applications/bin/linux64GccDPOpt/sonicDyMFoam"
Stopped

Any help you can provide will be kindly appreciated.

[challenger]

Hi abminternet,
maybe you can try running the debugger version.

[abminternet]

Hi Challenger,

I really don't think that recompiling OpenFOAM again is the best solution.

[nakul]

Hi,
Have you made any changes to the code?

Basically what the error implies is that the constructor of internal energy(e) is being called somewhere in the code. But the thermophysical model being used in sonicdymfoam never calls this constructor and hence 'e' is never implemented during the execution of the code.

You may be calling this function somewhere but the thermophysical model doesn't allow this.

You may provide some more details regarding the thermophysical setup you are using to get more precise advice.

[abminternet]

Hi Nakul,

I haven't made any changes to the code. In fact, in line 75 there is an "#include eEqn.h", could it be here where e() is being called?

Thanx for the reply,

[nakul]

Hi,
Right now I am not the pc with OF installed. The header file eEqn.H is where most probably the code for handling the 'energy equation' is written.

If you haven't made any change in the code then look where e() is being called. Most probably it would be called in createFields.H. Alternatively you may post your createFields.H here so that I can have a look.

[abminternet]

Hi Nakul,

sure thing, my createFields.H looks like this:

Info<< "Reading thermophysical properties\n" << endl;

autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();

volScalarField& p = thermo.p();
volScalarField& e = thermo.e();
const volScalarField& psi = thermo.psi();

volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh
),
thermo.rho()
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

# include "compressibleCreatePhi.H"


Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);

Info<< "Creating field DpDt\n" << endl;
volScalarField DpDt =
fvc:Dt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);


As you mentioned, e() is being called here, but this is the original code, i checked the code online and it has e() too, does that mean that there is an error in the original code?

[nakul]

Hi,

The code is correct. What I understand is that the "thermodynamics package" that you have chosen is not initialising the variable e. You have chosen a package which does all the calculations based on enthalpy and e is not required.

But due to the code e() is being called and hence the error. So you should decide whether internal energy is required by you or not.

Now there are two ways :
1) Change the thermodynamics package and involve some internal energy based calculation. I think eConstThermo instead of hConstThermo may do the trick.

or

2) You may comment out the line calling e() in createfields.H and recompile the solver. But do make sure that the object 'e' is not being used anywhere.

Alternatively a far lesser complicated way is to refer a tutorial of sonicDymFoam and see what thermodynamics package is being used there. As I am unable to access OF presently so I don't know if a tutorial for sonicDymFoam exists or not. You can however google it or can also get some hints by looking at sonicFoam's tutorials.

Whatever you do please post your findings here,

[SMesser]

I'm having the same issue as abminternet. I'll try playing with the code, but the original reason I went to hPsiThermo + hConstThermo is that ePsiThermo + eConstThermo was crashing on the calculation of e for my particular sim. (Non-convergence errors.... I was hoping the enthalpy could be reasonable even if the internal energy was getting close to a singularity.)

I haven't found any tutorials for sonicDyMFoam. Do you have a link? The tutorials for sonicFoam use eConstThermo, which is why I started with that.

I also tried hPsiThermo + eConstThermo, but that's not a valid pairing.

[SMesser]

So, it looks like e() does get used, but maybe it shouldn't be?

When I run wmake from the sonicDyMFoam directory, output is as follows:

olorin:/opt/openfoam171/applications/solvers/compressible/sonicDyMFoam # wmake
\SOURCE=sonicDyMFoam.C ; g++ -m32 -Dlinux -DWM_DP -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-40 -I../sonicFoam -I/opt/openfoam171/src/thermophysicalModels/basic/lnInclude -I/opt/openfoam171/src/turbulenceModels/compressible/turbulenceModel -I/opt/openfoam171/src/finiteVolume/lnInclude -I/opt/openfoam171/src/dynamicMesh/lnInclude -I/opt/openfoam171/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam171/src/OpenFOAM/lnInclude -I/opt/openfoam171/src/OSspecific/POSIX/lnInclude -fPIC -c $SOURCE -o Make/linuxGccDPOpt/sonicDyMFoam.o
In file included from sonicDyMFoam.C:74:0:
../sonicFoam/eEqn.H: In function ‘int main(int, char**)’:
../sonicFoam/eEqn.H:4:23: error: ‘e’ was not declared in this scope
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:8:10: warning: unused variable ‘momentumPredictor’
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:11:10: warning: unused variable ‘transonic’
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:14:9: warning: unused variable ‘nOuterCorr’
make: *** [Make/linuxGccDPOpt/sonicDyMFoam.o] Error 1

advice?

I'm going to try a couple other thermodynamics models...

[SMesser]

Someone mentioned that simulations tend to act funny if grid cells get aspect ratios > 4. That turns out to be a major problem with my setup. I changed up the mesh to turn down the spatial resolution in some spots, aiming for cells that don't change much in size (<1.2 ratio between widths of neighboring cells) and have small aspect ratios (<4)... sim was able to run a lot longer.

I'll eventually need the higher spatial resolution, but it looks like I'll be playing with overlapping patches and non-rectangular grids rather than tweaking the thermodynamics.

[nakul]

Can you provide some more details about the case and maybe then I might be able to help you with your problem?

[SMesser]

Hiyas, nakul.

Not sure if you were addressing me or abminternet or both...

The following is my blockMeshDict. Decreasing the resolution let me run out to the 80 us end time, which was my initial goal. The "poppet" moves toward negative x at a constant -25 m/s. I'm declaring this working for the moment. My next step is to learn some more advanced-meshing techniques (non-rectangular grids, overlapping patches, and so on) in hopes that I can resolve the changing geometry without going over the aspect-ratio limit I mentioned above. (My previous failing runs had aspect ratios greater than 11.)

convertToMeters 0.01;

vertices        //(x y z) of each vertex
(
    (-3 3 0)
    (-0.2 3 0)
    (-0.1 3 0)  //2; moves with poppet
    (0 3 0)     //3; fixed
    (-3 2 0)
    (-0.2 2 0)
    (-0.1 2 0)  //6; moves with poppet
    (0 2 0)     //7; fixed
    (-0.1 1 0)  //8; moves with poppet
    (0 1 0)     //9; fixed
    (8 1 0)
    (-3 0 0)
    (-0.2 0 0)
    (-0.1 0 0)  //13; moves with poppet
    (0 0 0)     //14; fixed
    (8 0 0)
    (-3 3 1)
    (-0.2 3 1)
    (-0.1 3 1)  //18; moves with poppet
    (0 3 1)     //19; fixed
    (-3 2 1)
    (-0.2 2 1)
    (-0.1 2 1)  //22; moves with poppet
    (0 2 1)     //23; fixed
    (-0.1 1 1)  //24; moves with poppet
    (0 1 1)     //25; fixed
    (8 1 1)
    (-3 0 1)
    (-0.2 0 1)
    (-0.1 0 1)  //29; moves with poppet
    (0 0 1)     //30; fixed
    (8 0 1)
);

blocks          //vertices should be listed here in cyclic order; hex: x1^ is vertex 0->1; x2^ is 1->2; x3^ is 0->4
(
    hex (4 5 1 0 20 21 17 16) (40 20 1) simpleGrading (0.2 1 1)
    hex (5 6 2 1 21 22 18 17) (4 20 1) simpleGrading (1 1 1)   //1; alongside poppet
    hex (6 7 3 2 22 23 19 18) (4 20 1) simpleGrading (1 1 1)
    hex (11 12 5 4 27 28 21 20) (40 60 1) simpleGrading (0.2 1 1)  //3; majority of plenum
    hex (8 9 7 6 24 25 23 22) (4 20 1) simpleGrading (1 1 1)
    hex (13 14 9 8 29 30 25 24) (4 20 1) simpleGrading (1 1 1)
    hex (14 15 10 9 30 31 26 25) (100 20 1) simpleGrading (7 1 1)  //6; capillary
);

edges           
(
);

patches      //when fingers curl in direction of vertices, thumb should point outward   
(
    wall sideWall
    (
        (0 16 17 1)     //0
        (1 17 18 2)
        (2 18 19 3)
        (9 25 26 10)
    )
    wall slamWall
    (
        (0 4 20 16)
        (3 19 23 7)     //1
        (4 11 27 20)
        (7 23 25 9)
    )
    patch orifice
    (
        (10 26 31 15)
    )
    wall poppet
    (
        (6 22 21 5)
        (5 21 28 12)    //2
        (13 29 24 8)
        (8 24 22 6)
    )
    symmetryPlane midplane
    (
        (11 12 28 27)
        (14 30 29 13)
        (15 31 30 14)
    )
    empty frontAndBack 
    (
        (0 1 5 4)
        (1 2 6 5)
        (2 3 7 6)
        (4 5 12 11)     //3
        (6 7 9 8)
        (8 9 14 13)
        (9 10 15 14)
        (16 20 21 17)   //7
        (17 21 22 18)
        (18 22 23 19)
        (20 27 28 21)
        (22 24 25 23)   //10
        (24 29 30 25)
        (25 30 31 26)
    )
);

mergePatchPairs 
(
);

// ************************************************************************* //

Cheers

[nakul]

Hi SMesser,

Right now I am not on my workstation with OF installed. But it would be good if you could provide a pic of your mesh.

Secondly, if you are getting accurate results then why you want to increase skewness in your current mesh? Your mesh is good enough if its giving you results!!

Finally if you actually want to mesh something complex where you can't avoid skewness then let me suggest that blockMesh utility may not be the best option for you. You may try using snappyHexMesh utility, which is more efficient or you can even try commercial meshing softwares like gridgen.They may generate good quality mesh even for very complex geometries.

If you have any problems with case setup do let me know!!

[SMesser]

An image should be attached.

One section of the mesh is only 4 cells across, from one wall to the opposite - and that doesn't seem like enough cells to get good spatial resolution. (In the image, this is the region which is at lowest pressure.) When I ran with 20 cells across this region, the simulation crashed as noted above. The mesh was the only thing that changed, and since the aspect ratio was well in excess of what others have told me is a reasonable upper limit, I'm assuming the aspect ratio is the problem.

I said it runs. I didn't say it was accurate.

The geometry isn't all that complex, but it changes - which means the aspect ratios of grid cells can change substantially. Ideally I'd start the simulation at zero distance between two sections, and run it until the moving part bounces off the far wall and returns to its original position. That means that in a full run, the aspect ratio of the grid cells will hit infinity on three separate occasions. Given that I've had trouble with aspect ratio 11, I suspect infinity may be overly large.

Thus I'm looking for tools which might allow a dynamic _number_ of cells, or which would otherwise do sneaky things with the cell geometry. (Maybe there's a way to have the grid flow around corners? Sort of a hybrid between Eulerian & Lagrangian dynamics?) My understanding of snappyHexMesh is that it deals with spatially complicated boundaries, but isn't really focused on dynamic geometry.

The other challenge I'll be facing is trying to decrease the pressure in the downstream portion of the simulation. Reality has it at vacuum. The pictured simulation has it at half the maximum density. When I've tried to decrease the pressure via setFields, I've gotten crashes - but I don't know if that's due to the choice of solvers or thermodynamic models. I tried using setFields to introduce stepwise changes to minimize the ratio across any one gridline, but my initial attempt there failed.... Do you know of any workarounds for large pressure differences?

Thanks

[deepsterblue]

Sarah,

Do you have changes in mesh connectivity? Or is it just mesh-motion?

[SMesser]

As it stands, it's just mesh-motion, but I need changes in connectivity.

In the image posted above, you can see the mesh is continuous from one end of the device to the other. When the valve is closed, there'll be a break in the middle. (I.e. the grid will not be simply connected as the narrowest channel in the image starts at zero width.)

Since the poppet is never far from the orifice and since its motion will influence gas flow, I don't think just changing face / patch boundary conditions will be sufficient... but I don't know of appropriate examples.

Thanks.

[deepsterblue]

Dynamic layering should work just fine in your case. You can also maintain cell aspect ratios / thickness. You would probably have to approximate opening / closure events using attach/detach sliding-interfaces. Take a look at the movingConeTopo tutorial and the engine classes in 1.6-ext.

[SMesser]

So I downloaded OpenFOAM-1.6-ext_kubuntu-10.04.1.iso from Sourceforge, but it doesn't seem to have a "movingConeTopo" case in it, just a bunch of papers which talk about things various folks have done with OpenFOAM. Only a few examples, though, of how they were done. (The various engine cases seem to be just pictures and text - no code.)

Google gives me several references to movingConeTopo, but none of them seem to include code proper.

Could you provide links to relevant code, so that I know what to do with the dynamicMeshDict, etc.?

Thank you

[SMesser]

Found OpenFOAM movingTopoCase at http://openfoam-extend.git.sourcefor...b77157;hb=HEAD

downloading & investigating now...