This is my first post on this site, since I only recently started practicing with OpenFOAM. So far, however, I noticed that for an open-source software suite that is not thát widely used, it is quite mature. However, not all settings are well documented, and unfortunately a lot of users of the cfd-online forums don’t post their solutions if they got their case finally working. That is also why I started this whole website, to show you how I approached certain cases that are relevant to my own work.

Download the case files here: duct_case.tar

A ducted propeller in OpenFOAM

In my PhD project, I focused on inland ships. These are self-propelled ships and due to limitations on their dimensions (locks, bridges, fairway dimensions), these ships are really full (that is, a high block coefficient ). Together with limited draft imposed by limited water depth, these ships can only fit a small propeller. To increase the amount of thrust and power per square meter of propeller area, and to improve efficiency, inland ships are equipped with ducted propellers.

A 19A duct

The fixedJump boundary condition

Looking around the cfd-online forums and OpenFOAM-related websites in general, it seemed like modeling an actuator disk as used in Maritime Hydrodynamics, a pressure jump, is not straightforward. It apparently is not possible to – without coding – add body forces in a certain region of the mesh (at least I did not manage to, if you know how, leave a comment!). Looking through the C++ guide of OpenFOAM 4.1, I came across the FixedJump boundary condition. It uses a cyclic boundary condition at its base, but allows a fixed difference for a certain field variable between the two connected boundaries. This is exactly what an actuator disk requires. All field variables should be continuous, except for a fixed jump in pressure. The value of the jump can be computed from the required total thrust force and the area of the disk.

Creating baffles in snappyHexMesh

In order to implement the fixedJump boundary condition, we require that there are two patches, with the same points and topology and preferably very, very close to each other. I tried multiple approaches – I did not find the right solution on the cfd-online forum, such as creating a closed solid inside the duct, and using the opposing sides of this solid to implement the cyclic and fixedJump boundary conditions. However, snappyHexMesh puts different, non-conformal vertices on both sides. The cyclicAMI boundary conditions – which should allow for a mapping between two patches with unequal topolgy, did not work either.

Then I found: OpenFOAM 2.2.0: snappyHexMesh. Within the refinementSurfaces sub-dictionary in snappyHexMeshDict, it is possible to define a faceZone. This faceZone can be defined as a baffle, i.e. a very thin surface. With this setting, snappyHexMesh will ensure that interior grid points are snapped to this surface, and finally create two patches (a master and a slave patch). In my case, the surface got the name ‘disk1’, and therefore two cyclic patches, ‘disk1’ and ‘disk1_slave’, will turn up in the resulting mesh.

My approach

I used the creation of baffles in snappyHexMesh as well as the fixedJump boundary condition to create the pressure jump across the ducted propeller. Furthermore, I used simpleFoam to solve the flow equations. simpleFoam is fast and thus suitable to find out how a ducted propeller should be modeled.

Adding the faceZone

First, I have to add the faceZone to the snappyHexMeshDict file. This is done as shown in the following code block. During meshing, snappyHexMesh will create two internal patches: disk and disk_slave. The reason it uses ‘disk’ is because it is defined with that name in the geometry section of the snappyHexMeshDict file.

Unfortunately, I did not manage to find how to directly set the internal patches as cyclic patches from snappyHexMesh, it seems to ignore additional patchInfo settings for faceZones. However, this can be overcome using createPatch (or createBaffles, if you like).

castellatedMeshControls
{

    features
    (
        { file "disk1.eMesh"; level 1; }
    )

    refinementSurfaces
    {
        // add the duct (and other surfaces) here
        
        disk1Object
        {
            level (1 1);
            faceZone diskFaces;
            faceType baffle;

        }
    }

    // Remaining settings
}

Some snapping and layer addition experience

In the case that can be downloaded at the bottom of the page, note that I use two separate dictionaries for snappyHexMesh: one for refinement and snapping, and a second one to do the layer addition. This allows me to first check whether the surface is snapped smoothly to the duct (which is fortunately is), before adding layers (which otherwise fails in places where the mesh is badly snapped and ‘bumps’ occur on the patches).

Apart from the separated snappyHexMeshDict files, I also subdivided the duct STL in four parts. This allows me to have the leading and trailing edge separated, and simply fix the refinement level of the mesh on both. The reason for this is that if snappyHexMesh is refining based on feature angles or local curvature, small ‘islands’ of higher refinement levels occur, and those seem to be locations where snappyHexMesh clearly has troubles getting a well-snapped mesh. Hence, I use equal values for minimum and maximum refinement per surface. The border between the lower and  higher refinement levels then becomes an almost straight line, something that snappyHexMesh seems to work well with.

