GGI in OpenFOAM
 

Dear All,

I am happy to announce that I have now finished the implementation of the General Grid Interface (GGI) in OpenFOAM, including massive parallelisation and algebraic multigrid solver (this was a struggle). The SVN snapshot, including the updated tutorial will be uploaded later today.

Please note that as a part of parallelisation I have changed the domain decomposition tools AND the definition of the GGI - you will have to (slightly) update the existing cases.

Since GGI is pretty popular and quite tricky to set up, I would like to propose a skype teleconference with a demonstration, some background, tutorial and a series of slides on how to use the feature. Please drop me a line if you are interested and we at Wikki will do the organisation.

Many thanks to people who contributed to the code - you know who you are. I also expect to present the details of numerics, implementation and parallelisation on the Workshop in Montreal

Enjoy,

Hrv

Hrv,

that sounds great to me. Thank you very much. I will give it a try soon.

Regards BastiL

Thanks for all the efforts, and what is your skype account?
I am insterested in.
Anyone else?

Hi Hrv,

This sounds good to me. Please put me on the list for the proposed teleconf.

Thanks,
Helmut

Hello Hrv,

Thanks for this great work. I like to know more about Ggi, and how to work with that. Would you add my name to the list for teleconf?

Thank you
Mojab

Hrv,

I have run first tests now and everything works fine so far. Thank you very much.

I would like to join you skype session, but have no webcam and never used it before. Is it possible to just listen and watch? Would be nice :-)

Fabian

I'd like to be in too

Hi Hrv

Thanks for the huge amount of effort that this all must have taken. Would it be possible to record and post the video conference as a tutorial file in OpenFOAM? Also I am unable to find ggi on the subversion repository - What is it called and which package is it in?

Thanks again

Nick

Hello,

The GGI is now part of the main OpenFOAM-1.5-dev libraries, and is not provided as a separate package.

Martin

OpenFOAM GGI webcast
 

Dear All,

Thank you for your responses, both through the Forum and directly. We have picked Thursday 14/May/2009 14.00 GMT for the webcast, which will be done as a teleconference + pdf slides.

To register, please E-mail enquiries@wikki.co.uk with your name, affiliation and contact details.

Hrv

Thanks for your great work!
 

Hi,

I would like to thank you for this very useful development.

I'm sure that Martin Beaudoin 'knows who he is', but I would like to acknowledge him in name, as one of the guys who did much of the work.

Håkan.

missed OpenFOAM GGI webcast
 

Hi Hrv

Is there any available documentation about the GGI implementation, as video recordings or slides? Unfortunately I missed the webcast!

Thanks again

Daniele

Thats really crazy, Daniele... Tell me how to live in future... ;)

Marico,

You are right, and fortunately I am living in the present with some look to the future :). I just have not read the month correctly.
In any case Hrv please add me into your list

Daniele

GGI user report
 

First, I'd like to applaud Martin Beaudoin and Hrvoje Jasak for their amazing efforts with the GGI. Truly inspiring work!

It's now time to expose the GGI to the real world and begin to accumulate the experiences. Here is my brief report which includes both success stories and problems. (The story will continue in Montreal, but hopefully without the problems section.)

Case:
A single-blade pump in a volute. Incompressible, turbulent flow.

Mesh:
Two separately generated GridPro meshes, one for the impeller (rotor) and the other for the volute (casing). Total number of cells appr. 1.15M and on the impeller avg. y+~=42 (wall functions). The GGI patches were completely non-matching: insideGGI/outsideGGI = 10080/9744 faces.

MRFSimpleFoam:
I had no difficulty setting up the problem and computing a quasi-steady (frozen-rotor) solution with both kEpsilon and kOmegaSST turbulence models. The error in the mass flux across the GGI was:
GGI pair (insideGGI, outsideGGI) : 0.0588090734 0.0588615586 Diff = -8.7289918e-05 or 0.148429338 %.

turbDyMFoam:
I managed to start a single-processor time-accurate run using the quasi-steady solution as an initial guess, but I had great difficulties keeping the pressure under control (now I know why). Once the pressure did settle - I had to use a very small time step - I managed to compute about 1.2 revolutions with the deltaT corresponding to 0.5deg per time step. So, the GGI seems to work! The pressure behavior continued to be 'noisy', but remained sensible as you can see in the picture (head.jpg).

