Hi, is there a general way to quantify numerical dissipation regardless of the spatial discretization scheme, if you're not using upwind? I'm working with FEM (2D incompressible, fractional step), and I'm getting a laminar pattern for a Re = 500 around a cylinder with about 10000 points in the mesh. Also, when running cavity flows, I can't get the top left vortex that appears after Re = 2500. For this Reynolds I'm using 165 X 165 points, but I've tried with 215x215 and the result was the same.

TIA

Márcio Ricardo

Marcio,

Check your scaling method - if you use one. In other words, if you draw in millimetres.

I scale firstly on residence time through the model & then on Reynold's number. The first scaling - if not performed will result in the same Re, but will not give the correct time-relationship.

Try it & see if it changes your result.

diaw...

Thanks, diaw. I'll try it. I'm running on a scale 1:10 ( my characteristic length is 10 times lower than the physical length: for a given Reynolds I use a velocity profile u = 1.0 and a length = 0.1, and then the viscosity is 10 times lower too, with a unity-density).

-Márcio Ricardo

Good question. There is certainly some information about numerical dissipation buried in your discrete equations, but to analyze them including all aspects of the method can be quite difficult. Of course, the obvious (and useless) answer to your question is this: If you know the exact solution, compare it with your numerical solution, and you can extract the numerical dissipation. Well, one thing you can do is increase the grid and time resolution (significantly... at least by a factor of 2). This will effectively reduce any numerical dissipation. If that doesn't solve your problem, I would tend to think that there is something else, aside from numerical dissipation, affecting your results.

About the cylinder, what do you mean by "laminar pattern"? On the wall, in the wake....? The wake (starting with the separated shear layer) should become turbulent first, before the boundary layer transition, as you increase Re above 180. One thing you should be aware of, though: 2D flow at Re=500 doesn't exist in reality. At that Reynolds number the wake is three-dimensional. That may or may not be linked to turbulence... And also: 10000 points seems quite low for this kind of flow. Are you using high-order discretization?

>>laminar pattern:

The flow past the cylinder looks like at Re=30, 40. There is only a static bubble downstream, which is less than 1D long. At Re = 500, even knowing that the flow is three-dimensional, I expected some kind of vortex shedding.

>>number of points and order of discretization:

The exact number of points I'm using is 16142, shared by 3978 quadrilateral elements. I'm using 2nd order lagrange elements for the velocity (9 nodes per element) and 1st order elements (4 nodes per element) for the pressure. I'm using constant weight functions (CVFEM), and the formulation is based on weighted residuals, not variational formulation.

If the problem is not numerical dissipation, could this be some kind of error in the spatial discretization of the advective term?

Thanks in advance

Márcio Ricardo

I see. Well, you're certainly right in expecting vortex shedding in your unsteady computation, and I can now see why you would blame numerical diffusion. It's possible that your scheme is so diffusive that the instability of the steady shear layers is suppressed. Considering your grid resolution, this seems indeed likely.

We have run this flow in the laminar regime quite extensively, because it's such a well-known benchmark case for unsteady viscous flow computations. On the grid resolution for a second-order accurate scheme (in space and time) I would agree with Alonso and Martinelli who showed that they needed on the order of 400 cells around the cylinder circumference to accurately and reliably predict the vortex-shedding frequency. For the laminar case! You will also want to have a large enough domain to be able to use a free-stream condition, and you need to resolve both the boundary layer and the wake region with a fairly fine grid, so I would not expect too much from a 16000 element grid. Especially not, if you're looking for turbulence.

Take a look at this paper before venturing any further:

http://aero-comlab.stanford.edu/Papers/reno95.pdf

Be skeptical about their results for the Strouhal number. Just like you, they ran Re=500 in 2D which is simply not realistic, so the Strouhal number does not match with experimental data. However, their arguments on the numerical issues are still valid, their resolution studies are very informative, and I suggest you read what they have to say on grid resolution and time resolution, even though your scheme is different.

