[EnricoDeFilippi]

Hi, I'm simulating a stage of an axial turbine. I'm making a comparison between CFX and OpenFoam and the velocity, pressure and density fieds are the same, while the temperature field is unuual in the rotor, thus affecting the Mach calculation. As it's an wxpansion, the temperature should decrease and so it does in the stator, but at the mixing plane it has an unusual increase and then it decrease again but obviously it is not the correct field. Maybe it's due to a energy non conservation across the mixing plane or the MRF zone introduce some error. I've had the same results with steadyCompressibleMRFFoam and steadyUniversalMRFFoam. I'm simulating a compressible real gas, but also the same case with ideal gas shows the same behaviour, even with different mesh. Did anyone already face this kind of problem or do you have an idea about this mistake and how to solve it?
Thanks in advance.

[calim_cfd]

Hi.
Just guessing here..

Are the boundary conditions properly set for the T file/variable in the 0 directory?

[EnricoDeFilippi]

Yes, they are the same as in the tutorial https://github.com/Unofficial-Extend...ineMixingPlane , so I guess they are correct.

[calim_cfd]

Hi.
That was my first guess. Second guess would be for you to check the other boundary conditions of every other variable. Check the reference pressure, the transport and thermophysical properties. Also, check for wrong boundary condition assignments. You have many patches, one may be wrong across the input files.. idk.. But, IF every other result matches the CFX one, this is really weird. You might wanna check the type of solver both codes are running and coupling physics among the momentum-pressure-density.
Also, always check every tutorial you run. They are written by experts, but they are still humans, sometimes it's a matter of an uploaded file with an older version and everything just fails.


It's hard to pin point these issues in these case w/o carefully inspecting it.

Keep digging and you'll surely find the nasty bug!
gl!

[EnricoDeFilippi]

Thank you very much for your suggestions. All the boundary conditions seem fine for me. Probably the error is inside the thermophysical properties since I'm simulating a real gas and real gas library is still under development in the extended version, so maybe ther's somethng tricky inside. I hope I can find a solution as soon as possible, even if I'm investigating it for a while.
Thanks for the encouragment!

[HenrikJohansson]

Hi Enrico,

I recently updated my version of foam-extend 4.0 to the latest version and as well installed 4.1.
With the old 4.0 the mixing plane works great and my solution converges.
But with the updated version almost all fields are incorrect by the mixing plane just like yours.
Looks like the might have made some changes in the code in general that affects the mixing plane
because I can not find any changes in the mp when checking the history on sourceforge.

When it comes to 4.1 my solution crashes. It seems like there were alot of changes within the
steadyUniversalMRFFoam solver. My p field takes forever to solve and crashes after a few itterations.

So in conclusion I think I will stay with the old version of 4.0 until an official release of 4.1 is done.

I could share that version with you if that might solve your problem.

[ncasari]

