[Naruto]

Hello,
I am currently using OpenFOAM for 3-D simulation of Vertical Axis Wind Turbine. Actually I moved on to 3-D simulation after 2-D simulation. My current problem is I am getting too large negative Cp. Infact sometimes it may go beyond 1. Here is the result of my CheckMesh utility:
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.2.2 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.2.2-9739c53ec43f
Exec : checkMesh
Date : Jan 06 2014
Time : 23:47:53
Host : "hpcmaster01"
PID : 17376
Case : /home/Turbine_BUET_0_exp.
nProcs : 1
sigFpe : Floating point exception trapping - not supported on this platform
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 1001176
faces: 6019344
internal faces: 5930454
cells: 2665263
faces per cell: 4.48353
boundary patches: 10
point zones: 0
face zones: 0
cell zones: 2

Overall number of cells of each type:
hexahedra: 603260
prisms: 0
wedges: 0
pyramids: 82226
tet wedges: 0
tetrahedra: 1979777
polyhedra: 0

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
*Number of regions: 2
The mesh has multiple regions which are not connected by any face.
<<Writing region information to "0/cellToRegion"

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
blade-2 31703 30935 ok (closed singly connected)
Rotor 5144 2692 ok (non-closed singly connected)
blade-1 31703 30935 ok (closed singly connected)
Rotor-AMI 5782 5900 ok (non-closed singly connected)
topAndBottom 2678 1447 ok (non-closed singly connected)
Stator 5144 2692 ok (non-closed singly connected)
inlet 171 200 ok (non-closed singly connected)
sides 612 700 ok (non-closed singly connected)
outlet 171 200 ok (non-closed singly connected)
Stator-AMI 5782 5900 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-15 -10 -2.5) (25 10 3)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (1.48176e-18 7.67635e-17 -1.5866e-18) OK.
Max cell openness = 3.39725e-16 OK.
Max aspect ratio = 61.3245 OK.
Minimum face area = 2.18162e-08. Maximum face area = 1.57259. Face area magnitudes OK.
Min volume = 1.75721e-12. Max volume = 0.535229. Total volume = 4400. Cell volumes OK.
Mesh non-orthogonality Max: 70.075 average: 19.2599
*Number of severely non-orthogonal (> 70 degrees) faces: 2.
Non-orthogonality check OK.
<<Writing 2 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 1.66154 OK.
Coupled point location match (average 0) OK.

Mesh OK.

End

As you may see from the checkMesh utility the mesh is fine. But there is only one thing that may be a reason for abnormality and that is:
"The mesh has multiple regions which are not connected by any face."

What does it mean? I got the same message while 2-D simulation too. But in case of 2-D I got good results. I have modelled my whole problem after the propeller tutorial. In my opinion the problem may arise from the fvSchemes or fvSolutions too. I am posting it here.

FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) cellLimited Gauss linear 1;
}

divSchemes
{
default none;
// div(phi,U) Gauss upwind;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,k) Gauss upwind;
div(phi,omega) Gauss upwind;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
pcorr ;
p ;
}


//
solvers
{
pcorr
{
solver GAMG;
tolerance 1e-7;
relTol 0.005;
smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration off;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 20;
mergeLevels 1;
maxIter 100;
minIter 1;
}
p
{
$pcorr;
tolerance 1e-5;
relTol 0.01;
}

pFinal
{
$p;
tolerance 1e-6;
relTol 0;
}

"(U|k|epsilon|omega|nuTilda)"
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 1;
tolerance 1e-07;
relTol 0.05;
maxIter 100;
minIter 1;
};
"(U|k|epsilon|omega|nuTilda)Final"
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.05;
maxIter 100;
minIter 1;
}
}

PISO
{
nCorrectors 2; // PISO loop
// not used in turbSimpleDyMFoam
nOuterCorrectors 2; // PIMPLE loop
nNonOrthogonalCorrectors 0;
correctPhi true;
}

PIMPLE
{
nOuterCorrectors 2;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
correctPhi true;

}

SIMPLE
{
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}

relaxationFactors
{

p 0.3;
U 0.7;
k 0.6;
omega 0.6;
epsilon 0.6;

}


I am at lost at the moment. I would really appreciate any help from you.

[Naruto]

For your convenience I am attaching some pictures. It is the picture of the overall computational domain:
Capture.jpg

These are the AMI PATCHes:
Capture2.JPG

Thanks

[waku2005]

Hi Naruto

I also simulate micro VAWT.
Are your simulations successful for stopped or very slow rotation cases ?
In addition, plz. check in your constant/controlDict for axis definition of moment and air density etc.

Best regards,
waku2005

[Naruto]

Thanks for your reply. Actually I am simulating a high speed turbine. The Reynolds number is quiet high. But I am getting negative power co-efficient. Although I was able to solve the problem a little. Previously I used to get power co-eff around -0.5.

