The lattice Boltzmann model (LBM) is a powerful fluid dynamics model which incorporates the physics of mesoscopic fluid “particles” propagating and colliding on a simple lattice. This propagating and collision of “particles” mimics the interaction of atoms in a real fluid, but on a much larger scale. The upshot of this technique is that the correct fluid behavior emerges from the model
FIGURE 7.11 Currentvoltage curve for a fuel cell calculated using a continuum mathematical model of a fuel cell. 
without us having to necessarily worry about solving mathematically complicated fluid dynamics models. In other words, the averaged longwavelength properties of the lattice Boltzmann system obeys the desired NavierStokes equations, but rather than solving the NavierStokes equations directly, these properties emerge quite naturally from the local interactions within the model.
As mentioned, the power of this technique lies in the simple steps of propagating fluid “particles” to neighboring lattice sites, and colliding these particles when they reach a site. Here, fluid particles are representative of mesoscopic (length scale between that of atomistic and continuum systems) portions of the fluid, and are described by a particle distribution function. The particle distribution function is a continuous variable that exists in discrete directions (either going nowhere, in the 4 {10} directions, or in the 4 {11} directions of a simple square lattice) such that “particles” flowing in a given direction will reach a neighboring lattice site in a single time step. Upon reaching a lattice site the particle distribution functions are collided with each other, such that the resultant particle distribution has lower hydrodynamic stresses. In particular, we can write the time evolution of the particle distribution function as
where f is a ninedimension vector where each component of the vector is the particle distribution function in a given direction. The component, fi, is the
ith component of this vector, whose directions are depicted in Figure 7.12. The zeroth component goes nowhere, the directions labeled from one through four go to nearest neighbors, and directions five through eight go to nextnearest neighbors. This scheme is referred to as a D2Q9 model, because it is a twodimensional model and there are 9 components to the particle distribution function. Other models exist, in two – and threedimensions, although we are obviously restricted here to twodimensions if we want to solve the model in a spreadsheet. r is the position on the lattice, and e) is a vector which points to the neighboring lattice site in the ith direction. At is the time step in the simulation (usually unity), and t is the current time (usually dimensionless). The first part of the above equation,
fi(r + eiAt, t + At) = f*(r, t) (7.45)
simply takes the postcollision particle distribution function, f*, and moves it over to the next lattice site (r + eiAt) at the next time step (t + At). The second part of the evolution equation,
is a more general statement that the particle distribution function is modified upon collision by a collision operator, Oi. This collision operator depends on all of the components of the particle distribution function at this location and at this time. However, relaxing the system such that the hydrostatic stresses go to zero at each iteration is computationally more stable (which is what we will do here). It is worth noting that this limits the viscosity that can be simulated, and if it is required to change the viscosity of the fluid then the hydrodynamic stresses must be more slowly relaxed towards zero. Furthermore, if it is required to vary the shear and bulk viscosities separately then one can relax the deviatoric (components not along the trace) and hydrostatic components (trace components) of the nonequilibrium stress tensor separately. However, as in fluid dynamics it is typical to simulate systems in dimensionless numbers and map the simulation on to real systems using characteristic quantities (speed of sound in the fluid, Reynolds number, Womersley number, etc…) this is only really an issue if the simulation has different viscosities in different locations within the same simulation (multicomponent fluids or turbulent fluid flow, for example). We can, therefore, write down the typical update of the particle distribution function in the form
fi(r + ejAt, t + At) = f*(r, t) = feq (r, t) (7.47)
where fieq is the distribution which, for a given density and velocity of fluid, minimizes the hydrostatic stresses. The equilibrium particle distribution function for the system considered here is given by
The weights, the w’s, are given by w0 = 9, wi = w2 = w3 = w4 = 1, and w5 = w6 = w7 = w8 = jig. These values are obtained when the model is mapped onto fluid dynamic models (such as the NavierStokes equations), but we are not concerned about the derivation of the model here. p is the density of the fluid, and v is the fluids velocity. These hydrodynamic quantities (mass density p, and momentum density, pv) are moments of the distribution function and can be calculated from the particle distribution function
P = E fi (7.49)
i
and
pv = E fie (7.50)
i
In other words, the particle distribution function can be used to calculate the local density and velocity of the fluid. This in turn is used to calculate the equilibrium particle distribution function (minimizing hydrodynamic stresses). The equilibrium particle distribution function is then propagated to neighboring lattice sites. The resulting particle distribution functions are again used to calculate the local density and velocity of the fluid, and the cycle continues with the fluid dynamics being updated with each cycle of events or iteration. This is depicted in Figure 7.13.
The twodimensional regular square lattice of the lattice Boltzmann is represented in Figure 7.13. In the first step, the particle distribution functions are represented by the circles on the lattice nodes and their directions are represented by the arrows. Only the particle distribution functions pointing towards the central node are depicted, although in the actual simulation there are particle distribution functions pointing in all directions and at all nodes. In the second step, these particle distribution functions are propagated to the next site over, which is the central node (why these particle distribution functions are the only ones depicted). Notice that the size of these particle distribution functions can be quite random at this point. The third step is to minimize the hydrostatic stresses by setting the particle distribution function equal to the equilibrium particle distribution function (which still maintains the density and fluid velocity). This would usually result in a more ordered looking set of particle distribution functions. The final, fourth, step depicted is exactly the same as that depicted in the second step. The particle distribution functions are propagated to their neighboring nodes. It is really quite amazing that the correct fluid dynamics behavior emerges from this simple system where the particle distribution functions are propagated and collided with each other.
The lattice Boltzmann model is a powerful technique which is especially useful when modeling complex boundaries, and is very popular for simulating flows through porous media. This is because the lattice Boltzmann model can include solid regions in an efficient and incredibly simple way. In particular, we need to add noslip boundary conditions around the edge of the fluid


