Turbulent demos

Asked by Marcos Samudio on 2012-04-17

I have my laminar NS solver running with SUPG and PSPG stabilizations, and to test it I used the turbulent channel problem, which is included in the CBC package. My idea is to couple this with the turbulent k-epsilon solver later on.

The turbulent channel is driven by a constant pressure gradient (applied as body force) equal to 0.05**2, and a viscosity of 1/7900. With this data, in a closed channel of width 2, the velocity profile should be V_x = 9.875 - 9.875*y**2.

The laminar solution reproduces this analytical solution very accurately, but I have noticed that the StandarKE implementation, which is in turbulen_channel.py does not.

Is there something missing in the turbulent model or am I doing something wrong?

Thanks,

Marcos

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
Last reply:

This question was reopened

Hi Marcos,

Perhaps I'm missing something, but the turbulent solution should not look
anything like the laminar. What does your solution look like?

Mikael

On Tuesday, 17 April 2012, Marcos Samudio wrote:

> New question #193953 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/193953
>
> I have my laminar NS solver running with SUPG and PSPG stabilizations, and
> to test it I used the turbulent channel problem, which is included in the
> CBC package. My idea is to couple this with the turbulent k-epsilon solver
> later on.
>
> The turbulent channel is driven by a constant pressure gradient (applied
> as body force) equal to 0.05**2, and a viscosity of 1/7900. With this data,
> in a closed channel of width 2, the velocity profile should be V_x = 9.875
> - 9.875*y**2.
>
> The laminar solution reproduces this analytical solution very accurately,
> but I have noticed that the StandarKE implementation, which is in
> turbulen_channel.py does not.
>
> Is there something missing in the turbulent model or am I doing something
> wrong?
>
> Thanks,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

Marcos Samudio (marcossamudio) said : #2

Hi Mikael,

I modified the original NS solver and coupled it with StandardKE_Coupled. The turbulent solution is the same you get when running the original turbulent_channel demo (the largest velocity is 0.698). The laminar solution yields the parabolic profile I already mentioned.

How do you validate your turbulent model then? I was thinking I should get the same profiles for some reason. Should I get some particular velocity deflect law or something similar? Now I'm confused!!

Thanks,

Marcos

Den Apr 17, 2012 kl. 8:25 PM skrev Marcos Samudio:

> Question #193953 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/193953
>
> Marcos Samudio posted a new comment:
> Hi Mikael,
>
> I modified the original NS solver and coupled it with
> StandardKE_Coupled. The turbulent solution is the same you get when
> running the original turbulent_channel demo (the largest velocity is
> 0.698). The laminar solution yields the parabolic profile I already
> mentioned.
>
> How do you validate your turbulent model then? I was thinking I should
> get the same profiles for some reason. Should I get some particular
> velocity deflect law or something similar? Now I'm confused!!

You probably have some reading up to do on turbulence modeling:-) The solution close to the wall should be u = utau and then there's a region where the log-law applies (http://en.wikipedia.org/wiki/Law_of_the_wall). Since the flow is turbulent there's no analytical solution and you have to compare your results with direct numerical simulations (3D transient and chaotic) of the Navier-Stokes equations. There's really a lot of data out there.

Best of luck, don't hesitate to ask (I have a PhD in turbulent reactive mixing).

Mikael

>
> Thanks,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Marcos Samudio (marcossamudio) said : #4

Fun! More reading! haha I have read some of Pope's book but I think I should refresh some concepts!

Thanks a lot for your help Mikael!

Marcos Samudio (marcossamudio) said : #5

Ok, I have realized how inappropriate my question was! The turbulent flow is actually simulating all small scales energy loss, so the profiles must indeed be very different.... sorry for that.

Anyway, I am comparing the calculated velocity profiles with the universal laws, and the results are not as good as one would expect. I have plotted the linear variation law (u_plus = y_plus) and the log-law (u_plus = ln(y_plus)/k +B) against the results of the simulations using the implemented solver (both the modified solver and the original one yield the same result). This is the result: https://docs.google.com/open?id=0B9qsTlUZxTpdWGtOY09fSVcyMXM with u_tau = 0.05 and Re_tau = 750.

The slope is actually very good, and I think that using a wall-function may be the fix for this problem, do you think so?

I started thinking of a way to build this wall-function in 3D. My first attempt would be setting u in the nodes closer to the boundary to the following:

