[mturcios777]

I was browsing through the old extend/dev code and found there used to be a lot more code in the $FOAM_SRC/engine folder than there is in the current 2.1.x version. The engineTopoChangerMesh/layerAR looks like it could be really useful to do engine simulation so I decided to try and port/adapt it to 2.1.x as an engineMesh subclass.

I've uploaded my first attempt which requires a few changes to the files in engineMesh to compile and link properly. In the Make/files, add

engineMesh/layeredEngineMeshAR.C

Additionally, we need a pure virtual bool function in engineMesh.H:

[CODE]
virtual bool deform() = 0;
[CODE]

And the function deform should be defined in all the other the subclasses (layeredEngineMesh, fvMotionSolverEngineMesh, etc). A simple one to use is

bool deform()
{
    this->move();
    return false;
}

Anyways, I compile, link and after adding the following entries to the engineGeometry dictionary:

engineMesh    layeredAR;
.
.
.
deformAngle 10;
delta 1e-4;
offSet 1e-4;
piston
{
        patch   piston;
        coordinateSystem
        {
                origin (0 0 0);
                coordinateRotation
                {
                        type    axes;
                        e1      (1 0 0);
                        e2      (0 1 0);
                        e3      (0 0 1);
                }
        }
        minLayer        1e-3;
        maxLayer        1e-2;
}

(the actual values are not so important, as I still need to test and see what the parameters do). Everything goes well until its time to do the layer addition/subtraction, when I get the following:

Courant Number mean: 0.000371093475 max: 0.00648838644
deltaT = 1.56621499e-06
Crank angle = 342.02585 CA-deg
deltaZ = 8.26588501e-06
pistonLayerID: 0
Piston layering mode
Cannot find point in pts1 matching point 0 coord:(0.000187666667 6.83333333e-05 0.000473727705) in pts0 when using tolerance 5.07608117e-08
Searching started from:0 in pts1
    Compared coord:(0.000187689683 6.82700891e-05 0.000473727705) with difference to point 6.73022297e-08
.
.
.
cyclicPolyPatch::order : Writing neighbour faces to OBJ file "/windows/D/dEFRectMeshLayerTest/side_half0_faces.obj"
cyclicPolyPatch::order : Writing my faces to OBJ file "/windows/D/dEFRectMeshLayerTest/side_half1_faces.obj"
cyclicPolyPatch::order : Dumping currently found cyclic match as lines between corresponding face centres to file "/windows/D/dEFRectMeshLayerTest/side_half1_faceCentres.obj"
--> FOAM Serious Error : 
    From function cyclicPolyPatch::order(const primitivePatch&, labelList&, labelList&) const
    in file meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C at line 1383
    Patch:side_half1 : Cannot match vectors to faces on both sides of patch
    Perhaps your faces do not match? The obj files written contain the current match.
    Continuing with incorrect face ordering from now on!
--> FOAM FATAL ERROR: 
Not implemented

    From function cloud::autoMap(const mapPolyMesh&)
    in file fields/cloud/cloud.C at line 65.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-2.1.x/src/OSspecific/POSIX/printStack.C:201
#1  Foam::error::abort() at ~/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/error.C:249
#2  Foam::Ostream& Foam::operator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at ~/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/errorManip.H:85
#3  Foam::cloud::autoMap(Foam::mapPolyMesh const&) at ~/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/fields/cloud/cloud.C:66
#4  Foam::mapClouds(Foam::objectRegistry const&, Foam::mapPolyMesh const&) at ~/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/mapClouds.H:52
#5  Foam::fvMesh::mapFields(Foam::mapPolyMesh const&) at ~/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/fvMesh/fvMesh.C:489
#6  Foam::fvMesh::updateMesh(Foam::mapPolyMesh const&) at ~/OpenFOAM/OpenFOAM-2.1.x/src/finiteVolume/fvMesh/fvMesh.C:674
#7  Foam::polyTopoChanger::changeMesh(bool, bool, bool, bool) at ~/OpenFOAM/OpenFOAM-2.1.x/src/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChanger.C:293
#8  Foam::layeredEngineMeshAR::deform() at ~/OpenFOAM/OpenFOAM-2.1.x/src/engine/engineMesh/layeredEngineMeshAR/layeredEngineMeshAR.C:257
#9  
 at ~/OpenFOAM/OpenFOAM-2.1.x/applications/solvers/combustion/myEngineFoam/myEngineFoam.C:93
#10  __libc_start_main in "/lib/libc.so.6"
#11  
 at /usr/src/packages/BUILD/glibc-2.11.3/csu/../sysdeps/i386/elf/start.S:122
