Propeller dynamics/body force model
 

Hello Folks,

i'm investigate some elementary dynamics of submerged cylindrical bodies including propulsion
units, i.e. ductet propellers. Therefore i have to approximate the propeller dynamics
using a body force model. A 'real' simulation (describing the propulsion unit in detail)
using moving meshes isn't intended for this thesis? Are there any experiences doing that
with OpenFOAM?
In my estimation i need to integrate a 'black box' in my meshed domain, within a dynamic
boundary condition have to implemented!? The 'black box' got to be an arc disk with a
diameter according to the propeller one (momentum disk). Thus i'll have a prop inlet/outlet
and a disk surface certainly.
On the prop inlet i need to specify the pressure and velocity (wake) vectors in cartesian
co-ordinates. With these parameters i'm able to calculate the dynamics using i.e. a simple
approximation like the 'actuator disc theory' or any further. On the prop outlet i need to
commit the new calculated parameters to the simulation.
Has anybody helpful suggestions (or references) to perform such procedures usually?

Thx,
/Stefan


PS:
I'm using OF-1.6 on a multicore node. (Perhaps other OF versions are needed?)

Hi Stefan

look here for starters http://www.cfd-online.com/Forums/ope...e-bc-of15.html

Our company has paid for development of a more advanced model which uses non-uniform pressure and angle distribution based on the radial coordinates.

This will be a part of the 1.5-dev version when we decide to release it as Open Source, don't ask me when this will happen as it's not up to me.

Hi again

Just talked it over with my company and it will released as a part of the 1.5-dev edition when Hrv has time to include it.

The thing it lacks is the non-uniform wake-field as we did not need it, but you can surely base you work on this and then add some form of wake-field code to it.

Some tutorials of the use will also be available when it is released.

This is a interesting information! I exactly need the non-uniform pressure and velocity distribution based on the radial coordinates of the disk area.
Currently i'm looking for a principal way of implenting a body force model.
Which model i choose isn't decided yet. A suitable model would be a standard propeller series to be approximated (i.e. Wageningen B-Series)!
My idea was to calculate the propeller dynamics in a uniform stream condition to validate the body force model with given propeller characteristics (K_T, K_Q, eta). Following that, i want to join the momentum disc with the hull to investigate the powering performance. Just like the standard proceeding on naval architecture.
I'm sure that will be a time-consuming work, but i see no alternatives!

/Stefan

Stefan,

There are two approaches to introducing momentum sources:

1. explicit introduction of a pressure jump as a boundary condition, such as fanFvPatchField.

2. addition of a localized momentum source (propeller) or sink (turbine) as a body force.

Personally, I prefer the latter. For complex geometries, I'd rather not have the constraint of gridding a disc in the middle of an already complex problem. Instead, we specify (in constant/bodyForceDict) the swept volume of the propeller (hub radius, tip radius, chord), and then search for all of the cells which lay within this volume. This is very general and easily permits modeling of multiple propellers, and even control surfaces. The final step is to decide on a loading distribution.

First-order approach for loading is to uniformly distribution the total thrust and torque over the volume. If you are doing self-propulsion simulations, you will need the open-water curves of Kt and Kq vs. J so that you can specify loading as a function of speed.

Second-order approach would be to use something like the Hough and Ordway circulation distribution which unloads the root and tips.

Want more? Then couple to either a blade-element model or a panel code. Just depends on your objectives.

I should note that in the process of developing this capability, it is a good idea to use basic momentum analysis to verify your implementation. You can derive the velocity jump for a specified body force, and then grid up a simple propeller-in-a-box test case as shown in the image below http://dl.dropbox.com/u/2182201/Scre...51.33%20AM.png

To modify simpleFoam, do the following:

1. create the file createBodyForce.H. For the simple case of a single prescribed and uniform body force, it might look something like:

Info<< "Creating body force" << endl;