domain such that the velocity of the fluid at the boundary is equal to the velocity of the boundary (here we’ll take this to be zero, but moving solid interfaces can be included too). In practice, when the particle distribution function propagates from the fluid domain into a region which is solid, then the particle distribution function simply reflects, or bounces back, into the fluid domain. This is depicted in Figure 7.14.
The crucial step in the bounceback boundary conditions is the implementation. The particle distribution functions that are propagating towards solid nodes have to bounced back, such that after propagating the particle distribution functions, the particle distribution function is on the same node but pointing in the opposite direction. The first step is to take the particle distribution function and turn it around so it points in the opposite direction, and place it on the solid node from which it is bouncing. Then we can perform the regular propagation step, moving the particle distribution functions to their neighboring nodes, and move the particle distribution functions from the solids back into the fluid domain. Mathematically, the bounceback boundary conditions can be written in the form
fk (r, t + At) = f*(r, t)
FIGURE 7.15 The serpentine flow channels often seen in proton exchange membrane fuel cells. 
where k represents a direction which is opposite to the direction i. In other words, if the particle distribution function is traveling in the direction i towards the wall then this particle distribution function will be bounced back by the wall to where it came from, but it will now be traveling in the opposite direction, k.
The lattice Boltzmann model allows us to simulate the hydrodynamics of Newtonian fluids in a relatively easy manner, and consider complex geometries by simply assigning some nodes as solid areas and some nodes as fluid areas. Now let’s turn our attention to implementing this model within a spreadsheet to simulate the flow of gases in fuel cells.
In the fuel cells the gas enters and circulates through a serpentine arrangement of channels (at least in some designs) and this is depicted in Figure 7.15. The flow channels allow the gas to come into contact with the permeable electrodes. We can therefore apply the lattice Boltzmann model to the simulation of the fluid flow in these channels. While the real channels are threedimensional and the flow can change as the air navigates through the serpentine channels, here we will simulate a simpler configuration. Because we are solving this model using a spreadsheet we will limit the simulation to twodimensions and, just to see how the model works, we will only model the fluid flow as it goes around a corner. This is shown in Figure 7.16. Notice that we have to assign the density and velocity of the fluid at the boundary of the model. The density and velocity of the fluid within the simulation domain is then solved in response to these boundary conditions. Where the boundary is solid material, we can set the density and velocity to zero as these regions won’t be solved in the model. For the boundary sites at the entrance and exit of the section of the channel to be simulated we have to set the density (here we set it to unity) and velocity of the fluid. The fluid velocity is assumed to obey Poiseuille flow. In other words, the velocity of flow in the channel is of
where vmax is a constant that we can change, W is the width of the channel, and x is the distance from the center of the channel. In response to these boundary conditions we can now use a spreadsheet model to solve the fluid dynamics within the simulation domain.
The layout of the spreadsheet model is depicted in Figure 7.17. The simulation domain is replicated multiple times across the simulation (the boxes), and different aspects of the simulations are stored in these regions. For example, on the far left of the spreadsheet the top box represents an area of the spreadsheet which contains either 0’s or 1’s depending on whether the lattice site is fluid or solid, respectively. Below this are three boxes, each representing the same simulation domain, which store the fluid density and the velocity components. Even though the simulation size is relatively small (here we only capture 22 x 22 lattice sites within the simulation domain) we need multiple simulation domains in the spreadsheet, each storing different parameters within the model.
As we move across the spreadsheet, particle distribution functions are stored. The particle distribution function is a 9 component vector and so the different components are stored as we go down the spreadsheet. The first column of boxes representing the particle distribution function is the current particle distribution function in the simulation, just after the particle distribution functions have been propagated. The next column of boxes represent

