[antar]

Dear Foamers

Greetings!

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)

where

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:
actuationDiskSourceCoeffs
{
   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::
addRadialActuationDiskAxialInertialResistance
(
    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_[0]
      + 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] =
            T
           *(radialCoeffs_[0] + 
              radialCoeffs_[1]*r2 + 
              radialCoeffs_[2]*sqr(r2))
            /intCoeffs;

        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

Thanks
Antar

[lebc]

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,
Luis

[SGWANG]

hello, guys

Is there someone still working on this?

I think it should be

PNG image 2019-10-25 10_01_34.jpg

[vahidmazidi]

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