To check on the spatial discretization you could actually start with steady computations. Run laminar 2D flow at low Reynolds numbers (below 47) and compare the size of the separation region and the position of the separation point with experiments, for various Reynolds numbers.

You are in great company with this problem, because it's such a standard test case that a lot of people have used it to validate their codes and done it successfully. It's safe to predict that you will eventually succeed. Simply because you must.

.... and some more questions: What initial conditions do you use for the unsteady flow computation, what time-step did you choose, and how many time-steps did you run?

The domain dimensions were: a rectangle with L=5, H=2.5, a circular cylinder with D=0.1 located at x=1.65 from the west boundary at y=1.25. The time step was dt=0.0003 and I've run 162000 time iterations. I always use density = 1.0 and set the Reynolds number only by varying the viscosity.

The initial conditions were u, v, p = 0 in the whole domain. As boundary conditions, a unity velocity profile on the west face, and on the south, east and north faces I used homogeneous neumann bcs for vel, with p=0 on them. Of course, u = v = 0. on the wall, dp/dn = 0 on it too. Also, in my fractional step scheme I did not include the pressure in the intermediate velocity equation. This leads to a 1st order scheme for pressure, but the velocity is still 2nd order in time, isn't it? I'm using explicit time advance with Adams-Bashforth for convective terms, with a lumped mass formulation for the transient term (diag(Mij) = sum_j(Mi), off_diag_i = 0);

Thanks a lot. I'll check the paper right now.

Márcio Ricardo

Have you got a convergent solution, for Re 500, starting up your problem from a stationary velocity field withou any kind of Reynolds ramping strategy?

" ...By the formulation you are working (CVFEM) I guess you are from UFSC in Brazil, am I right? ... "

Renato.

Yes, I've got a convegent solution, but not a satisfactory result, as I've posted. I'm brazilian, too, from Federal University of Uberlândia - MG(UFU). Actually, I'm with 2 advisors, one from here and the other from UNESP at Ilha Solteira. We look at CVFEM just like any other finite element method, by calculating the element matrices, assembling them on a global matrix and so on. I'm not sure, but I think in UFSC they use the method as an unstructured finite volume method, the shape functions and derivatives being used just to calculate variable and flux values at the control volumes faces, but no element matrices are assembled. Indeed, professor Maliska call their method "Element-based Finite Volume Method".(Anyone from UFSC please correct me if I'm mistaken). You work with FEM too, don't you?

-Márcio

Yes, I'm from COPPE/UFRJ and I've been working with edge-based stabilized finite element formulation (SUPG/PSPG to be more specific) with linear tetrahedral elements. I've just finished to implement the parallel transient part of my code and ran this problem for Reynolds 100 to validate my implementation. The result seems to be fine up to now and I haven't ran this problem with Re 500 because I intend to implement some kind of LES-based turbulence before (maybe a standard Smagorinky model before something more sophisticated).

It's good to find another Brazilian here... We could keep in contact (rnelias@gmail.com)

Regards

Renato.

I'm thinking about using SUPG too, but in a CVFEM fashion. I've seen it in the paper

