[aspera]

Originally Posted by sbaffini

I have to say that i overlooked at this method of computing just because it had another small parameter () which was, again, much smaller than what i typically need (order 1 without pressure gradient, order 0.1 in problematic cases when i use the pressure gradient term). What is your experience on this?

In this methods in all Mach number then , but then we need solve and if we set and solve the preconditioned equations corresponding (approximately) to . Therefore, we choose the as little as possible, but so as not to receive a loss of precision in the division.

I do not have experience in the calculation of the implicit scheme. In the explicit schemes without preconditioning I do not get convergence for low speed flows, when with preconditioning i get it. I did not compare improvements the convergence in explicit schemes whit different .

Originally Posted by sbaffini

This is something i didn't really got (in the paper) and i suspect i'm really missing something here. While i agree that using the square roots for each flux might slow down things a little bit (and especially so in the 90's), i can't see this as my bottleneck in the code performance nor can any profiling tool. So i simply decided to stick to the original Roe formulation. Actually i couldn't see any major difference in my tests (still, i didn't try this on the 5 Toro shock tube tests, some of which were failed by Fluent with both Roe and AUSM 1st order ). Is there any particular reason for which you decided to use averages (instead of Roe averages), besides trying to follow closely the paper?

So, i guess, you also use simple averages for on the faces when computing the fluxes.

Yes, you're right, I used the average by reason of Weiss paper and because it is easier than the Roe (+ code has been written on the GPU, and I was paranoid about performance ))) ). I think it is possible to average as you'd like when .

I also noticed that the Fluent strange working on some simple tasks (for example hypersonic).

[youmengtear]

Dear Paolo.
I am trying to develop a implicit density-based solver with time derivative preconditioning method in framework of foam-extend 4.0. But when I use a CFL number beyond 4.0, it diverged. I am not sure the system matrix I used is correct. Could you in convinent to explain how to discrete it in a cell with one or more face is boundary part? Thank you very much

[sbaffini]

Originally Posted by youmengtear

Dear Paolo.
I am trying to develop a implicit density-based solver with time derivative preconditioning method in framework of foam-extend 4.0. But when I use a CFL number beyond 4.0, it diverged. I am not sure the system matrix I used is correct. Could you in convinent to explain how to discrete it in a cell with one or more face is boundary part? Thank you very much

Dear WangYi,

first of all I want to stress that I am a practitioner more than an expert with a specific study in the field. So, take everything with a grain of salt.

In my code I have two kinds of boundaries, body fitted solid ones (walls and symmetry) and non solid ones (body fitted open boundaries and immersed solid boundaries).

The way I treat non solid boundaries is via ghost cells or sort of, in the sense that I apply there the convective scheme as if they were internal boundaries. This means that the implicit part comes from the convective scheme, but I simply add its matrix to the internal cell row-column, as the ghost cell is just fictitious for me (just some space in the variables array where I conveniently store what I need to apply bcs straightforwardly).

On body fitted solid boundaries I directly build the fluxes and apply them. However, for the convective part, I add (still, only to the adjacent fluid cell row-column) the jacobian of the flux to the implicit part, as I am solving in primitive variables.

For the diffusive fluxes I simply use the TSL approximation for the implicit part. What changes on the boundary, thus, is just the length at the denominator (and, again, the fact that I only add the implicit part to the row-column relative to the adjacent fluid cell).

In addition to this, note that on the boundary I don't even care about the exact scaling of these terms (even if, in the future, I will certainly do it). Thus, for example, the jacobian I add for the convective part on body fitted solid boundaries is the same I use for internal faces (so it has a wrong 2 factor in it).

In practice, except for the dissipation of the convective scheme (which is only added congruently), my implicit part on the boundary is identical to that in the interior, using for the geometric part a fictitious cell on the other side of the boundary that doesn't actually exist. In some cases I get the exact jacobian, in some others not. But for the second order scheme that wouldn't be exact in any case, so I don't care that much (yet).

When using this approach I can reach a CFL as high as 1E6.

Obviously, that's not always possible, and there are cases which need to run on much lower values. The range 10-100 is more typical of practical cases. Very few ones need to start (or pass trough) a 5 (default value in Fluent), but when I need to go below that, that's typically the signal for other problems.

In case it might help, I want to also add how I solved some of the issues I raised in my original question at the beginning of this thread.

First of all, I finally switched to a second order scheme that, given the gradient extrapolated values on the faces for p and T, would then compute with them all the other fields (instead of independently extrapolating density and enthalpy with their own gradients). Shock tube tests confirmed that this is the way to go (especially for near vacuum cases). This complicates the handling of terms that typically end up into the pressure (e.g., gravity potential) or enthalpy (e.g., gravity potential energy and turb. kinetic energy), but it is definitely the way to go.

For what concerns the pressure in the preconditioning velocity. I tried computing it in the classical way suggested in some papers (with loops on faces and using max over neighbors), but that produced worst results and code with respect to my original approach based on the pressure gradient already available.

In the Darmofal paper, cited in the first post, are also presented two scalings of the preconditioning velocity (you first compute it classically, then scale it with a function of M). One is that also used in Fluent; the other one is an original proposal from Darmofal. I finally found that the Darmofal proposal is better in my case. It is less aggressive of the Fluent one, but much better than the unscaled case over a large M range (but reduces to that for M = 0). The key here seems to be uniformity.

My convective scheme still uses a Roe average of the adjacent preconditioning velocities, but on boundaries now I revert to the max between the one in the adjacent fluid cell and the hypothetical one arising on the face for its current values of the variables.

Another aspect I underestimated is the unsteady preconditioning, which needs consideration if you go unsteady. There aresome paper around describing how they solved it in OVERFLOW.

[youmengtear]

Dear Paolo,
Thanks to your quick reply and sorry for my poor English. From your explaination, the diagonal part maybe wrong in my system matrix, as the boundary contribution is not added.
I'm only dealing with Euler equation, so the diffusive fluxes are not noted. The second order scheme comes with the same way as you describe. I have not read the Darmofal paper you mensioned, and I will read it later.
Thank you very much!