// Allocate the body force vector field and initialize it to zero
volVectorField bodyForce
(
IOobject
(
"bodyForce",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
mesh,
dimensionedVector("zero",
dimForce/dimVolume/dimDensity,
vector::zero)
);

// Read a dictionary to configure the body force
IOdictionary bodyForceDict
(
IOobject
(
"bodyForceDict",
runTime.time().constant(),
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);

// Define the (cylindrical) region for the body force from dictionary values
const scalar xo (readScalar(bodyForceDict.lookup("xo")));
const scalar yo (readScalar(bodyForceDict.lookup("yo")));
const scalar zmin (readScalar(bodyForceDict.lookup("zmin")));
const scalar zmax (readScalar(bodyForceDict.lookup("zmax")));
const scalar tipRadius (readScalar(bodyForceDict.lookup("tipRadius")));
const scalar hubRadius (readScalar(bodyForceDict.lookup("hubRadius")));

// Define the body force magnitude and direction
vector unitVector(bodyForceDict.lookup("directionVector") );
unitVector /= mag(unitVector);
const scalar magnitude(readScalar(bodyForceDict.lookup("magnitu de")));

label numCells(0),numInside(0);
scalar totalVolume(0.0);

forAll(mesh.cells(),cellI)
{
// Loop over all the cells in the mesh
++numCells;

// Check containment of each cell center
const vector& cellCenter = mesh.C()[cellI];
const scalar dx = cellCenter[0] - xo;
const scalar dy = cellCenter[1] - yo;
const vector r(dx,dy,0);

if ((cellCenter[2] < zmax) &&
(cellCenter[2] > zmin) &&
(mag(r) < tipRadius) &&
(mag(r) > hubRadius)) {
++numInside;
totalVolume += mesh.V()[cellI];
bodyForce[cellI] = unitVector;
}
}

// Reduce values for parallel execution
reduce(numInside,sumOp<label>());
reduce(numCells,sumOp<label>());
reduce(totalVolume,sumOp<scalar>());

Info << "total volume = " << totalVolume << endl;
Info << numInside << " of " << numCells << " cells are inside" << endl;

// Set the body force magnitude
if (totalVolume > 0) {
bodyForce *= magnitude/totalVolume;
}

// Write the resulting field to a file
bodyForce.write();


2. Add the bodyforce to UEqn.H

tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
- bodyForce
);

3. make sure to include "createBodyForce.H" in simpleFoam.C


Good luck

These are quite helpful lines. I'm still new of implenting code in solver sources. I'll take up my employment next days (probably for weeks, indeed) and report about the progress.
Your suggested way of proceeding seems to be rational, thanks for that!

/Stefan

Hi everybody,

Very interesting your post Eric, I chose the same approach to model an actuator disk inside simpleFoam, and it works very well. However I don't use the part where you look for the cells which are inside the "propeler" volume, I get those cells directly from my grid generator, create a cellZone file, then I use funkySetFields to impose a specific value to those cells to a variable force that I add the same way you do in the equation.

Each selection method has its advantages I think. For example, with your method, how do you select cells if you want to model a disk (propeler, turbine...) which is not aligned with one of the directions? If you only use the method you described here with Zmin and Zmax, you select too many cells, right? So I guess you are doing something else?

Regards,

Vincent

@Eric
What is the reason for reducing these values in case for parallel execution?
Are the equations solved in simpleFoam cell centre based?


@Vincent
Which grid generaor do you use and what is the structure of the extracted cellZone file? Is 'funkySetFields' a standard function included in OF?
I use 'snappyHexMesh' for meshing my domain, so i'll need to specify my designated cells located at propeller position in your case, afterwards. Which utility is suggested for that (i.e setFields, cellSet)?

Hi Stephan,

I use GridPro to generate my grids for the actuator disk. The cellZone file is simply a list of the cells adjacent to the disk. FunkySetFields is not directly included in the OF distribution but search for the term on the forum or on google you will find the wiki page dealing with it where you can download the sources and find some tutorials.
Hope this helps.

Vincent

Hi Vincent,

I agree that the posted code snippet is not general. Note the wording in my post, "...it might look something like..." I wanted to give a simple example.

Of course, this can be generalized for any direction and loading distribution.

Eric

Hi Eric,

Yes, don't worry I understood what you meant.
But I am still wondering, don't you get a kind of "ladder" effect when you impose a value to cells which are crossed by your "virtual" disk? Or you use a very fine mesh so you can simply say that this effect is negligeable?
Or how do you avoid that?
Would you happend to have an image of what I am talking about?

Thanks for the information.

Vincent

Vincent,

If the grid is very poor, then this stair-stepping effect can be seen. However, in practice, the grid is typically clustered near the outer radius of the turbomachine (e.g., due to geometry or meshing technique). Therefore, it is not a big deal.

Eric

Ok I see. Thanks for your reply.

Vincent

Hello again,

last days i did some work on the body force modell. First i have to say, that the posted lines by Eric above still work with any major issue. Now i'm facing some tricky things and got some questions?

For what stands the number in
(cartesian co-ordinate directions ?)

I'm looking for a viable possibility of indicating the collected cells due to e.g. visualize them in paraFoam. I'm reading some details of cellSet via dictionary. Is this wokable for that?

Further issues are reading some data of interest. In the same way of marking' cells i want to read the velocity vectors in each of the collected cells. Same routine i'll use to implement a non-uniform thrust loading. Therefor i've got a function of thrust distribution. (BTW: The work of 'Hough and Ordway' seems to be a fundamental investigation of approximate propellers. Does anybody hold this paper, 'cause it is not provided online!?)

