I recently tested a cavity flow case on PyFR with ac-navier-stokes. It is an OpenFOAM tutorial case which simulates a cavity with three non-slip fixed walls and a moving walls moving at speed of 1m/s. In my test simulation,
Since it is a 2d flow, I applied one no-slp-wall with u = 1.0, v = 0.0 and the other three fixed walls with no-slp-wall u = v = 0.0.
“ac-navier-stokes” solver was used. Different schemes and pseudo time step schemes were tried.
Artificial compressibility beta “ac-zeta” was varied.
Mesh contains 100*100 structured grids.
2nd, 3rd, 4th and 6th order of solver were tested.
However, it is quite surprising that none of those tests above on PyFR could maintain the boundary conditions (I also tried ac-in-fv, no helpful) and express a symmetric velocity profile in that cavity. The results simulated by AC method are more like compressible flow but incompressible one.
The best solution showed in 6th-order ac-NS solver, while the boundary conditions could still not be held, even with very large ac-zeta (=100).
For that higher order solver, time elapse on a GPU is twice than IcoFOAM on one CPU. Actually, it is because that a much less Delta t was required for PyFR, but after all, longer time spent.
Additionally, the variation of artificial compressibility factor has a huge effect on the solutions, which makes me suspect the reliability of my results. (The maximum velocity could varied between 1.059 -1.066 when ac-zeta is tuned between 2.0 - 6.0)
The comparison between simulation solution on PyFR with 6th order accuracy (left) and OpenFOAM icoFOAM (right) is shown below.
Hence, I am just confused about the solution of AC method in PyFR. Frankly, I am confused about whether AC-method is reliable for incompressible flow simulations.
Can you send me your simulation files (.pyfrm, .ini) and I will investigate. The artificial compressibility method for internal flows is more tricky because of the pseudo-pressure waves are trapped and interact with the viscous flow field. We have applied the AC solver to several external incompressible flows and got excellent agreement against projection based solvers and very low Mach experimental data.
Here are the simulation files (attached). So have you tried internal pipe flow for both steady and unsteady conditions, which are the cases I really concern about.
You are comparing the solutions at two different Reynolds numbers, 100 on PyFR and 10 on OpenFOAM. The example case of OpenFOAM uses a side length of 0.1 whereas in your PyFR mesh it is 1. I have attached a file showing the results for P=4 at Re=10 (left) and Re=100 (right) using the artificial compressibility solver with the p-multigrid. You can see that Re=10 resembles your OpenFOAM solution.
About the speed. PyFR targets unsteady turbulent flows and icoFoam is tailored for steady laminar incompressible flows. You cannot compete with an unsteady solver against a pure steady state solver in a steady problem. Moreover, the problem size is very small and does not fully leverage the parallelism of GPUs.
I have not simulated any internal flows using the artificial compressibility solver. You can search for older work by Dochan Kwak on internal flows using AC.
Thanks a lot. But still I am wondering why the boundary velocity is 1.075 instead of 1.000 on the solutions from PyFR. And does it change much with the AC factor?
I dont know whether my understanding is right. To approximate incompressibility through AC method, the AC factor should be relatively very large than pressure so that the continuity eqn would be very close to the incompressible continuity eqn.
I tried a test in 2d internal pipe flow with constant pressure gradient of 4. There seemed to be something like normal shock happened when ac-factor/P = 0.5, while if ac-factor/P reached 10,000, the flow is quite incompressible.
Does ac-factor/P matter the incompressibility of the flow?
Artificial compressibility factor sets the pseudo speed of sound and therefore adjusts the characteristics of the system. Typical values for ac-zeta are (1 – 10)*maximum velocity magnitude. If you set ac-zeta very large, it limits your explicit pseudo-time step and kills the performance.
A strong pseudo pressure wave occurs always in the beginning of the simulation. With a large ac-zeta value, it seems weaker than with lower ac-zeta because it has travelled further with respect to the flow field (pseudo-speed of sound/flow speed). In an ideal case, this initial wave would be killed during the first physical time step by just performing a lot of iterations, but it is not feasible. Therefore, the residual of this initial wave tends to propagate in physical time during the initial transient phase of the simulation. With external flows the wave can be easily dissipated by a sponge region at the far-field boundaries. In internal flows you just need let it bounce during the initial transient phase and wait until it is naturally dissipated by viscosity. Please not that this wave occurs only in the very beginning when you start the simulation from initial conditions that are “far away” from the actual solution. You only need to develop the flow once and then you can keep restarting from the developed solution to gather statistics.
About the boundary conditions. The values in the solution file are domain-sided values from the polynomial representation of the solution. Remember that the solution is discontinuous in FR! The actual boundary velocity would be
BCvelocity = 0.5*(u_domainside + u_ghostside),
which would be exactly one for the lid in the cavity flow.
Thanks for your advice. I have got good alignment on cavity flow now.
I am now trying to simulate pressure driven laminar 2d pipe flow with periodic inlet/outlet BC. The velocity profile at roughly 50 seconds is still not parabolic. Instead, the flow generates some vortex at ReD less than 200 when ac-zeta = 20. The result is attached also the .pyfrm and .ini.
The pressure gradient is constant 4 and nu = 0.01.
If you look at the pressure field early in the simulation you can clearly see that there is a periodicity issue. The source term you impose is perfectly fine, but it does work for the gradient between the periodic back face and the front face.
p(x) = 0 - 4x,
p’(x)= -4
but derivative at the discontinuity reads p’(0) = [p(0) - p(8)]/dx = [0 - (-32)]/dx.
Currently, we do not have a functionality to force the constant pressure gradient in a periodic domain. However, for this case you can add the following line to pyfr/solvers/baseadvec/kernels/negdivconf.mako, which adds a “pressure derivative” source term (-4.0) to the x-momentum equation (index [1]).
% for i, ex in enumerate(srcex):
tdivtconf[{i}] = -rcpdjac*tdivtconf[{i}] + ${ex};
% endfor
I did not properly test it, but I have attached a picture what I got after running for a minute. Please note that the pressure looks constant if you add the source term this way, because p = du/dx + dv/dx = 0 + 0. I also suggest using BDF2, rk4 or tvd-rk4 and multigrid.
Can you tell me more about the meaning of this solver code? Actually, I need not only add const pressure gradient but also pulsatile dp/dz as well. So I hope I can understand a little bit of it so that I can arbitrarily apply different pressure condition.
By the way, PyFR in my machine is installed through pip3. I dont know whether it matters, since it won’t change anything after the modification of that .mako file. The speed of my flow is zero in the whole pipe. Actually, even if I delete .mako file pyfr can still run.
Additionally, I saw that, in the negdivconf.mako file, it shows ‘ndim = 2’. Does that mean that solver is for 2D case. What file should I modify If I am going to run 3D pipe with pressure gradient on z-direction?
Sorry to bother you again. Just want to know what “srcex” is and how the code imports the constant parameters, such as nu, in .ini file into the kernel.
How do I refresh pyfr after modifying the kernel .mako file?
Dear Niki,
similar to Will I’m trying to use PyFR to simulate incompressible, turbulent pipe flow. As for the nature of pipe geomtries the simulations would greatly profit from a streamwise cyclic boundary condition. Hence my question: Would your workaround of enforcing a constant pressure source term in the negdivconf.mako file work for a 3d case as well? Would a pressure drop in z-direction be achvied by adding for example “tdivtconf[3] = tdivtconf[3] - 4.0;”?
Do you plan to implement such a streamwise cyclic boundary condition in PyFR in general?