I stopped to update my dev installation to this latest version and followed the mixerGgi-tutorial's example in changing the setup. (Btw, the zone treatment is now much better than before.) I think I got everything right, but when I restarted the time-accurate computation, the pressure problems came back two fold. Now that I had to dig a bit deeper I noticed the following:

1. Mass flux (phi) across the impeller wall is, clearly, not zero.
2. Uz on the impeller wall is not zero. (The rotation axis is the z-axis.) See the picture (Uz_impeller.jpg). In fact, they had not been zero since I started the turbDyMFoam run.

I've checked my case setup many times, but I cannot find anything wrong with it. The wall velocities worked beautifully with the MRFSimpleFoam computations. So, I guess there is a bug in the mixerGgiFvMesh, but I'm unable to identify it. Perhaps this is a enough for Hrv and Martin to get to the bottom of this. The tutorial is 2D so it doesn't yield any help on this. Has anyone else seen something like this?

This story will continue in Montreal, where I hope to present something good. I hope this problem has an easy fix -- perhaps at my end.

Best regards,
- mikko

Hi Mikko

We are doing similar simulations except that I stayed in 2D for now. I did use MRFSimpleFoam and turbDyMFoam as well, using the ggi, and I had the same problem than you, having trouble with the pressure. I came to the same conclusion than you, having a small time step solves the problem, and that is quite logical, as the mesh is physically moving, the time step should be so that you rotate the mesh not more than one or two cells each step.

However, as I am in 2D I didn't see the problem with the walls, but I have seen a case quite recently where the wall boundaries were not correctly set, but it still worked in MRFSimpleFoam, as nothing really turned.
Did you put a zeroGradient for the pressure on your walls, and a fixedValue (0 0 0) for the velocity at the wall?

Good luck!

Olivier

Hi Olivier,

The tangential velocities on the rotating wall computed by MRFZone (for MRFSimpleFoam) and mixerGgiFvMesh (for turbDyMFoam) agree nicely. It's the z-direction that is faulty. It can be easily seen how the wall velocities are computed by the MRFZone, but I cannot say the same about mixerGgiFvMesh. Hopefully someone with greater code fluency can bring further light onto this issue.

I added a new picture with a tighter scale to show how the Uz follows the z-component of the face normals, with maximum/minimum values around 45deg. If you're working with 2D geometries, you won't see this at all. You wouldn't have any 3D geometries lying around?

Perhaps the cylindrical coordinate system is at the root of this problem.

Yours,
- mikko

Hi Mikko

I see what you mean. As you said with my 2D geometry, I don't have this problem. I do however have a 3D geometry of the ercoftac pump, meshed already, maybe I could share it with you if you want. I need to clean it a bit though and to work on it a bit first.

Olivier

Hi Mikko,

What you need is is called a movingWallVelocity boundary condition: this will make sure the velocity on the surface of the impeller (which MUST be in global Cartesian coordinates) is automatically updated to make sure the mass flux is exactly zero.

Looking forward to some updates + pretty pictures to show everybody what an open source effort can lead to.

Best,

Hrv

further data
 

Hello Hrvoje,

I am using the movingWallVelocity BC. Let me provide some hard ascii data to be more specific.

When I stopped the GGI run, with the previous dev version, U and phi demonstrated the following:

