
[Sponsors] 
October 28, 2018, 05:19 
Unphysical temperature with mixing plane

#1 
New Member
Enrico De Filippi
Join Date: Jul 2018
Location: Brescia, Italy
Posts: 14
Rep Power: 6 
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. 

October 29, 2018, 14:05 

#2 
Senior Member
mauricio
Join Date: Jun 2011
Posts: 167
Rep Power: 16 
Hi.
Just guessing here.. Are the boundary conditions properly set for the T file/variable in the 0 directory?
__________________
Best Regards /calim "Elune will grant us the strength" 

October 29, 2018, 15:14 

#3 
New Member
Enrico De Filippi
Join Date: Jul 2018
Location: Brescia, Italy
Posts: 14
Rep Power: 6 
Yes, they are the same as in the tutorial https://github.com/UnofficialExtend...ineMixingPlane , so I guess they are correct.


October 30, 2018, 05:33 

#4 
Senior Member
mauricio
Join Date: Jun 2011
Posts: 167
Rep Power: 16 
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 momentumpressuredensity. 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!
__________________
Best Regards /calim "Elune will grant us the strength" 

October 30, 2018, 11:53 

#5 
New Member
Enrico De Filippi
Join Date: Jul 2018
Location: Brescia, Italy
Posts: 14
Rep Power: 6 
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! 

February 15, 2019, 06:31 

#6 
Member
Henrik Johansson
Join Date: Oct 2017
Location: Gothenburg
Posts: 38
Rep Power: 7 
Hi Enrico,
I recently updated my version of foamextend 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.
__________________
/ Henrik Johansson 

March 31, 2020, 04:42 

#7 
New Member
Nicola Casari
Join Date: Mar 2020
Posts: 3
Rep Power: 5 
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 nonconservation 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/openfoame...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 NASAGE). Feel free of writing to me if you need further details. Best regards, Nicola 

April 8, 2020, 11:34 

#8 
New Member
Join Date: Mar 2017
Posts: 25
Rep Power: 8 
Ciao Nicola,
could you please describe your evaluation of UTheta more detailed? Last edited by .bastian; April 9, 2020 at 05:39. 

April 9, 2020, 03:12 

#9 
New Member
Nicola Casari
Join Date: Mar 2020
Posts: 3
Rep Power: 5 
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: Code:
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; } Nicola 

April 9, 2020, 07:55 

#10 
New Member
Join Date: Mar 2017
Posts: 25
Rep Power: 8 
Ciao Nicola,
thank you for your fast reply. I will look if i get it implemented. 

July 9, 2020, 05:33 
More details on the implementation

#11 
New Member
Nicola Casari
Join Date: Mar 2020
Posts: 3
Rep Power: 5 
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 Code:
volVectorField UTheta ( "UTheta", U ); At this point, one needs to actually calculate UTheta. I did this in MRFZone.C, including the following function: Code:
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; } } } 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 Code:
Info << "Calculate velocity projection in tangential direction" << endl; UTheta == U; mrfZones.UTheta(UTheta) ; Code:
// Create rotational velocity (= omega x r) URot == U  Urel; Code:
Velocity fields URot or UTheta not found. 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 Code:
src/thermophysicalModels/basic/derivedFvPatchFields/mixingPlaneEnthalpyJump/mixingPlaneEnthalpyJumpFvPatchFields.C Last edited by ncasari; July 28, 2020 at 04:14. 

July 9, 2020, 06:05 

#12 
New Member
Duran Martin
Join Date: May 2018
Posts: 1
Rep Power: 0 
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 

July 29, 2020, 03:55 

#13 
Member
Join Date: Dec 2015
Posts: 74
Rep Power: 9 
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: Code:
volVectorField UThetaV ( "UThetaV", U ); volScalarField UTheta ( "UTheta", mag(U) ); volVectorField URot ( "URot", U ); 2) At this point, one needs to actually calculate UTheta. Include the following function in MRFZone.C: Code:
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; } } } Code:
void UTheta(volVectorField& UTheta ) const; 3) The function has also to be present in MRFZones.C, the file that manages the simultaneus presence of more MRF Zones: Code:
void Foam::MRFZones::UTheta(volVectorField& U) const { forAll (*this, i) { operator[](i).UTheta(U); } } Code:
// Compute UTheta (needed for mixing plane) void UTheta(volVectorField& U) const; Code:
~/foam/foamextend4.1/src Code:
UThetaV == U; mrfZones.UTheta(UThetaV); UTheta = mag(UThetaV); URot == Urot; 6) Compile the Solvers using the Allwmake in: Code:
~/foam/foamextend4.1/applications 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
whats the cause of error?  immortality  OpenFOAM Running, Solving & CFD  13  March 24, 2021 07:15 
mixing plane causes fatal error  amin.z  FLUENT  0  March 8, 2015 04:26 
mixing plane model: too much backflow and does not converge  sbhi  FLUENT  1  December 4, 2014 00:31 
mixing plane vs frozen rotor  feafan  STARCCM+  5  August 25, 2014 06:14 
MFR and mixing plane boundary .  Mica  FLUENT  0  January 26, 2007 09:39 