[SD@TUB]

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]

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]

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]

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]

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

[SD@TUB]

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]

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]

@Eric

// 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

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]

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]

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

Originally Posted by vinz

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]

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]

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]

Ok I see. Thanks for your reply.

Vincent

[SD@TUB]

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

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

fm.first().first()

or

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]

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]

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]

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]

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]

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]

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