[CFD-Lover]

Hi all,

I have been reading the following code about introducing a volume source term in the simpleFoam to solve the flow over HAWT. Snipped of the code is as follows;

    if(debug >= 1) {
        Info << "Calculating volume force from actuator disk.\n";
    }

    ReadGeometry(iMesh);

    scalar RadialDist2;
    vector LineTangent;
    vector CircumferentialDirection;
    
    vector TotalForce(0.0,0.0,0.0);
    scalar TotalTorque = 0.0;
    
    scalar DiskVolume = 0;

    //Loop over all cells and check if the cell center is in the actuator disk region
    for(label i = 0; i < iMesh.C().size(); i++) {

        if(PointIsInDisk(mPointStartCenterLine,mPointEndCenterLine,iMesh.C()[i],RadialDist2, LineTangent,CircumferentialDirection)) {

            if(debug >= 3) {
                Info << "Point: " << i << " is in the actuator disk. Coordinates: " << iMesh.C()[i] << "\n";
            }
            
            vector axialForce = LineTangent*CalcAxialForce(sqrt(RadialDist2),mRho)/mRho;
            ioVolumeForce[i] += axialForce;
            
            //compute the total force added to the actuator disk, this is just for control
            TotalForce += axialForce*iMesh.V()[i];


            vector circForce = CircumferentialDirection*CalcCircForce(sqrt(RadialDist2),mRho)/mRho;
            ioVolumeForce[i] += circForce;

            TotalTorque += (CalcCircForce(sqrt(RadialDist2),mRho)/mRho)*sqrt(RadialDist2)*iMesh.V()[i];
            DiskVolume += iMesh.V()[i];
        }
    }

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//    Check if a given point is located in the actuator disk region.
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    vector VecLineTangent(iPointEndCenterLine - iPointStartCenterLine);
    scalar LineTangentLength = sqrt(VecLineTangent.x()*VecLineTangent.x() + VecLineTangent.y()*VecLineTangent.y() + VecLineTangent.z()*VecLineTangent.z());

    if(LineTangentLength != 0.0) {
        VecLineTangent /= LineTangentLength;
    }
    else {
        Info << "Warning: The centerline tangent has zero length.\n";
        return false;
    }

    oLineTangent = VecLineTangent;
    
    vector VecStartLineToPoint(iPoint - iPointStartCenterLine);
    scalar PointProjOnLine = VecStartLineToPoint & VecLineTangent;

    //Check if the point is inside the actuator disk in the axial direction
    if(!(PointProjOnLine >= 0.0 && PointProjOnLine <= LineTangentLength)) {
        return false;
    }
    
    vector VecLineToPoint(VecStartLineToPoint - (VecLineTangent*PointProjOnLine));
    scalar RadialDist2 = VecLineToPoint.x()*VecLineToPoint.x() + VecLineToPoint.y()*VecLineToPoint.y() + VecLineToPoint.z()*VecLineToPoint.z();
    oDist2 = RadialDist2;
    
    oCircumferentialDirection = VecLineTangent ^ VecLineToPoint;
    oCircumferentialDirection /= mag(oCircumferentialDirection);

    //Check if the point is inside the actuator disk in the radial direction
    return(RadialDist2 <= mExtRadius*mExtRadius && RadialDist2 >= mIntRadius*mIntRadius);
    
}

vector VecLineTangent(iPointEndCenterLine - iPointStartCenterLine);
    scalar LineTangentLength = sqrt(VecLineTangent.x()*VecLineTangent.x() + VecLineTangent.y()*VecLineTangent.y() + VecLineTangent.z()*VecLineTangent.z());

    if(LineTangentLength != 0.0) {
        VecLineTangent /= LineTangentLength;
    }
    else {
        Info << "Warning: The centerline tangent has zero length.\n";
        return false;
    }

    vector VecStartLineToPoint(iPoint - iPointStartCenterLine);
    scalar PointProjOnLine = VecStartLineToPoint & VecLineTangent;

    //Check if the point is inside the actuator disk in the axial direction
    if(!(PointProjOnLine >= 0.0 && PointProjOnLine <= LineTangentLength)) {
        return false;
    }
    
    vector VecLineToPoint(VecStartLineToPoint - (VecLineTangent*PointProjOnLine));
    scalar RadialDist2 = VecLineToPoint.x()*VecLineToPoint.x() + VecLineToPoint.y()*VecLineToPoint.y() + VecLineToPoint.z()*VecLineToPoint.z();
    
    //Check if the point is inside the actuator disk in the radial direction
    return(RadialDist2 < mIntRadius*mIntRadius);

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//    Compute the force component in the axial direction. The force is computed from a simple equation resulting in a force
//    that varies with the radial distance.
//    If you have a better model of a rotor, comment the four lines below and add your own calculation of the axial force.
//    Do not forget to also change the calculation of the tangential force (CalcCircForce()) below.
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    scalar axialForce = 0.0;
    scalar radiusScaled = (iRadialDist/mExtRadius - mIntRadius/mExtRadius)/(1.0 - mIntRadius/mExtRadius);
    scalar Ax = (105.0/8.0)*mThrust/(CalcDiskThickness()*mPI*(3.0*mIntRadius+4.0*mExtRadius)*(mExtRadius-mIntRadius));
    axialForce = Ax*radiusScaled*sqrt(1.0 - radiusScaled);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    
    if(debug >= 2) {
        Info << "Axial force: " << axialForce << "\n";
    }
    
    return axialForce;
}

scalar actuatorDiskExplicitForce::CalcCircForce(const scalar &iRadialDist, const scalar &iRho) {
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//    Compute the force component in the tangential direction. The force is computed from a simple equation resulting in a force
//    that varies with the radial distance.
//    Change the four lines below if you have a better model.
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    scalar tangentialForce = 0.0;
    scalar radiusScaled = (iRadialDist/mExtRadius - mIntRadius/mExtRadius)/(1.0 - mIntRadius/mExtRadius);    
    scalar At = (105.0/8.0)*mTorque/(CalcDiskThickness()*mPI*mExtRadius*(mExtRadius-mIntRadius)*(3.0*mExtRadius+4.0*mIntRadius));
    tangentialForce = (At*radiusScaled*sqrt(1.0 - radiusScaled)/(radiusScaled*(1.0 - mIntRadius/mExtRadius) + mIntRadius/mExtRadius));
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    if(debug >= 2) {
        Info << "Tangential force: " << tangentialForce << "\n";
    }
    
    return tangentialForce;
}

The code is written to solve forces only in 3D but would like to rewrite to solve forces as axisymmetric. I don't really know what to change in the upper snipped code. The C++ file is also attached.

Many thanks,