the postcollision particle distribution functions. Recall, that in the current model the postcollision particle distribution functions are set to the equilibrium particle distribution functions, which only depends on the current density and velocity within the fluid. The third column of boxes representing the particle distribution functions modify the particle distribution functions to account for the bounceback boundary conditions. In particular, the particle distributions which are propagating towards a solid site are placed on that solid site, but directed back in the opposite direction towards the fluid. In the next column of boxes, the particle distribution functions are propagated to the neighboring sites in the direction of the particle distribution functions. This concludes the iteration in the model and these values have to be copied back to the first column of boxes which represents the current particle distribution functions, so that the iteration can be completed and the spreadsheet will be repopulated with calculations for the next iteration. There are two additional columns of particle distribution functions in the spreadsheet to the far right. The first is simply the initial values of the particle distribution function assuming no fluid velocity (these can be copied to the first column to reset the simulation back to the beginning). The last is a copy of the propagated particle distribution functions. Rather than copy the propagated particle distribution functions, one at a time, to the beginning of the spreadsheet to complete the iteration, the values are first stored here and then copied to the beginning. In this way all of the components of the particle distribution functions are copied across all at once. Without this step, the lower number components might be copied across and the later components would be recalculated based on newer lower components of the particle distribution functions, before being copied across. At the top of the simulation we will place two buttons which will either copy the initial values of the particle distribution functions across to the lefthand side of the simulation and reset the simulation, or copy the latest particle distribution functions across and iterate the simulation. Let’s look at the actual spreadsheet implementation, with snapshots in Figures 7.17 to 7.27 taken from the spreadsheet.
Figure 7.18 depicts the top of the simulation. The maximum velocity in the center of the channels for the boundary conditions (fixed velocity and density) is assigned in cell D1. Interestingly, we need no other constants in this dimensionless model. Two buttons are included at the top of the spreadsheet. The first button runs a macro which initializes the system, and copies particle distribution functions from another section of the spreadsheet (corresponding to zero velocity and unit density) to the active part of the spreadsheet, the region where the current particle distribution function is stored (see Figure 7.22). This allows us to initialize the system and start the simulation from scratch. The second button runs a macro which copies updated particle distribution functions (subsequent to propagation and collision) back to the beginning of the iteration. In other words, the current particle distribution functions are contained within the spreadsheet, and used to obtain an updated particle distribution function. By copying this updated particle distribution
FIGURE 7.19 Spreadsheet implementation of the lattice Boltzmann model. The fluid density. 
function back to overwrite the current particle distribution functions, then the simulation has progressed an iteration, and the spreadsheet is automatically populated with the next updated particle distribution functions. At the bottom of Figure 7.18 is 22 x 22 size region corresponding with the simulation domain size, where cells that contain the number 1 are designated solid sites and cells that contain the number 0 are designated fluid sites. Recall that we solve the fluid dynamics in the fluid regions and have implemented fluid bounce back boundary conditions (noslip boundary conditions) at the fluid solid interface.
The density in the simulation is depicted in Figure 7.19. Around the edge of the simulation domain the density is assigned as part of the boundary conditions, as is either 0 if the site is solid or 1 if the site is fluid. Inside the domain, the density is determined from the current particle distribution function. In particular, the code in cell C41 is
(AA13+AA39+AA65+AA91+AA117+AA143+AA169+AA195+AA221))
In other words, if the site is solid then just put the density to zero. However,
Velocity in xdirection is set to a given value in the fluid regions of the boundary, around the edge of the simulation domain. This establishes Poiseuille flow at the boundaries. Г
A 
В 
C 
D 
E b F 
G 
H 
1 
J 