Additionally as input value for the iteration, i have to put the net thrust according to the resistance of investigated body. In the source of force function (forces.C) the summarized pressure and viscous forces are stored by or What would be the best way of reading the forces every time step. When i put the fm.* in my createBodyForce.H i get the message, that fm ist not declared.

/Stefan

Hough and Ordway paper: http://dl.dropbox.com/u/2182201/1029317.pdf
Nice to have for historical reasons, but it take some work to read and then create a model for body forces.

CFDSHIP-IOWA report:
http://www.iihr.uiowa.edu/~shiphydro...Report_TOC.pdf
http://www.iihr.uiowa.edu/~shiphydro...3_Modeling.pdf
Good example if you want to see the equations for implementation (equations 36 and 37 of the second pdf)

Thank you for the papers!

Three things that were coming in my mind reviewing the equations (36) and (37) (knowing these are technical questions, not a problem of OF):
The functions of distribution in eq. (36) seems to be non-dimensionalized. But with division by DXPROP in eq. (37) there will be a dimension [1/m] left!?
The terms of A_x and A_phi are functions of T or C_T respectively Q or K_Q. Isn't it a conflict, if the results towards integration in eq. (38) over distribution volume are coupled with T an Q (by rho*L^2*U^2 and rho*L^2*U^3) as well?
The maximum of f_bphi(r/R) in eq. (36) is near to the hub in opposite to the known maximum of f_bx or T with r/R=0.7!? I guessed equal positions (r/R) of maxima for T and Q.

Back to OF:
In your screenshot you visualize the included cells collected by the loop. Could you give me a hint how to do that? I'm not sure if cellSet works for that, cause the cellSetDict do not provide cylindrical boundings.
Additionally i want to use such construct to read out velocities (i.e. input value for iteration and to produce isotach lines to expose typical wake fractions afterwards).

Thanks in advance!

/Stefan

Eric,

I've added a force to the pisoFoam solver based on the above procedure. I'm getting some oscillations in the flow velocity through the center of the turbine, not at the edges of the turbine. I've refined the mesh 4 levels and am still getting oscillations. I've experimented with the differencing schemes for the div(phi,U) term in the fvSchemes files. The Gauss upwind scheme seems to work the best, while any linear scheme produces oscillations that extend from the disc to the inlet of the domain. What schemes are you using in your test simulation? From the picture above, it doesn't look like you have any of these oscillations.