Aborted

I've attached a picture of a periodic face and the corresponding dumped *.obj files. It looks like the solver dies after a cell is removed. Could someone familiar with these classes give me a hand? Again, I think this would be VERY useful for lots of people.

[mturcios777]

Alright, seeing if I can get some more help with this particular problem

My code crashes like this:

--> FOAM FATAL ERROR: 
object is not allocated

    From function Foam::autoPtr<T>::operator->()
    in file /home/mturcios/OpenFOAM/OpenFOAM-2.1.x/src/OpenFOAM/lnInclude/autoPtrI.H at line 153.

FOAM aborting

Due to code that looks like this:

autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(false);
...
       if (topoChangeMap->hasMotionPoints())
...

Clues?

[deepsterblue]

There's a reason why more of the engine code is in OF-ext. There's a fair amount of variation in the base-level code (the primitiveMesh/polyMesh/fvMesh level) between different flavors of OpenFOAM, particularly when topo-changes are involved.

Any particular reason why you want to stick to OF-2.1.x?

[mturcios777]

We are doing a lot of work with sprays and are impressed by the new Lagrangian spray classes, as well as a lot of the improvements to thermo/chemistry libraries. I'm going to attempt to compile the current version of extend to see how the engine mesh works and if it is worth it to continue this line of thought.

[mturcios777]

I've found what you mean about the different versions of topology changes. It appears the old version of polyTopoChanger's function changeMesh is more like the current polyTopoChange member function changeMesh, with the mapping included.

I was able to determine this by taking a step back and compiled the old 1.6-ext version under 2.1.x (compiling only those parts relevant to the layerAR/layerARGambit subclasses). engineTopoChangerMesh is implemented with dynamicMesh and polyTopoChanger. Wouldn't it make more sense to use topoChangerFvMesh for this (topoChangerFvMesh exists in both 1.6-ext and 2.1.x, though I haven't looked at implementation and it would probably be much different)?

[kmpang]

Hi Marco,

Did you manage to compile and run the solver with engineMesh of layeredAR successfully?

I'm interested in this too.

Please advise.

Thanks.
Pang

[mturcios777]

Yes, have a look at the following bug report:

http://www.openfoam.org/mantisbt/view.php?id=653

With this fix, layerAR works as expected (assuming you deal with updating the phi, meshPhi, the turbulence and thermo models).

[kmpang]

Many thanks for the link.

[kalle]

Marco, do you know if this works in parallel? If so, do you have constraints in that the domain decomposition has to be "columns", such that no processor attempts to remove cells from other processors?

Regards,
Kalle

[mturcios777]

Hi Kalle,

I've never tried it in parallel before, so I can't say for sure. The layerAdditionRemoval works by checking the thickness of the layer of cells directly above the faceZone pistonLayerFaces. As long as the layer remains the same across all the processors it should work. I would argue against trying to use this modifier in parallel and an motionSolver; better to explicitly move the piston boundary.

Good luck, and if you obtain good results please let me know!

[kalle]

Thanks for you input. Some pure mesh motion testing suggest that you do need to have domain decomposition as "columns" (like engineScotch decomposer in 1.6-ext does), plus all processors need a share of the layerAR modifier, like back in 1.5-dev. (i.e. you cannot have some processor only calculating flow outside the cylinder, which is a serious restriction if in-outflow is considered)

I would really prefer to have non-topological mesh motion for these reasons.

K

[mturcios777]

It makes sense that each processor would need the modifier, but since the modification happens on the faceZone you may be able to get away with it. Thanks for the update!

[lucasal]

Hi Marco,

I am currently running a simulation with sprayEngineFoam by using a simple box as engine mesh with layered addition removal as described by you here.

At the moment, I am simulating only compression/expansion, starting from -180°CA (BDC) to +180°CA (next BDC).

What I noticed is that temperature and pressure are not specular respect to TDC, e.g. at -50°CA during compression temperature(640K) and pressure are substantially higher than at +50°CA during expansion (T=285K!!).
From theory they should be the same, since neither combustion nor fuel injection occur.

The thermophysicalProperties used are:

thermoType
{
type hePsiThermo;
mixture reactingMixture;
transport sutherland;
thermo janaf;
energy sensibleEnthalpy;
equationOfState perfectGas;
specie specie;
}

And the fvSolution file is the following:

solvers
{
rho
{
solver PCG;
preconditioner DIC;
tolerance 1e-07;
relTol 0.1;
}

rhoFinal
{
$rho;
tolerance 1e-07;
relTol 0;
}

"(U|k|epsilon)"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-07;
relTol 0.1;
}

