CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Propeller dynamics/body force model (http://www.cfd-online.com/Forums/openfoam-solving/71705-propeller-dynamics-body-force-model.html)

SD@TUB January 12, 2010 13:06

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?)

linnemann January 12, 2010 13:21

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.

linnemann January 13, 2010 04:36

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.

SD@TUB January 13, 2010 06:07

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

egp January 21, 2010 07:01

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

SD@TUB January 21, 2010 16:52

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

vinz January 22, 2010 08:24

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

SD@TUB February 3, 2010 07:23

@Eric
Quote:

// Reduce values for parallel execution
reduce(numInside,sumOp<label>());
reduce(numCells,sumOp<label>());
reduce(totalVolume,sumOp<scalar>());
What is the reason for reducing these values in case for parallel execution?
Are the equations solved in simpleFoam cell centre based?


@Vincent
Quote:

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
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)?

vinz February 3, 2010 10:41

Quote:

Which grid generaor do you use and what is the structure of the extracted cellZone file? Is 'funkySetFields' a standard function included in OF?
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

egp February 4, 2010 07:42

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


Quote:

Originally Posted by vinz (Post 243446)
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


vinz February 4, 2010 07:55

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

egp February 4, 2010 08:05

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

vinz February 4, 2010 08:08

Ok I see. Thanks for your reply.

Vincent

SD@TUB February 25, 2010 18:19

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
Code:

cellCenter[?]
(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
Code:

fm.first().first()
or
Code:

fm.first().second()
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

egp February 26, 2010 07:10

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)

SD@TUB February 26, 2010 18:22

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

tstovall March 15, 2010 12:50

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

vinz March 15, 2010 16:39

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

egp March 15, 2010 17:02

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.

tstovall March 16, 2010 12:26

1 Attachment(s)
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


All times are GMT -4. The time now is 14:49.