# 3D volume force term to axisymmetric!

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

February 1, 2017, 18:30
3D volume force term to axisymmetric!
#1
Member

OpenFoam
Join Date: Jun 2016
Posts: 79
Rep Power: 3
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";
}

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(debug >= 3) {
Info << "Point: " << i << " is in the actuator disk. Coordinates: " << iMesh.C()[i] << "\n";
}

ioVolumeForce[i] += axialForce;

//compute the total force added to the actuator disk, this is just for control
TotalForce += axialForce*iMesh.V()[i];

ioVolumeForce[i] += circForce;

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();

oCircumferentialDirection = VecLineTangent ^ VecLineToPoint;
oCircumferentialDirection /= mag(oCircumferentialDirection);

//Check if the point is inside the actuator disk in the radial direction

}```
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
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;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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
 actuatorDiskExplicitForce.cpp.zip (3.1 KB, 3 views)

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Frenk_T OpenFOAM 5 November 24, 2016 04:23 MBttR OpenFOAM Native Meshers: snappyHexMesh and Others 4 August 15, 2016 03:21 Ganesh FLUENT 13 January 22, 2014 05:11 colopolo CFX 13 October 4, 2011 22:03 Rick FLUENT 0 March 9, 2006 11:28

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