64 
VX 

65 
1 
2 
3 
4 V 
6 
7 
8 
9 

66 
1 
0 
0 
0 
0 0 
0 
0 
0 
0 
67 
2 
0 
0 
0 
0 0 
0 
0 
0 
0 
68 
3 
0 
0 
0 
0 0 
0 
0 
0 
0 
69 
4 
0 
0 
0 
0 0 
0 
0 
0 
0 
70 
5 
0 
0 
0 
0 0 
0 
0 
0 
0 
71 
6 
0 
0 
0 
0 0 
0 
0 
0 
0 
72 
7 
0 
0 
0 
0 0 
0 
0 
0 
0 
73 
8 
0 
0 
0 
0 0 
0 
0 
0 
0 
74 
9 
0 
0 
0 
0 0 
0 
0 
0 
0 
75 
10 
0 
0 
0 
0 0 
0 
0 
0 
0 
76 
11 
0 
0 
0 
0 0 
0 
0 
0 
0 
77 
12 
0.000475 
0.000280487 
0.000327523 
0.000352163 0 000365797 
0.000370323 
0.000377376 
0.000386017 
0.000417597 
78 
13 
0.001275 
0.000940408 
0.000959451 
0.000978272 0.000993261 
0 001001876 
0.001016402 
0 001038644 
0.00109166 
79 
14 
0.001875 
0 001502338 
0.001488616 
0.001483939 0.00148668 
0.001490377 
0.001503776 
0.001524241 
0 001565293 
80 
15 
0.002275 
0 001887601 
0.001861127 
0 001839251 0.001828771 
0 001823375 
0 001828635 
0.001836099 
0.00185098 
81 
16 
0.002475 
0.002082959 
0 002052479 
0.002023175 0.002004986 
0.001992106 
0 00198813 
0.001980662 
0.001971573 
82 
17 
0.002475 
0 002084376 
0 002054935 
0 002026454 0 002008235 
0 001993287 
0.001983982 
0.001966216 
0 001941566 
83 
18 
0.002275 
0.001891621 
0 001868156 
0 001848814 0 001838743 
0 001828346 
0 001819646 
0.001798643 
0.001768323 
84 85 86 
19 19 21 
0.001875 0.001275 0.000475 
0 001508166 0.00094646 0.000283881 
0.001498905 0.000969986 0.000333047 
0.001498386 Аз.001503025 0.000993609/ 1001012282 0.0003604161 r000377099 
0 001501858 0.001019222 0.000382627 
0.001498478 0.00102193 0.000387147 
0 001480583 0.001010694 0.000382477 
0.001453494 0.000993393 0.000378441 
87 
22 
0 
0 
0 
0 
0 
0 
0 
0 
However, inside the simulation domain the velocity is obtained from the second
moment of the particle distribution function (see text). The code in cell C67 is
=IF(C13=1, О, (AA39AA65+AA143AA169+AA195AA221))
FIGURE 7.20
Spreadsheet implementation of the lattice Boltzmann model. The velocity in the xdirection.
if the site is fluid then the density is given by
P = E fi (7.53)
i
the first moment of the particle distribution function. As we relax the particle distribution function directly to its equilibrium state, the postcollision particle distribution function is calculated directly from this density and the fluid velocity, which is depicted in Figure 7.20 and 7.21.
Figure 7.20 depicts the section of the spreadsheet which calculates the velocity in the xdirection, while Figure 7.21 depicts the section of the spreadsheet which calculates the velocity in the ydirection. Similar to the density, the velocity at the boundary conditions is assigned at the beginning of the simulation and is fixed for the entire duration of the simulation. In particular, if the site is solid the velocity of the fluid is set to zero but if the site is fluid then the fluid is set according to Poiseuille’s flow. In other words, the velocity in the direction parallel to the channel is given by
u = umax (R2 – x2) (7.54)
where umax is the maximum velocity in the center of the channel (cell D1), R is
Boundary conditions dictate the velocity m the уdirection m exactly the same way
as in the xdirection.
In response to the prescribed velocities at the boundary, the fluids velocity m the simulation domain are obtained dynamically. The current velocity, calculated from the particle distribution function is obtained. The code m cell C93 is
=IF(C13=1, 0, (AA91 AA117+AA143AA169AA195+AA221))
the radius of the channel, and x is the distance from the center of the channel. The velocity in the active region of the simulation domain is calculated from the second moment of the particle distribution function
pv = fi*ei (7.55)
if the site is fluid (otherwise the velocity is set to zero). The code in cell C67 is
=IF(C13=1, 0, (AA39AA65+AA143AA169+AA195AA221))
where only the components of the particle distribution function whose representative vectors ei have components in the xdirection play a role in determining the velocity in the xdirection. In particle if ei is oriented with a component to the right then the particle distribution function is added, while if ei is oriented with a component to the left then the particle distribution function is subtracted. Similarly, in the ydirection, the code in cell C93 is
=IF(C13=1, 0, (AA91AA117+AA143AA169AA195+AA221))
where only the components of the particle distribution function whose representative vectors ei have components in the ydirection play a role in determining the velocity in the ydirection.
The current particle distribution function, within the spreadsheet model, is shown in Figure 7.22 (or at least the first component of the particle distribution function, fo, with the remaining eight components in the cells below this). The particle distribution function for the boundary nodes do not change during the simulation and are calculated directly from the density and velocity. In particular, in cell Z12 the particle distribution function f0 is calculated with the code
=IF(B12=1, 0, (0.44444*B40*(11.5*(B66*B66+B92*B92))))
Similarly, for the particle distribution function f1 (shown in Figure 7.23), the code in cell Z38 is
4.5*(B66*B66) – 1.5 *(B66*B66+B92*B92))
Here, we are setting the particle distribution function at the boundary condition to the equilibrium distribution based on the assigned boundary conditions (constant density and velocity). Recall that the equation for calculating the equilibrium particle distribution function, for the current D2Q9 model, is
where the weights, the w’s, are given by wo = 9, wi = W2 = W3 = W4 = 1
The boundary values of the particle distribution function are obtained from the density and velocity values, which are the fixed boundary conditions. The code in cell Z12 is
=IF(B12=1, 0, (0.44444*B40*(11.5*(B66*B66+B92*B92))))
The values of the particle distribution function m the simulation domain are later calculated in the spreadsheet, and then their values are copy and pasted back here as part of an iteration.
The particle distribution function for the central simulation domain is again copied from later calculations in the spreadsheet.
Y 
Z 
AA 
AB 
AC 
AD 
AE 
AF 
AG 
AH 

