CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Big problem predicting adiabatic wall temperature (wing, mach 0.33, rhoSimplecFoam) (https://www.cfd-online.com/Forums/openfoam-solving/116578-big-problem-predicting-adiabatic-wall-temperature-wing-mach-0-33-rhosimplecfoam.html)

fredo490 April 22, 2013 13:30
Big problem predicting adiabatic wall temperature (wing, mach 0.33, rhoSimplecFoam)
 
Dear all,
I'm working on rhoSimplecFoam / rhoLTSPimpleFoam (OF 2.2) to simulate wing at mach 0.33. The simulations run smoothly and give good results (and convergence) except for the wall temperature ! All the other variables are pretty good, the pressure distribution is acceptable, ...

I run the case with k-omega SST, a structured mesh of y+ = 0.1, and I tried a large range of settings.

Every time, when I check the results, I get a wall temperature varying of 10°C while theoretical results and Fluent simulations give a variation of about 2°C maximum. I am 100% sure that the Fluent result is accurate with experimental results and other software. The Fluent simulation was made with a mesh of y+ 0.2 , k-omega SST (compressibility effects + viscous heating), sutherland viscosity and kinetic-theory for the thermal conductivity.

Here is a rough plot to compare the wall temperature of the 2 cases:
http://www.fredo490.fr/public/temperature-compare.png

When I plot the temperature field into the domain, they are very comparable:
http://www.fredo490.fr/public/temper...re-general.png


However a tiny detail is missing. There is a kind of "hot" boundary layer in Fluent that flows all around the airfoil. It turns out that this characteristic is missing in my OF result. Here is a plot of a region near the upper leading edge (but it is also true everywhere else).
http://www.fredo490.fr/public/temper...pare-zooml.png


To make the things simple, here is my full case including mesh (9MB):
http://www.fredo490.fr/public/Naca00...malProblem.zip


Do you have any idea of where the problem is located ? Is it the viscous heating in the boundary layer ?

fredo490 April 23, 2013 04:36
I have ran extra test with Fluent and it turns out that the difference come from the Viscous Heating. If I uncheck the viscous heating box in k-omega sst, the Fluent results look similar to my OpenFoam data.

Now the question is: How to implement the Viscous Heating in rhoSimplecFoam ?!

fredo490 April 23, 2013 10:07
It seems that the answer to my problem is in this topic: http://www.cfd-online.com/Forums/ope...-enthalpy.html

I'm currently running a case with the viscous heating to try the code.

fredo490 May 21, 2013 10:41
2 Attachment(s)
I come back to my topic to give the solution of my problem.

To make things simple, OpenFoam neglects lots of terms and we need to edit the source code of rhoSimplecFoam and rhoPimplecFoam. Enclosed is the source of my equation (Jasaks book's). The code now include: the traditional energy equation with the Viscous Heating + Pressure Work + Kinetic Energy.

For OpenFoam 2.2, you need to edit EEqn.H to get:
Code:
    // Viscous Heating Source terms
        volTensorField gradU = fvc::grad(U);
        volTensorField kin1 = turbulence->muEff() * (gradU + gradU.T());
        volVectorField kin4 = kin1 & U;               
        volScalarField kin5 = fvc::div(kin4);       
               
        volScalarField divU = fvc::div(U);
        volVectorField kin2 = 2.0/3.0 *turbulence->muEff() * divU * U;
        volScalarField kin3 = fvc::div(kin2);


    volScalarField& he = thermo.he();

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K) // Kinetic Energy
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt // Pressure Work
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
    ==
        fvOptions(rho, he)
      + kin5
      - kin3
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    thermo.correct();
For rhoSimplecFoam, the terms are exactly the same so edit the EEqn.H and adapt it to add kin3 and kin5 to the equation.

This new solver doesn't request any special things except the integration of the new schemes. The gradient scheme can be linear and the solution is stable with divergence schemes like upwind and linearUpwind.

To get an accurate result of the viscous heating, you need to have a y+ < 1 (because it is generated by the gradient of velocity in the sub-layer). In my case, I have used a y+ of 0.2 of average and the data looks very good.

I have put enclosed a picture of the heated viscous layer obtained by OpenFoam.

Tushar@cfd May 22, 2013 00:03
Quote:
Originally Posted by fredo490 (Post 428926)
I come back to my topic to give the solution of my problem.

To make things simple, OpenFoam neglects lots of terms and we need to edit the source code of rhoSimplecFoam and rhoPimplecFoam. Enclosed is the source of my equation (Jasaks book's). The code now include: the traditional energy equation with the Viscous Heating + Pressure Work + Kinetic Energy.

For OpenFoam 2.2, you need to edit EEqn.H to get:
Code:
    // Viscous Heating Source terms
        volTensorField gradU = fvc::grad(U);
        volTensorField kin1 = turbulence->muEff() * (gradU + gradU.T());
        volVectorField kin4 = kin1 & U;       
        volScalarField kin5 = fvc::div(kin4);   
       
        volScalarField divU = fvc::div(U);
        volVectorField kin2 = 2.0/3.0 *turbulence->muEff() * divU * U;
        volScalarField kin3 = fvc::div(kin2);


    volScalarField& he = thermo.he();

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K) // Kinetic Energy
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt // Pressure Work
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
    ==
        fvOptions(rho, he)
      + kin5
      - kin3
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    thermo.correct();
For rhoSimplecFoam, the terms are exactly the same so edit the EEqn.H and adapt it to add kin3 and kin5 to the equation.

This new solver doesn't request any special things except the integration of the new schemes. The gradient scheme can be linear and the solution is stable with divergence schemes like upwind and linearUpwind.

To get an accurate result of the viscous heating, you need to have a y+ < 1 (because it is generated by the gradient of velocity in the sub-layer). In my case, I have used a y+ of 0.2 of average and the data looks very good.

I have put enclosed a picture of the heated viscous layer obtained by OpenFoam.

Nice work.. :)

fredo490 May 22, 2013 00:15
Thx but I have to say that the topic of Chrisi1984 help me a lot.

fredo490 May 22, 2013 00:28
I forgot to say : the heat prediction highly depends of the free steam turbulence parameters! A high turbulence intensity (let say 1% such as a wind tunne) will not give the same rresults as a low intensity (0,01% such as a high flight condition). You need to get good k and omega values. So I highly recommend to use the tool given by this website to estimate the two values.

AMitgM95 May 12, 2022 02:24
Hello,

Could you please help in writing the viscous dissipation term in OpenFOAM v9? I tried but was unable to figure out how to write the expression for shear-stress symmetric tensor ((gradU + gradU.T())).


All times are GMT -4. The time now is 04:17.