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 coordinates. 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 OF1.6 on a multicore node. (Perhaps other OF versions are needed?) 
Hi Stefan
look here for starters http://www.cfdonline.com/Forums/ope...ebcof15.html Our company has paid for development of a more advanced model which uses nonuniform pressure and angle distribution based on the radial coordinates. This will be a part of the 1.5dev 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.5dev edition when Hrv has time to include it. The thing it lacks is the nonuniform wakefield as we did not need it, but you can surely base you work on this and then add some form of wakefield 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 nonuniform 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 BSeries)! 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 timeconsuming 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. Firstorder approach for loading is to uniformly distribution the total thrust and torque over the volume. If you are doing selfpropulsion simulations, you will need the openwater curves of Kt and Kq vs. J so that you can specify loading as a function of speed. Secondorder 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 bladeelement 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 propellerinabox 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
Quote:
Are the equations solved in simpleFoam cell centre based? @Vincent Quote:
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)? 
Quote:
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 Quote:

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 stairstepping 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 Code:
cellCenter[?] 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 nonuniform 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() Code:
fm.first().second() /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. CFDSHIPIOWA 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 nondimensionalized. 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. 
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.cfdonline.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 10:20. 