Hex mesh
 

I want to create a hexahedral mesh of a channel with a parabolic profile..

Any suggestions on how to do this?

Salome geometry of this profile and extrusion is generated in two ways, spline for parabolic profile and also a number of vertices lying on the parabola joined by straight lines to approximate the profile.

I have been working on this for the last 2 weeks and have not figured out how to do it yet!!!!!

Please make a sketch of what you are trying to do, it will help us to help you.

Hi Niels, thank you very much for your reply!

I have attached 2 screenshots of the geometry basically im trying to mesh a 3D extrusion of this profile..

I have tried different ways of generating the profile; using a series of points connected by straight lines, or using a spline through all of those points to create the parabolic shape.

If you upload the hdf I can do that for you and give you the hdf back, so you can see what have been done.

I will post the hdf as soon as I can.

I have meshed this geometry in hexahedral cells using the body fitting algorithm + create boundary faces, hex (i,j,k) + quadrangle mapping, also automatic hexahedralisation methods of meshing. The problem I have is that I need to export to unv and then convert using unvToFoam.. The problem is that the export to export from Salome to unv doesn't handle certain type of cells; or the conversion from unv to openfoam does not work it give an error message as follows:

ideasUnvToFoam:*ideasUnvToFoam.C:881: int main(int, char**): Assertion `nrFaceCells[faceI] == 1 || nrFaceCells[faceI] == 2' failed.
Aborted

Hi Niels,

I have attached a dump of the hdf file. Kind Regards,
Matthew Skaras

Hello

It was not that difficult, see the attached dump and a little teaser screen.

The only thing missing was to split every face into something that has only 4 edges.

I can then create a 2D sketch and make segments of blocks that has 6 faces like a cube.

Then in the mesh you can use hypothesis automatic hex mesh.

Best

Hello Niels,

Thank you for your help, I have tried this method also, the problem is I need further refinement in some areas in order to control the aspect ratio of the cells.

The problem here is in the conversion process- I can convert the mesh to openfoam and view it in paraview but when I run checkMesh I get a long series of the following errors:

Checking topology...
Boundary definition OK.
--> FOAM Serious Error :
From function bool zone::checkDefinition(const label maxSize, const bool report) const
in file meshes/polyMesh/zones/zone/zone.C at line 207
Zone cellZone_0 contains invalid index label 32000
Valid index labels are 0..31999
--> FOAM Serious Error :
From function bool zone::checkDefinition(const label maxSize, const bool report) const
in file meshes/polyMesh/zones/zone/zone.C at line 207
Zone cellZone_0 contains invalid index label 32001
Valid index labels are 0..31999

and ends as follows:

--> FOAM Serious Error :

--> FOAM FATAL ERROR:
Too many errors

From function messageStream::operator OSstream&()
in file lnInclude/messageStream.C at line 200.

FOAM aborting

I want to run cyclic patches and I cannot because of these errors I believe

Kind Regards,

Matt

I think you need to give a sketch or something of what and where you want to refine.

You can create groups of edges in Salome Geo and use those to create submeshes of the overall mesh.

This is how I usually do it when making graded hex mesh in Salome.

First of all thank you very much for your help so far! It is much appreciated.

Sorry I think I may have been a bit unclear. I am not having trouble in creating submeshes for refinement.. The issue I have is when I do create submeshes (and export to unv then convert to OF) I get errors with cell zone references (see error message above) when running checkMesh. The same mesh exported to unv prior to creation of submeshes converts from unv to OF and running checkMesh gives no errors.

I am wondering if you could perhaps explain why this may be happening?

Can you post the py dump of the whole shebang with mesh?

Attached is a dump of the hdf,

Cheers

Matt

The steps I then take are to export the mesh to unv format then run ideasUnvToFoam <filename>.unv

Running checkMesh gives the above mentioned output.

Also note I have not named any of the boundary faces as I generated this mesh purely to test if it converted to OF without any errors. For the actual mesh I would use I would create groups of faces in the mesh (create group and select using selection filter, leaving the walls along the z-axis undefined and renaming the defaultFaces patch to walls) and then export as above.

Kind Regards

Matt

Hi

I get no problems at all.

Could you try and delete your constant/polyMesh folder before converting.

If that does not fix it something must be wrong in your OF install.

Here is the two output files.

Thank you for such a quick response!

I also get no errors after deleting constant/polyMesh- my output look the same as what you posted.

However after converting and changing constant/polyMesh/boundary to add cyclic patches and wall patches and rename defaultFaces patch I get the following output when running setFields and checkMesh:

matt@matt-desktop:~/Desktop/parabolicHex3D$ setFields
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.0.1 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.0.1-51f1de99a4bc
Exec : setFields
Date : Jul 17 2012
Time : 19:55:35
Host : matt-desktop
PID : 3581
Case : /home/matt/Desktop/parabolicHex3D
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0



--> FOAM FATAL ERROR:
face 0 area does not match neighbour by 91.7559% -- possible face ordering problem.
patch:front my area:0.000244545 neighbour area:0.000659134 matching tolerance:0.001
Mesh face:143500 fc:(-1.09222 0.697452 0.15)
Neighbour fc:(-0.28762 1.34188 0)
If you are certain your matching is correct you can increase the 'matchTolerance' setting in the patch dictionary in the boundary file.
Rerun with cyclic debug flag set for more information.

From function cyclicPolyPatch::calcTransforms()
in file meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C at line 220.

FOAM exiting

matt@matt-desktop:~/Desktop/parabolicHex3D$ checkMesh
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.0.1 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.0.1-51f1de99a4bc
Exec : checkMesh
Date : Jul 17 2012
Time : 19:55:46
Host : matt-desktop
PID : 3583
Case : /home/matt/Desktop/parabolicHex3D
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0



--> FOAM FATAL ERROR:
face 0 area does not match neighbour by 91.7559% -- possible face ordering problem.
patch:front my area:0.000244545 neighbour area:0.000659134 matching tolerance:0.001
Mesh face:143500 fc:(-1.09222 0.697452 0.15)
Neighbour fc:(-0.28762 1.34188 0)
If you are certain your matching is correct you can increase the 'matchTolerance' setting in the patch dictionary in the boundary file.
Rerun with cyclic debug flag set for more information.

From function cyclicPolyPatch::calcTransforms()
in file meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C at line 220.

FOAM exiting


I cannot figure out why this error is occuring!!!

Can you give med the case without the mesh and then a list of commands/edits you go through to obtain your final mesh now?

I think its in your general setup some minor thing is preventing it from working properly.

Attached is my case (below)


To convert to the final mesh I first run ideasUnvToFoam
I then run checkMesh- no errors
I then view with paraFoam
I then edit to constant/polyMesh/boundary file to what is given in the case
I then run setFields and get the error above
So I re-run checkMesh and get the error above again

I also tried changing the matchTolerance to 1 and it worked however the initial field is not correct..

Thank you for all your help so far!

Kind Regards

Matt

Hi

I see no files

Attached is the case directory

removed the polyMesh directory since it was too large

cheers

here is the boundary file again

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.0.1 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

5
(
front
{
type cyclic;
nFaces 5000;
startFace 143500;
matchTolerance 1e-3;
neighbourPatch back;
}
back
{
type cyclic;
nFaces 5000;
startFace 148500;
matchTolerance 1e-3;
neighbourPatch front;
}
outlet
{
type patch;
nFaces 1000;
startFace 153500;
}
inlet
{
type patch;
nFaces 4;
startFace 154500;
}
walls
{
type wall;
nFaces 1996;
startFace 154504;
}
)

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

i JUst realised that you might need this

cheers

Hi

Your setup is really not good. Looking at the boundaries in paraview from the file you have provided give me the following view.


If I use autoPatch 5 on the other hand I get the last picture.
Then you can always rename the patch names in the boundary file in the new time folder.

So its no wonder why checkMesh and setFields complains.

To avoid autoPatch you need to export the BC names from Salome using "create groups from geometry" option in the mesh module.

You need to define your groups in the GEO module. Right click the geo and select "Create group" use faces, select each face for the individual BC's and then in Mesh module use create groups from geometry and use the groups in the geo module in the tree structure.

Hello,

Given your last post, I must divide the parabolic curve into 2 parts (2 half parabolas) with a box in the middle to allow me to create a group on the geometry for the inlet see images for clarification (underside shows the inlet and middle box more clearly)..

This makes meshing the entire domain more difficult as I need to create many submeshes to ensure relatively uniform cell sizes when using the nb segments 1d hypothesis for wire dicretisation algorithm.

As you can understand this is very tedious.. so what I would like to know is can I mesh the 2 volumes either side of the box in the middle with salome and use blockmesh to generate the middle box and then use mergeMesh to create the overall mesh?

OR is there a more efficient way of doing this problem in Salome?

Also a question regarding export of salome mesh to unv.. there are certain cell type the export utility in salome does not handle, is this correct? If so what are they?- Mainly concerned with triangular prisms, since quadrangle mapping the geometry shown creates some triangles

I am struggling!! I have been working on this for weeks now and have been unable to find a good solution.. There must be an easier way to do this in salome!!!

Kind Regards

Matt

Hi

You can also use the mesh as is and use surfaceToPatch command.

This takes an stl/obj as input and will convert the faces that matches the surface file into a patch.

You might have to play with the -tol option so that it doesnt pick faces that are part of the cyclic patches.

By use the mesh "as is" do you mean the mesh I sent you (.py dump) or my new mesh without any groups of faces added to the mesh?

I have tried surfaceToPatch with no success im not sure if I am using it correctly?- I just exported my mesh to both unv and stl then converted the unv mesh and ran surfaceToPatch with the stl file I created.

I also tried creating groups from geometry in the mesh module- this does not work either I still have issues with face match tolerance for front and back faces!!

Yes to the first.

You need to make stl files that only cover the BC's, so for the inlet you would export only the face to an stl that has the inlet faces.

You can export directly to stl from the geo module. So just explode the geometry in your previous post into faces and export the small strip which I assume is the inlet to an stl.

If you just export one stl for the whole geometry it will create a patch with all the boundary faces.

So just to clarify..

I should use the geometry with no groups of faces and then create stl files for each boundary and then use surfaceToPatch to name the boundaries?
I have tried this but every time I run surfaceToPatch i get a new time directory with a new polyMesh dir, can I merge them all into 1?

Or do use the original geometry (ie the a single parabolic curve for the curved profile) i sent you with named groups for all boundaries except for inlet (which does not yet exist) and then use surfaceToPatch for the inlet only?

I would go for the last method.

And you just move the new polyMesh from the new time directory back into constant. Just remember to rename constant/polyMesh to constant/polyMesh.bak that way you have a backup if something gets screwed up.

You can also use paraview to select the faces you want to be in the inlet bc and then triangulate those cells using a filter then go to files save data and you can save to stl.

This way you wont have to use Salome to create the stl and you are sure that the stl lies exactly on the faces you want.

http://www.vgtc.org/PDF/slides/2008/..._selection.pdf

So I have used Paraview to generate the stl file for the inlet and everything works well except for one problem I use setFields to initialize alpha I get the following error:

--> FOAM FATAL ERROR:
More than six unsigned transforms detected:
6(((2.18444 0 -0.15) (1 0 0 0 1 0 0 0 1) 0) ((2.18702 0 -0.15) (1 0 0 0 1 0 0 0 1) 0) ((2.18952 0 -0.15) (1 0 0 0 1 0 0 0 1) 0) ((2.19195 0 -0.15) (1 0 0 0 1 0 0 0 1) 0) ((2.19431 0 -0.15) (1 0 0 0 1 0 0 0 1) 0) ((2.1966 0 -0.15) (1 0 0 0 1 0 0 0 1) 0))

From function void Foam::globalIndexAndTransform::determineTransforms ()
in file primitives/globalIndexAndTransform/globalIndexAndTransform.C at line 185.

FOAM exiting

I just set up the case using the autoPatch approach and renamed the BC's and ran bubbleFoam on it without an inlet and it runs fine for a couple of iterations without any cyclic errors.

I think your patches have just been defined wrong so that front patch does not match back patch because faces from the wall/outlet etc is also in the cyclic patches.

sorry.. I edited the post above. It now works only problem is with paraview as mentioned above after running setFields

I dont get any errors after running setFields.

Here is the final case. I've disabled turbulence as it make it run very unstable.

https://dl.dropbox.com/u/15968063/case.tar.gz

I think something is wrong with your specification of the densities as is seems alpha1 mixes in a strange way, see the pics.

EDIT: Its not the densities its the BC specification in the 0 folder, look to the bubbleColumn tutorial case for the right setup.

EDIT2: I would also consider making the mesh more "straight" and refine around the interface. Having these coarse cells which are angled to the alpha interface is bound to give you stability issues.
Also is the mesh in m or mm, because checkMesh says it has the extend of roughly 2.24m wide and 1.5m high, and the mesh seems quite coarse for such large dimensions.

Hello,

Thank you for supplying the mesh,

I have a customized version of bubbleFoam to account for turbulence when the phases are separated ie where the phase fraction approaches 1. I have run this solver on a column geometry meshed with blockMesh and it is quite stable with turbulence ON. When I come to running on this parabolic geometry it is not so stable and I have tried on tet meshes and the hex mesh you have posted.

In regards to EDIT 2, yes I agree, I have been having trouble with this for the last few weeks and the mesh you got working for me was only an initial first run... However it proved unstable also. Yes the mesh is very coarse- the domain is in m.

Can you suggest a more appropriate solution for meshing this domain?
It apparently more difficult than it seems!

Ps I can send you the modified bubbleFoam solver if you would like to have a look


Kind Regards,

Matt

Hi

After using a laminar approach for the case with modified BC's and ran it for about 40s I can see that the domain is filling with alpha0. Since you have alpha0 set as 0.6 at the inlet you will start filling the domain.

Don't know if this is what you intended?

To get this case being more stable you can try the following things.

• Refine mesh at interface
• Make the mesh to be more parallel to the interface, at the interface
• Change the fvSchemes to upwind for k/epsilon and maybe even for Ua and Ub
• Change the cyclics to symmetryPlanes, just to see if it influences stability.
• Make the height of the alpha1 portion bigger to help numerics at the outlet

Hi, tried changing all but the cyclic option in my list and no change in convergence (maybe a little but still crashes).

After changing the cyclics to symmetryPlanes there is absolutely no issues.
So somehow the cyclics are screwing with the convergence.

We have discussed a similar behaviour here.

http://www.cfd-online.com/Forums/ope...lgorithms.html

By alpha0 I assume you mean the phase fraction alpha (dispersed/gas phase) to be 0?

The reason for setting to 0.6 at the inlet was to increase the flow rate of air into the domain without having to increase the velocity.. Should I leave this as 0.5 and then alter velocity only to achieve the desire inlet flow rate?

REgarding refinement at the interface would you suggest having "straight/parallel" cells below and above the interface or will it suffice to have this refinement only above the interface?


So it appears that the cyclic boundaries were causing the issues..

I looked at your mesh and would like to reproduce it using the autoPatch method. Every time I run autoPatch I am unable to determine which patch has been allocated to which face since paraview crashes every time I run paraFoam..
Could you outline the way you do it?

Just an aside, I also tried running with symmetryPlane patches on a relatively coarse tet mesh without any improvements in stability, (note I had also changed gradschemes to cellLimited leastSquares in fvSchemes- to improve stability on the unstructured grid)

Thanks for your help with everything I very much appreciate it!

Kind Regards

Matt

Hi

Well I do it the way you were trying to do it with paraview, if that don't work your back at guessing or you need to make a stl for every surface and run surfaceToPatch on all of them.

I would highly recommend getting paraview to work though.
just Install paraview manually and skip the whole paraFoam stuff.

Just open paraview and open system/controlDict and if it asks about file type just choose openfoam, as paraview > 3.12 has native support for OF.

Setting alpha to 0.6 is fine, after thinking a little about the physics its clear that the domain has to start filling since alpha0 has no where to go until it hits the outlet so probably adding outlets as shown in the pics (black bars) will be a good idea.

I've also made the case available to you and the study dump for Salome.
I had to create a 2D mesh and extrude it. If I didn't I got the OF errors about the cyclics not matching, but if you are to use symmetryPlanes it does not matter.

https://dl.dropbox.com/u/15968063/Case2.tar.gz

How do you create a 2d mesh and then extrude? I have been trying to do this for a while but the method I have used instead (due to lack of knowledge on extruding a mesh) is to use quadrangle mapping and then Hexahedron(i,j,k) algorithm.

Thanks for helping me out with this! It look like I'm going to have to play around with salome and ill most likely use autoPatch to name the boundaries.

One other point; I downloaded the latest Paraview (3.14.1) and tried to open the first case you sent me.. it seemed to have trouble with the fact that there were multiple entries for walls in the boundary file, however using paraFoam to open the case from the case directory had no issues. This wont really matter because after running autoPatch there wont be any duplicated entries in the boundary file, and after I name them I will be able to view them using paraFoam.


Kind Regards,

Matt

Go to the mesh module.

Mesh -> create mesh -> assign a set of hypotheses -> 2d: Automatic quadrangulation -> specify number of segments -> apply and close

Mesh -> Compute

Modification -> extrusion -> select triangle at top -> tick select whole mesh, submesh or group (select mesh if not already selected in the tree) -> Extrusion to distance -> Z = 0.01 -> Number of steps = 15 -> Apply and close

Thank you!

That makes meshing so much easier!!

I have a question regarding naming patches- Using paraFoam to view the mesh allows the mesh to be read even though there are multiple patches with the same name however it only displays the first one in the list in the boundary file. Using the latest version of paraview and opening the controlDict will not display the mesh since there are multiple patches with the same name.

My question is when solving does OpenFOAM recognise multiple wall patches with the same name?

Kind Regards

Matt

It can for serial run.

For parallel runs you need to use wildcards and name each wall patch by itself.

You can see how I have used wildcards in the Case2.tar.gz case.

So look in the boundary file and f.ex the 0/U file you will see how wildcards work.

If you want them all in one singular partch you can use the createPatch utility.

There is a sample dict for that in the OF install folder just do this under the OF install folder

Should get you the location of it