[MUSIC] Hi, welcome again to the fifth week of our class, modeling and simulation of natural processes. As I said before, boundary conditions in free dynamics are just important as the actual equation which we are solving for. Not only is it very important, it's also quite tricky to implement. And this is the topic of the present module. As a reminder, if we take an example of flow around an obstacle, there are different types of boundaries here. One type of boundary is the physical boundary around the surface of the obstacle. At this place, we need to tell the simulation that while there is no fluid inside the obstacle, there's only fluid everywhere else. And there's some condition the fluid needs to respect along the surface of the obstacle in order to represent the surface properties of the physical object which is placed there. We also have, additionally to the obstacle, which is the physical boundary, other types of boundaries which have no physical meaning. And which are just there to implement a constraint of numerical system like on the last in flow condition, on the right an outflow condition, and on the top and the bottom some other condition. We need, in summary, to have a way numerically speaking, to implement different types of boundary conditions. On the inflow we need a boundary condition which is able to impose a certain velocity profile, which is a way to tell the system at which velocity we want the fluid to flow around the obstacle. On the outflow, we need something which makes the outflow behave as discretely as possible so that the simulation can pretend that the actual simulation domain is much larger. And we don't see the effect of the cut restricted boundaries on the flow field. The same is also about for the upper and lower boundary. Finally, we need a boundary condition which has a physical meaning to implement the physics of the physical obstacle. It should be pointed out, that we already have about the condition implemented in the code you learned so far. If we do nothing, our code is periodic because it's periodicity what we implemented when we take the streaming step. As a reminder, periodicity means all fluid portions which are leaving the domain on the right will be reinjected on the left. And if they are leaving domain on the lower part, they are reinjected on the upper part. Periodicity can be a use of bound recognition. We're going to use it for the lower and the upper boundary. But you are not going to exploit this property for the inflow or the outflow, for which we are going to do something special. Let's start with the obstacle. Most obstacles, which we have in physics, have a so-called no-slip wall. iacroscopically speaking, this means that the flow velocity Is zero around the contour of the obstacle. If you look at the physics from a molecular point of view, it means that the molecule, in average, adheres somehow to the wall. It means that microscopically speaking, if you look at the wall under a microscope, you will see it has a certain surface roughness. And because of this roughness, molecules at average, will adhere to the wall and at average, adopt a zero velocity in the very small area surrounding either wall. The no-slip condition is not the right condition for all walls in physics. There are cases, like for example, an air bubble in water, which have a different condition. Because if you look at an air bubble under microscope, its surface is very smooth even at the molecular level. In this case, another condition applies. But today, we will just look at how to implement the most common condition for obstacles, no-slip walls in lattice Boltzmann. One numerical scheme to achieve this is the so called bounce-back bound recognition. The bounce-back bound recognition is quite easy to achieve. Look at what we do, if we take a cell, which is closed to an obstacle. On the image here we have one lattice cell, which has an obstacle placed right on top of it. During this training step, three of the population will leave the fluid domain and will stream into the off-cycle and need to be treated in a special way. What you're going to do is to following. A population which leaves the cell [INAUDIBLE] the wall, will just be reflected back and get back to the cell it came from. It's value will not change during this training step. We'll just take the same value, copy it back to the cell and invert its direction. That's all. This is how the bounce-back boundary condition works, for both orthogonal directions and diagonal directions. Once more, let's look at the population, which leaves the cell [INAUDIBLE]. It moves towards the wall. It is inverted. It is reflected back and gets back to the same cell, precolation state green, without any modification of the value of the population. The same thing is done for populations which move along diagonal directions. They move to the wall during streaming, are half-way through their motion, reflected back, inverted, and copied back to the cell they came from. Be careful, because this is not intrinsic behavior. You would expect that something which hits a wall at a certain angle, is reflected like on the mirror. This is not what happens here. We don't want reflection. We reflecting back to where they come from, and copy them to the same cell. This is what we needs to be done to implement a no-slip wall in lattice Boltzman. The bounce-back boundary condition is quite easy to implement. All you need to do is to identify all cells which are close to an obstacle and identify all the lattice directions, which link these fluid nodes to the nodes which are inside the obstacle. It's not difficult, but it involves quite some bookkeeping. You need to identify to find all the cells, which are closed in obstacle, and then go through all populations and see which ones are affected and need to undergo special treatment during the streaming step. It is not difficult, but it is a bit awkward. Fortunately, there is an easy solution to this. Instead of actually modifying this streaming step, we just stream all the populations into a solid node. That means we do allocate a memory for all cells which are inside the obstacle as well, and accept the populations are streamed from the fluid into solid nodes. Once they are in the solid node, instead of doing further correlations inside the solid, because solid is not a fluid, we just revert them inside the node. [INAUDIBLE] location, and that the next streaming step, we stream them back to where they came from. Let's see how this works in detail. Here we have two cells. One on the lower left is a fluid cell. The one of the upper right is a non-fluid cell, it's placed inside the obstacle. We are going to call this cell a bounce-back cell. Bounce-back cells, they don't do any collision sets. They don't do a BGK in collision. In collision, all they do is, take all the population and copy them to your pocket location. During streaming, they do the same thing as all our other cells. They take the population and copy them to the neighbors. So let's see in detail what happens if we focus on a population, which is the population with index eight, which is streamed from the lower left cell, the fluid cell into the solid cell that we're streaming,. It would just be copied without modification, onto the bounce back node, where it will reach the bounce back node in pre correlation green state. Then during correlation the bounce back node does nothing all populations, including the one it just received, and copy them to the reposit location. So collision takes the incoming populations, the green ones, maps them onto out populations, the red ones which are the value of the opposite population and the incoming state, they just revert all the populations back in opposite directions. And then in the next streaming set would happens naturally, is that population which [INAUDIBLE] as a came in from the lower left cell goes back to the cell it came from in the opposite direction. The advantage of doing so is that we don't need to do any book keeping, we don't need to remember which directions connect a fluid to a solid and which don't. We just go until the bounce-back nodes, the one inside the solid and divert all the populations we have there to collision without any distinction, everything is automatic. In this way, we have a neat and natural way to implement noted condition in the method. This is the formula which describes what happens on a bounce-back node instead of the collision term. During collision term, during collision bounce back nodes just take the outgoing population at a given direction, which we call j, and copy them into the opposite direction, which we call the direction with syntax i. So, Fi incoming X times T plus 1 is equal to Fj outgoing X times T. And I and J are linked together by saying that I is the opposite direction of J, or other words Vi is equal to minus Vj. As a Python code, this mechanism works as follows. We focus first on all nodes which are to implement a bounce-back approach. So we loop over all nine lattice directions, take the in-going populations, In opposite directions and copy them back to the algorithm populations in the given direction. We enumerated the lattice velocities in this example, in the In this way, you'll find formulating lattice algorithm, in such a way that if you go through the lattice velocities you will see that the opposite of a given velocity is always given by 8 minus the given velocity. Paid for example direction 0 your opposite direction is 8, take 1, your opposite direction is 7, your opposite of is 8 minus i. We have to make sure that these mechanisms only apply to bounce-back nodes and not to the other nodes. For this we use another matrix which we call the obstacle matrix. It is, this matrix, consists of zeros and ones. Zero is to label all positions where there's fluid, and one for all positions which are bounce-back nodes. We create this matrix beforehand and use it as a mask to restrict the operation which can formulate on this line. Two nodes which have the badly wanted the obstacle, which are obstacle nodes. NumPy accepts an array based on syntax based on matrices, on masks, to restrict dominance to which applies an operation. This ends the treatment of the obstacle of the conditions on the obstacle, and we now turn to the inflow condition, which is very important as well. The purpose of the inflow condition, being to impose, again, velocity or a velocity profile on all the cells which are placed along the left boundary of the domain. On the left boundary, we have a sort of an issue of the streaming, because after streaming, some of the populations are unknown. The populations f0, f1, and f2, and the incoming state, in the green state, are unknown, because they were supposed to be streamed from neighboring cells on the left, which do not exist because they are outside the numerical domain. Therefore, the boundary condition needs to find a way to assign them a value. We already know the velocity of these cells because the velocity is the value we want to impose on the boundary condition. But we don't know the density rho and we don't the three populations f0, f1, f2. Let's find the completions key to assign a reasonable value. Calculating density is quite easy, to do this, let's remember the formulas for calculating the macroscopic density and the velocity from the incoming populations. Density is just the sum of all populations. Velocity rather than momentum which is velocity times density is again the sum of all populations but awaited by a constant, which comes from the lattice velocities. We will now neglect the values of delta x and delta t we'll put them to 1 to switch to so called lattice units to make their equation simpler. All the coach are going to deadlock, will be not dis en genio. Let's split the sum calculate the density into three partial sums. The first partial sum is rho one, it is the sum of all three populations on the left column of f zero plus f one plus f two. Similarly, rho two with the partial sum in the middle, it is equal to f three + f4 + f5 and rho 3 is the partial sum on the right, it is equal to f6 + f7 + f8. Rho 1 is unknown, because the sum of the three populations, we don't know, whereas rho 2 and rho 3 are known. It is now easy to calculate the density from these three partial sums. The density is just the sum of row one plus row two plus row three. Momentum, which is density times velocity is not much more difficult. Let's focus on the x component on momentum, roe times UX. It is equal to the sum of the populations times the x component of the corresponding lattice velocity. For the right row for F6, F7 and F8, this component is just 1, because the x component of the lattice velocities. Excuse me, it is minus 1, because they are all pointing to the left. The x component of lattice velocities is negative 1. For the middle sum, for f3, f4, and f5, the x component of is just 0. So we neglect them. And for f0, f1, and f2, which are all pointing to the right, the x component of the velocity is plus 1, so the x component of momentum Rho times 2x = rho 1- rho 3. We can use these two equations, put them together to eliminate the unknown rho 1, and find an explicit value for the density rho. By subtracting these two equations, we get rho being equal to rho 2 + 2 rho 3 / 1- ux. By purely arithmetic terms, we find the right value of the density on the inflow boundary. Now, we also need to find a reasonable boundary to assign to the three populations, f0, f1, and f2. Remember that populations are just equal to equilibrium term plus a small perturbation. Therefore, as a first attempt, let's just initialize the equilibrium and say that f0 incoming is equal to equilibrium term for the zero direction. By now, we know the density and the velocity equilibrium only depends on macroscopic variables density and velocity. So we can compute the equilibrium term explicitly and assign it to f0. And the same for f1 and for f2. But let's get even more precise and correct the equilibrium term by the monitor by an appropriate perturbation. To get a good guess for this perturbation, we look at the opposite directions, see how much they deviate from their equilibrium, and copy the value of this deviation to the current populations. Example, f0 is first initialized to the equilibrium in 0 direction, then we go to the opposite direction f8. See how much it deviates from equilibrium, we take the difference between f8in minus the equilibrium term in direction 8. Take this difference and add it up to the f0 direction. In summary, the incoming of 0 population is assigned a value of equilibrium in direction 0, plus a difference between the incoming population in direction i8 minus equilibrium in direction 8. Do the same procedure for f1 and for f2, and we are done. Let's see what this looks like in terms of Python code, written with based syntax. To calculate the density, we need to first define the indices of the three columns. Column one having the indices 0, 1, 2. Column two, the indices 3, 4, 5. And column three, the indices 6, 7, 8. The density on the left wall, on the inflow, is equal is to 1 over 1- ux on the same row, on the same column times the sum of f3 plus f4 plus f5, which we call rho 2. In Python with the non-Python syntax we used the keyword sum to compute a sum, and I indicate with indices col2 to which population we want to restrict ourselves and indicators to parameter access equals zero. That the sum must be performed over the latest populations over the first index of the three-dimensional [INAUDIBLE] matrix. The same is done to calculate the sum over the column three. And the final formula is density being equal to 1 / 1- ux times rho 2 + 2 x rho 3. Once we have calculated density, we can use our completion scheme to assign a value to the three unknown populations f0, f1, and f2. They are equal to the equivalent term 0, 1 and 2. Plus the difference between the incoming populations 8, 7, 6, which are the opposite ones, minus the equilibrium at 8, 7 and 6. To end with, let's find an appropriate scheme] to implement the outflow condition. Remember that the outflow condition should behave as if the domain didn't end. We want to implement it in such a way that the outflow perturbs the system as little as possible. We'll do this with a very simple trick. On the outflow condition we have, again, three unknown populations, f6, f7, and f8. In this case we will not do any complicated calculations, we'll just take the same populations on the previous cell, which is located right on the left inside the fluid domain. Take its value, copy it to the unknown values on the outflow boundary, and we are done. By doing this simple copy, we will pretend that the fluid flow does not vary a lot any longer. Far from the obstacle on the outer boundary condition. It is like a zero gradient on the three populations close to the outflow column. In terms of Python, what that means is we take the incoming populations, along column 3, which is the indices 6, 7, and 8 on the previous column. Take their values, copy them again to the indices 6, 7, 8, which is column 3 along the outflow column, and we are done. By doing this very simple scheme, we get a visually very pleasing and also physically interesting way to implement an outflow condition. This ends our module on boundary conditions using the lattice Boltzmann method. Stay tuned for the next module. [MUSIC]