"A flow-condition-based interpolation finite element procedure for incompressible fluid flows" [ Bathe and Zhang, Computers and Structures 80 (2002) 1267â€"1277]

They use SUPG for the advective term only, but their formulation is based on a discontinuous finite element method, which I'm not familiar with. For a 9 noded quad element, they take patches of 4 noded sub element and use linear shape functions in it. The weight functions are unity-constant too, and split each patch in four sub-control volumes. Do you know this method?

Regards

Márcio

My email: pivello@gmail.com

With H=1.25 and D=0.1 your upper and lower boundary are about 0.575 away from the cylinder wall, correct? That's only about 6 times the cylinder diameter and I don't think it's distant enough to use free-stream conditions (unless you use some kind of free-stream correction). As a rule of thumb, keep your free-stream boundary about 30-50 diameters (or chord lengths) away from the object (cylinder, airfoil...).

For the time step, you have to think in non-dimensional terms. Why 0.0003? How do you know that's going to be small enough, or maybe it's excessively small? If you are using explicit time stepping your time step requirement is most likely due to stability reasons (and your time step may be small enough for accuracy), but you should still check on it. The Strouhal number you expect is roughly 0.2, which means that the vortex shedding frequency is

f = 0.2 x free-stream speed / cylinder diameter.

Depending on your time stepping scheme (let's say second-order), you will need in the order of 50-100 time steps per oscillation period. The period is given by T=1/f, so your time step should be no larger than

dt = T/50 = 1/(50f)

I don't know if that's ok with your dt=0.0003 because you didn't tell what the free-stream velocity is. If you use a non-dimensional speed of 1, for example, it turns out that your time step is extremely small. That will give you good accuracy, but also means trouble: It will take an extremely long time before you see vortex shedding (maybe run your code for a couple of weeks).

The best is to use implicit time-stepping (e.g. dual time), so you can choose the time-step required for accuracy without having to worry about stability.

The reason why I asked for initial conditions: It may take a very long time (running through dozens of oscillation periods T), before you will see vortex shedding develop from an initially uniform flow field. Here's a trick to speed up the process: Start your unsteady computation from an asymmetric solution. First obtain a steady-state solution at some moderate angle of attack, say 5-10 degrees, and then start the unsteady computation (angle = 0 degrees) from that solution.

Mani wrote: As a rule of thumb, keep your free-stream boundary about 30-50 diameters (or chord lengths) away from the object (cylinder, airfoil...).

---------

Excellent point. Can I ask how you specify the 'free stream boundary' - assuming you are refering to the top & bottom lateral boundaries? Do you set a velocity & pressure condition - or 'do nothing'?

I would heartily agree to set any lateral boundary imposition as far away as possible, in order to minimise its interference.

--------- Mani wrote: Depending on your time stepping scheme (let's say second-order), you will need in the order of 50-100 time steps per oscillation period. The period is given by T=1/f,...

----------

Very good advice. What is the '50-100 time steps per oscillation' based on? I guess experience?

Thanks for that shared wisdom, Mani...

diaw...

>Can I ask how you specify the 'free stream boundary'

Diaw,

I suppose it depends on the equations you are solving. With a compressible code I am typically using a quasi-1D Riemann approach as a non-reflecting farfield boundary condition.There must have been early studies on the distance requirements, depending on the BC, but don't ask me about references. These rules of thumb are passed down from researcher to researcher to become common knowledge.

In case of the cylinder I actually did check on convergence as the distance to the boundary is increased. It's not a clear-cut case, because the wake is never going to be "undisturbed farfield", even far downstream. The question there is: How far downstream does the wake still induce an influence on the cylinder? And: How far downstream do you actually get an accurate solution near the wake, considering the numerical diffusion of the wake?

>What is the '50-100 time steps per oscillation' based on? I guess experience?

Yes, others' experience as well as my own, specifically for the laminar cylinder. You need a certain temporal resolution to minimize numerical dissipation, which would alter the vortex-shedding frequency or even suppress shedding altogether (damping effect of numerical "viscosity"). In that sense, Marcio's initial question on numerical dissipation really hits a nerve, but a solution to this problem is relatively straight-forward: increase both spatial and temporal resolution.

Mani,

Thanks so much for the excellent comments.

I have noticed that the proximity of boundaries, coupled with their influence on the 'wave nature' of the flow field certainly has a very significant effect on 'stiffening' the equivalent bulk modulus of the system. I call this 'container effect'. This effect thus alters the 'second viscosity' of the fluid 'system'.

The temporal & spatial discretization requirements again seem to correspond to those required to resolve a very low-speed event in a 'wave field'. I term this the 'viscous shock, or bump'. I have seen effects in 1D at speeds as low as 0.5 mm/s. This phenomenon is the thing we have been smudging out with upwinding & convection-stabilisation. In corresponds to the 'onset of unstable flow', but really only becomes dominant when large-enough to overcome the viscous damping of the fluid & system.

diaw...