36 
F 1 

37 
1 
2 
4_______ з 
4 
5 
6 
7 
8 
9 

38 
1 
0 
0 
Д 
7___ о 
0 
0 
0 
0 
0 
0 
39 
2 
0 
0 
0 
0 
0 
0 
0 
0 
0 

40 
3 
0 
0 
0 
0 
0 
0 
0 
0 
0 

41 
4 
0 
0 
0 
0 
0 
0 
0 
0 
0 

42 
5 
0 
0 
0 
0 
0 
0 
0 
0 
o 

43 
6 
0 
0 
0 
0 
0 
0 
0 
0 
0 

44 
7 
0 
0 
0 
0 
0 
0 
0 
0 
0 

45 
8 
0 
0 
0 
0 
0 
0 
0 
0 
0 

46 
9 
0 
0 
0 
0 
0 
0 
0 
0 
0 

47 
10 
0 
0 
0 
0 
0 
0 
0 
0 
0 

48 
11 
0 
0 
0 
0 
0 
0 
0 
0 
0 

49 
12 
0.111269509 
0.111269509 
0.111301741 
0.111317662 
0.111321606 
0.111318493 
0.111313681 
0.111307575 
0.111305046 

50 
13 
0.111536642 
0.111536642 
0.11152904 
0.111530039 
0.111530594 
0.111527604 
0.111523924 
0.111519851 
0.111520886 

51 
14 
0.111737272 
0.111737272 
0.111722321 
0.11170916 
0.11170033 
0.111692488 
0.111686702 
0.111681757 
0.111681623 

52 
15 
0.111871158 
0.111871158 
0.111854345 
0.111835431 
0.111819842 
0.11180692 
0.111797676 
0.111789616 
0.111784851 