Thanks,

Tim

Hi,

I don't know what Eric will say, but for me it is normal that you make the oscillations Disappear when you refine the mesh, it should be even worse actually. The coarsest your mesh, the more dissipation you get which is damping your solution. So what I usually do, I start with a coarse mesh, converge on this one and then refine, converge, and refine again.
Normally I use three levels of refinement, coarse medium fine.

For the schemes, I noticed the same behavior than you did. The upwind scheme is a first order scheme from what I understood in the user guide. So it is easier to get convergence with it. So with the same idea than before, I do the same, I start with the upwind scheme to get a converged solution and then switch to another linear scheme which is more precise.

Hope this helps.

Regards,

Vincent

For the test problem, I used

div(phi,U) Gauss limitedLinearV 0.5;

It looked pretty smooth, so I'm not sure why you are seeing oscillations. Did you use

div(phi,U) Gauss linear;

If so, this would definitely lead to oscillations in the flow field for convection dominated flows.

I have been using div(phi,U) upwind; but I switched and tried div(phi,U) LimitedlinearV 0.5; and got the same oscillations. I've included a contour and line plot of the flow. Are you able to get runs without this level of oscillations?

http://www.cfd-online.com/Forums/att...1&d=1268756621

Would you be willing to post your test case so I can verify if the issue is in my solver implimentation or case setup?

Thanks,

Tim

Hello everybody,
i'm interesting to body force analysis.
What is the dimension "magnitude" in the code?
Is it possible add the rotational term to the thrust? Maybe could be correct using MRF + body force to the same volume of cells?

Thanks for any reply

Andrea

Hello Andrea,

the dimension of magnitude is L^4/T^2 (force per density) as you can guess from:
I added tangential forces on every involved cell to approximate the rotational effect of the propeller. What you intend with adding a 'rotational term to the thrust'?

~Stefan

Hi Stefan,
thank you for you reply!
This is what I mean!
In the code posted dy Eric there is just thrust, is it?
I changed t to put a propeller in whatever direction I want, can you help me to add also tangential force to the involved cells?

Thanks

Andrea

Hello Andrea,

maybe i could...

The code above is an example by Eric to implement a uniform force distribution in direction of a selectable directionVector. For testing you could also calculate a uniform tangential force across the propeller radius with a given magnitude and an tangential vector which point for every cell in tangential direction of the local position vector (from propeller centre to each cell centre). It is a little bit tricky, but with some careful considerations you will calculate the appropriate direction with basic trigonemtrical functions.
The local position vector to cellI ist described by r and phi according to the local point of origin (propeller centre). In case of a directionVector = (1 0 0), the tangential vector, say tan' schould be something like this:
tan'(x_i,z_i)=(cos[phi],sin[phi])^T with phi=arctan(y_i/z_i)

It is a good help to plot your developed functions and cecking the rightness. Lately in paraView you see, if the tangential vectors are correct.
With a factor, say f you are able to define the rotation direction of the propeller:
tanrot'=f*tan'(x_i,z_i) with f={-1,1}

x_i,z_i are the global cell centre positions, if the global point of origin is (0 0 0) and y_i=0 (one central screw propeller).
If the above things will working, you can determine the ratio from axial and tangential forces by a propeller nomogram (e.g. standard series). Additionaly a radial distribution of the magnitudes could be implement (see PNA Vol.II for qualitative radial circulation distribution of an open propeller).

I hope that helps a bit!

~Stefan

Broken links for the literature provided by Eric
 

Hi Eric
Your post dated :

February 26, 2010,

have link to some articles. These links seem to dead now. If possible please
provide with the new links to these articles.

Thanks a lot
jaswi