U
impellerWall
{
type movingWallVelocity;
value nonuniform List<vector>
18488
(
(4.9647339 -5.2771553 -0.0018872856)
(4.7601108 -5.3303373 -0.0020557974)
(4.5670492 -5.3605736 -0.0019166547)
(4.3854936 -5.3684121 -0.0014655788)
(4.2141501 -5.3542863 -0.00090790046)
(4.0513943 -5.3195429 -0.00031373207)
(3.8939149 -5.2681244 0.00025266462)
...

phi
impellerWall
{
type calculated;
value nonuniform List<scalar>
18488
(
-1.6940659e-21
-1.6940659e-21
1.6940659e-21
0
0
0
0
...

So, we have strange non-zero z-velocities (comparatively small however) on the impeller wall, but the mass flux across it is practically zero. Up until now I had no severe problems - except at the very beginning - but then I went ahead and updated the dev version to the current one and the same results, after a very small time step, depict this:

U
impellerWall
{
type movingWallVelocity;
value nonuniform List<vector>
18488
(
(4.9682308 -5.2911397 -2.0632922e-05)
(4.762548 -5.346407 -2.2850246e-05)
(4.5680695 -5.3785635 -1.9691256e-05)
(4.3844853 -5.388396 -1.4211187e-05)
(4.2105156 -5.3759217 -8.5769172e-06)
(4.0448415 -5.3421346 -2.9240217e-06)
(3.8852911 -5.2907393 2.8557031e-06)
...

phi
impellerWall
{
type calculated;
value nonuniform List<scalar>
18488
(
-3.432816e-12
-2.0036404e-11
6.8836892e-11
5.7797203e-11
1.764198e-10
-2.6405148e-11
7.8682583e-11
-1.7541699e-10
...


The velocity values are now little different and the z-components still non-zero (still comparatively small, but a little scrolling revealed values upto 0.012). The problem is the mass flux values, which can no longer be labeled as 'practically zero'. Little scrolling revealed values upto 1.0e-7 which lie in the same ball park with the internal face values. Within this same time step the pressure had gone wild and so far I've been unable to tame it with my trials. This is where I need a second opinion.

Yours,
- mikko

Hallo, foamer,
Can anybody confirm me that turbDyMFoam working on GGI is fully paralleled? I mean whether the movingCells zone can be cut through.
And in which code module is this machenism located? I would just look it through to get a deep feeling.

Regards.

Hai

Skype failed...
 

Did anyone record the conference call? My skype bandwidth was too low to hear anything useful.

Example/Tutorial
 

Dear Foamers,
Is there any example/tutorial case working with the current 1.5-dev version? My old cases don't run anymore and I missed the skype presentation...
Any help much appreciated,
Pal

Hello all,

I missed this tutorial conference, can anyone please point me to some tutorial or examples of how to set GGI run on OF?

Thanks,
-Pablo

Some pdf were mentionned. Is it possibly to gain access to them ?

Hello,

If you want a working case with the GGI and some hints on how to set up a case with a GGI, you can go to the OpenFOAM Wiki, to the turbomachinery wiki page, there is a tutorial on how to set up a case with a GGI in it. To go there, just follow the following link, http://openfoamwiki.net/index.php/Si...vaned_diffuser.

And the tutorials are available on the svn. Just follow the steps written on the wiki page.

Good luck.

Moving part of GGI-mesh cutting a boundary
 

Hello,

allow me to bother you with a question. Before I pose it, please see the attached image

Attachment 1605.
red = static part of the GGI-mesh, blue = moving part of the GGI-mesh, black = boundaryPatch ("wheel"), gray = boundaryPatch ("road").

As you can see the moving part of the mesh is cutting the boundaryPatch ("road").

Can the present GGI implementation (svn-revision 1503 of /src) treat this behavior correctly? (Simulation is run with turbDyMFoam [svn-revision 1197])

Yes, it will work, but I am not sure if this is what you want. Do you want to impose a boundary condition on the exposed part of the rotating mesh below the road or you simply don't care what happens there?

Hrv

What I actually want is to quantitatively compare (especially Cl, besides Cd) an MRF approach and a GGI approach for automotive simulations with turning wheels and moving floor.

On the first shot I did not care about the boundary condidtion on the exposed rotating part but ran into difficulties.

Do you have any recommendation on how to treat this kind of simulations in a proper way?

We are getting too close to answering support questions on an open forum: I would need to understand much better what you are trying to do and what kind of results you expect to see.

Sorry, we should really deal with this directly - I am sure you understand.

Hrv

The goal is to use the GGI approach for automotive simulations with turning wheels and moving floor. Therefore I already ran a siumulation with MRF approach (MRFSimpleFoam) and now wanted to run the same car with GGI approach and compare coefficients. [I know MRFSimpleFoam is steady-state and turbDyMFoam not, but I wanted to give it a try.]

Sure

Hello,

I am trying to use ggi for a noncomformal interface without any motion as an alternative to stitchMesh. I created the two meshes with snappyHexMesh and merged them with mergeMeshes. I set the typ the interfacing boundaries to ggi type with:
GGI_OUT_BLOCK_maxX_1
{
type ggi;
// type patch;
nFaces 6108;
startFace 10116664;
shadowPatch GGI_IN_BLOCK_minX;
zone GGI_OUT_BLOCK_maxX_1_Zone;
bridgeOverlap true;
}

GGI_IN_BLOCK_minX
{
// type patch;
type ggi;
nFaces 5626;
startFace 10991895;
shadowPatch GGI_OUT_BLOCK_maxX_1;
zone GGI_IN_BLOCK_minX_Zone;
bridgeOverlap true;
}

Unfortunately, I get this error, when running with simpleFoam:
--> FOAM Warning :
From function GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file /opt/OpenFOAM/OpenFOAM-1.5-dev/src/OpenFOAM/lnInclude/GGIInterpolationWeights.C at line 412
polygonIntersection is returning a zero surface area between
Master face: 6092 and Neighbour face: 1197 intersection area = 0
Please check the two quick-check algorithms for GGIInterpolation. Something is missing.

From function void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
in file /opt/OpenFOAM/OpenFOAM-1.5-dev/src/OpenFOAM/lnInclude/GGIInterpolationWeights.C at line 533
Uncovered faces found. On master: 23 on slave: 0
Largest slave weighting factor correction : 0.9994067 average: 0.001887828
Largest master weighting factor correction: 0.9991873 average: 0.00420317

Floating point exception



Somehow, the 'intersection area = 0' is a bit strange to me...
Does anyone have a suggestion, what I am doing wrong!? Both patches should have the the same size (but might vary slightly due to snappyHexMesh).

Regards!
Fabian

Fabian,

it should work but setup is little tricky. Did you create faceZones from patches? Do you run in parallel? If yes did you ensure correct parallelisation (globalFaceZones)? Do the Zones match? Did you check this in Paraview?
Btw. Running it in parallel is quite slow for parallel cases. Are you happy with stitchMesh? I have lots of problems with it.

Regards BastiL

Hello,
based on the ercoftac pump tutorial, I am trying to set a turbine computation but I have few difficulties.
- I would like to simulate only one runner channel (instead of the full 360°), how should I set the GGI parameters (is there any "pitch angle" to specify ?).
I guess that I have 2 cyclicGGI bondary condition, a wall for the blade, what about the interface rotor/stator ?
Here are my BC:
10
(
SP_HIGHP
{
type patch;
nFaces 1017;
startFace 10749226;
}
SP_LOWP // this is the interface stator side (360°)
{
type ggi;
nFaces 17350;
startFace 10750243;
shadowPatch INTERFACEGGI;
bridgeOverlap false;
zone SP_LOWP_ZONE;
}
SP_WALL
{
type wall;
nFaces 143896;
startFace 10767593;
}
STAY
{
type wall;
nFaces 49852;
startFace 10911489;
}
GUIDE
{
type wall;
nFaces 45575;
startFace 10961341;
}
BLADE
{
type wall;
nFaces 24016;
startFace 11006916;
}
INTERFACEGGI // this is the rotor/stator interface - rotor side (51°)
{
type ggi;
nFaces 3136;
startFace 11030932;
shadowPatch SP_LOWP;
bridgeOverlap false;
zone INTERFACEGGI_ZONE;
}
PERIODSIDE1
{
type cyclicGgi;
nFaces 5488;
startFace 11034068;
shadowPatch PERIODSIDE2;
zone PERIODSIDE1_ZONE;
bridgeOverlap false;
rotationAxis (0 0 1);
rotationAngle 51.42857143;
separationOffset (0 0 0);
}
PERIODSIDE2
{
type cyclicGgi;
nFaces 5488;
startFace 11039556;
shadowPatch PERIODSIDE1;
zone PERIODSIDE1_ZONE;
bridgeOverlap false;
rotationAxis (0 0 1);
rotationAngle -51.42857143;
separationOffset (0 0 0);
}
INLET
{
type patch;
nFaces 3360;
startFace 11045044;
}
)

- Secondly, I think I don't totally understand the difference between MRFSimpleFoam and simpleTurboMRFFoam. The first one is frozen-rotor and is the second one what CFX call "stage" ?

Is there any tutorial using simpleTurboMRFFoam or turbDyMFoam ?

Thanks a lot !
Flo

Hello Fabian,

Check your mesh for concave faces on the GGI patches.

The cutting algorithm being used by the actual GGI implementation cannot handle concave faces on the GGI patches.

snappyHexMesh is known to generate mesh with concave faces.

You can check your mesh for the presence of concave faces using checkMesh. But watch out, checkMesh concave angle threshold is 10 degree, so all concave angles from 10 degree and less will not be reported.

You can change this default setting in your main controlDict by changing the tolerance parameter primitiveMeshFaceAngleThreshold.

Martin

Hi Learned Ones

I have a similar problem that works/rotates with icoDyMFoam but gives a floating point error after the first step with turbDyMFoam. I'm using the movingWallVelocity boundary condition in the U file for the blade which I assume is also necessary for the icoDyMFoam to work.

The GGI is working for the icoDyMFoam case so I don't think it's a mesh interface problem.

I've tried reducing the time and changing the number of p correctors enough that I don't get any changes with further iterations but am clueless as to what to try next.

Any suggestions would be gratefully received

Thanks

Nick

This is the output

Initializing the GGI interpolator between master/shadow patches: InterT/InterR
Evaluation of GGI weighting factors:
Largest slave weighting factor correction : 2.85271e-05 average: 2.77807e-06
Largest master weighting factor correction: 4.47993e-05 average: 7.11812e-07

BiCGStab: Solving for Ux, Initial residual = 1, Final residual = 2.71511e-08, No Iterations 7
BiCGStab: Solving for Uy, Initial residual = 1, Final residual = 4.66166e-08, No Iterations 8
BiCGStab: Solving for Uz, Initial residual = 1, Final residual = 4.14282e-09, No Iterations 6
BiCGStab: Solving for p, Initial residual = 1, Final residual = 9.96992e-08, No Iterations 267
BiCGStab: Solving for p, Initial residual = 0.705496, Final residual = 9.4242e-08, No Iterations 211
BiCGStab: Solving for p, Initial residual = 0.0302484, Final residual = 3.65547e-08, No Iterations 199
time step continuity errors : sum local = 9.62158e-14, global = 2.47271e-15, cumulative = 2.47271e-15
BiCGStab: Solving for p, Initial residual = 0.046603, Final residual = 5.78336e-08, No Iterations 206
BiCGStab: Solving for p, Initial residual = 0.00832937, Final residual = 4.69163e-08, No Iterations 147
BiCGStab: Solving for p, Initial residual = 0.000697974, Final residual = 9.73038e-08, No Iterations 78
time step continuity errors : sum local = 2.45325e-13, global = 5.62079e-14, cumulative = 5.86806e-14
BiCGStab: Solving for p, Initial residual = 0.00268834, Final residual = 6.2455e-08, No Iterations 178
BiCGStab: Solving for p, Initial residual = 0.000414399, Final residual = 4.79598e-08, No Iterations 125
BiCGStab: Solving for p, Initial residual = 3.51465e-05, Final residual = 9.8068e-08, No Iterations 34
time step continuity errors : sum local = 2.47464e-13, global = -1.67022e-13, cumulative = -1.08342e-13
BiCGStab: Solving for p, Initial residual = 0.000847403, Final residual = 2.87118e-08, No Iterations 154
BiCGStab: Solving for p, Initial residual = 0.000142032, Final residual = 9.4255e-08, No Iterations 127
BiCGStab: Solving for p, Initial residual = 1.09702e-05, Final residual = 9.90043e-08, No Iterations 20
time step continuity errors : sum local = 2.4977e-13, global = 1.42343e-13, cumulative = 3.40009e-14
Floating point exception


Thanks again

Hey Nick, could you upload your case file so we could have a look at what you're trying to simulate?

Hi James

These are the case files less the mesh which was too big to upload with this lot. If you need it just say. I've got about three iterations running now before a floating point error again.

Thanks Nick