53 
16 
0.111938142 
0.111938142 
0.111921024 
0.11190018 
0.111881618 
0.111865719 
0.111853714 
0.111842303 
0.111832439 

54 
17 
0.111938142 
0.111938142 
0.111921232 
0.111900708 
0.111882384 
0.111866367 
0.111853595 
0.111840258 
0.111826934 

55 
18 
0.1118, 
71158 
0.111871158 
0.111854938 
0.111836941 
0.111822082 
0.111808963 
0.111797795 
0.111784596 
0.111770192 
56 
19 
0.1117 V7272 
0.111737272 
0.111723171 
0.111711338 
0.111703683 
0.111695893 
0.111687983 
0.111676514 
0.111663107 

57 
20 
0.11V V642 
0.111536642 
0.11152983 
0.111532083 
0.11153392 
0.111531325 
0.111526247 
0.11151669 0.111505352 

58 
21 
0.1111 
9509 
0.111269509 
0.111301754 
0.111318087 
0.111322553 
0.11131941 
0.111313709 
0.111304432 
0.111294972 
59 
22 
0 
0 
0 
0 
0 
0 
0 
0 0 

The particle distribution function for the boundary is again calculated from the fixed density and velocity profile at the boundaries. The calculation here is for fi which is propagating to the right. The code in cell Z38 is
=0.1111111*B40*(1 + 3*(B66) + 4.5*(B66*B66) – 1.5 *(B66*B66+B92*B92))
and w5 = w6 = w7 = w8 = 36 .In this way the current boundary values of the particle distribution functions ensure that the velocity and density are fixed at the boundaries.
The current values of the particle distribution function in the active region of the simulation, as depicted in Figures 7.22 and 7.23, are not calculated here. The current iteration particle distribution functions are copied and pasted from the calculated values from the previous iteration (the values only are copied and pasted, so there is no code in these cells). As mentioned earlier, the current iteration particle distribution functions are used to calculate the next iteration elsewhere in the spreadsheet, and when the updated values are copied and pasted here the simulation is iterated. Let’s turn our attention to how we calculate the updated particle distribution functions for the next iteration.
Figure 7.24 depicts the part of the spreadsheet calculating the postcollision particle distribution function for the component f0. Other components are calculated in the spreadsheet below this. Recall that in the current model we are relaxing the stresses directly to local equilibrium (for simplicity and stability) and the postcollision particle distribution function is simply the equilibrium particle distribution function. The code in cell AY13 is
=IF(C13=1, 0, (0.444444*C41*(1 – 1.5*(C67*C67 + C93*C93))))
where if the site is solid then the particle distribution function is zero, else it is calculated from the equation above for the equilibrium particle distribution function. The boundary condition particle distribution functions do not change and, therefore, the values are simply copied across from the current particle distribution function. In other words, the code in cell AX12 is
and these values remain unchanged.
Now that we have calculated the postcollision particle distribution function, we can calculate the propagation of the particle distribution functions to the next lattice site in the direction they are oriented. This process of particles propagating along in space and then colliding with each other mimics the microscopic fluid dynamics and is at the heart of why the lattice Boltzmann is such a powerful model. However, the lattice Boltzmann’s main strength is the ease at which complex boundaries can be easily implemented, and before we propagate the particle distribution functions we must prepare the bounce – back boundary conditions (the noslip boundary conditions at the solidfluid interface).
The idea behind this step of setting the particle distribution functions for bounceback boundary conditions is to take the particle distribution function propagating towards a solid site, and turn it around so that when we propagate the particle distribution functions it will be propagated back to where it came from, but be facing in the opposite direction. Obviously the particle distribution function components f0 don’t propagate as the zeroth direction
The postcollision particle distribution function is calculated from the density and velocity. The density and velocity are, of course, calculated from the particle distribution function. This relaxes the stresses in the fluid toward (or actually to in the current model) local equilibrium. The code in AY13 is =IF(C13=1, 0, (0.444444*C41 *(1 – 1.5*(C67*C67 + C93*C93)))) The distributions at the boundaries remain unchanged. The code in cell AX 12 is simply =Z 12 
is the rest direction (or no direction). Figure 7.25, therefore, shows the section of the spreadsheet which sets up the bounce back boundary conditions for the /1 components. The code in cell BV38 is
=IF(AND(B12 = 1, C12 = 0), AY64, AX38)
Note that the same code is used for sites regardless of whether they are boundary or active regions, as we want the sites at the boundaries (which are solid) to store the particle distribution functions that will be propagated back into the fluid domain in the next step. However, if we are looking in the right direction, or [10] direction, then the last sites on the right of the simulation won’t need to be calculated (there is obviously nothing to the right of these sites which could propagate from the left). The code in cell BV38 starts with
which is saying if the current site is solid and the site to the right is fluid, then there will be a particle distribution function propagating from the right (in the left direction) which needs to be turned around and propagated back out to the right. If this is the case then find the particle distribution function component to the left /2 on the lattice site to the right of the current site, and set the particle distribution function /1 at the current site equal to this. This way when we propagate it back into the fluid it will end up on the lattice site to the right of the current lattice site but in the opposite direction (the e1 direction from the e2 direction). We do a similar operation, always switching the directions of the particle distribution functions propagating into the solid sites, for all components of the particle distribution function. In other words, at the end of this step, we have particle distribution functions at all other solid sites adjacent to fluid sites, ready to propagate the required particle distribution functions back into the fluid.
The final step in the iteration is to propagate the particle distribution functions in the directions towards the neighboring lattice sites. (In other words, propagate the /’s in the directions of the e vectors to the next lattice site.) Figure 7.27 depicts the sections of the spreadsheet which calculate the propagation of the particle distribution function in the e1 direction. Note that only the particle distribution functions in the active region are required here, as this is what we will be copying and pasting back to the current particle distribution functions as depicted in Figure 7.22. The propagation step is really quite simple, as all we have to do is move the particle distribution functions over. So in the e1 direction the particle distribution functions are copied to the right. The code in cell CU39 is
which is the particle distribution function in the same direction, /1 , but from the lattice site to its immediate left (from the bounceback boundary condition step). Similarly, particle distribution functions are propagated in the
The bounceback boundary conditions are implemented at the interface between the fluid and the solid. Before we propagate the particle distribution functions we modify the distributions at the interface. In particular, if the distribution is going to be propagated into a solid, then place the distribution on the solid node facing back towards the fluid. The code in cell BV38 is
=IF(AND(B12 = 1, C12 = 0), AY64, AX38)
This sets the distribution up, ready to be propagated at the next step back into the fluid. /ч
other directions to complete the iteration. We now have to copy the values of the propagated particle distribution functions back to the current particle distribution function (see Figure 7.22). However, if we simply copy each component of the particle distribution function in turn, then some components will be automatically updating while the first components to be copied across are copied across. For this reason, we copy the values for all of the particle distribution functions to a separate location in the spreadsheet (all the way over to the right) and then copy these values one at a time into the current particle distribution functions. The act of copying and pasting all of these values within the spreadsheet is stored in a macro that can be run at the push of a button, and each time the button is pushed the simulation iterates. The lattice Boltzmann model is a dynamics model and so we can see how, in response to the boundary conditions, the flow field within the simulation domain evolves. Figure 7.26 shows the fluid flow within the lattice Boltzmann simulation described above. The fluid velocity is depicted as arrows, where the size of the arrows indicate magnitude and the direction of the arrows the direction of fluid flow. The fluid velocity is assigned at the boundaries and then the lattice Boltzmann model is used to calculate the fluid velocity inside the fluid domain of the active region of the simulations (the region where the fluid is updated).
Typically, the lattice Boltzmann model might be solved in three dimensions and incorporate a larger simulation domain. The small twodimensional model implemented here is for instructional purposes, to demonstrate how simple and powerful this technique is. Furthermore, the lattice – Boltzmann can be used to capture multicomponent fluids, turbulence, or flow through porous media and is a wonderful model for capturing the fluid dynamics associated with fuel cells and many other alternative energy systems.
1. The photovoltaic and electrochemical equivalent circuit model couples a solar cell with an electrochemical cell. However, there are a number of parameters which can depend on the the exact system. Vary the parameters (within realistic ranges) and maximize the efficiency of these devices.
2. Implement the simple electrochemical model of a fuel cell. Determine how sensitive the currentvoltage curve is to ohmic losses and limit current.
3. The current voltage curve can be obtained using the continuum electrochemical model of a fuel cell. Obtain such a currentvoltage curve and fit this to the simple electrochemical model of a fuel cell. For example, to get a good
FIGURE 7.26
Fluid flow around a corner in a serpentine fuel cell channel, calculated using the lattice Boltzmann model.
The particle distribution function is now propagated to the next lattice site over. In particular, the distribution is propagated to the right, and the new distribution in cell CU39 is simply
=BV39
which is the same distribution fl but from the lattice site to its immediate left.
CS 
CT 
CU 
CV 
CW 
CX 
CY 
CZ 
DA 
DB 