normal = grad(y)/||grad(y)|| where y is the distance to the wall
u = (1/||u||)(u - dot(u,normal)*normal)*u_wall where u_wall depends on y(wall-functions)

This way u would satisfy the wall function and would is some sense conserve streamwise direction.

To implement this, I would use a MeshFunction containing vertices_inside_boundary, mark them with some bid and then do

L = (1/FacetArea(mesh))*(1/||u||)(u - dot(u,normal)*normal)*u_wall*v*ds(bid)

and then assemble, to get a vector V_wf, right?.

The next step would be writing this values in the just-calculated vector, and this is where I have my doubts... I have seen the question https://answers.launchpad.net/dolfin/+question/95871

Sparsity patterns might be different in this case, since I am only integrating over bid. I could modify this by adding, perhaps, -DOLFIN_EPS*v*ds, and then do V = V +V_wf. Is there a way to add these vectors directly?

Well, I hope I was clear in my explanations, and thanks a lot for your time Mikael!

Marcos

Den Apr 18, 2012 kl. 11:40 PM skrev Marcos Samudio:

> Question #193953 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/193953
>
> Status: Solved => Open
>
> Marcos Samudio is still having a problem:
> Ok, I have realized how inappropriate my question was! The turbulent
> flow is actually simulating all small scales energy loss, so the
> profiles must indeed be very different.... sorry for that.
>
> Anyway, I am comparing the calculated velocity profiles with the
> universal laws, and the results are not as good as one would expect. I
> have plotted the linear variation law (u_plus = y_plus) and the log-law
> (u_plus = ln(y_plus)/k +B) against the results of the simulations using
> the implemented solver (both the modified solver and the original one
> yield the same result). This is the result:
> https://docs.google.com/open?id=0B9qsTlUZxTpdWGtOY09fSVcyMXM with u_tau
> = 0.05 and Re_tau = 750.
>

The k-epsilon model is not very good, at least not for low Reynolds numbers like this where you try to resolve the flow all the way up to the wall. The results should be more or less exact for the V2F model.

Your yplus is more than 10 at the node closest to the wall. You need this to be about 1 to fully resolve the flow. See problems/channel.py and uncheck the bit with

x = m.coordinates()
x[:, 1] = arctan(pi*(x[:, 1]))/arctan(pi)

That will give you a mesh that is denser close to the walls.

> The slope is actually very good, and I think that using a wall-function
> may be the fix for this problem, do you think so?

If your problem has a high Reynolds number and you're not really interested in what goes on close to the wall, then wall functions are a good option. With wall functions you don't have to resolve all that much close to walls, yplus = 30 is ok

>
> I started thinking of a way to build this wall-function in 3D. My first
> attempt would be setting u in the nodes closer to the boundary to the
> following:
>
> normal = grad(y)/||grad(y)|| where y is the distance to the wall
> u = (1/||u||)(u - dot(u,normal)*normal)*u_wall where u_wall depends on y(wall-functions)

I think wall functions can be implemented much like the boundary conditions that are already there.

>
> This way u would satisfy the wall function and would is some sense
> conserve streamwise direction.
>
> To implement this, I would use a MeshFunction containing
> vertices_inside_boundary, mark them with some bid and then do
>
> L = (1/FacetArea(mesh))*(1/||u||)(u -
> dot(u,normal)*normal)*u_wall*v*ds(bid)
>
> and then assemble, to get a vector V_wf, right?.

Look at the Yplus class in Walls.py. It is generic and can be used for 3D. It probably does a lot of what you're thinking about here.

Not sure I follow the rest, but best of luck

Mikael

>
> The next step would be writing this values in the just-calculated
> vector, and this is where I have my doubts... I have seen the question
> https://answers.launchpad.net/dolfin/+question/95871
>
> Sparsity patterns might be different in this case, since I am only
> integrating over bid. I could modify this by adding, perhaps,
> -DOLFIN_EPS*v*ds, and then do V = V +V_wf. Is there a way to add these
> vectors directly?
>
> Well, I hope I was clear in my explanations, and thanks a lot for your
> time Mikael!
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Marcos Samudio (marcossamudio) said : #7

Thanks for the fast response Mikael. Will provide feedback on my results!

Marcos

Marcos Samudio (marcossamudio) said : #8

Thanks Mikael Mortensen, that solved my question.