[brammekeuh987]

Dear Foamers

Currently I'm trying to create a mesh by using SHM on a very detailed blockmesh. All curves and edges are already predetermined in the blockmesh (see blockMesh.png) .
The only reason I use SHM is to refine the mesh near the fixed objects, in such a matter that the overal mesh quality remains good. No stretching or turning or whatever is needed, just cut up the blocks (as in castellation.png) and optimize the edges a bit where the refinement changes level in order to reduce non-orthogonality and off you go.

Unfortunately, the mesh buckles in the Z-direction, which shouldn't move at all. In leadingEdge.png and trailingEdge.png you can see the buckling near the leading and trailing edge of the wing together with a featureEdge (in red). The snapping also doesn't happen to the featureEdge of the surrounding cylinder as seen in snapping0.png

So I already tried following techniques:

1. Searachble disks in the geometry input of my SHM-dict. These disks are used in the refinement surface (level 0) in order for the snapping to take place on the disk.
2. featureEdges, which are generated from an stl, which correspond with the most outward edges of the geometry.
3. Playing around with the quality of my stl. As Tobi pointed out in his thread (Creating waterproof STL using snappyHexMesh or salome ) this could influence the meshing process. Different quality types of stl have been used.

None of the above seemed to work. The SHM-log file gives me following warnings:
-FOAM Warning : Displacement [....] at mesh point [....] coord [....] points through the surrounding patch faces
-FOAM Warning: [...] Did not successfully snap mesh. Continuing to snap to resolve easy surfaces but the resulting mesh will not satisfy your quality constraints

which looks like a logical warning given the situation.

I looked around in the 'flange' tutorial, but it didn't seem to help. The settings of my snapControls are rather 'normal':

snapControls
{
    nSmoothPatch 5; 
    nSmoothInternal $nSmoothPatch;
   
    tolerance 4.0; 
    nSolveIter 100;     
    nRelaxIter 5;    
    nFeatureSnapIter 10; 
    implicitFeatureSnap false; 
    explicitFeatureSnap true ;
    multiRegionFeatureSnap false ;
}

Some imperfections:

• the extractFeatureEdge extracts per wing 4 'strange' strait lines inside the stl. Since these are located not near the edge, I would doubt this would influence the buckling on the sides. This only happens with the 'high' quality stl, the low quality doesn't have this'.
• the trailing edge in the stl is made separately, so the featureEdge extracted from it is not 'perfect' but interrupted on the upper and bottom side of the wing. Since the buckling happens all the way across the wing, this won't do any harm (I think)

Anyone any tip/hint to get this straightened out?

