Week 3: Progress And Pretty Pictures

Hello everyone, and welcome to week 3! In this blog post I’ll continue discussing some of the details of LBM and include some pictures from my simulation. 

~~~ 

Cylinder in a channel: Last week, I introduced how the lattice grid worked in LBM. In the channel, I’ve also placed a cylindrical obstacle, and projected onto our 2d lattice grid, it’s just a circle. How this is represented in code is I have a 2d list of boolean values (the same size as my original lattice grid), and values are True where the cylinder is and False everywhere else. The code to create this is rather simple because we can use the equation of a circle and a less than symbol to denote we want everything inside the cylinder to have the value True.  

Equations in more detail: (I have not figured out how to type nice equation text in WordPress yet, but just bear with me here) 

The most popular Lattice Boltzmann Model uses the Bhatnagar-Gross-Kook (BGK) approximation, which describes the relaxation of a non-equilibrium distribution to an equilibrium distribution. 

The pressure-corrected equation governing LBM, a discretized Boltzmann equation, proposed by Chen et al. is

where ƒ(x,t) is the pre-collision particle distribution, Ω[f(x,t)] is the collision operator using the BGK approximation, and ƒ (with a tilde) is the post-collision distribution function. ƒ tilde is streamed to its neighboring lattices and thus becomes ƒ(x+∂e,t+∂).  

Ω is the collision operator, given by the Bhatnagar-Gross-Kook approximation

Ω is a function of the current distribution function and τ (tau), the relaxation time, which describes the rate at which the distribution progresses to equilibrium. tau relates to the kinematic viscosity, ν, of the fluid as ν =  c2(2τ − 1)/2, where c is the speed of a propagating wave in the fluid, c2 = 1/3. ƒeq is the equilibrium distribution function, and for a D2Q9 lattice, is given by the equation

where w_i are the predetermined weights for the nine directions and ρ is the macroscopic density of the fluid. D2Q9 uses the following weights for the nine directions described in last week’s blog post. 

From the particle distribution, the fluid’s macroscopic velocity and density can be calculated and used in the equation for ƒeq. These are the densities and velocities we’re more used to seeing in physics classes. 

~~~ 

Back to the Python implementation: Phew that was a lot of equations… But it all looks a lot more manageable when put into code. 

A little refresher from last week:

From the Stream to Boundary Conditions step, it’s all in a for loop that runs for however many iterations you choose. First thing in the loop is Stream. Here, I used numpy’s roll function, which makes this step a lot easier. The next step is calculating the macroscopic fluid variables, made easy with np’s sum. With these macroscopic fluid variables, we can calculate feq in a for loop that loops through the nine directions and calculates values across the entire lattice grid at once with : as list indexes. The next step is Collision. Here, we use the feq that was just calculated to find the f distribution for the next time step. And finally, there is an if statement that will be carried out if we have set wall_boundary to True in the parameters section up top. This just means we want to simulate the presence of walls at the top and bottom of the channel. If not, it means we’re looking at a section of fluid from a larger expanse. The cylinder boundaries are treated the same way, with the bounce-back method. Essentially what this method does is it reflects the incoming velocities at the boundary. If a particle comes in with a velocity to the right, it will be “bounced back” and leave with a velocity to the left. It’s basically elastic collisions with the wall at a particle level. In the code, you’ll see the lists [0,2,1,4,3,7,8,5,6] which is the order of the reversed velocities (original [0,1,2,3,4,5,6,7,8].) Directions per the lattice definition in Blog 2. After setting reflective boundaries, the LBM part of the loop is done, and it’s time for graphing! 

Graphing: Vorticity is defined as the curl of the velocity, and it’s essentially a measure of the tendency of a fluid to “spin”. When you graph vorticity with different values mapped to a color gradient in matplotlib, it makes it really easy to see which parts of the fluid are turning clockwise and counter-clockwise. 

This section of code is placed in the for loop, and for every time step, it graphs the values of vorticity across every lattice point. ~~~ 

Progress: 

So… the first time I graphed my results, this happened! That does not look anything like a fluid. I checked and double checked my equations, and everything seemed alright. After quite a while of attempting to debug, I realized the order of my code was not the same as the algorithm flowchart. I had placed the reflective boundary conditions in between the stream and collision step. Turns out, the order of this in the for loop matters a lot!!! Yikes. 

After changing the order in the loop, the results looked alright! 

Re number 60: laminar flow Re number 250: turbulent flow with vortices Once the flow gets to higher Re numbers, it starts looking a bit messy, but vortices are still visible. 

~~~ 

Next Steps: 

Now that the basic bounce-back method is working at showing the presence of walls at the top and bottom of the channel, I’ll look into other boundary conditions that I can implement that improve the accuracy of my results. So far, I’ve seen that the extrapolation scheme or half-way bounce back method is pretty popular, so I’ll probably read more into that.   

Thank you for reading! See you next week with some more updates