Now it is around -0.2. I achieved this change by changing fvSolution.
solvers
{
pcorr
{
solver GAMG;
tolerance 1e-7;
relTol 0.005;
smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration off;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
maxIter 100;
minIter 1;
}
p
{
$pcorr;
tolerance 1e-5;
relTol 0.01;
}

pFinal
{
$p;
tolerance 1e-6;
relTol 0;
}

"(U|k|epsilon|omega|nuTilda)"
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 1;
tolerance 1e-07;
relTol 0.05;
maxIter 100;
minIter 1;
};
"(U|k|epsilon|omega|nuTilda)Final"
{
solver PBiCG;
preconditioner DILU;
nSweeps 2;
tolerance 1e-07;
relTol 0.05;
maxIter 100;
minIter 1;
}
}

PISO
{
nCorrectors 2; // PISO loop
// not used in turbSimpleDyMFoam
nOuterCorrectors 2; // PIMPLE loop
nNonOrthogonalCorrectors 0;
correctPhi true;
}

PIMPLE
{
nOuterCorrectors 2;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
correctPhi no;

}



relaxationFactors
{

p 0.3;
U 0.7;
k 0.6;
omega 0.6;
epsilon 0.6;

}

cache
{
grad(U);
}
It is my changed fvSolution. My rotational axis definition is ok. Could you suggest me anything else?

[waku2005]

Hi Naruto

Well, I think fvSolution is not the reason of negative Cp.

I'd like to know your TSR and free stream velocity.
If you set too high TSR, it may bring the negative Cp.
How about set a relatively slow TSR such as around 1.

Best regards.
waku2005

[Naruto]

Hi,
Thanks for your interaction. My freestream velocity is 7m/s. The TSR is only 3.0. I really could not figure out the problem. For your conveinence I am attaching my case files. If you could manage time please take a look. I would be very glad
https://dl.dropboxusercontent.com/u/...T.S.R.-3.0.rar

[waku2005]

Hi Naruto

Thanks for your case files.
I made a run of the case without any changes on my desktop (OF 2.2.x/pimpleDyMFoam),
and got Cp and Torque variations until t=0.6 as attached figures.

Cp is relatively small and negative in some period,
but I think that is due to its wing section and solidity( NACA0024 ?) ,
and these results are seemed to be reasonable.
# Cp will vary and some averaging process would be required.

For the above Cp calculation, I used　below values.

• V(free stream)=7.0m/s
• TSR=3.0
• A(area)=H*D=0.1m*0.784m
• omega=-54.1401134 rad/s (=517rpm)
• rho(air density)=1.25 kg/m3

Best regards,
waku2005

[Naruto]

Dear Wakku,

I think, I have partially solved the problem. Actually I am trying to validate with an experimental work. According to the experiment, the average power co-efficient is 0.15.

Anyways, actually my boundary layer meshing was incorrect. After increasing the boundary layer thickness a little bit, I got a little better result. Although, It is still not sufficient. So far I had been able to get an average of 0.08.
That's horrible.

I have found one very interesting thing though. If you use upwind scheme for div(phi,omega) and div(phi,k), you would get better result. From my very limited knowledge of CFD, I am guessing that at the current configuration the flow is fully turbulent. That is why it is more stable with upwind scheme instead of linear upwind one which is more suitable in viscous sublayer.

I have two questions. Like for setting the value of k I use the following formula k=1.5(U*I)*(UI) where I is the turbulent intesity. As I am using SST k-omega model for omega my equation is omega=sqrt(k)/((turbulent length scale)*Cmu^0.25). I got this value after studying the code of OpenFOAM. Am I correct? I am normally choosing a turbulence intensity of 1% and turbulent length scale of like 0.25 at the moment.

Could you suggest anything else to improve my result? Thanks.

[waku2005]

Hi naruto

I think turb. length of 0.25 is somewhat long. How about 0.1 or 0.05?
And if you use more refined boundary layers, please consider to use komegaSST_LowRe for low Re effect.

Sincerery yours
waku2005

[Naruto]

Dear wakku,

I heard about low Re k-omega SST model. But my current RE_no is aroudn 1.1e05. Would this turbulence model be efficient for this high Re number?
And I will use length of 0.05 and 0.1.
Thank you

[waku2005]

Hi Naruto

RASmodels such as komegaSST_lowRe
includes low Re number effects near wall.
Those models are not for low Re cases, and also appricable for hi Re cases.

It might improve numerical results
but you will require more finer boundary mesh
and change wall treatments of wings.
I hope this will help you.

Best regards
waku2005

[Naruto]

Dear wakku,

Thank you for your suggestion. Could you please give me a little hint on how to give boundary condition using kklOmega model?

Thank you.

[waku2005]

Hi Naruto

Sorry, I've never used kklOmega model for low Re RANS simulation.
Please refer the below thread:
http://www.cfd-online.com/Forums/ope...el-low-re.html

Regards
waku2005