36 
F 1 

37 
1 
2 
3 
4 
5 
6 
7 
8 
9 

38 
1 
V 

39 
2 
0 
0 
0 
0 
0 
0 
0 
0 

40 
3 
0 
0 
0 
0 
0 
0 
0 
0 

41 
4 
0 
0 
0 
0 
0 
□ 
0 
0 

42 
5 
0 
0 
0 
0 
0 
0 
0 
0 

43 
6 
0 
0 
0 
0 
0 
0 
0 
0 

44 
7 
0 
0 
0 
0 
0 
0 
0 
0 

43 
6 
0 
0 
0 
0 
0 
0 
0 
0 

46 
9 
0 
0 
0 
0 
0 
0 
0 
0 

47 
10 
0 
0 
0 
0 
0 
0 
0 
0 

48 
11 
0 
0 
0 
0 
0 
0 
0 
0 

49 
12 
0.111269509 
0.111301606 
0.111317863 
0.111321207 
0.111318891 
0.111313049 
0.11130815 
0.111304215 

50 
13 
0.111536642 
0.111528933 
0.111530251 
0.111530235 
0.111528016 
0.111523329 
0.111520439 
0.111520079 

51 
14 
0.111737272 
0.111722208 
0.111709393 
0.111699985 
0.111692933 
0.111686128 
0.111682392 
0.111680841 

