CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Formulation of the Radial Acutation Disk Source

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

LinkBack Thread Tools Search this Thread Display Modes
Old   September 30, 2013, 09:48
Default Formulation of the Radial Acutation Disk Source
New Member
Antar Netra
Join Date: Sep 2013
Location: India
Posts: 2
Rep Power: 0
antar is on a distinguished road
Dear Foamers


I am scratching my head to understand a very small piece of code related to fvOption radialActuationSource.

The function calculates the thrust for a wind turbine and then distributes the thrust over the wind turbine disk using a radial function. In the radialActuationSource.H, the description says the following and I am quoting :

thrust = 2 \rho A U_{o}^2 a (1-a)


a = 1- Cp/Ct 
A = disk area
U_o = upstream velocity

Thrust is distributed by a radial function

Thrust(r) = T (C_0 + C_1 r^2 + C_2 r^4)

Example Usage:
   fieldName       U;          // name of field to apply source
   diskDir         (-1 0 0);   // disk direction
   Cp              0.1;        // power coefficient
   Ct              0.5;        // thrust coefficient
   diskArea        5.0;        // disk area
   coeffs          (0.1 0.5 0.01); // radial distribution coefficients
   upstreamPoint   (0 0 0);    // upstream point
Now the radialActuationSourceTemplates.C with comments regarding the point i do not understand is given next:

template<class RhoFieldType> 
void Foam::fv::radialActuationDiskSource::
    vectorField& Usource,
    const labelList& cells,
    const scalarField& Vcells,
    const RhoFieldType& rho,
    const vectorField& U
) const
    scalar a = 1.0 - Cp_/Ct_;
    scalarField Tr(cells.size());
    const vector uniDiskDir = diskDir_/mag(diskDir_);

    tensor E(tensor::zero);
    E.xx() = uniDiskDir.x();
    E.yy() = uniDiskDir.y();
    E.zz() = uniDiskDir.z();

    const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
    const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);

    const vector avgCentre = gSum(zoneCellVolumes*zoneCellCentres)/V();
    const scalar maxR = gMax(mag(zoneCellCentres - avgCentre));

   // Question : how do we get this expression. by the variable name 
   //                it seems that we are integrating the radial function
   //                 (C_0 + C_1 r^2 + C_2 r^4). 
   //                Is that so ?
    scalar intCoeffs =
      + radialCoeffs_[1]*sqr(maxR)/2.0
      + radialCoeffs_[2]*pow4(maxR)/3.0;

    vector upU = vector(VGREAT, VGREAT, VGREAT);
    scalar upRho = VGREAT;
    if (upstreamCellId_ != -1)
        upU =  U[upstreamCellId_];
        upRho = rho[upstreamCellId_];
    reduce(upU, minOp<vector>());
    reduce(upRho, minOp<scalar>());

    scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1.0 - a);
    forAll(cells, i)
        scalar r2 = magSqr(mesh().cellCentres()[cells[i]] - avgCentre);

        // Question 2 : Here we are distributing the thrust as per the 
        //                  radial function. Can somebody please explain
        //                  how does the expression below ensures that
        //                  the total thrust T = sum of (Ti)s when this 
        //                  approach is used.

        Tr[i] =
           *(radialCoeffs_[0] + 
              radialCoeffs_[1]*r2 + 

        Usource[cells[i]] += ((Vcells[cells[i]]/V_)*Tr[i]*E) & upU;

    if (debug)
        Info<< "Source name: " << name() << nl
            << "Average centre: " << avgCentre << nl
            << "Maximum radius: " << maxR << endl;
Any help is highly appreciated

antar is offline   Reply With Quote

Old   March 4, 2017, 15:27
Luis Eduardo
Join Date: Jan 2011
Posts: 85
Rep Power: 15
lebc is on a distinguished road
Hi Antar,

Reopening an old post...

I'm also struggling to find out what is done within radialActuationDiskSource, did you find out the answer to the questions you posed?

For Question 1, the expression
scalar intCoeffs = radialCoeffs_[0] + radialCoeffs_[1]*sqr(maxR)/2.0 + radialCoeffs_[2]*pow4(maxR)/3.0
seems to be the integration of
radialCoeffs_[1]*maxR + radialCoeffs_[2]*4*pow3(maxR)/3.0
and for me it doesn't make sense...

Also, maybe the intCoeffs is used to assure that the total Thrust is applied, but I couldn't find a mathematical proof of it.

Did you keep working with it?

Best Regards,
lebc is offline   Reply With Quote

Old   October 24, 2019, 07:40
New Member
Shuguang Wang
Join Date: Feb 2018
Posts: 13
Rep Power: 8
SGWANG is on a distinguished road
hello, guys

Is there someone still working on this?

I think it should be

PNG image 2019-10-25 10_01_34.jpg

Last edited by SGWANG; October 25, 2019 at 06:07.
SGWANG is offline   Reply With Quote

Old   November 3, 2019, 02:20
New Member
Join Date: Oct 2015
Posts: 1
Rep Power: 0
vahidmazidi is on a distinguished road
this comment is very very goog

thank you ...

but i have a question :
Where did this formula(Tr=fx=A*(c0+c1*r^2+c2*r^4)) come from?
vahidmazidi is offline   Reply With Quote


Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
GPU Linear Solvers for OpenFOAM gocarts OpenFOAM Announcements from Other Sources 37 August 17, 2022 14:22
wmake compiling new solver mksca OpenFOAM Programming & Development 14 June 22, 2018 06:29
blocking an annular disk with circumferential and radial channels user0314 Main CFD Forum 0 June 9, 2011 11:10
Version 15 on Mac OS X gschaider OpenFOAM Installation 113 December 2, 2009 10:23
no enthalpy change across the momentum source Atit Koonsrisuk CFX 2 December 19, 2005 02:33

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