CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Problem with secondary flow in a square duct (

cfdmarkus October 29, 2008 07:10

Dear all I am trying to sim
Dear all

I am trying to simulate the secondary flow in the corners of a fully developed square duct. To capture the secondary flow I am using an EASM turbulence model.

Instead of getting a nice smooth flow field I get some sort of oscillations in the flow field, see figure. The figure is obtained using upwind schemes on the convection terms and shows the shear stress Rxy in the lower left corner.

I don't quite understand where these artificial structures come form.
Can someone point me in the right direction to understand the origin of these structures.


Thomas Baumann April 15, 2010 04:17

Hallo Markus,

I'm trying to simulate a fully developed quare duct, too. I implemented some different EASM models.
But solving a square duct I get the secondary flow near the corners wrong spinning. I tested with existing nonlinear turbulence models like the NonlinearKEShih and I get the same wrong directions of the flow. Do you have a similar problem or know what's wrong implemented?

Thank you very much, Thomas

cfdmarkus April 16, 2010 09:36

Hi Thomas,

I think I cannot help you much here - sounds a little bit like you have swaped a +/- sign somewhere in the nonlinear terms.
Have you validated your model on a case like flat plate or channel flow in order to check whether the stress components are correctly predicted.


Thomas Baumann April 20, 2010 10:47

Hi Markus,

I am testing the models at a 2d channel flow now. Have you solved square ducts with ASM's and the secondary flows right-spinning?

I know that changing
-(S_ik*Skj...-0.33*S_kl S_kl delta_ij) to +(S_ik*Skj...-0.33*S_kl S_kl delta_ij), the secondary flows spin right, but I'm sure that the code has no minus/plus error here.

But after testing the 2D-cases I will be able to tell you more.

Thank you and regards


cfdmarkus April 21, 2010 03:59

Hi Thomas,

when I did the sqaure duct the secondary flow was in the correct direction.
However, the EASM I am using has a (+) sign in front of the (SikSkj -...) term.
Which EASM are you trying to use?


Thomas Baumann April 21, 2010 06:35


I've implemented the basic ASM model from Gatski ansd Speziale and the normalized one ,too. Papers: On explicit algebraic stress models for complex turbulent flows by Gatski and Speziale (normalized model equation 52), Towards the development of second-order closure models for nonequilibrium turbulent flows by Speziale and Xu (Standard-model equation 13).
But I've seen many ASM-models and they are all implemented in this kind of way:

tau_ij=2/3 K delta_ij -constant(+a1 k²/epsilon S_ij + a2 k³/epsilon²(S_ij w_kj + S_jk*w_ki) -a3 k³/epsilon² (S_ik J_kj - 1/3 S_kl S_kl delta_ij))

with a1, a2, a3 > 0. So I think my implementation is right. Which ASM model have you implemented?

I've problems because of the velocity of the 'nearest cell to wall' has a velocity in orthogonal-direction to the wall, see the pic on the top. Do you have such problems too? Do you use damping functions or corrections for the near wall cells in the tau equation or something in this way?

Regards Thomas

cfdmarkus April 21, 2010 08:15

I am using the EASM of Wall&Johansson with the SSG pressure strain model.
I never had any problems with that model and it also works without damping functions (I am using a omega based formulation).

It still would be interesting to see your normal Reynolds stress components for a simple channel flow - I have a feeling these are not correct, other wise you should get the correct orientation of the secondary flow.


Thomas Baumann April 21, 2010 10:15

Hi Markus,

at the moment I'm on holiday and I don't have access to OF. Next week I'll be back and be able to simulate the 2d channel flow and will post the post processing.

Regards Thomas

Thomas Baumann April 27, 2010 10:33

Hi Markus,

I have implemented the basic 2d ASM from "An explicit algebraic Reynolds stress model for incompressible and compressible turbulent flows" by Wallin and Johansson. (Equation 1.18, 1.19, 2.3, 2.4, 2.5, 2.7) The secondary flow rotates correct :).
I'm implementing damping functions to get a low reynolds turbulence model, that the flow near the wall is parallel to the wall... Is it okay to start with the high-Reynolds-kEpsilon-model, or should I use a low-Reynolds-model? Because of problems using damping functions in epsilon, nut and in anisotropy tensor, too.

After all I will compare the results to experimental data of a channel flow.

Regards Thomas

cfdmarkus April 27, 2010 10:45

HI Thomas,

I am glad to hear the rotation is correct now.

Why dont you just use the same kEps formulation which Wallin&Johansson have used. If I remember correctly, it should allow for an integration to the wall.


Thomas Baumann June 16, 2010 08:13

Hi Markus,

sorry for my late report. It took me a long time to implement some new models and share their results with DNS data.

