Hi dear all;
I am modelling fluid flow arround a circular cylinder to find the drag forces. But I have some problems with boundary conditions. I implement Zou and He method for velocity boundary condition. To initialize the problem I set the same velocity for both inlet and outlet (left and right boundary respectively) and periodic boundary condition for up and down the domain. The ratio of domain width to circular diameter is 25. After simlulation the profile of velocity seems to be correct; but pressure of all nodes do nor reach a steady state and increasing with a relative high rate. Can anyone help me?
Hello,
I had the same problem with the pressure increase. The reason is a global mass increase of the system due to the velocity boundary conditions. Have a look at this thread here.
Reducing your Mach number slows down the mass increase very efficiently.
Maybe this helps you a bit,
Timm
Hi Dear Ghassemi,
I am also working to develop code for flows around a circular cylinder. I try to check the bounceback, periodic and symmetry boundary conditions and different types of inlet/outlet boundary conditions. But i observed some numerical insatbilities in my results. You also do the same work but using Zou and He boundary conditions for inlet/outlet and wall boundary conditions. If you don’t mind can you post here just your inlet/outlet and wall boundary conditions may be it help me and many others to know betterly of using such types of boundary conditions.
Good luck and bye for now.
Regards
Khan
I think that the problem of mass increase inside a channel can be avoided by using a pressure boundary condition on the outlet, instead of a velocity boundary condition. For the sake of illustration, take the example program “cylinder2d” in the current OpenLB release (0.5), which uses a pressure outlet. Set the reference velocity uIn to 0.1 to exaggerate the effect of mass increase, and observe the time evolution of the average density in the channel. Although the density initially oscillates, due to the initial condition, it remains bounded. If you modify the outlet boundary condition from pressure to velocity (this requires only changes to one line of code), you will see that the total mass gradually increases. You can watch a figure with these two time-evolution curves and download the source code on the community Wiki page corresponding to this post:
http://www.lbmethod.org/community:avoid_mass_increase_with_pressure_outlet
Here is an intuitive explanation of why this works, based on the discussion provided by Timm on this forum thread. Consider a 1D version of the channel, for the sake of simplicity. If the inlet has the velocity uIn, and the outlet has the velocity uOut=uIn, then the mass flux at the inlet is rhoInuIn, and the mass flux at the outlet is rhoOutuIn. Given that rhoOut<rhoIn, because pressure decreases along the channel, and density is proportional to pressure, the mass increases inside the channel. As a conclusion, it is not right to set uIn=uOut, because this presumes that the fluid is fully incompressible, which is never the case in lattice Boltzmann. As Timm suggested in the above-mentioned thread, one could fix this by computing an appropriate value of uOut which guarantees mass conservation. But although this is simple to do in a Poiseuille flow, it is more difficult when there are obstacles in the channel which affect the exact value of the pressure gradient. Now, if instead you impose a constant density on the outlet, the flow is allowed to do by itself what it wants, and will naturally evolve to a regime in which there is a mass balance. There is of course no proof that it always works, but it seemed to be working in all cases I tried out.
Hi Jonas,
this is of course a nice way to avoid mass increase. But you introduce an additional error source to the solution. Since the mass in conserved, the inlet and outlet flux are the same (at least averaged over some time). But the velocity profiles differ at the inlet and outlet, because
\rho_in u_in = \rho_out u_out
I have calculated that this error in the velocity is larger than the LB error itself. What are your results there?
Timm
I forgot one question:
When using inlet velocity and outlet pressure boundaries: If I impose a velocity profile at the inlet, how do I know the pressure value at the outlet?
From a theoretical standpoint, I would think that the error introduced by using a slightly different velocity on the outlet than on the inlet is consistent with the finite-Mach-number error of the lattice Boltzmann scheme. This is all part of the seemingly contradicting phenomenon of a “constant yet fluctuating” density when simulating an incompressible flow with a quasi-compressible fluid solver like lattice Boltzmann. More precisely, when one pushes a compressible fluid solver towards the limit of fluid incompressibility, one can assume that the density has two distinct components:
rho = rho_0 + Ma^2 rho_{hydro},
where rho_0 is a constant, “thermodynamic” component, and rho_{hydro} a “hydrodynamic” contribution, which merely reflects pressure variations (p_{hydro} = c_s^2 rho_{hydro}), and which is proportional to the square Mach number (Ma). For more details, and if you are at ease with German, you may want to read section 9.5.2 in the book by Dieter Haenel. The pressure drop in a channel is a hydrodynamic effect. Therefore, this term should scale like Ma^2. The thing is, the compressibility error in lattice Boltzmann sales like Ma^2 anyway. Therefore, multiplying the velocity by a factor proportional to a fluctuation of rho_{hydro} introduces an error which is compatible with the overall discretization error (at least asymptotically).
It is possible that you measure something different in a Poiseuille flow. But be aware that Poiseuille has so many symmetries that you cannot draw for sure any conclusion from it. For example, the paper on the Inamuro boundary condition shows that Poiseuille flow with this boundary condition is implemented with machine accuracy, independent of the channel width. This does certainly not mean that lattice Boltzmann is able, in general, to resolve flows with machine accuracy. It just means that Poiseuille has symmetries which are in favor of the lattice Boltzmann scheme and the Inamuro boundary condition. I would be quite surprised if, in case you solve a flow in a physically meaningful geometry, the artifact from introducing a pressure outlet instead of a velocity profile turned out to be anything but negligible.
As for the density to choose on the outlet, it seems to me that any positive value should be fine (you may want to create a false sense of security by choosing rho=1). The argument is the following. While on the outlet with fixed density the flow automatically constructs an appropriate velocity profile, on the inlet with fixed velocity it automatically leads to the appropriate density value. You know that the density should drop by delta_rho along a channel. If you fix the outlet at a density value rho_out, run the simulation, and wait until the initial transient flow fades out (you can of course shorten this transient period by choosing a good initial condition), the density at the inlet will reach the value rho_in = rho_out + delta_rho.
Thanks for your fast reply.
I think that I must contradict. In fact, you are right that both errors scale with Ma^2, but the prefactors are quite different.
For example: If I take a channel of dimension 25x50x200 with Ma = 0.04, Re = 10 (tau = 0.86), then I see that the LB L2-error using both inlet and outlet velocity BC is appr. 0.0006. A rough estimate of the error I make when I enforce equal fluxes on inlet and outlet is however 0.01. Both errors scale with Ma^2, but the “artificial” error, which is created in order to compensate the mass increase, is much more important.
Nevertheless I realize the charm of the velocity inlet / pressure outlet BC and I will try this approach.
Dear Timm and Jlatt;
Thank you very much. Your discussions became very usefull for me. I applied inlet velocity and outlet pressure boundary in my simulations and got good results for drag coefficient for low reynolds flow.
Now I have two more questions. My simulations take long time to reach steady state (drag force becomes relatively constant), for example with a 301*301 lattice dimension, running time by a high performance computer is about 10 hours (300,000 computation cycles). Also in the first cycles of time history, there are high fluctuations for the drag force and pressure (or density). This causes the program to stop particularly for higher reynolds flows (the problem arises when rho beomes zero). So I want to now:
- I want sure whether my algorithm is efficient or I can modify it to run faster.
- How I can solve the problem with shorter transient period or at least with flauctuations with lower domains so that rho (density) never becomes zero.
Hello ghassemi,
your code is able to process about 7.5e5 site updates / second. If you use a computer with a single CPU, this is okay. Maybe there is a factor of 2 possible by optimizing the code.
But I wonder why it takes 300,000 time steps. What is your Mach number and what is your kinematic viscosity in lattice units?
Timm
Dear Timm;
I use the following parameters in my simulation for flow around a circular cylinder (Re=20)
lattice= 301*301
lattice space=0.005
physical viscosity= 5e4, Lattice viscosity: 0.1
velocity=0.25
Ma=0.13
tau=0.8
I think that if I want to reduce Mach number with the same viscosity and velocity, The only way is to apply finer lattice space. In this case I should enlarge the domain to gain a suitable domain width to circular diameter ratio. Is it correct?
Sorry for mistake, physical viscosity=5e-4
For Timm: Yes, I agree with your point. Sometimes, a complexity analysis doesn’t say it all, and the prefactor matters. I was just suspecting that the low L2 error of 0.0006 was a “feature” of the Poiseuille flow, and that such good values cannot be expected in more challenging geometries. But on the other hand, more challenging geometries don’t have analytical solutions, which is why we can’t run an analytical-profile-vs-constant-pressure comparison on them, and my remark doesn’t make much sense anyway…
For ghassemi: By choosing your initial condition carefully, you can avoid instabilities and converge faster. My guess is that it is OK to initialize the particle populations at their equilibrium with a given value of rho(x) and u(x). A first thing to try is to choose the analytical Poiseuille solution for rho and u. This is much better than choosing (rho=1, u=0) and is mostly problematic close to the cylinder. If you are still not satisfied with the result, you can multiply the velocity by a function which drops to zero at the cylinder boundaries, and tends to 1 far from the cylinder.
@Jonas: I see that for more complex geometries the LB errors are larger in the first place, whereas the error due to fixed mass stays more or less constant. So you are right to assume that in some geometries keeping the mass constant does not heavily influence the accuracy. I will keep your precious advise in mind, because it may save me a lot of trouble in the future. After I have finished my current project I will deal with the question of keeping the mass constant in more detail, including velocity inlet / pressure outlet conditions.
Thanks again for your tireless advices!
Dear all,
I want to compare the results of the cylinder code in LBM with finite element method.
to solving the same problem at both method, I have some problems with inlet and outlet boundary conditions.
In FEA , I use a poiseuille velocity profile at the inlet completely similar to LBM code.
For the outlet I use an Outflow/Pressure at FEA which set the outflow presure equall to 0. (relative presure)
For LBM code cylinder.m , how should I define the outlet boundary condition to be equivalent to FEA problem? I used rho=1 in FEA. so regarding to the Inflow velocity/Outflow pressure boundaries which are mentioned in this topic how should I define my boundary conditions in LBM?
Thanks for your help,
Hello,
the cylinder.m simulation uses Zou-He boundary conditions and the outlet boundary condition is already a constants pressure BC. In fact in lattice boltzmann since pressure and density are closely related you specify the density at the outlet. Here you have rho=1. To get the pressure you hav to multiply this by the speed of sound squared (c_s^2=1/3). Then since the pressure is defined up to an additive constant you simply subtract the outlet pressure (p=1/3) on all nodes and you can compare the pressure profiles with the one you get with FEM. All this to say that you don’t need to change a thing in cylinder.m.
Thank you very much Orestis
So, you mean that if I set Rho(outlet)=1 at all nodes in Cylinder.m it is equivalent to put P(outlet)=p0 in FEM, Is it true?
The problem is just about the outlet,
Please consider the FEM solution as base which I set P(outlet)=p0, then what will be the outlet boundary in LBM regarding this? Is it p=1/3?
Thanks a lot for your Guides
As long as you impose something constant for the pressure on the outlet it does not really matter whether it is 1, 0 or 982374879233, physically you should get the same answer. In this case one just imposes p(outlet)=1/3 to avoid any numerical stability problems (if you use 982374879233 as p(outlet) then the code would probably become unstable) but theoretically you can do whatever you want.
To compare your simulations you then just have to take the same reference pressure. If you want it to be p(outlet)=p0 you run the simulation with p(outlet)=1/3 and then when you do your comparison you just do
p(whole domain)=p(whole domain)-1/3+p0
and you are in business.
Thanks a lot Orestis,
So it’s not really matter that the outflow pressure of one of the soulutions is p0 and the other is 1/3. and the solutions are equivalent, Is it true?
anyway, can you please tell me that p=1/3 is in Pascall or Bar and what is the relation?
Thanks a lot
Casy,
also here I refer to this thread.
The computer does not know any units, it does only know numbers. For this reason you have to specify a relation between the physical units (m, s, kg, Pa etc.) and the dimensionless LB units.
If the pressure in LB units for example is 1/3 then the density is 1, since p = c_s^2 rho and c_s^2 = 1/3. This means that the density 1 in LB units corresponds to the physical censity of, say, 1000 kg/m^3 in physical units. But be careful: In fluid dynamics only the pressure gradient is physical, you can choose the pressure reference value to be 0, 1000 Pa, 3.14159 MPa or any other value.
Timm