CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Propeller dynamics/body force model

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   January 12, 2010, 13:06
Default Propeller dynamics/body force model
  #1
Member
 
Stefan
Join Date: Jan 2010
Location: Kiel, Germany
Posts: 72
Rep Power: 7
SD@TUB is on a distinguished road
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?)
SD@TUB is offline   Reply With Quote

Old   January 12, 2010, 13:21
Default
  #2
Senior Member
 
linnemann's Avatar
 
Niels Nielsen
Join Date: Mar 2009
Location: NJ - Denmark
Posts: 445
Rep Power: 14
linnemann will become famous soon enough
Hi Stefan

look here for starters Fan type BC in OF15

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

PS. I do not do personal support, so please post in the forums.
linnemann is offline   Reply With Quote

Old   January 13, 2010, 04:36
Default
  #3
Senior Member
 
linnemann's Avatar
 
Niels Nielsen
Join Date: Mar 2009
Location: NJ - Denmark
Posts: 445
Rep Power: 14
linnemann will become famous soon enough
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.
__________________
Linnemann

PS. I do not do personal support, so please post in the forums.
linnemann is offline   Reply With Quote

Old   January 13, 2010, 06:07
Default
  #4
Member
 
Stefan
Join Date: Jan 2010
Location: Kiel, Germany
Posts: 72
Rep Power: 7
SD@TUB is on a distinguished road
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
SD@TUB is offline   Reply With Quote

Old   January 21, 2010, 07:01
Default
  #5
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
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

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
egp is offline   Reply With Quote

Old   January 21, 2010, 16:52
Default
  #6
Member
 
Stefan
Join Date: Jan 2010
Location: Kiel, Germany
Posts: 72
Rep Power: 7
SD@TUB is on a distinguished road
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
SD@TUB is offline   Reply With Quote

Old   January 22, 2010, 08:24
Default
  #7
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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 is offline   Reply With Quote

Old   February 3, 2010, 07:23
Default
  #8
Member
 
Stefan
Join Date: Jan 2010
Location: Kiel, Germany
Posts: 72
Rep Power: 7
SD@TUB is on a distinguished road
@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)?
SD@TUB is offline   Reply With Quote

Old   February 3, 2010, 10:41
Default
  #9
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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
vinz is offline   Reply With Quote

Old   February 4, 2010, 07:42
Default
  #10
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
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 View Post
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
egp is offline   Reply With Quote

Old   February 4, 2010, 07:55
Default
  #11
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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
vinz is offline   Reply With Quote

Old   February 4, 2010, 08:05
Default
  #12
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
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
egp is offline   Reply With Quote

Old   February 4, 2010, 08:08
Default
  #13
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Ok I see. Thanks for your reply.

Vincent
vinz is offline   Reply With Quote

Old   February 25, 2010, 18:19
Default
  #14
Member
 
Stefan
Join Date: Jan 2010
Location: Kiel, Germany
Posts: 72
Rep Power: 7
SD@TUB is on a distinguished road
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
SD@TUB is offline   Reply With Quote

Old   February 26, 2010, 07:10
Default
  #15
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
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)
egp is offline   Reply With Quote

Old   February 26, 2010, 18:22
Default
  #16
Member
 
Stefan
Join Date: Jan 2010
Location: Kiel, Germany
Posts: 72
Rep Power: 7
SD@TUB is on a distinguished road
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

Last edited by SD@TUB; February 28, 2010 at 13:27.
SD@TUB is offline   Reply With Quote

Old   March 15, 2010, 12:50
Default
  #17
New Member
 
Tim Stovall
Join Date: Mar 2009
Posts: 12
Rep Power: 8
tstovall is on a distinguished road
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
tstovall is offline   Reply With Quote

Old   March 15, 2010, 16:39
Default
  #18
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
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
vinz is offline   Reply With Quote

Old   March 15, 2010, 17:02
Default
  #19
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
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.
egp is offline   Reply With Quote

Old   March 16, 2010, 12:26
Default
  #20
New Member
 
Tim Stovall
Join Date: Mar 2009
Posts: 12
Rep Power: 8
tstovall is on a distinguished road
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?



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
Attached Images
File Type: jpg T428R03wake.jpg (72.6 KB, 1692 views)
tstovall is offline   Reply With Quote

Reply

Tags
momentum disk, propeller, propulsion

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Two-phase air water flow problems by activating Wall Lubrication Force challenger85 CFX 5 November 5, 2009 06:44
Drag force in massless particles in discrete phase model payam_IUST FLUENT 2 October 18, 2009 23:25
DEFINE_CG_MOTION and pressure force Teo Fumagalli FLUENT 0 April 11, 2008 10:25
DPM model w/ Wave model - errors in documentation HS FLUENT 0 April 12, 2006 04:37
Propeller model (CFX-5.7.1) Jesper CFX 6 April 18, 2005 19:31


All times are GMT -4. The time now is 19:33.