Fluid Dynamics and the Lattice Boltzmann Model

The lattice Boltzmann model (LBM) is a powerful fluid dynamics model which incorporates the physics of mesoscopic fluid “particles” propagating and col­liding 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 Current-voltage curve for a fuel cell calculated using a continuum mathemat­ical model of a fuel cell.

without us having to necessarily worry about solving mathematically compli­cated fluid dynamics models. In other words, the averaged long-wavelength properties of the lattice Boltzmann system obeys the desired Navier-Stokes equations, but rather than solving the Navier-Stokes 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 propa­gating 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 par­ticular, we can write the time evolution of the particle distribution function as

where f is a nine-dimension 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 next-nearest neighbors. This scheme is referred to as a D2Q9 model, because it is a two­dimensional model and there are 9 components to the particle distribution function. Other models exist, in two – and three-dimensions, although we are obviously restricted here to two-dimensions 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 post-collision 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 hy­drodynamic 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 com­ponents (trace components) of the nonequilibrium stress tensor separately. However, as in fluid dynamics it is typical to simulate systems in dimension­less numbers and map the simulation on to real systems using characteristic quantities (speed of sound in the fluid, Reynolds number, Womersley num­ber, 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 func­tion 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 Navier-Stokes 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 calcu­late 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 two-dimensional 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 no-slip boundary conditions around the edge of the fluid

FIGURE 7.13 The propagation and collision steps in the lattice Boltzmann model are de­picted. It is through the iteration of these propagation and collision steps that the fluid dynamics are captured in the model.

FIGURE 7.14 The bounce-back boundary conditions in the lattice Boltzmann model. The particle distribution function is propagating towards a solid boundary, and we want to bounce this back into the fluid. Therefore, we can turn it around and place it on the solid node. That way, when we propagate the particle distribution function it ends up back in the fluid domain.

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 bounce-back boundary conditions is the imple­mentation. 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 direc­tion, and place it on the solid node from which it is bouncing. Then we can perform the regular propagation step, moving the particle distribution func­tions to their neighboring nodes, and move the particle distribution functions from the solids back into the fluid domain. Mathematically, the bounce-back 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 to­wards 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 arrange­ment 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 perme­able electrodes. We can therefore apply the lattice Boltzmann model to the simulation of the fluid flow in these channels. While the real channels are three-dimensional 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 two-dimensions 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 sim­ulation domain is replicated multiple times across the simulation (the boxes), and different aspects of the simulations are stored in these regions. For exam­ple, 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 represent­ing 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 multi­ple 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 distri­bution functions have been propagated. The next column of boxes represent

FIGURE 7.17 The layout of the spreadsheet implementation of the lattice Boltzmann model.

the post-collision particle distribution functions. Recall, that in the current model the post-collision particle distribution functions are set to the equi­librium 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 bounce-back boundary conditions. In particular, the parti­cle 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 as­suming 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 left-hand 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 den­sity) 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 (correspond­ing to zero velocity and unit density) to the active part of the spreadsheet, the region where the current particle distribution function is stored (see Fig­ure 7.22). This allows us to initialize the system and start the simulation from scratch. The second button runs a macro which copies updated parti­cle 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 bot­tom 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 (no-slip 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

=IF(C13=1, 0,

(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 x-direction 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.000360416-1 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, О, (AA39-AA65+AA143-AA169+AA195-AA221))

FIGURE 7.20

Spreadsheet implementation of the lattice Boltzmann model. The velocity in the x-direction.

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 post-collision parti­cle 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 x-direction, while Figure 7.21 depicts the section of the spread­sheet which calculates the velocity in the y-direction. 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 x-direction.

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+AA143-AA169-AA195+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)

i

if the site is fluid (otherwise the velocity is set to zero). The code in cell C67 is

=IF(C13=1, 0, (AA39-AA65+AA143-AA169+AA195-AA221))

where only the components of the particle distribution function whose rep­resentative vectors ei have components in the x-direction play a role in de­termining the velocity in the x-direction. 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 y-direction, the code in cell C93 is

=IF(C13=1, 0, (AA91-AA117+AA143-AA169-AA195+AA221))

where only the components of the particle distribution function whose repre­sentative vectors ei have components in the y-direction play a role in deter­mining the velocity in the y-direction.

The current particle distribution function, within the spreadsheet model, is shown in Figure 7.22 (or at least the first component of the particle distri­bution 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 veloc­ity. In particular, in cell Z12 the particle distribution function f0 is calculated with the code

=IF(B12=1, 0, (0.44444*B40*(1-1.5*(B66*B66+B92*B92))))

Similarly, for the particle distribution function f1 (shown in Figure 7.23), the code in cell Z38 is

=0.1111111*B40*(1 + 3*(B66) +

4.5*(B66*B66) – 1.5 *(B66*B66+B92*B92))

Here, we are setting the particle distribution function at the boundary condi­tion 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*(1-1.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 post­collision particle distribution function for the component f0. Other compo­nents 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 post-collision 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 distribu­tion 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

=Z12

and these values remain unchanged.

Now that we have calculated the post-collision particle distribution func­tion, we can calculate the propagation of the particle distribution functions to the next lattice site in the direction they are oriented. This process of parti­cles 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 no-slip boundary conditions at the solid-fluid interface).

The idea behind this step of setting the particle distribution functions for bounce-back boundary conditions is to take the particle distribution function propagating towards a solid site, and turn it around so that when we prop­agate 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 post-collision 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 bound­ary 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

IF(AND(B12 =1, C12 = 0)

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

=BV39

which is the particle distribution function in the same direction, /1 , but from the lattice site to its immediate left (from the bounce-back boundary con­dition step). Similarly, particle distribution functions are propagated in the

The bounce-back 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 com­ponent 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 two­dimensional model implemented here is for instructional purposes, to demon­strate how simple and powerful this technique is. Furthermore, the lattice – Boltzmann can be used to capture multi-component fluids, turbulence, or flow through porous media and is a wonderful model for capturing the fluid dy­namics associated with fuel cells and many other alternative energy systems.

7.5 Exercises

1. The photovoltaic and electrochemical equivalent circuit model couples a solar cell with an electrochemical cell. However, there are a number of param­eters 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 current-voltage curve is to ohmic losses and limit current.

3. The current voltage curve can be obtained using the continuum electro­chemical model of a fuel cell. Obtain such a current-voltage 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 poros­ity). 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 three-dimensional modeling which would be very impractical using a spreadsheet.