CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

3D volume force term to axisymmetric!

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

Reply
 
LinkBack Thread Tools Display Modes
Old   February 1, 2017, 18:30
Default 3D volume force term to axisymmetric!
  #1
Member
 
OpenFoam
Join Date: Jun 2016
Posts: 79
Rep Power: 3
CFD-Lover is on a distinguished road
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;

Code:
    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];
        }
    }
Code:
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//    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);
    
}
Code:
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);
Code:
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//    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,
Attached Files
File Type: zip actuatorDiskExplicitForce.cpp.zip (3.1 KB, 3 views)
CFD-Lover is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
multiphaseEulerFoam high Courant number Frenk_T OpenFOAM 5 November 24, 2016 04:23
sHM layer process keeps getting killed MBttR OpenFOAM Native Meshers: snappyHexMesh and Others 4 August 15, 2016 03:21
how to set periodic boundary conditions Ganesh FLUENT 13 January 22, 2014 05:11
Force can not converge colopolo CFX 13 October 4, 2011 22:03
Modelling a force as a momentum source term - HELP Rick FLUENT 0 March 9, 2006 11:28


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