This document describes the steps needed to implement a simple body force propeller model in OpenFOAM. Body forces
are a way to represent objects interfering with a fluid without discretizing said object. Consider a body moving in a fluid.
This body will disturb the surrounding fluid by transferring some of its momentum through viscous and normal stresses.
This can be approximated by introducing a forcing term 
in the Navier-Stokes momentum equation as shown in
Eq. (1). A body force model of the moving object is a formulation for the distribution of 
in the fluid; such that the
generated momentum is accurately represented.
 | (1) |
In this case, let’s consider the momentum introduced by a rotating propeller. Approximating this using a body force model presents a major advantage in ship performance calculations. Discretizing a complex rotating geometry behind an already complex geometry (ship’s hull + appendages) is both challenging and computationally intensive.
There are many approaches to describing a propeller using a distribution of 
. For example, the open source self
propulsion framework SHORTCUt has several formulations implemented, based on Blade Element Momentum theory and
Lifting Line/Surface theory. See Windén (2021a,b) for more information. In this report, a very simple model will be
created to demonstrate the principle of how to modify OpenFOAM solvers for this kind of purpose. These are the goals
that will be addressed in the following sections:
First, let’s create a directory where you will keep the source code of the new solver.
Now, let’s get a copy of simleFoam and call it my_simpleFoam. Note the path to the root installation of OpenFOAM; if yours is different then adjust accordingly.
Remove derivative solvers we don’t need.
Finally, before we start modifying the source code, let’s modify the make files. We need to change the name and location of the target executable to avoid overwriting anything.
Next, we need to have the solver create/read the body force field. In createFields.H, add the following just under the
intitialization of 
.
The boundary conditions of the volumeForce field are defined by adding a file called volumeForce in the 0-directory in any case that is to use the my_simpleFoam solver. The volumeForce file should contain the following:
which sets the boundary condition of 
to zeroGradient for all boundaries.
my_simpleFoam.C is the source code for the solver. We need to add an expression for the volume force here. Add the following above #include "Ueqn.H"
Also in my_simpleFoam.C, add the following below #include "initContinuityErrs.H"
Finally, we need to add the body force field to the Navier Stokes momentum equation. In Ueqn.H, Add the volume force to the right hand side of the equation, modifying the definition of tmp<fvVectorMatrix> tUEqn as:
In Section 4, initProp.H was added to the simpleFoam source code to initialize the propeller model. Here, we define how the solver should read propeller properties from a dictionary. Let’s define a new dictionary, it should be called propellerDict and be placed in the system-folder of a case where you would like to use my_simpleFoam. In this example, the propeller thrust and torque will be obtained by curve fit to an experimental open water curve. The dictionary therefore contains 5th order polynomial coefficients that describes the relationship between advance ratio and thrust/torque for a MARIN 7967 propeller. The dictionary can be modified with different coefficients to fit any other available open water curves.
propellerDict should contain the following information:
Now we need to set up the solver to read this dictionary. Create a new file called initProp.H in the same directory as my_simpleFoam.C and initialize the propeller model as follows:
Let ![[U ]
xpp](Eqstyles/inline_eqn6.png)
be the total velocity 
at a probe location 
, a distance 
upstream of the propeller origin
at a radius of 
. Also, let the propeller orientation be defined by the propeller disk normal vector 
and the
propeller vertical direction be defined by the vector 
. Then
 | (2) |
and the propeller inflow velocity 
can be calculated as
![U0 = [U]xpp ⋅P](Eqstyles/eqn3.png) | (3) |
From this, an approximate propeller advance ratio 
can be calculated as
 | (4) |
where 
is the propeller rotation rate and 
is the propeller radius. Using a 5th order polynomial fit to
experimental open water curves, the coefficient of thrust 
and torque 
can be calculated from this advance ratio
as
 | (5) |
 | (6) |
Here, the propeller origin 
, the propeller orientation 
, the distance to the probe point 
, the propeller
radius 
, the rotation rate 
as well as the 5th order coefficients 
and 
are all read from the propeller
dictionary at run-time (as shown in the source code for initProp.H) and do not need to be hard-coded into the
solver.
From Eq. (5) and Eq. (6), the total thrust 
and torque 
of the propeller are given as
 | (7) |
 | (8) |
