CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Heat Transfer from a Rough Cylinder in Tunnel RE=2.2E5 M=0.07 (

aerothermal December 6, 2010 08:56

Heat Transfer from a Rough Cylinder in Tunnel RE=2.2E5 M=0.07
Dear All,

I am running in rhoSimpleFoam (OpenFoam 1.7.x), realizableKE, rough wall functions, a cylinder (D=0.15 m and length=0.5 m) with roughness (equivalent sand-grain height, Ks=0.00135 m) in a tunnel of cross section (0.9 by 0.5 meters) operating with flow of Re=2.2E5 and inlet Tu=0.45%.

some assumptions for simulation:
- upstream tunnel length = 10*D
-downstreal tunnel length = 20*D
-atmospheric back pressure imposition = 101325 Pa
-inlet X velocity = 23.45 m/s
-intlet Tu = 1%
-inlet mut/mu = 50
-first prysm layer at 1e-3 m
-last prysm layer at 6.645e-2 m
-number of layers 26
-prysm growth rate 1.05
-Nominal Ks = 0.00135 m
-despite the laminar-turbulent transition is located at 1.2 degree (found by integral method and Rek=600 criteria for nominal Ks), OpenFoam does not have a laminar-transition yet available.
-inlet air temperature = 303.15 K
-surface temperature DeltaT = 6 K
-mutRoughWallFunction used for momentum
-mutalphatWallFunction used for heat

It has been hard to get convergence with kEpsilon, realizableKE and even kOmegaSST. I have tried several configurations and followed recommendations of several forums threads.

Sometimes the residuals stall or oscillates, which is common in CFD. As I am interested in forces (by analogy with heat transfer), this may not as critical as it looks. However, I do not see an expected symmetry in the U or P solution (contours) as I got in other softwares like CFD++ or even SolidWorks Flow Simulation. Other issue is that p residuals stall at very high value, well above 1e-2.


    default        steadyState;
    default        Gauss linear;
    grad(U)        cellMDLimited Gauss linear 1;
    grad(p)        cellMDLimited Gauss linear 1;
    default              none;
    div(phi,U)      Gauss linearUpwindV cellMDLimited Gauss linear 1;
    div((muEff*dev2(grad(U).T()))) Gauss linear;
    div(phi,h)      Gauss linearUpwind cellMDLimited Gauss linear 1;
    div(phi,epsilon) Gauss upwind;
    //div(phi,omega)  Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phid,p)      Gauss upwind;
    div(U,p)        Gauss upwind;
    default                        none;
    laplacian(muEff,U)                Gauss linear limited 0.5;
    laplacian(alphaEff,h)        Gauss linear limited 0.5;
    laplacian((rho|A(U)),p)        Gauss linear limited 0.5;
    laplacian((rho*rAU),p)        Gauss linear limited 0.5;
    laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5;
    //laplacian(DomegaEff,omega)        Gauss linear limited 0.5;
    laplacian(DkEff,k)                Gauss linear limited 0.5;
    laplacian(1,p)                  Gauss linear limited 0.5;
    laplacian((rho*(1|A(U))),p)        Gauss linear limited 0.5;
    default                  linear;
    interpolate(U)        linear;
    default        limited 0.5;
    default        no;
    p              ;


        solver          GAMG;
        tolerance      1e-10;
        relTol          0.001;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels    1;
        nSweeps  2;
        nPreSweeps  0;
        nPostSweeps  2;
        solver          PBiCG;
        preconditioner  DILU;
        tolerance      1e-08;
        relTol          0.005;
        solver          smoothSolver;
        smoother        GaussSeidel;
            nSweeps                2;
        tolerance      1e-08;
        relTol          0.005;
        solver          smoothSolver;
        smoother        GaussSeidel;
            nSweeps                2;
        tolerance      1e-08;
        relTol          0.005;;
        solver          smoothSolver;
        smoother        GaussSeidel;
            nSweeps                2;
        tolerance      1e-08;
        relTol          0.005;

    nCorrectors 3;
    nNonOrthogonalCorrectors 5;
    pRefCell 0;
    pRefValue 0;
    pMin            pMin [ 1 -1 -2 0 0 0 0 ] 10000;
    rhoMax          rhoMax [ 1 -3 0 0 0 0 0 ] 4;
    rhoMin          rhoMin [ 1 -3 0 0 0 0 0 ] 0.5;
    p              0.2;
    rho            0.1;
    U              0.5;
    k              0.5;
    epsilon        0.5;
    h              0.2;


application    rhoSimpleFoam;
startFrom    startTime;
startTime      0;
stopAt    endTime;
endTime    2500;
deltaT    1;
writeControl    timeStep;
writeInterval    100;
purgeWrite      0;
writeFormat    binary;
writePrecision  6;
writeCompression uncompressed;
timeFormat      general;
timePrecision  6;
graphFormat    raw;
runTimeModifiable yes;


Mesh stats
    points:          242333
    faces:            1613006
    internal faces:  1584830
    cells:            720237
    boundary patches: 7
    point zones:      0
    face zones:      0
    cell zones:      0

Overall number of cells of each type:
    hexahedra:    0
    prisms:        316888
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    403349
    polyhedra:    0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch              Faces    Points  Surface topology                 
    cylinder            12188    6158    ok (non-closed singly connected) 
    topbottom          6802    4340    ok (non-closed singly connected) 
    topbottom2          6820    4349    ok (non-closed singly connected) 
    side1              964      552      ok (non-closed singly connected) 
    inlet              220      131      ok (non-closed singly connected) 
    outlet              220      131      ok (non-closed singly connected) 
    side2              962      551      ok (non-closed singly connected) 

Checking geometry...
    Overall domain bounding box (-1.47689 -0.45 -0.25) (2.97689 0.45 0.25)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (1.20477e-17 4.90409e-16 7.44393e-18) OK.
    Max cell openness = 3.20177e-16 OK.
    Max aspect ratio = 7.90905 OK.
    Minumum face area = 3.77495e-06. Maximum face area = 0.0165682.  Face area magnitudes OK.
    Min volume = 3.95177e-09. Max volume = 0.000473042.  Total volume = 1.99538.  Cell volumes OK.
    Mesh non-orthogonality Max: 70.0352 average: 22.7072
  *Number of severely non-orthogonal faces: 1.
    Non-orthogonality check OK.
  <<Writing 1 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 1.0217 OK.

Mesh OK.



PS: the main reference is the paper below...
[1] Achenbach, E., The effect of surface roughness on the heat transfer from a circular cylinder to the cross flow of air, Int. J. Heat and Mass Transfer, Vol. 20, 1977, pp. 359-69. paper link

aerothermal December 6, 2010 16:21

boundary conditions
1 Attachment(s)
Dear All,

See attached the boundary conditions folder "0".



aerothermal December 8, 2010 08:49

Cylinder Cp data
1 Attachment(s)
Dear All,

I have not posted yet the mesh but I will do in next few days.
See attached the experimental data for Cp over a rough and smooth cylinder regarding the case above.

There is a correlation for Cp over the rough cylinder fitted in range 2.2E05 < Re_D < 4.0E06 by Stefanini et. al[1] that may help the analysis of the problem:
a <- 2.646
b <- 1.194
cpCorr <- 1-a*sin(b*thetaRad)^2

where 0<thetaRad<pi/2.



[1] Stefanini, L. M., Silvares, O. M., Silva, G. A. L., and Zerbini, E. J. G. J., Boundary , Heat Transfer on Iced Cylinders, AIAA Paper 2010-7672, AIAA Atmospheric and Space Environments Conference, Toronto, Ontario, Aug. 2-5, 2010, American Institute of Aeronautics and Astronautics, Reston, 2010.



aerothermal December 18, 2010 11:12

5 Attachment(s)
Dear All,

Maybe this subject is not so interesting to everybody but I will continue posting on that because I think that some day some people may need information about the issue. So I am describing what I did to get some (very) preliminary convergence despite not so satisfactory for my purposes (heat transfer).

See attached the residuals, continuity, forces and heat transfer along the iterations (time in steadyState). Only the main residuals are outside the zip, for the rest see zip file.

Also I am posting the the "/0" and "/system" folders (see zip). As the "constant" folder is too big (approx 10 MB) to be placed here, I uploaded it to a virtual folder at address below:

It appears that OF goes well until iteration 200, when it starts diverging and oscillating about at higher level of residuals than the level reached in 200 time. This is caused mainly by the wake effect. Despite I could filter most of it, there is some residual oscillation that could not be killed. Is it my wrong impression that OF171 is not truly steady state solver so the wake cannot be averaged or filtered out?

See attached also, the results of the velMag countours and Cp for OF171 compared with those of CFD++ ( The OF171 results are in blue.

Both software use the realizableKE model, assume wall functions to simulate momentum and heat transfer over the rough surface and have the same mesh. However, CFD++ has the laminar-turbulent transition and OF171 does not have since it is full turbulent only. The transition in CFD++ was triggered at 1.5 degree (of the cylinder) that is actually happened in experimental runs.

The mesh was developed in mesh generator MIME (, exported as CGNS by MIME, and converted from CGNS to FOAM by OF1.5-dev tool cgnsToFoam. Due to some errors in conversion from CGNS format to FOAM, I had to autoPatch to find the patch faces of the mesh since the converter only recognized the internal mesh.

As I have the heat transfer and Cf results validated I will post them all.


Guilherme da Silva

aerothermal May 12, 2011 12:41

5 Attachment(s)
Dear All,

I managed to get convergence in OF1.6-ext but not in OF1.7.1 with same configuration files in folders 0/ constant/ system/. I will post them shortly.

In order to make it work, I had to modify the alphatWallFunction to a new alphatRoughWallFunction proposed by our team. I changed also the mutRoughWallFunction to be compliant with models used in icing literature. I hope to post them in a repository soon and make it available to everybody.

Meanwhile, I am posting the results.

For those who want details, we are publishing the paper in next SAE Aircraft and Engine Icing International Conference in Chicago.

aerothermal July 17, 2011 09:34

I will organize data and post here soon.

aerothermal December 20, 2011 09:47

validation of new alphatWallFunction for heat transfer over rough walls
Dear colleagues,

The recent results of validation performed in OpenFoam 1.6-ext and 1.7.1 can be found at our last paper:

Comments and questions are welcome in order to improve the model and the code.



fredo490 May 21, 2013 07:11


Originally Posted by aerothermal (Post 287702)
Dear All,

Maybe this subject is not so interesting to everybody but I will continue posting on that because I think that some day some people may need information about the issue.

Dear Guilherme,
Your sentence perfectly matchs my case and I'm really interestet in your work.
I have read your paper and I have to say that you made a great job. I went to take a look on your website and on the links you give at the end of your paper. As appendix, you say that your code is available on but it seems that the files have been deleted.

May I ask you to share it with me ? Did you try to implement it into OF 2.x ?

Regards, Fred

aerothermal May 22, 2013 11:02

Dear Frederick,

Yes. All files are there but not in binary download section.
The source is in the repository.

To close that you will need Mercurial package. Please google it to check how it works.

To know more about Mercurial:

After installing mercurial, type the command:
hg clone


Guilherme - aerothermal

fredo490 May 23, 2013 04:12

Thank you very much for your answer.
For those who also want to get the source code with Ubuntu:

1) Install Mercurial

sudo apt-get install mercurial
2) Open a new terminal where you want to save the source code and write:

hg clone
After a while, all the source code will be downloaded ;)

All times are GMT -4. The time now is 16:13.