Hi Enrico and Henrik,
I know this thread has not been updated for a while but I would like to let you know that I run into the same issue of yours. By analyzing the problem, I found that it is related to the non-conservation of the total enthalpy due to an error in the calculation of the rothalpy jump (there is also a ticket open at https://sourceforge.net/p/openfoam-e...ndrelease/316/). The fixing of the problem involves to add UTheta in the createFieldsDict of the solver you are using (e.g. steadyCompressibleMRFFoam) and the updating of its value in the iEqn.H. To evaluate UTheta, I created a new function in MRFZone. The function is then called from within the solver.

The solver is validated on an HPT stage of a compressible axial turbine (EEE by NASA-GE).
Feel free of writing to me if you need further details.

Best regards,
Nicola

[.bastian]

Ciao Nicola,

could you please describe your evaluation of UTheta more detailed?

[ncasari]

Dear Bastian,

Basically you calculate the tangential unity vector by calculating the cross product of the rotational speed and the distance between the cell center and the axis. Then you can easily project the velocity vector in that direction. Something like:

   
    forAll (cells, i)
    {
          label celli = cells[i];
          const scalar magUTheta = mag(rotVel ^ (C[celli] - origin));
          const vector eTheta = (rotVel ^ (C[celli] - origin))/magUTheta;
          UTheta[celli] = (U[celli] & eTheta)*eTheta;
     
    }

Regards,
Nicola

[.bastian]

Ciao Nicola,



thank you for your fast reply.
I will look if i get it implemented.

[ncasari]

As a follow up of some private messages I received, I am posting down here the solution I implemented. Please, note that this might not be the best way to proceed, but the solution seems to work so I am quite confident it might be helpful.

The issue I found is that the rothalpy jump is not evaluated since the field UTheta is not found by the mixingPlaneEnthalpyJumpFvPatchFields.C boundary condition.
The fix requires to work on three different files: createFields.H, iEqn.H and MRFZone.C.

The first thing needed is to add UTheta as a field in createFields.H. This can be done for example with

volVectorField UTheta
    (
            "UTheta",
            U
    );

after U has been defined.

At this point, one needs to actually calculate UTheta. I did this in MRFZone.C, including the following function:

void Foam::MRFZone::UTheta(volVectorField& UTheta ) const
{

    const volVectorField& C = mesh_.C();

    const vector& origin = origin_.value();
    const vector rotVel = Omega();

    const labelList& cells = mesh_.cellZones()[cellZoneID_];

    forAll (cells, i)
    {
          label celli = cells[i];
          const scalar magUTheta = mag(rotVel ^ (C[celli] - origin));
          const vector eTheta = (rotVel ^ (C[celli] - origin))/magUTheta;
          UTheta[celli] = (UTheta[celli] & eTheta)*eTheta;

    }

   // Included faces
    forAll (includedFaces_, patchi)
    {
        forAll (includedFaces_[patchi], i)
        {
            label patchFaceI = includedFaces_[patchi][i];
            UTheta.boundaryField()[patchi][patchFaceI] = vector::zero;
        }
    }
    // Excluded faces
    forAll (excludedFaces_, patchi)
    {
        forAll (excludedFaces_[patchi], i)
        {
            label patchFaceI = excludedFaces_[patchi][i];

            const scalar magUTheta = mag(rotVel ^
                  (C.boundaryField()[patchi][patchFaceI] - origin));
            const vector eTheta = (rotVel ^
                  (C.boundaryField()[patchi][patchFaceI] - origin))/magUTheta;

            UTheta.boundaryField()[patchi][patchFaceI] -= (UTheta.boundaryField()[patchi][patchFaceI] &
                 (rotVel ^ eTheta)) *eTheta;

        }
    }
}

(keep in mind of declaring it into the relative .H file).
This projects the U velocity in the theta direction. The last step required is the update of the UTheta value for each iteration. The modification proposed is to include

 

    Info << "Calculate velocity projection in tangential direction" << endl;
    UTheta == U;
    mrfZones.UTheta(UTheta) ;

After

    // Create rotational velocity (= omega x r)
    URot == U - Urel;

In this way UTheta is updated with the value. You should not see anymore the warning

Velocity fields URot or UTheta not found.

I hope this can be of help.

Regards,
Nicola

P.S. I was said that there are at least two versions of fe41 that are circulating, and the difference lies in the UTheta type in

src/thermophysicalModels/basic/derivedFvPatchFields/mixingPlaneEnthalpyJump/mixingPlaneEnthalpyJumpFvPatchFields.C

In one version it is declared as volScalarField and in the other as volVectorField. The description so far applies in the case of volVectorField. If this is not the case, what I would do, to change only one file, is to create two variables in the createFields Dict. One will be called UThetaV (volVectorField) and the other UTheta (volScalarField). This way you can pass UThetaV to mrfZones.UTheta and then UTheta is simply mag UThetaV.

[Duran_Martin]

Thank you Nicola for taking the time to detail and share your solution!

I'll be implementing this and hope it solves my problems.

Kind regards,
Duran Martin

[WhiteW]

Thanks to Nicola we also identified a solution for the case when UTheta is defined as volScalarField in:
src/thermophysicalModels/basic/derivedFvPatchFields/mixingPlaneEnthalpyJump/mixingPlaneEnthalpyJumpFvPatchFields.C

Furthermore, running the tutorials of the steadyCompressibleMRFFoam, for Fe4.1, I also encountered the following warning:

From function void gradientEnthalpyFvPatchScalarField::updateCoeffs(c onst vectorField& Up)
in file derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C at line 135
Velocity fields U or URot or UTheta not found. Performing enthalpy value update for field i and patch 0

From function void gradientEnthalpyFvPatchScalarField::updateCoeffs(c onst vectorField& Up)
in file derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C at line 141
Velocity fields U or URot or UTheta not found. Performing enthalpy value update for field i and patch 1

This solution also calculates URot as a copy of Urot to remove the previous warning.


1) The first thing needed is to add 3 fields in createFields.H:
- tangential velocity vector UThetaV
- tangential velocity scalar UTheta
- rotational velocity vector URot
This can be done for example with:

 volVectorField UThetaV
    (
            "UThetaV",
            U
    );

    volScalarField UTheta
    (
        "UTheta",
        mag(U)
    );

    volVectorField URot
    (
            "URot",
            U
    );