I implemented a ASM from Baglietto and Ninokata "Anisotropic Eddy Viscosity Modeling to Industrial Engineering Internal Flows" from 2006. It's with a low-reynolds formulation. I implemented it and simulated some basic flows (pipe flow, channel flow, a heated bundle and a square duct). The secondary flows at the square duct and bundle are good. The velocity-profil in the channel and pipe flow is good, too, the Reynolds-stresses are okay (far from the wall). I simulated a DNS test case from

The problem is, the fully developed channel flow in x-direction (umean ca. 18 m/s has velocity changes in y-direction, so vmean isn't zero). Regard the plot. I plotted the pressure profil, which should be constant in y-direction ,too. The problem is maybe the near-wall pressure condition (which is in my case the standard zeroGradient condition).

Have you had such problems near the wall? Here is the plot of the Reynolds-stresses for a channel flow with Re_tau = 395.

Or has anybody an idea about a better bc for the pressure?

Thanks a lot,


cfdmarkus June 18, 2010 05:51

Hi Thomas,

how is the convergence of the EASM model? How many mangnitudes do the residuals drop? I found that EASM sometimes dont convergence very well which may also results in some kind of oscillations in the flow field.

Does the blue curve correspond to the vv component and why does it have such a strange behaviour close to the wall?


Thomas Baumann June 18, 2010 10:14

Hi Markus,

the residuals are good (all of them < 10e-7).
the simulation of the fully developed channel flow showed velocity changes in y-direction (see the pics above). And because of that, S_22=2*(du2/dx2) isn't zero... and I think that's why there are changes within the Reynolds-stresses near the wall. Maybe it's the definition of the divdevReff, I took the one from the LRR-model:
+ fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
- fvm::laplacian(nuEff(), U)

You're right, the blue one is the one in y-direction (so R_vv).

How have you defined the divdevReff? Do you update nut after determine k and epsilon but before calculatin R or afterwards?



cfdmarkus June 21, 2010 05:59


I am using the same appraoch as it is used in the other non-linear models, i.e.
divDevReff =
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

I update nu_t after computing k and epsilon!!

Hope this helps.

Thomas Baumann June 25, 2010 14:18

Hi Markus,

I have implemented the turbulence model like the other non-linear models in OpenFOAM. I have still the same problems with the velocities to the walls and the pressure.
But simulating a channel flow with these models like the nonlinearShih or the LienCubickLowRe, I get here the same problems with the non-constant pressure over the height and the normal velocities, too.
Have you simulations with these models without this effect?

So I don't know what to change at the model, now. Which model did you implement?

Best regards,


cfdmarkus June 30, 2010 03:58

Hi Thomas,

I have results for a planar channel flow using an EASM.
In my simulation the V-velocity is zero and the pressure is constant across the channel.

Can you confirm that the initial Residuals for the pressure p_0 does drop below 10^-7; as you have mentioned before. If this is the case, I can't think of anything else that could go wrong.


Thomas Baumann July 1, 2010 09:59

Hi Markus,

the residuals of the pressure are below 10^-7.
I found out that using divDevReff without the nonlinearStress_-part,

divDevReff =
without: fvc::div(nonLinearStress_)
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))

the static pressure is constant over the whole height, but of course, there is no anisotropy of the turbulence for calcutating the NS-equation. But from the velocity-gradient, the k, epsilon I get quite good Reynolds-stresses.

But using the fvc::div(nonLinearStress_) the pressure isn't constant.
Could you please show me an example of one of your implemented EASM? I have no more idea at the moment how to solve the problem.

Are you using mapped boundary conditions?

Regards Thomas

Gearb0x July 1, 2010 12:03

Dear Thomas and Markus,

sorry to step in the discussion but I'm working on a U-Bend with a square section too.
Could you tell me if you used the RSM model from OpenFoam or one you implement yourself?

Thanks for the help,

cfdmarkus July 2, 2010 04:18

1 Attachment(s)
Jérémy: Unfortunately, EASM models are not currently available in OF. But the RSM seems a good idea for your case.

Thomas: This is what I get for a channel flow (cyclic BC) (U,V,p).

If you want I can take a quick look at your EASM model implementation.
You can also send me the test case you are running if you want.


ancsa May 10, 2012 10:50

Hi, it's been a while that someone posted here but maybe you made some progress since then.

I was checking the NonlinearKEShih model for a simple boundary layer flow, so that the only velocity gradient is dU/dZ. For this the analytical result can be derived so I could compare to that.

What I see is that with the current implementation of the model the resulting trace of the Re-stress tensor does not give the correct k value. I think this is because some dirac terms of the Shih model are not implemented.

Also if adding those terms the R_11 and R_33 terms are changed compared to the analitical solution, which is similar to me to the wrongly rotating square duct secondary flow.

Did you finally work out how to use the NonlinearKEShih model or use models implemented by you?

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