The case can be found at: https://we.tl/t-tqdzILyacA
The total size is 150Mb (because of the stl's), pictures in attachment.

OF-version: v1906

Thanks!
Bram

[jherb]

I have downloaded your case and used the Python code given at https://codereview.stackexchange.com/a/101298 to look at your mesh data. (I first had to convert the case into ascii file format)

What you could do, is using this Python class to find the points of your mesh beloging to the back and front and set the z-coordinate to a value you want. Then you need to save the full point file in constant/polyMesh.

Perhaps, there is also a way to do this with OpenFOAM internal utilities, but I can't remember, if there was and which.

Here is the code to extract the min/max values of the z-coordiantes of your front and back:

back_points = np.array(pm.points)[np.array(pm.faces[int(pm.boundary['back']['startFace']):int(pm.boundary['back']['startFace'])+int(pm.boundary['back']['nFaces'])]).flatten()]
front_points = np.array(pm.points)[np.array(pm.faces[int(pm.boundary['front']['startFace']):int(pm.boundary['front']['startFace'])+int(pm.boundary['front']['nFaces'])]).flatten()]


print(np.min(back_points[:, 2]), np.max(back_points[:, 2]))
print(np.min(front_points[:, 2]), np.max(front_points[:, 2]))

[brammekeuh987]

Hi Joachim

You're coming back to my rescue, thanks ! I see your answer a bit late, since I already made some progress, which I was planning on reporting now, but I'll try your approach later on. Much appreciated!

So I found out some more things and I have corrected some fundamental errors. It was this site from Elia Agnani which put me in the right direction (https://sites.google.com/site/snappywiki/ )

-First Error:
Searchable objects / user defined regions, are ONLY used for refinement regions and NOT for refining surfaces or snapping! You can only use the imported objects (stl's, obj's, ... ) for snapping.
This was a biggy and already improved the surface a lot.

-Second Error:
I made the cylinder and the disks separately in an stl-format. I needed to make them separately because I use FreeCAD and it can't make a multi region stl, but further this doesn't have any implications. But wat was more important was that the holes where the objects go through are also in the disks (see cylinderSTL_front.png) and the objects that protrude the disk didn't have any sides (see objectSTL.png)! It's really important that you only have the stuff in the STL which should attract points, no more (!!) and no less (obviously ). Otherwise it still wrinkles and buckles a lot.

In attachment the improved result (improved_result.png)

But it isn't perfect yet.
One edge isn't snapping yet (no_inner_snapping.png) and it still buckles a bit (slightBuckling.png). But by the looks of Joachim's python code, this could be easily resolved I think.

I'll come back when I tried Joachim's magic and/or find something more!

Cheers
Bram

[brammekeuh987]

So played a bit more around and I have following remarks:

1. try to keep the nSmoothPatch as low as possible. This smoothing algorithm enhances the mesh when you work with a non-perfect mesh. But since your patch is already in an almost perfect state it just distorts things and it causes non-orthogonality and so on. I put it on 0 (but this gives other issues);
2. put the nSmoothInternal firmly large (20). The nSmoothInternal doesn't get activated when the nSmoothPatch is 0.
3. tolerance low (1). Because you don't distort the mesh too much, so snapping should go fairly easy. The only snapping that needs to be done is the better approximation of the cells to the surface after the castellation step.
4. nFeatureSnapIter around 20, just to make sure nothing goes wrong
5. don't do any refining of feature Edges. The refinement does help with the snapping, but you're always messing with the 'transition' cells in between levels. In my case this means that I have these transition cell's with a greater non-orthogonality in my boundary layer around my objects, which is not desired.

When putting the nSmoothPatch 1, it messes around with the second row of points away from my featureEdge and causes a non-orthogonality of 60 (which in my case normally should be 18-ish) in the face between the first and second row of cells (see nSmoothPatch1.png). I'm still trying to get it out, because otherwise I need an external smoothener for the mesh. And talking about smootheners is a bit of topic here, but as a starter I tried extBlockmesh, but it doesn't work with the not completely hexahedral mesh that you get after SHM. So, I'll try to look a bit both ways (getting it out SHM or with an additional smoothener)

If anyone has other experiences or other ideas with SHM, let me know!

I also tried Joachim's approach: this does exactly what it should do: you can set all coordinates to whatever you want, the only problem is that with a large nSmoothPatch, the cells tend to change from patch. So when you set all the z-coordinates of a certain patch to a fixed coordinate, the mesh gets distorted. Underneath the code I used in python:

from testFolder import parsMesh

pm=parsMesh.ParsedMesh(meshDir)

back_points_indices=np.hstack(pm.faces[int(pm.boundary['back']['startFace']):int(pm.boundary['back']['startFace'])+int(pm.boundary['back']['nFaces'])])
front_points_indices=np.hstack(pm.faces[int(pm.boundary['front']['startFace']):int(pm.boundary['front']['startFace'])+int(pm.boundary['front']['nFaces'])])


for index in back_points_indices:
    pm.points[index][2]=47.100
for index in front_points_indices:
    pm.points[index][2]=2*47.100

    
back_points = np.array(pm.points)[back_points_indices]
front_points = np.array(pm.points)[front_points_indices]
print(np.min(back_points[:, 2]), np.max(back_points[:, 2]))
print(np.min(front_points[:, 2]), np.max(front_points[:, 2]))

writeFile=open(meshDir+"points","w")
writeFile.write("//########################### HEADER FILE BIT**ES ########################### \n\n")
writeFile.write("FoamFile\n{")
for x in pm.points_hdr:
    strr=str(x)+"\t\t"+str(pm.points_hdr[x])
    if str(x) == "location":
        strr=str(x)+"\t\t\""+str(pm.points_hdr[x])+"\""
    #print (strr)
    writeFile.write("\t"+strr+";\n")
writeFile.write("}\n")
writeFile.write(str(len(pm.points))+"\n( \n" )

for punt in pm.points:
    strr="("
    for el in punt:
        strr+=str(el)+" "
    strr=strr[:-1]+")\n"
    writeFile.write(strr)


writeFile.write(")")
writeFile.close()

print("Done")