# Formulation of the Radial Acutation Disk Source

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

 September 30, 2013, 09:48 Formulation of the Radial Acutation Disk Source #1 New Member   Antar Netra Join Date: Sep 2013 Location: India Posts: 2 Rep Power: 0 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 : Code: 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: Code: template 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 zoneCellCentres(mesh().cellCentres(), cells); const Field 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()); reduce(upRho, minOp()); 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