p
{
solver GAMG;
tolerance 1e-07;
relTol 0.1;
smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;

}

pFinal
{
$p;
tolerance 1e-07;
relTol 0;
}

"(U|k|epsilon)Final"
{
$U;
tolerance 1e-07;
relTol 0;
}

"(Yi|O2|N2|H2O)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0;
}

h
{
$Yi;
relTol 0.01;
}
hFinal
{
$Yi;
}
}

PIMPLE
{
transonic no;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
momentumPredictor yes;
}

relaxationFactors
{
fields
{
".*Final" 1;
}
equations
{
".*Final" 1;
}
}


Have you noticed these P and T behaviours with layered A/R and engineMesh motion?
I was wondering if there could be an issue in the meshPhi flux that is not correctly computed expecially when new cells are added during expansion.

Thank you very much for any advice.

Regards,

Luca

[mturcios777]

I haven't run into any issues before, after making some adjustments to the pEqn. When you say you are using sprayEngineFoam, how are you adjusting for the flux on new faces? How have you modified the solved to account for topology changes?

[lucasal]

Thanks for the reply. I didn't make any change to the solver yet.

Due to a mismatch of units, I noticed that there was a rho missing in the pEqn.H when computing the surfaceScalarField phiHbyA, so I added it in this way:

(fvc::interpolate(rho*HbyA) & mesh.Sf())

Which other adjustments to the pEqn.H need to be done?

Have you also tried some simulations in parallel?

In my case, the simulation runs well until fuel injection. When there are fuel droplets in the layer that is going to be removed, the simulation crashes without telling anything relevant about the cause.
This does not happen in serial, and the lagrangiang cloud is handled correctly with layered A/R.

Regards,

Luca

[mturcios777]

Have a look at rhoPimpleDyMFoam on the official release and sonicTurbDyMEngineFoam from the extend project. You need to adjust the flux when the topology changes (its in a file called compressibleCorrectPhi.H or correctPhi.H).

[lucasal]

Hi Marco,

Thank you for the advice. In the last week I had a look at the files you mentioned, and I built in the correctPhi.H from rhoPimpleDyMFoam into sprayEngineFoam, but the correction of pressure had a weird behaviour.

So, yesterday I found the thread "0001195: dynamicRefinement update misses meshPhi" made by you, where you were using sprayDyMFoam with the compressibleCorrectPhi.H. In my opinion, this file could suit my case too, because layer A/R or mesh refinement are both topological changes. Am I right?

Therefore, I built in this file in my simulation and the results seem much better.

I would like to ask you something about this change.
Why are you updating the variable "phi" in this way? And why is the "pcorrEqn" identic to the "pEqn" but with the variable "pcorr" instead?

tmp<fvVectorMatrix> UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
);

UEqn().relax();

Info << "Creating phi" << endl;

volScalarField rAU(1.0/UEqn().A());

phi = fvc::interpolate(rho)
*((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho,U));

while(pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::ddt(psi, pcorr)
+ fvc::div(phi)
- fvm::laplacian(rho*rAU, pcorr)
);

pcorrEqn.solve();

if (pimple.finalNonOrthogonalIter())
{
phi += pcorrEqn.flux();
}
}

If you compare with the file correctPhi.H from rhoPimpleDyMFoam, "phi" is not updated before "pcorrEqn" and the "pcorrEqn" itself is a bit different. Moreover, the "phi" at the end coming out of the "pcorrEqn.flux" has a - instead of a + as in your case.

dimensionedScalar rAUf("rAUf", dimTime, 1.0);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divrhoU
);

Info<< "I am in Pcorr.H" << endl;

pcorrEqn.solve();

if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();

I would be very grateful if you could explain me a bit more these differences.
Thank you very much.

Regards,

Luca

[mturcios777]

When the official version of OpenFOAM came out with rhoPimpleDyMFoam, I compared it to the correctCompressiblePhi I was using from the ext version, and noticed that the ddt(psi,pcorr) term made a big difference with the results. In general when I compared the official version with the ext version, I found the modified ext-version gave better results overall. One thing I haven't tried in a while is to use the official rhoPimpleDyMFoam on a layerAR case and see if things have gotten any better. If you end up doing that, could you report back on how things went? One thing to be aware of is that the latest updates in 2.3.x cause any toplogy modifier that adds boundary faces (layerAR, attachDetach) incompatible with RAS models, as there is a change in the way near wall distance is calculated.

[lucasal]

Hi Marco,

