I have been stuck with this question for a while and I have not obtained the answer that will enable me move forward. I am seriously at a cross road as I cannot move forward without this solved.

I simulated flow, driven by body force, in a porous media. After a particular value of the body force, I get nan. However, I need to simulate at a high body force value (High Reynolds number) to obtain the necessary results to carryout further analysis.

I like to split the manifold reasons for which a simulation can be unstable into three categories: (1) You are simulating “wrong” physics, which means, a physical setup which is inadequate for your numerical method, or (2) there’s an error in the way you use the numerical method, or you ran into a limitation of the numerical method, or (3) there’s a bug in your code.

To eliminate class (1) errors I’d suggest you go through the following checklist:

Is the Mach number small enough? Look for the maximum velocity value in your simulation (“computeMax(computeVelocityNorm(lattice, lattice.getBoundingBox()))”), and make sure it is smaller than, say, 0.05. If you want to increase the Reynolds number further, decrease the viscosity instead of the body force.

Is your simulation laminar? If the Reynolds number exceeds, say, 100, you should be aware of the fact that the simulation might be turbulent and require that you are knowledgeable about fluid turbulence (the simulation will almost never spontaneously do the right thing when the flow is turbulent, unless the resolution is very large).

As to class (2) issues, consider the following measures:

Numerical instabilities are very often related to boundary conditions. If your porous media is in contact with the inlet/outlet, then this is something you may want to avoid by adding a porous-media-free zone close to inlet and outlet.

If the simulation is unstable with a pressure condition at the inlet, and a pressure and the outlet, consider using a flow condition instead (fixed velocity at the inlet, and Neumann condition (zero velocity-gradient) at the outlet).

Try out both the local and the non-local version of the boundary condition (“createInterpBoundaryCondition” vs. “createLocalBoundaryCondition”).

If you are desperate, try to trade BGK for one of the schemes that are reported to be sometimes more stable (RegularizedBGKdynamics and MRTdynamics). However, it is good to rule out the other items mentioned in this post first, because it is always better to heal a sickness than covering up the symptoms.

For class (3) you’re on your own, because the number of possible bugs is infinite. As a hint, you should try taking snapshots of the simulation right before the instability occurs, in order to understand where in the geometry of the system the instability starts.

Good luck with this. If you’re able to fix your problem, it would be nice if you could tell us what the problem was, so we gain experience on typical issues encountered by end-users.

Thank you Jonas for your reply a I understood it very well. I want to start by stating the conditions of my simulations. I used a bod force and periodic boundary conditions in my simulation

First, I want to rule out class (3) because I already simulated at lower body forces and my simulations are awesome, Results came out close to analytical values and I have visualized by taking snapshots and creating videos and they are all good.

For class(2), for the reasons that I already have good results from my earlier simulations with periodic boundary conditions, I think boundary condition is out of the way. Also, I have the porous media free layers around my simulation domain. However, I only worked with the local boundary condition (CreateLocalBoundaryCondition). I would run the simulation using the non-local boundary conditions. But I do not understand the difference between the two of them(would appreciate an explanation) , and also, is the non-local applicable to periodic boundary conditions?

Now, I think the failure of the simulation at the high body force is due to a high Mach number. The velocity for calculating the Mach number is it the maximum in the media, or the average?

I know that two ways to increase Re is to either reduce viscosity or increase resolution. I cannot increase resolution because I work with an image file that was presented to me. I am left with the option of reducing Reynolds number. However, for the BGK method that I used, simulation result is highly dependent on viscosity (defined by the relaxation time). Thus, I have already started thinking of the MRT method whose result is independent of the fluid viscosity. Thus, before going into it, I want to ask, does it also suffer from failure at high Mach number.

In all, I know that LBM is accurate at low Mach numbers. However, even if the Mach number is large, should it not still give an output, even though it may be wrong?
Thanks

Also the MRT is only valid in the small Mach number regime. Unfortunately, there is no other way for you as to increase resolution. Basically, what you could do: Increase the number of grid points along each axis to a factor of 2. If you are in 2D, then each obstacle node becomes four obstacle nodes (2x2). In 3D, you get eight from one single node (2x2x2). This does not increase your physical resolution but the numerical resolution.