after U has been defined.

2) At this point, one needs to actually calculate UTheta. Include the following function in MRFZone.C:

void Foam::MRFZone::UTheta(volVectorField& UTheta ) const
{

    const volVectorField& C = mesh_.C();

    const vector& origin = origin_.value();
    const vector rotVel = Omega();

    const labelList& cells = mesh_.cellZones()[cellZoneID_];

    forAll (cells, i)
    {
          label celli = cells[i];
          const scalar magUTheta = mag(rotVel ^ (C[celli] - origin));
          const vector eTheta = (rotVel ^ (C[celli] - origin))/magUTheta;
          UTheta[celli] = (UTheta[celli] & eTheta)*eTheta;

    }

   // Included faces
    forAll (includedFaces_, patchi)
    {
        forAll (includedFaces_[patchi], i)
        {
            label patchFaceI = includedFaces_[patchi][i];
            UTheta.boundaryField()[patchi][patchFaceI] = vector::zero;
        }
    }
    // Excluded faces
    forAll (excludedFaces_, patchi)
    {
        forAll (excludedFaces_[patchi], i)
        {
            label patchFaceI = excludedFaces_[patchi][i];

            const scalar magUTheta = mag(rotVel ^
                  (C.boundaryField()[patchi][patchFaceI] - origin));
            const vector eTheta = (rotVel ^
                  (C.boundaryField()[patchi][patchFaceI] - origin))/magUTheta;

            UTheta.boundaryField()[patchi][patchFaceI] -= (UTheta.boundaryField()[patchi][patchFaceI] &
                 (rotVel ^ eTheta)) *eTheta;

        }
    }
 }

and in MRFZone.H

void UTheta(volVectorField& UTheta ) const;

This projects the U velocity in the theta direction.

3) The function has also to be present in MRFZones.C, the file that manages the simultaneus presence of more MRF Zones:

void Foam::MRFZones::UTheta(volVectorField& U) const
{
forAll (*this, i)
{
operator[](i).UTheta(U);
}
}

and in MRFZones.H:

//- Compute UTheta (needed for mixing plane)
  void UTheta(volVectorField& U) const;

4) Compile the MRF Library using the Allwmake in:

~/foam/foam-extend-4.1/src

5) The last step required is the update of the UTheta value for each iteration. The modification of iEqn.H includes:

 UThetaV == U;
mrfZones.UTheta(UThetaV);
UTheta = mag(UThetaV);
 URot == Urot;

after: "// Create rotational velocity (= omega x r) Urot == U - Urel;"

6) Compile the Solvers using the Allwmake in:

~/foam/foam-extend-4.1/applications