Reynolds doesn't have an impact on the simulation

I’m trying to solve the following equation (incompressible NS) over a 2D domain (Tesla Valve):
CodeCogsEqn (2)
S(u) = (grad(u) + grad(u)’)/2 is the symmetric part of the gradient of u, with ’ being the transpose operator.

I’m forcing a velocity profile on the inlet (Poiseuille flow) and zero velocity on the walls.

To handle the non-linear term (u*grad(u)) I’m using Newton’s method. Subbing u and p with u+du and p+dp everywhere and solving in du and dp. I’m using an initial condition for u and p to be 0 everywhere, and iteratively set u := u+du and p := p+dp.

Finally, I’m solving for Reynolds = 50, then use that solution as starting condition for Reynolds = 100, then for 150 etc… until 500.

Everything runs pretty smoothly, converging with ease in a couple iterations per Reynolds number, but the results I’m getting are off.


Normalized Velocity at various Reynolds. A pressure plot snuck in the gif, my bad.

A tesla valve should cause a steep pressure drop in the backwards direction. At low Reynolds, the flow is basically a Stokes flow, and as such has pretty much the same pressure drop in both directions. As Reynolds increases I should begin to see its effects. Instead, the plots for velocity at the various reynolds are VERY similar, only scaled up as Reynolds increases. The plots for pressure DO change with varying Reynolds, as they should to allow those velocity plots.


Normalized Pressure at various Reynolds. The pressure plot from the previous plot snuck out from here and is missing.

Comparing to other visualizations I found on the internet, it doesn’t look like I’m getting close to the real solution, and either way by simulating the forward flow (which should have a much smaller pressure drop) I’m getting approximately the same results.

Do you have any pointers on anything I might have gotten wrong?

If you have the time (and will), I also uploaded the whole project here. It’s definitely not the best code I’ve written, but if you want to give it a look and maybe find out what I did wrong you would absolutely make my month (that’s, more or less, how long I’ve been trying).

Your linearization does not depend on the Reynolds number, so why would you expect your flow to do so?

Reynolds changes based on the inlet speed, wouldn’t it be reasonable to guess that Reynolds has an impact on the flow?

Also, you can non-dimensionalize the equation to:
CodeCogsEqn (3)
making it depend explicitly on Reynolds.

You just change the inlet speed. I wouldn’t expect much of a difference in the flow far from the inlet, e.g., near the valve. If your linearization depends on the Reynolds number, then you’re much likely to see changes. You could just use this script to see if this is indeed the case. You’d have to change the mesh and put the appropriate BC, but it’s solving steady-state NS with a continuation on the Reynolds number (exactly what you are doing, I guess?).

Could you expand on what you mean when you say a linearization that depends on the Reynolds number? Do you mean that the linearized equation needs to depend on Reynolds or something else?

Also, I can’t run the code because I don’t have mpi installed (yet, I might have to do so to test out the code later), but I’m sure the results should be different since I’m basically trying to imitate this paper. The only thing the paper does not provide is the governing equations used, but I’m fairly confident in mine, since they are basically your average incompressible NS with stress-free natural neumann condition.

Do you mean that the linearized equation needs to depend on Reynolds or something else?

Yes, if you look at the aforementioned script, you’ll see that the varf depends on Re, and thus, as the Reynolds number increases, you move from something very laminar to something slightly more “turbulent”. What you have looks like something very laminar, albeit with different forcing terms on the inlet as you do your continuation.

I don’t have mpi installed

This script should run in sequential, you can see that the same applies with respect to the Reynolds number/varf throughout the continuation on Re.

You need to update your variable mu throughout the continuation.

I’ve been trying to understand the latest example and I couldn’t get it to work in my scenario: I haven’t really learned how to use the most advanced features, and while I could mostly understand how the varf was being used (though I’m still not familiar with it), I couldn’t manage to figure out what code was specific for that example and which is needed in the general scenario.

Either way, I rewrote my code as to mimic the one in the example:

  • Removed Reynolds from the dirichlet condition and made it appear into the equation through non-dimensionalization
  • Generated an initial state through Stokes, by imposing flow at the inlet, non-penetration and no-slip at the walls and no conditions at the outlet
  • Solved iteratively for various Reynolds forcing du and dp to be 0 on every boundary, since the Stokes flow already verifies all the BCs

The backwards flow (the one that should have a high pressure drop) stops converging at Re > 350 (the iterative steps don’t converge anymore), while the forward flow stops at Re > 450.

This kind of asymmetry could be a good sign, but I’m still getting very close pressure drops, so I’m still having the same problem as in the OP.

You don’t start from the Stokes solution for the Re > 350 solves, right? You do perform a continuation, by reusing the solution from a lower Reynolds solve?

Stokes as the basis for Re=50, Re=50 for Re=100 and so on, so yeah, I’d guess it’s what you call continuation.

You can find the updated full code here for reference.

I’m not really interested in what happens after Re=250, to be honest, so my main concern is why the solutions at Re=250 don’t have a large enough pressure drop.

Just a remark, Your pressure are complete wrong

my the solution

the scriptNS-Dipilato.edp (1.9 KB)
and the zip of the mesh file bb.msh.zip (108.8 KB)

The foward direction looks VERY promising and pretty similar qualitatively to what it should be.

I’m trying to do the backwards direction, too, but to do so I need some geometrical information on the .msh you are using. Can you share them? Or do you know of a software that can let me open .msh files to see the geometric properties?

In particular, I need the coordinates of the inlet to correctly force the inlet velocity.

Dear LEONARDO,

I just build le mesh in 10 mn with emc2 software (very very old cf. emc2-v2.17.tar.gz )
work under mac and linux with X11.

the inlet segment is :

[15., -.5] to [15., .5 ] label 2

the outlet segment is
[ 7.27817459 , -6.57106781] [ -6.57106781, -7.27817459] label 3

In inlet, I just put a Poisseile flow with max velocity = 1

the emc2 data file is

bb.msh_emc2.zip (117 KB)

I tried it out and it worked awesome! Thank you (both) for all the support!