The snapped mesh around the duct. Note the (almost) sharp border between different refinement levels. This seems to be important for snappyHexMesh to snap the surface well enough. I did need to have 2 buffer layers (nCellsBetweenLevels) though. With this setting at 1, there would still be some bumps.

Another thing I found out for this specific case is that in the snapControls sub-dictionary, it is useful to keep nSmoothPatch low, but increasing nSolveIter doesn’t hurt and may increase the quality of the final mesh. Futhermore, although I added a specific feature edge for the connection between the duct and the disk patch (that will later be used to put the pressure jump on), the setting explicitFeatureSnap does not make things better. In fact, it was better for me to use implicitFeatureSnap.

When adding layers, I feel that snappyHexMesh prefers to have each of the surfaces on which you want layers specified separately rather than just using the patch group ‘duct’ (which I in fact use to define the no-slip BC later on). Hence, like this:

layers
{
  leadingedge	{ nSurfaceLayers 3; }
  trailingedge	{ nSurfaceLayers 3; }
  innerring	{ nSurfaceLayers 3; }
  outerring	{ nSurfaceLayers 3; }
}

With regard to settings, increasing nSmoothNormals – which smoothens the normals of the interior grid, not those on the patch – works pretty well. Also in cases with a more complex geometry, I found that it is worthwhile to wait a bit for this process. With nSmoothSurfaceNormals, however, you should be careful. For example, the normals close to the trailing edge (but still on the ‘flat’ part) will start pointing backward if too many surface smoothing iterations are used, which creates a distorted layer mesh.

Adding the fixedJump boundary condition

With a nicely snapped mesh including good-looking layers (y+ is 100, which is OK when using wall functions), the boundary conditions on the patches can be defined. In the case you can download, you’ll see that createPatch is used to ensure that all patch names are correct, that all parts of the duct are grouped under ‘duct’, and – most important – that the internal baffle faces are set to be cyclic patches. I also give them somewhat better names, namely disk_in and disk_out.

Next, the boundary conditions need to be defined. This is done in the ‘0’ folder. We need to implement a pressure jump, so I’ll start with what is required for the ‘p’ file. I do assume, however, that other boundary conditions for ‘p’ are known to  you already.

disk_in
{
    type	fixedJump;
    patchType	cyclic;
    jump	uniform -0.8;
}
disk_out
{
    type	fixedJump;
    patchType	cyclic;
}

Note that in the above code, only one of the patches gets assigned the jump! This is actually the master patch. Remember that snappyHexMesh created both ‘disk’ and ‘disk_slave’. Using createPatch, I renamed the first to ‘disk_in’ and the latter to ‘disk_out’. Hence, the jump should be defined on ‘disk_in’, which is the master patch.

As for the value of the jump, I simply wanted a thrust coefficient (of the propeller only!) between 1.0 and 2.0:

I used a thrust coefficient of 1.6, and since the used duct has a radius of 1.0 meter, and I use a speed of 1.0 meter as well, the total thrust is approximately 2500 newtons. Dividing this by the are of the disk (which conveniently is exactly ), the pressure jump is around 800 newtons per square meter. However, simpleFoam, being an incompressible, single-phase solver, uses the pressure divided by density for the ‘p’ variable! Thus, to have a 800 pascal pressure jump, the jump in the above code must be set to 0.8. Finally, it is negative because of the orientation I used: the flow is directed from large to small x.

For other variables than pressure, the boundary condition on the disk patches is simply cyclic, as displayed in the following:

disk_in 
{ 
    type cyclic; 
} 
disk_out 
{ 
    type cyclic;
}

After setting the boundary conditions, I used a – turbulence model with wall functions, so that a very fine grid around the duct is not required. Further settings are pretty standard and are taken from the tutorials that are shipped with openFoam.

Note, however, that the grid downstream of the duct should be finer in order to really capture the flow phenomena occuring there. At this moment, the (circular) wake of the duct is quickly smoothened out across the coarse grid cells downstream of the duct. However, this is just a test case and me practicing with OpenFOAM.

Download the case files

I hope this post helped for those needing a ducted actuator disk in OpenFOAM. The files can be downloaded through this link: duct_case.tar