What do you mean with the last sentence? In my simulations I use RAS compressible with kEpsilon and layerAR and I didn't experience any crash when layers are added. Are k and epsilon not computed correctly due to the change in near wall distance calculations?

I had a look at the results with and without the compressibleCorrectPhi.H taken from sprayDyMFoam. I ran an engine simulation with sprayEngineFoam.
As you can see from the pictures during the compression phase, this correction causes a non-physical T distribution at the bottom layer of the combustion chamber (where layers are removed), where T is higher than internal (~600K @ 60°CA bTDC e.g.) and piston boundary field (433K).
This does not happen without the correction.

This is the compressibleCorrectPhi.H I am using.

if (meshChanged)
{

Info << "Starting pressure correction" << endl;

U.correctBoundaryConditions();

U.oldTime().correctBoundaryConditions();

    forAll(U.boundaryField(), patchi)
    {
        if (U.boundaryField()[patchi].fixesValue())
        {
            U.boundaryField()[patchi].initEvaluate();
        }
    }

    forAll(U.boundaryField(), patchi)
    {
        if (U.boundaryField()[patchi].fixesValue())
        {
            U.boundaryField()[patchi].evaluate();

            phi.boundaryField()[patchi] =
                rho.boundaryField()[patchi]
               *(
                   U.boundaryField()[patchi]
                 & mesh.Sf().boundaryField()[patchi]
               );
        }
    }

    forAll(p.boundaryField(), patchI)
    {
        if (p.boundaryField()[patchI].fixesValue())
        {
            pcorrTypes[patchI] = fixedValueFvPatchScalarField::typeName;
        }
    }
}


{
    Info << "pcorr time" << endl;

    volScalarField pcorr
    (
        IOobject
        (
            "pcorr",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("pcorr", p.dimensions(), 0.0),
        pcorrTypes
    );

  pcorr == p;
  pcorr.oldTime() == p.oldTime();

  tmp<fvVectorMatrix> UEqn
  (
     fvm::ddt(rho, U)
     + fvm::div(phi, U)
     + turbulence->divDevRhoReff(U)
  );
  
  UEqn().relax();

  Info << "Creating phi" << endl;

  volScalarField rAU(1.0/UEqn().A());

  phi = fvc::interpolate(rho)
    *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho,U)); 

  while(pimple.correctNonOrthogonal())
  {
      fvScalarMatrix pcorrEqn
      (
       fvm::ddt(psi, pcorr)
       + fvc::div(phi) 
       - fvm::laplacian(rho*rAU, pcorr)
      );
      
      pcorrEqn.solve();
      
      if (pimple.finalNonOrthogonalIter())
      {
	phi += pcorrEqn.flux(); 
      }
  }
}

The central part of the sprayEngineFoam.C is the following:

    while (runTime.run())
    {
        #include "readEngineTimeControls.H"
        #include "compressibleCourantNo.H"
        #include "setDeltaT.H"

        phi = mesh.Sf() & rhoUf;

        runTime++;

        Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
        
        bool meshChanged=mesh.deform();

        if (meshChanged)
        {
            Info << "Mesh topology has changed!" << endl;
	    if (correctThermo)
               {
                   Info << "Correcting thermodynamics for changed mesh" << endl;	
                   thermo.correct();
                   rho = thermo.rho();
                   rho.oldTime() = p.oldTime()*psi.oldTime();

                   forAll(Y, i)
                   {
                       Y[i].correctBoundaryConditions();
                   }
                   combustion->correct();
               }
            if (correctPhi)
               {
                   Info << "Correcting convective flux for changed mesh" << endl;

                   #include "compressibleCorrectPhi2.H"

                   if (correctTurb)
                    {
                       Info << "Correcting turbulence model" << endl;
                       turbulence->correct();
                    }
               }
    
        }
        // Make the fluxes relative to the mesh-motion
        fvc::makeRelative(phi, rho, U);

        if (meshChanged && checkMeshCourantNo)
        {
            #include "meshCourantNo.H"
        }

        parcels.evolve();

        #include "rhoEqn.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "YEqn.H"
            #include "EEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

        #include "logSummary.H"

        rho = thermo.rho();

        if (runTime.write())
        {
            combustion->dQ()().write();
        }

Do you have any suggestions where I could investigate a bit more, or why the correction behaves like that?

Thank you very much,

Regards,

Luca

[mturcios777]

The temperature seems very high indeed. I don't remember seeing anything that high on layer removal. What does your pEqn.H look like? I remember having issues because there was an issue with the dpdt term for a bit. Might have something to do with it...