I was trying to run a few cases at a higher Re on PyFR.1.12.1. The first case was a flow over 2d square cylinder at Re=2.2e4 using AC solver. As a start, I initiated the cases with p=1 and a sponge field for 50s, dt = 5e-3 and pseudo dt=2e-4. Everything seemed fine except the pressure residual was quite stiff and remained to be about 0.3~0.5. Then I restarted the case for more 50s with p=2 and [solver-dual-time-integrator-multip], the pressure residual dropped down to about 0.1~0.2, but the forces on the square were strange and had no periodicity. Thus, I extended the calculation for another 200s, but in vain. The forces still remains unstable, and the pressure residual was over 0.1. The profile of force-time and the contours of pressure and velocity (t=120s and 240s) are shown below.
t=120s
t=240s
However, It’s very encouraging that every time I reviewed the topics posted on the forum, I was able to capture some new tips.
The first tip is using the anti-alias, due to the aliasing error that may occur on non-linear problems, it’s advised to adopt anti-alias = (flux, surf-flux) option to stabilize the result. But once I activated the option, no matter how small dt and pseudo dt I set, the Nan problem appeared in the first period.
The second tip is to refine the mesh and decrease the dt and pseudo dt, I generated a finer mesh (almost 2 times the first mesh and kept the yplus < 3) but still got the same unstable problem even for a long run.
The third tip is using a higher polynomial order like setting p = 3/4. But due to the explicit pseudo time stepper and the multi-p method, the cost of calculation was quite high (nearly 100hs), especially for a long run. I know it’s hard to maintain balance between the efficiency and accuracy, but stability, in my opinion, would be supposed to be guaranteed in the first place.
Can someone give me any hints on how to stabilize the case and advice on improving the performance of the calculation? The mesh file and ini file are attached below.
There is quite a lot to unpack here, but I’ll try to address some of the points.
You say that the forces are ‘unstable’, from what you are saying I’m assuming you mean that there isn’t some nice regular pattern? I don’t have a lot of experience from square cylinders, but with round cylinders, as you increase the Reynolds number things do become more ‘chaotic’. A quick search comes up with this paper. Their flows look not too dissimilar from what you have.
The upshot is that if you aren’t getting NaN then I don’t think you need stabilisation, as this dynamics isn’t unstable in that sense.
Addressing the pressure convergence, pressure is almost always the slowest quantity to converge. It will also benefit the most from p-multigrid due to its elliptic character. Some really good configuration files for the SD7003 case can be found in Niki’s paper. These are broadly applicable to other flows.
Finally, I would question the physicality of this case. It is known that above a Re of about 300 the flow for circular cylinder starts to become three dimensional. Again, I haven’t done much work on square cylinders, but I imagine this carries over to them. Sharp features can make the physics less dependent on Re, but you’ll still have significant shear to induce three-dimensional effects. So when using a scale-resolving method on this case in 2D, I really think you’re forcing the solver to do something unphysical. This may be at the root of some of your issues.
Thanks for your kind reply. I have to admit that the case is actually a 3d problem and the loss of accuracy is inevitable. In fact, one reason I set the case in 2D is that I’d like to check my mesh was ok for the problem, so then I could extrude it in 3D with the least calculation cost. The other is I could compare with the result to weight how much the three-dimensional effect can be. So I will try a 3d simulation and see if any improvement can be achieved.
As for the regular pattern of the time history curve, I would be expecting a periodic sine-like curve rather than a random one. The paper you attached is quite impressive and has a certain similarity to what I’ve got. But from the time history curve of Cl, It’s a bit strange to see the average of Cl is non-zero, which may happen in the initial unstable period I think. I’ve looked through a few theses, and found the result showed a periodic curve at different degree, though most of them were 3d cases. Just weeks ago, I tried running circular cylinder cases with a relatively lower Reynolds number, and the results were pretty good. And that’s why I was surprised when not find the regular shedding of large eddies.
Last but not least, I did gain a lot from the p-multigrid and Niki’s replies from other posts. But It seemed that the pressure residual was choked and cannot reduce to below 0.1 even for a long period. Or maybe the time is not long enough? On the contrary, from the reply in the other post, with 3 cycles, the level of convergence of the turbulent Jet can be below 1e-3. I will try with longer time and a higher ac-zeta(if it’s necessary).
P.S. I found that sometimes the anti-alias = (flux, surf-flux) option can raise a Keyerror : Kernel mul has no providers, I would open a new topic if there exist any problem. Sorry for my being tediously long, I just want to explain everything I think is relevant.
Things are looking up a bit when I ran a 3d simulation with the pressure residual fluctuating around 0.01, smaller than 2d simulations. It should also be noted that the polynomial order in this case was only p1 with a P-multigrids cycle [(1,1),(0,2),(1,3)], the latter showed a strong control over the residual drops. However, the biggest weakness is the efficiency using such cycle, almost 20hrs to spend to step 20s.
Interestingly, I altered a few settings like ac-zeta, total time, dt/pseudo dt in 2d case, all of those helped little to the reduction of pressure residual. But when I set the number of pseudo-niters to a larger one (like 20 to 300), the residual dropped slowly but impressively from 0.2 to 0.01. After a few more rounds, the residual even dropped to 0.002. As an exchange, the cost of the calculation rose dramatically. I can explain this as the origin niter was set for laminar flow, which would certainly converged in a short period while for the turbulent flow, a longer period is required.
Additionally, in terms of the periodic curve I expected, I think I found an explanation that solved my question. See this paper. For large-vortex simulations, the sampling duration needs to be greater than 30-40 vortex shedding cycles (about over 300s) to obtain stable data. This picture also illustrates the problem.
I’m glad moving to 3D helped. I have two further points to make about this.
The first is that you are running p1. (order=1). DG, and by extension FR, does not have particularly good numerical characteristics at p1. You might want to try upping the order and maybe reducing the mesh resolution a little. If you are only going to run p1, you are probably better off using an FV scheme. (Mildly controversial viewpoint for the PyFR forum I know).
On the convergence of pressure, I wouldn’t get hung up on trying to get the pressure residual to machine zero, rather focus on making sure pressure is sufficiently converged to ensure the results are correct. I.e. if a pressure residual of 0.1 is fine then use that. If you need a pressure residual of 1e-6, then use that. Generally, I set the pressure residual tolerance to be fairly large and then set the velocity residual to be something sensible based on machine precision and Reynolds number. If you are going to investigate TKE budgets though, you probably wouldn’t want to use this tactic.
Once you have the dual time-stepping optimised, I typically end up with 1 to 5 pseudo iterations per physical time-step. Bear in mind that due to some deficiencies of linear multistep methods, you will typically see high residuals initially or after a restart. The simulation will need to drive out these transients before the dual time stepping will get into its stride. (i.e you might see 100 iterations per step for the first 30 steps or so and then it should decrease). When we replace BDF with DIRK implicit methods this will fix this.
But from the time history curve of Cl, It’s a bit strange to see the average of Cl is non-zero, which may happen in the initial unstable period I think. I’ve looked through a few theses, and found the result showed a periodic curve at different degree, though most of them were 3d cases. Just weeks ago, I tried running circular cylinder cases with a relatively lower Reynolds number, and the results were pretty good. And that’s why I was surprised when not find the regular shedding of large eddies.