52 
15 
0.111871155 
0.111854234 
0.111835677 
0.111819504 
0.111807396 
0.111797116 
0.111790298 
0.111784085 

53 
16 
0.111938142 
0.111920911 
0.111900438 
0.111881283 
0.111866223 
0.111853157 
0.111843025 
0.111831672 

54 
17 
0.111938142 
0.111921118 
0.111900976 
0.111882045 
0.111866891 
0.111853028 
0.111841007 
0.111826147 

55 
16 
0.111871158 
0.111854821 
0.111837216 
0.111821732 
0.111809496 
0.111797206 
0.111785357 
0.111769371 

56 
19 
0.111737272 
0.111723046 
0.111711614 
0.111703314 
0.111696425 
0.111687361 
0.111677271 
0.111662239 

57 
20 
0.111536642 
0.111529708 
0.11153234 
0.111533523 
0.111531834 
0.111525581 
0.111517424 
0.111504429 

58 
21 
0.111269509 
0.111301589 
0.111318323 
0.111322088 
0.111319895 
0.111312971 
0.111305147 
0.111293981 

59 
22 
FIGURE 7.27
Spreadsheet implementation of the lattice Boltzmann model. Propagation of the particle distribution function.
fit you could change the phenomonological constants (£1, £2, £3, and £4), the ohmic losses (RS), or the limit current (1L).
4. The nature of fluid flow through porous membranes might also be modeled using a lattice Boltzmann model. For example, when fluid travels towards a solid wall in the lattice Boltzmann model, it is entirely bounced back to where it came from by the solid wall. However, in porous media some of the fluid would continue on its journey. In particular, fluid flow through porous media can be represented using Darcy’s law which states that the flow rate of a fluid depends on the porousity of the structure it is flowing through. In the lattice Boltzmann model this can be incorporated by allowing some of the particle distribution function to propagate through the porous media (equal to the porosity) while some of the fluid bounces back (equal to one minus the porosity). Simulate fluid flow though a porous wall. How easily can the fluid move through this porous wall?
5. Could the lattice Boltzmann model be used to look at the aerodynamics of air flow in wind turbine studies, either as a way of looking at the air flow over land, such that the best locations could be identified for wind turbines or through the actual simulation of air flow over the wind turbine blades? This would presumably require threedimensional modeling which would be very impractical using a spreadsheet.