This methodology yields an approximation of the total forcing to apply. It does not however, model the distribution of
the force on the propeller disk. To do this, and thereby obtaining the distribution of 
on the Finite
Volume mesh, a Goldstein (1929) optimal distribution can be used. A non-dimensional radius 
is defined
as
 | (9) |
where 
is the propeller hub radius and 
is a vector from the propeller origin to the cell center of cell I,
projected onto the propeller plane. Two shape functions, 
and 
, describe the distribution of thrust and torque
respectively along the blade radius. These are defined, based on 
, as:
 | (10) |
 | (11) |
The distribution of 
on the Finite Volume mesh can be obtained by applying the shape functions 
and 
to
each cell of the mesh, with 
being applied outside of the propeller radius and inside the hub radius. The shape
functions are based on the cell center location from Eq. (9).
Finally, to ensure consistensy with the calculated thrust and torque from Eq. (7) and Eq. (8), the distribution of 
is first normalized by the sum of the forces over all cells inside the propeller disk as:
 | (12) |
 | (13) |
where 
is the volume of cell 
. Finally, the strength of 
in cell 
, denoted 
, is calculated by
multiplying the normalized force distribution by the thrust, torque and the appropriate directional vectors.
 | (14) |
This approach is very crude since only one velocity is probed and no consideration is given to the blade geometry
(Goldstein optimum used instead.) This means that the method is unable to detect variations in both the wake and
propeller geometry and therefore, is not suitable for studying propeller-hull interaction. It is provided here only as a simple
example of how a method like this can be implemented. The normalization done in Eq. (14) is consistent with what would
be done for more advanced models where 
, 
, 
and 
are calculated from the propeller geometry rather than
a fixed distribution.
In Section 4, forceEqn.H was added to the simpleFoam source code to calculate the distribution of 
. Create
a new file called forceEqn.H in the same directory as my_simpleFoam.C and calculate the body force as
follows:
After following the steps in the previous sections, to compile the solver run the following command in the my_simpleFoam-directory
After successful compilation, my_simpleFoam is now available as a solver to use in any OpenFOAM case. Bear in mind
that the case must contain the propellerDict dictionary, in the form shown in Section 5 and the 
boundary conditions,
in the form shown in Section 3 for the solver to work.
This section shows a very simple usage example for the my_simpleFoam solver that was created in the previous sections.
The propeller is modeled inside a cylindrical domain with a uniform inflow speed of 
. A shaft extends from
the inlet to the propeller position to mimic the setup of a cavitation tunnel. The dimensions of the domain are shown in
Fig. 1. The propeller radius is 
and the rotation rate is 
giving an advance ratio of 
. The
OpenFOAM case used here is attached to this document.
This example serves only to illustrate that the solver serves its purpose of accelerating the flow and does not go into detail regarding the flow-field accuracy with regards to the target propeller. For more advanced propeller models, such an analysis would be necessary to determine their fitness for purpose. Because the inflow velocity in this solver is probed at a single point upstream, it is unsuitable for studying propeller/hull interaction. It would however be possible to use it for limited propeller/rudder interaction studies by placing a rudder downstream.
The pressure distribution on the tunnel center plane is shown in Fig. 2. The pressure jump stemming from the propeller
blades can clearly be seen. Note that the pressure jump is not explicitly modeled or enforced but is a result of the
pressure-velocity coupling and the applied momentum source (
, see Eq. (1)). Three-dimensional streamlines are also
shown to illustrate how the flow is rotated at the propeller plane.
on the center plane of the tunnel. Streamlines (in 3D) passing through the
propeller disk are show in order to illustrate how the flow is rotated.The axial velocity on the tunnel center plane is shown in Fig. 3, with a detailed view near the propeller in Fig. 4. Three-dimensional streamlines are shown in these figures to illustrate how the flow is rotated.
on the center plane of the tunnel. Streamlines (in 3D) passing through
the propeller disk are show in order to illustrate how the flow is rotated.
on the center plane of the tunnel.
Streamlines (in 3D) passing through the propeller disk are show in order to illustrate how the flow is rotated.
Goldstein, S. (1929). On the vortex theory of screw propellers. Proc. R. Soc. Lond., 123.
Windén, B. (2021a). An open-source framework for ship performance cfd. In Proceedings of the 26th SNAME Offshore Symposium.
Windén, B. (2021b). Predicting the powering performance of different vessel types using an open-source cfd propulsion framework. In Proceedings of the SNAME Maritime Convention 2021.