Hi Tim,

I think your problem is the same that mine and the wiggles happen "always"
http://www.cfd-online.com/Forums/ope...l-wiggles.html

"Discrete body forces are used in
the present context to model the influence of wind turbines on
the flow. In
order to overcome the pressure wiggles introduced by discre
te body forces,
one approach is to smooth out the body forces by using a Gaussi
an distri-
bution instead of a Dirac delta distribution "
]
Using rotorDiskSource the problem is that wind speed is not calculated correctly in each cell so the angle of attack is incorrect, the force calculation, the wake....
Other option to solve the problem is to modify the Rhie-Chow algorithm.....

Any help?

Hi,

I'm new to OpenFoam and I'm wondering if there is a way to implement the localised momentum sources in OpenFoam 2.3.0, is the code any different to the OF 1.6 Version. Also, how would you create the BodyForceDict file?

Many Thanks

Matt

Dear all, please tell me how to solve this.
 

Hello everyone.

What I am currently trying to do is to study the flow around the hull taking into account the rotation of the propeller.

I think that calculating the propeller and hull as a single unit is very costly and unrealistic.

Therefore, we are considering alternative methods. The methods we are currently considering are

1. First, we calculate the propeller alone in a uniform flow.
2. Next, together with the hull model, we will place a simple propeller model (disk shape?) in place of the complex propeller model. Instead of the complex propeller model, place a simple propeller model (disk shape?) in its original location.
3. Based on the result of step 1, a volume force is added to the simple propeller model (disk shape?). I would like to add the volume force to the simple propeller model (disk shape?) and calculate around the hull.

Currently, I have completed the calculations in step 1 using both AMI and MRF methods, and I am wondering at what point to add the volume force.
What methods exist for approaching this in openfoam?

I'm using OpenFOAM-4.x.


I feel that this is a difficult problem that I cannot solve without your help.
Thank you very much for your help.

Hi,

its been a while that I implemented body forces to approximate a real propeller. In modern OF versions one choice would be using a function object. Or you may use the native API rotorDiskSource, see
https://www.openfoam.com/documentati...iskSource.html

There should be an example available online from Chalmers university on how to use it for a ship propeller.

Hope that points into the right direction.


/Stefan

Hi Stefan.

Thanks for your reply.
First of all, I'm a newbie when it comes to dealing with body forces, so I apologize for any inconvenience this may cause.

There are a few things I would like to confirm.
What has been discussed by stefan and others is
Are you saying that you have rewritten the code so that the acutuator-disk generates the general volume force that the propeller would generate without the pre-calculated results of the propeller?
Is this a little different from my way of doing things, which is to pre-calculate the propeller only and then use the result to give the volume force to the actuator-disk?

I would appreciate it if you could enlighten me in my ignorance.

By using a function object, do you mean that there is a way to give the result of the calculation as a volume force to the specified cell zone without rewriting the code?

I looked at some papers from Chalmers University, but I found that it requires a lot of rewriting of the code, which is very scary for me as a beginner in rewriting code.
If I have to rewrite the code, I am prepared to do so, but I am an inexperienced person who would like to avoid it if possible.

I would appreciate it if you could help me with my question.

Hi again,

No problem! Such a task takes time.
Unfortunately, I am not up-to-date what latest options you might have with a most recent version of OF, but I guess a little coding is - at least - still required.

Before the technique of 'function objects' was introduced in OF the only option was to adapt a given solver i.e. pimpleFoam to add source terms for body forces and to define where these forces shall apply.
With 'function objects' you can use a given OF solver and add own features as a function object that is compiled during run-time.

Using a cell set to pre-select cells where body forces shall be active is the right way.
The complexity of the propeller model depends on your needs. I used a weekly coupling model by determing the advance speed and J at the propeller plane from a resistance case and then used the KT at this J of a selected propeller to get the thrust to be used as a body force. This is maybe similar to your way.


/Stefan