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/)
-   -   Using k-omega SST to simulate simple channel flow (https://www.cfd-online.com/Forums/openfoam-solving/130805-using-k-omega-sst-simulate-simple-channel-flow.html)

MaLa March 4, 2014 13:40

Using k-omega SST to simulate simple channel flow
 
3 Attachment(s)
Hello people,

I have been tring to use the kOmegaSST model implemented in OF 2.2x, to simulate channel flow. I have defined the inlet and outlet as cyclic BCs, and am only evaluating for the half channel width. Moreover, depth in the z direction is one cell thick, so it is practically a 2 D channel. yPlus at the first node from wall ~ 0.75 (128 cells in the y direction).

Following are the boundary conditions for my turbulent parameters :

p

Code:


dimensions      [0 2 -2 0 0 0 0];

internalField  uniform 0;

boundaryField
{
  inlet
    {
      type            fixedValue;
      value          uniform 1.0;

    }

  outlet
    {
      type            fixedValue;
      value          uniform 0.0;
    }

  fixedWalls
    {
      type            zeroGradient;
    }

  frontAndBack
    {
      type            empty;
    }

  middle
    {
      type          symmetryPlane;
    }
}

U

Code:


dimensions      [0 1 -1 0 0 0 0];

internalField  uniform (0 0 0);

boundaryField
{
  inlet
    {
      type            zeroGradient;
        }

  outlet
    {
      type            zeroGradient;
        }

  fixedWalls
    {
      type            fixedValue;
      value          uniform (0 0 0);

    }

  frontAndBack
    {
      type            empty;
    }

  middle
    {
      type            symmetryPlane;
    }
}

k

Code:


dimensions      [0 2 -2 0 0 0 0];

internalField  uniform 0.325;

boundaryField
{
    inlet
      {
        type            zeroGradient;
      }

    outlet
      {
        type            zeroGradient;
      }

    fixedWalls
      {
        type            fixedValue;
        value          uniform 1e-10;
      }

    frontAndBack
      {
        type            empty;
      }
    middle
    {
      type            symmetryPlane;
    }


}

omega

Code:

dimensions      [0 0 -1 0 0 0 0];

internalField  uniform 400;

boundaryField
{
  inlet
    {
      type            zeroGradient;
        }

  outlet
    {
      type            zeroGradient;
        }

  fixedWalls
    {
      type            omegaWF;
      value          uniform 400;

    }

  frontAndBack
    {
      type            empty;
    }
  middle
    {
      type          symmetryPlane;
    }


}

Here, omegaWF is a modified omegaWallFunction which disables the re-evaluation of the production term at the first node next to the wall. Also, omega is calculated only from omegaVis (i.e. omega= omegaVis, in omegaWF) and not from omegaLog.

nut

Code:

dimensions      [0 2 -1 0 0 0 0];

internalField  uniform 0;

boundaryField
{
    inlet
      {
        type            zeroGradient;
      }

    outlet
      {
        type            zeroGradient;
      }

    fixedWalls
      {
        type            nutLowReWallFunction;
        Cmu            0.09;
        kappa          0.41;
        E              9.8;
      }

    frontAndBack
      {
        type            empty;
      }

    middle
      {
        type          symmetryPlane;
      }

}

The Reynolds_tau has been chosen as 395, hence the laminar viscosity (nu) is defined as 1/395 in constant/transportProperties.

When I compared my results with those obtained by CFDTM (http://cfd.mace.manchester.ac.uk/twi...tCase001Res001), my graphs looked as follows:

Does anybody know the reason for this? Has anybody tried validating kOmegaSST in the past? Can somebody looking at the figures and suggest why nut starts dipping at a certain y+?

Would be grateful for any insights. Thanks

ArathoN March 4, 2014 17:44

First why did you change omegaWallFunction? It works both for low-Re and High-Re meshes.
Second nutLowReWallFunction is meant for yplus between 2-5 and not below 1 and for kOmegaSST use nutUSpaldingWallFunction it is a continuous wall function and it workd in high and low re meshes.

MaLa March 5, 2014 05:02

Ok Thanks ArathoN. I will try this and let you know what I get.

I changed the omegaWallFunction to omegaWF because my yPlus for first node was so less. Hence I did not see why omega at the first node should be calculated as a mean of both omegaVis (omega at the first node in viscous sublayer) and omegaLog (omega at the first node in the log layer).

Code:

omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
What do you think about this?

Also, do you know the reference for the re-evaluation of the production term at the first node, in the omegaWallFunction?

Code:

G[faceCellI] =
            (mutw[faceI] + muw[faceI])
          *magGradUw[faceI]
          *Cmu25*sqrt(k[faceCellI])
          /(kappa_*y[faceI]);

In general, for other turbulence models (like say v2f), if the yPlus is below 1, what is the best nutWallFunction to use?

MaLa March 5, 2014 05:09

Hey ArathoN

I tried with your suggestions but the result did not change much! :(

Have you by chance tried running kOmegaSST for channel flow and validating it yourself?

ArathoN March 5, 2014 05:15

Generally if yplus i below 1 you need to set nut as calculated, but for the kOmagaSST it's better to set nut as nutUSpaldingWF (this was specified by the developer who implemented the model).

For omega i think it will check the yplus then it will consider the appropriate value (omegavis or omegalog), for this you better read the thread "compressible k omega SST" there are some insight on how the model was implemented by the user henry.

Edit: I'll check the source code too.

MaLa March 5, 2014 05:20

Do you know the name/tag/webpage of the developer who implemented kOmegaSST?

Would like to have a look at his comments because the kOmegaSST.H file does not specify any wall functions for nut as such! :!

ArathoN March 5, 2014 05:58

here the thread : http://www.cfd-online.com/Forums/ope...komegasst.html

And the developer is henry, if you read all the thread you'll see how he advice to use the nutUspalding as a WF for the komegasst and when i tried it i had better results, with calculated it underestimated the dissipation effects.

By the way now that i'm reading your BCs they are very strange:

First you imposed the pressure at inlet and outlet and you set that the velocity gradient must be zero there so it should be a fully developed flow in both patches, but in this case the variation of the pressure must be congruent with the dissipation ratio of the pressure otherwise the results will be really weird, it's better to set:
inlet: U fixed value / p zerogradient
outlet: U zerogradient / p fixed value 0
with these you'll have the region whre the flow will develop, then the fully developed flw so the length must be appropriate (it depends if the flow is laminar or turbulent).

Second you said you are using the cyclic BCs but i see no cyclyc type in you BCs, in you have connected the inlet with the outlet in blockmeshDict you need to set every variable as type cyclic in the patches inlet and outlet, otherwise OF will not recognize the cyclic condition and will switch to the default settings, like this
Code:

boundaryField
{
    inlet
      {
        type            cyclic;
      }

    outlet
      {
        type            cyclic;
      }

    fixedWalls
      {
        type            nutUSpaldingWallFunction;
      }

    frontAndBack
      {
        type            empty;
      }

    middle
      {
        type          symmetryPlane;
      }

}

Third what solver are you using? remember that simplefoam doesn't work well with the cyclic condition, in this case you better use the mappedpatch condition (this like the cyclic it will map the result of the desired patch to the target patch so u can map the values at the outlet to the inlet) or install the 1.6_extent version where there is the channalFoam solver now incorporated into pimpleFoam.

Forth even the BCs for k and omega are wrong, the problem as you stated is ill conditioned, you didn't give any initial value to these fields, change the inlet condition to fixed value where you compute the values by the known relations of k and omega (see this wiki page http://www.cfd-online.com/Wiki/Turbu...ary_conditions and take the relations with the mixing length or see the code and you'll find their expressions). And the value of Omega at the wall is given by this relation \omega_{wall} = 10 \frac{6 \nu}{\beta_1 (\Delta d_1)^2} where d1 is the distance of the first node to the wall.

MaLa March 5, 2014 14:22

Hi ArathoN

Thanks for your feedback. It made me think about quite a bit of stuff:

I initially set off with exploring channel flow. I assumed that since the flow would be developed, d/dx of all fields other than pressure would be zero. That is why I enforced a fixedValue for pressure at the inlet and outlet to create a pressure gradient, but kept the other quantities as zeroGradient at those patches.

On further thought, I suppose that a cyclic boundary condition is what might be more suited for channel flow. (Especially since one doesn't know before hand, if the flow is developed) However, declaring the inlet and outlet as neighboring cyclic patches does not allow me to have separate values for pressure at the inlet and the outlet.

Are you aware of some strategy to enfore a cyclic condition on ONLY a few fields and not all? In this case, for example, I would like to keep pressure as fixedValue at the inlet and outlet, but other fields as cyclic.


Quote:

Third what solver are you using? remember that simplefoam doesn't work well with the cyclic condition, in this case you better use the mappedpatch condition (this like the cyclic it will map the result of the desired patch to the target patch so u can map the values at the outlet to the inlet) or install the 1.6_extent version where there is the channalFoam solver now incorporated into pimpleFoam.
I am using simpleFoam and was not aware of its problems with cyclic condition. Regarding the usage of mapped Patch, I run into the same problem again: it maps, say the inlet to the outlet patch for all fields, which is not what I would like to have. (I must have a pressure gradient !)

What do you think?

Thanks for the thread!

MaLa March 5, 2014 15:07

Quote:

Fourth, even the BCs for k and omega are wrong, the problem as you stated is ill conditioned, you didn't give any initial value to these fields, change the inlet condition to fixed value where you compute the values by the known relations of k and omega
I have enforced Neumann conditions here (gradient is 0), so the problem is not ill-conditioned. Don't you think so?

But even if I were to declare the patches as cyclic, then also it should be okay, because the solver gets more equations by the relation between the cyclic patches and the number of variables is still equal to the number of equations. I think!

huangxianbei March 8, 2014 09:07

Hello,MaLa
I'm not sure why nu=1/395, for me this number is quiet large because in my case,Ret=194, the nu=1e-6(the fluid is water), so what kind of fluid are you focusing on?

I'm doing channel calculation too, while the results are poor. I use DNS, the process of the flow being quasi-steady turbulent is time consuming:(

MaLa March 10, 2014 02:25

Hello huangxianbei,

I assumed nu =1/395 for an imaginary fluid. I wanted to achieve a Reynolds_tau =395, so by chosing the dimensions of the channel in such a way that u_tau (wall shear velocity) = 1 and y_0.5 (half channel width) = 1, leads to setting nu =1/395. (Re_tau = u_tau * y_0.5 / nu).

Do you think this approach is ok?

huangxianbei March 10, 2014 02:42

Quote:

Originally Posted by MaLa (Post 479070)
Hello huangxianbei,

I assumed nu =1/395 for an imaginary fluid. I wanted to achieve a Reynolds_tau =395, so by chosing the dimensions of the channel in such a way that u_tau (wall shear velocity) = 1 and y_0.5 (half channel width) = 1, leads to setting nu =1/395. (Re_tau = u_tau * y_0.5 / nu).

Do you think this approach is ok?

In theory domain, I think it's ok to do such a imagination. I'm just interested in the fluid whose nu=1/395:),for I haven't got to known this kind of fluid.

By the way, I'd like to ask a question about the distribution of Ux(mean) in the y direction. In my result, the profile is not so smooth (some tiny drops of value somewhere), is your profile looks good?

huangxianbei March 10, 2014 02:50

http://www.cfd-online.com/Forums/mem...cture540-u.png
look at this picture, you can see obviously that the curve is not smooth

ArathoN March 10, 2014 03:19

Mmmmmm what value of the velocity did you set in the inlet? Because you need to correlate properly U_inf with u_tau.
The relation if I remember correctly is U_inf/u_tau=√(2/Cf).
So you have to evaluate the Cf of your case then you can calculate u_tau, particularly OF has an utility to calculate the wallshearstress and from there you'll have u_tau.

huangxianbei March 10, 2014 03:33

Well, I just run the channel395 tutorial to see if everything goes well. There are no inlet and outlet, using periodic in the x direction instead. And a perturbation is generated to make the flow goes turbulent quickly.

When the streamwise velocity oscillating around the Ubar, the flow is supposed to be statical turbulent.(I'm wondering there should be a more accurate method to judge when it reaches this situation, how do you judge this situation in your simulation?)

sivakumar June 19, 2014 04:48

Hi There,
I think you all are more clear about KOmegaSST wall treatment for low Re case.

In my case the y+ is 3 to 5,

So I am planning to use the following setup for KOmegaSST (low Re):

for K:

inlet
{
type fixedValue;
value uniform 0.39;
}

outlet
{
type zeroGradient;
}

top0 //Wall
{
type fixedValue;
value uniform 1e-10;
}

Omega:

inlet
{
type fixedValue;
value uniform 0.99;
}

outlet
{
type zeroGradient;
}

top0 //Wall
{
type omegaWallFunction;
value uniform 0.99; // ???????
}

nut:

inlet
{
type calculated;
value uniform 0;
}

outlet
{
type calculated;
value uniform 0;
}

top0 // Wall
{
type nutUSpladingWallFunction; Or fixedValue?
value uniform 0;
}

Please give me your suggestions,

Thanks,

Sivakumar


All times are GMT -4. The time now is 19:58.