#10 Issue with highly compressible flow in APU Testcase

[Lorenz Hammerschmidt]

Hello,

i am trying to validate the mixing plane respository for a highly compressible flowin a radial tubrine. t(ransonic/supersonic regime).
I am using the APU test case, for each geometry data is provided by Sauret ( https://doi.org/10.1115/IMECE2012-88315), validated using CFX, SU2.

Similar to the test case i have set up a total pressure / total tempreature inlet and a static pressure outlet (total pressure & pressureInletOutletVelocity). I have written a small wrapper around the pressureInletOutletVelocity, for the inlet as it is cyclic and i required a non normal but cyclic directed inlet direction. The BC is tested and validated.

I am using:
* ESI 2212 Version
* Mixing Plane Branch for the 2212 Version
* Running on a linux ubunutu 22.04 server with 30 cores

Simulation Setup
* 3 Domains - Stator|Rotor|Diffuser
* Between each domain there is a mixing plane
* MRF Zone with nonrotatingPatches for: Periodics, MixIn, MixOut (unshrouded domain)
* Each domain has cylic bounday conditions
* Inlet: Total Tempreature, Total pressure, cyclicPressureInletOuletVelocity
* Outlet: zeroGradient, Total Pressure, pressureInletOutletVelocity
* Ramping up pressure and omega for the rotating domain
* Using SST k-omega
* Parallisation was done using constrains

Mesh generation was done using TurboGrid and the same mesh was tested in CFX.
Post Processing was done using a self written TurboPost utility
Geoemetry generation and generation of iso surfaces was done with a self written CAD programm, special tailored for turbomachinery (will be made partly public in the future)

Issue:
Simulation runs fine for a few thousands iterations, but then the fields starts to act weird and accelerates further and further leading to a fail of the simulation. I am not sure if it is an issue regarding the mixing plane implementation (as its only an interpolation between two patches) or origins from somewhere else.

I have attached an isospan at 0.5 for the domain at a timestep 8000 (just after ramping is complete) and 13600. The latter shows already a disorted field.

What have i tried?
* Tried to under relax all fields and equations and increase ramping time with no success
* Tried simple to start of and tried pimple loop with LTS without success
* Modified boundary layers, increased mesh size wihtout any success
* Started to implement rothalpy (still in progress) for the energy equation including jump coditions for mixing plane.

Open to any suggestions or help. Looking forward to share any needed information and validation for the mixining plane and turbomachinery implementations in openfoam.

greetings,
Lorenz

[Martin Beaudoin]

Hello Lorenz,
Thank you for testing the mixingPlane with this very interesting case. I am no expert in compressible flows, but I can probably help out with some basic validations with the mixingPlane.

Obviously, having access to your test case would be easier, but there are still some remarks I can share with you based on the information supplied in your ticket report:

• The fact that your case ran for over 8000 timesteps without your simulation blowing up is a first indication that the mixingPlane is probably doing its interpolation job properly. When the mixingPlane is not properly set up, or has a software bug, the mixingPlane simulations tend to crash and burn much sooner than after 8000 iterations.
• Furthermore, we can see from the 2 snapshots you've produced for t=13600 that the instabilities are located rather far from the mixingPlane patches, and very close to the runner wall. When the mixingPlane is at fault, the instabilities usually appears very close to the mixingPlane interfaces. At least from my experience with incompressible flows...
• How large are the flux differences across your mixingPlane and AMI interfaces along the whole simulation? If you could plot the curves of the values reported by the 'mixingPlaneAndGGICheck' function object, using values taken from entries like: "Scaled flux: 0.7852258 0.7852258 Diff = 5.7797622e-10 or 7.3606372e-08 %" for the mixingPlane and "Flux: 6.04931e-08 6.04931e-08 Diff = -5.25853e-17 or 8.69278e-08 %" for the AMI interfaces, that would help to detect mass flow problems at the AMI or mixingPlane interfaces.
• Some comments about your mesh generation: you probably exported your TurboGrid meshes in CGNS format. For CFX, the import converter probably kept the mesh coordinates information in binary format, with full precision. For SU2, the file format being ASCII, how many decimals did the CGNS t to SU2 converter gave you for the mesh vertex coordinates? For OpenFOAM, how many decimals did the CGNS to FOAM gave you for the mesh vertex coordinates? Make sure you keep as much precision as your CGNS converter can give you.
• Do you have the computing budget (memory and time) to run this simulation with a full 360-degree geometry? If so, just run a full geometry case, hence removing the mixingPlane from the simulation and see if your OpenFOAM simulation still behaves properly compared to SU2 and CFX.

[Lorenz Hammerschmidt]

Hello Martin,

i have tried to rerun the whole process again. Some heads on what i have done because i have suspected something with the mesh or the energy equation for quite a while.

Mesh Generation*
* Geometry was generated in my own CAD programm with an export option for TurboGrid
* Generated there a mesh and exported it as CGNS
* Tried two things from there:
* Open Mesh in ICEM -> export as .msh file
* Open Mesh through in Fluent -> export .cas -> open in Fluent Meshing -> export as .msh
* In OpenFoam i used fluent3DMeshToFoam

I have not yet tried to the CGNS->OpenFoam converter.
I have posted about it in CFD Online before: https://www.cfd-online.com/Forums/openfoam-journal/239623-salehi-nilsson-2021-openfoam-francis-turbine-transients.html

i have prepared the case and mesh. I can't share it here but can provide you a link. There included is the mesh, the turbogrid input files, and the case with some results and the log file as well.

Regarding the SU2 case, i have not run it bymyself. I only compared it to literature.

[Martin Beaudoin]

Looking at your various case's solution fields, I would argue that the mixingPlane interfaces are behaving properly.

• The mixingPlane master patches are showing proper axisymmetric field patterns. The shadow patches fields are not axisymmetric because they have not been updated prior to saving on disk. This is a mere design choice; I wanted to display the state of the fields on the mixingPlane patches as the solver did evaluate them after the last iteration prior to saving on disk. I could add a debug switch to force the re-evaluation of the mixingPlane patch fields on write to file; this might be useful for display or report purposes, but that would be merely cosmetic.

• You can visualize the mixingPlane ribbon distribution from VTK files by adding the following dictionnary in your controlDict file:

// ************* //
libs
(
libMixingPlane.so libTurboTools.so
);

InfoSwitches
{
writeMixingPlaneRibbonsProfileToFiles 1;
}

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

The files you will be mostly interested in are:
VTK/mixingPlane_ROTORMIXOUT_DIFFUSERMIXIN_globalCoords_ribbonPatchFeatureEdges.vtk
VTK/mixingPlane_STATORMIXOUT_ROTORMIXIN_globalCoords_ribbonPatchFeatureEdges.vtk

• The patches for the mixingPlane ROTORMIXOUT_DIFFUSERMIXIN are completely conformal, or at least both patches cover the exact same angle span and are perfectly overlapping. I don't think you need a mixingPlane there. A simple AMI interface would do the job, unless of course you absolutely need your fields to be circumferentially averaged at this location.

• The patches for the other mixingPlane STATORMIXOUT_ROTORMIXIN are conformal in the Z direction, so you can observe that the ribbons distribution is following the vertical patches discretization you have enforced. You might want to play with another kind of ribbons discretization (like LINEAR_ALONG_STACK_AXIS), and see if a "coarser" ribbon averaging gives you more stable results. I would doubt that though.

• The flux imbalance across the mixingPlane, as reported by the mixingPlaneCheck functionObject, is mainly a reflection that the fields are not converged (at a mere 1e-3 after more that 15000 iterations).

• The mixingPlane schemes, as defined in system/fvSchemes, will apply "fluxAveraging" to k and omega (and epsilon if present), and select "areaAveraging" for the other fields. I am not experienced enough with your case, selected solver, and with compressible flows in general, but just make sure some of your other fields don't need to use the "fluxAveraging" option as well.

• Your mesh point coordinates are stored in ASCII format, with limited precision. Given the fact that you are chasing a stability problem with your overall solutions that starts developping after many thousands of iterations, I would advise keeping the maximum precision with your mesh point coordinates, hence keeping the double precision binary format all along your mesh generation, conversion and storage processes. This way, you will minimize the digital noise coming from the mesh discretization. Gradients are computed using those basic x,y,z coordinates. You have some very small zones with very large field variations (around the stator blade trailing edge for instance). I am not sure if this will have any impact on your problem, but it won't hurt either, and, as a bonus, you will save some disk space.

• As a debug strategy, you might want to simplify your case by simulating the stator and rotor alone, or the rotor and diffuser alone, or all three components individually, and compare to the results of a similar configuration using CFX or SU2. The idea is to make sure the parameters selected from your fvSchemes and fvSolutions files are correct for your type of simulation. You might also want to explore with other solvers, disable MFR, etc.

• One last comment: Using ParaView, I have been looking at the fields over various slices of your mesh. The most interesting features appear along the rotor blade. For the pressure field, we can almost observe a wave-like distribution of pressure very close to the surface of the blade. I am joining a few snapshots where we can see a somewhat structured distribution of high and low pressure points along the rotor blade. As I said earlier, I am not experienced enough with compressible flows (even less with supersonic flows), so could you comment on the possibility of a real structured pressure wave building up close to the rotor blade under the selected simulation configuration. Or are you convinced this is simply non-physical?

• One more comment: you are using a steady-state solver (rhoSimpleFoam) in a MRF configuration, meaning you are hoping to reach a steady-state solutions while using a rotating flow (ROTOR).

The fact that you are having difficulty converging your fields might be an indication that you are far from a steady-state solution, instead your solution is non-stationary due to the presence of this pressure wave acting up close to the rotor blade.

Hope this helps.

[Martin Beaudoin]

I can see from your log file that you were using rhoPimpleFoam instead of rhoSimpleFoam as specified in your system/controlDict. So you are computing a transient solution after all.

In the paper from Sauret, what kind of solver is "PushButton CFD" ? Steady-state or transient solver?

The paper is also using the Spalart-Allmaras turbulence model instead of kOmega SST. Have you tried using this turbulence model instead?

[Lorenz Hammerschmidt]

Hello, after trying rhoSimpleFoam, but having stability issues (similar results as with the ones presented) i have tried using rhoPimpleFoam as a quasi stationary solver (with localEuler) to stabilise the simulation. Changed the schemes and add PIMPLE Controls. In general it did work well. "Only" the problem around the blade persisted. I haven't tried Spalart-Allmaras turbulence yet, I will give it a try as well. In the paper a stationary solver was chosen. Parallel i am working on the tutorial case provided by foam-extend "axialStatorRotor" to set up different testing cases for the Mixing Plane implementation (in OpenFOAM-ESI). Idea is to use a geometry that can be easily reconstructed using the blockMeshfile. Planned is: Interface Treatment: MixingPlane / PeriodicAMI Different Number blades in each section Different BC (fixed velocity inlet with fixed pressure outlet & total pressure with pressueInletVelocity) greetings, Lorenz Martin Beaudoin mbeaudoin@users.sourceforge.net schrieb am 22. Apr. 2024 um 4:25: I can see from your log file that you were using rhoPimpleFoam instead of rhoSimpleFoam as specified in your system/controlDict. So you are computing a transient solution after all. In the paper from Sauret, what kind of solver is "PushButton CFD" ? Steady-state or transient solver? The paper is also using the Spalart-Allmaras turbulence model instead of kOmega SST. Have you tried using this turbulence model instead? [tickets:#10] Issue with highly compressible flow in APU Testcase Status: open Milestone: 1.0 Labels: MixingPlane ESI 2212 Validation Created: Thu Apr 04, 2024 08:07 AM UTC by Lorenz Hammerschmidt Last Updated: Sun Apr 21, 2024 02:58 PM UTC Owner: Håkan Nilsson Attachments: T (1.8 kB; application/octet-stream) U (2.3 kB; application/octet-stream) boundary (7.3 kB; application/octet-stream) fvSchemes_start (2.6 kB; application/octet-stream) p (2.1 kB; application/octet-stream) t13600_iso05_pressure.png (77.8 kB; image/png) t13600_iso05_velocity.png (95.4 kB; image/png) t8000_iso05_pressure.png (70.2 kB; image/png) t8000_iso05_velocity.png (90.9 kB; image/png) Hello, i am trying to validate the mixing plane respository for a highly compressible flowin a radial tubrine. t(ransonic/supersonic regime). I am using the APU test case, for each geometry data is provided by Sauret ( https://doi.org/10.1115/IMECE2012-88315 ), validated using CFX, SU2. Similar to the test case i have set up a total pressure / total tempreature inlet and a static pressure outlet (total pressure & pressureInletOutletVelocity). I have written a small wrapper around the pressureInletOutletVelocity, for the inlet as it is cyclic and i required a non normal but cyclic directed inlet direction. The BC is tested and validated. I am using: * ESI 2212 Version * Mixing Plane Branch for the 2212 Version * Running on a linux ubunutu 22.04 server with 30 cores Simulation Setup * 3 Domains - Stator|Rotor|Diffuser * Between each domain there is a mixing plane * MRF Zone with nonrotatingPatches for: Periodics, MixIn, MixOut (unshrouded domain) * Each domain has cylic bounday conditions * Inlet: Total Tempreature, Total pressure, cyclicPressureInletOuletVelocity * Outlet: zeroGradient, Total Pressure, pressureInletOutletVelocity * Ramping up pressure and omega for the rotating domain * Using SST k-omega * Parallisation was done using constrains Mesh generation was done using TurboGrid and the same mesh was tested in CFX. Post Processing was done using a self written TurboPost utility Geoemetry generation and generation of iso surfaces was done with a self written CAD programm, special tailored for turbomachinery (will be made partly public in the future) Issue: Simulation runs fine for a few thousands iterations, but then the fields starts to act weird and accelerates further and further leading to a fail of the simulation. I am not sure if it is an issue regarding the mixing plane implementation (as its only an interpolation between two patches) or origins from somewhere else. I have attached an isospan at 0.5 for the domain at a timestep 8000 (just after ramping is complete) and 13600. The latter shows already a disorted field. What have i tried? * Tried to under relax all fields and equations and increase ramping time with no success * Tried simple to start of and tried pimple loop with LTS without success * Modified boundary layers, increased mesh size wihtout any success * Started to implement rothalpy (still in progress) for the energy equation including jump coditions for mixing plane. Open to any suggestions or help. Looking forward to share any needed information and validation for the mixining plane and turbomachinery implementations in openfoam. greetings, Lorenz Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/turbowg/tickets/10/ To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/

[Martin Beaudoin]

I commented too quickly about the mixingPlane fields from the master patches being axisymmetric. I only checked with a couple of the fields. My mistake. The mixingPlane fields are not recomputed prior to being written to disk, so what you are seeing are the values computed and modified from the solver after the circumferential averaging constraint is imposed and the solver iterations completed for a given timestep. And basically, it depends on the solver implementation deciding when to evaluate the the fields on those patches. It is useful to leave the fields unmodified like this because one can evaluate how good or how badly the solver can respect the averaging constraint. But being able to see the averaged fields would be useful too. I'll see how to add this option, probably as a debug option.

[Martin Beaudoin]

Another possible simplification for debugging your simulation would be to replace the stator by a profile1DFixedValue BC with similar T, rho, U, p, k, etc values, and only simulate the flow through the rotor and the diffuser. By using an AMI interface between the rotor and diffuser, you would have a case without any mixingPlane interfaces. See if the instabilities around the rotor blade still persist under that setup.

[Martin Beaudoin]

Hey Lorenz,
Have you made any progress with your test case?
Do we still need to keep this ticket open?

[Lorenz Hammerschmidt]

Hello Martin,

i have made recently progress with my case and just about the validate the last few bits. As suspected, it is/was a problem with the energy equation. Implementing the rothalpy was fairly simple, getting the mixingplane jump running took me a little longer though.
I have switched back to the SIMPLE loop.

Temperature is dropping now nicely along the flow patch, while pressure and velocity fields do look quite good too. Still have some work to do on the wall boundary conditions for rothalpy.
Comparing with the foam-extend library, i think there is a mistake in the UTheta calculation needed for the rothalpy. Refering to the post: https://www.cfd-online.com/Forums/openfoam/210386-unphysical-temperature-mixing-plane.html. So far i have used zeroGradient for rothalpy, which gives wrong values for the rotating walls. Maybe Håkan or Jasak can be of assistant on this topic.

---

**Snapshot from the terminal of the current simulation: **

Time = 12629

velocityDampingConstraint dampU damped 0 (0%) of cells, with max limit 700
DILUPBiCGStab: Solving for Ux, Initial residual = 0.000111808, Final residual = 8.8678e-09, No Iterations 2
DILUPBiCGStab: Solving for Uy, Initial residual = 0.000330974, Final residual = 4.08012e-09, No Iterations 3
DILUPBiCGStab: Solving for Uz, Initial residual = 0.000255981, Final residual = 5.26444e-09, No Iterations 4
limitVelocity limitU Limited 0 (0%) of cells, 0 (0%) of faces, with max limit 800
DILUPBiCGStab: Solving for i, Initial residual = 0.000170392, Final residual = 2.55078e-09, No Iterations 4
limitTemperature limitT Lower limited 0 (0%) of cells, with min limit 230
limitTemperature limitT Upper limited 0 (0%) of cells, with max limit 800
limitTemperature limitT Unlimited Tmin 281.372
limitTemperature limitT Unlimited Tmax 598.382
DILUPBiCGStab: Solving for p, Initial residual = 0.0019947, Final residual = 8.26796e-05, No Iterations 18
DILUPBiCGStab: Solving for p, Initial residual = 0.000619989, Final residual = 7.15821e-05, No Iterations 17
time step continuity errors : sum local = 1.24058, global = -0.328815, cumulative = -565.972
limitVelocity limitU Limited 0 (0%) of cells, 0 (0%) of faces, with max limit 800
DILUPBiCGStab: Solving for omega, Initial residual = 2.81192e-05, Final residual = 6.65136e-09, No Iterations 2
DILUPBiCGStab: Solving for k, Initial residual = 0.0015513, Final residual = 2.07312e-09, No Iterations 4
ExecutionTime = 38616 s ClockTime = 38621 s

Checking flux 'phi' balance across non-conformal boundary patches:
Interface pair (ROTORPER1, ROTORPER2) : Interface type: cyclicAMI
Area: 0.00109392 0.00109392 Diff = 9.72366e-11 or 8.88884e-06 %
Flux: 0.0329633 0.0329633 Diff = 2.98678e-08 or 9.06091e-05 %
Interface pair (STATORPER1, STATORPER2) : Interface type: cyclicAMI
Area: 0.000148424 0.000148424 Diff = -5.6034e-11 or 3.77525e-05 %
Flux: 0.0312491 0.0312491 Diff = 2.88244e-08 or 9.22407e-05 %
Interface pair (DIFFUSERPER1, DIFFUSERPER2) : Interface type: cyclicAMI
Area: 0.00446935 0.00446935 Diff = 9.86112e-10 or 2.20639e-05 %
Flux: 0.0474891 0.047489 Diff = -6.11749e-08 or 0.000128819 %
Interface pair (STATORMIXOUT, ROTORMIXIN) : Interface type: mixingPlane
Area: 0.000127884 0.000151862 Diff = -2.39781e-05 or 15.7894 %
Flux: 0.0208437 0.0248669 Diff = -0.00402321 or 16.179 %
Scaled flux: 0.0208437 0.0209405 Diff = -9.68795e-05 or 0.462641 %
Interface pair (ROTORMIXOUT, DIFFUSERMIXIN) : Interface type: mixingPlane
Area: 0.000237219 0.000237219 Diff = -4.87891e-19 or 2.05671e-13 %
Flux: 0.0249191 0.024869 Diff = 3.11404e-05 or 0.124966 %
Scaled flux: 0.0249191 0.024869 Diff = 3.11404e-05 or 0.124966 %

fieldMinMax fminmax write:
min/max(mag(U)) = 0 582.54
min/max(rho) = 0.414343 3.0333
min/max(k) = 2.65823e-08 5770.42
min/max(p) = 56596.5 414016
min/max(T) = 281.374 598.559
min/max(omega) = 1776.4 1.26477e+09

surfaceFieldValue STATORMIXOUT_avg write:
weightedAverage(STATORMIXOUT) of T = 405.604
weightedAverage(STATORMIXOUT) of i = 183357
weightedAverage(STATORMIXOUT) of p = 219530

surfaceFieldValue ROTORMIXIN_avg write:
weightedAverage(ROTORMIXIN) of T = 405.549
weightedAverage(ROTORMIXIN) of i = 10674.2
weightedAverage(ROTORMIXIN) of p = 219113

surfaceFieldValue ROTORMIXOUT_avg write:
weightedAverage(ROTORMIXOUT) of T = 300.359
weightedAverage(ROTORMIXOUT) of i = 14836.6
weightedAverage(ROTORMIXOUT) of p = 68795.7

surfaceFieldValue DIFFUSERMIXIN_avg write:
weightedAverage(DIFFUSERMIXIN) of T = 304.291
weightedAverage(DIFFUSERMIXIN) of i = 17096.6
weightedAverage(DIFFUSERMIXIN) of p = 68747.6

Interestingly is the weightedAverage(by phi) good at the stator<>rotor interface but has some issues on the interface between rotor <>diffuser. I have used areaAveraging for rothalpy.
Attached some more screenshots from the simulation. As the temperature field shows, it cools down along the flow path. Not sure if there is any issue with the outlet condition, as can be seen in the temperature field, or if the simulation still needs a few more iteration. I am using total pressure/pressureInletOutletVelocity and inletOutlet for the tempreature as a outlet condition.

---

I will keep you posted on the progress. Planned is to upload the modified files to this working group, with a working validation case - ultimately merging with the main ESI branch.

[Martin Beaudoin]

No new comments/feedback since a few months ago. Let's close this.