5.5. Lees Edwards BCs

5.5.1. Background: Sheared molecular dynamics

Lees and Edwards [LeesEdwards1972] introduced the idea of sliding periodic boundary conditions in the context of molcular dynamics as a mechanism to impart shear. The algorithm has a number of advantages over other methods: it allows, in principle, high shear rates and it does not introduce any solid boundary which would disrupt the linear shear regime with boundary effects.

The computational picture is illustrated schematically below. A central computational system \(CS\) (shaded) is surrounded by periodic, or cyclic, images \(CS'\) (copies of which extend infinitely in all directions; only the immediate neighbours of \(CS\) are shown). One imagines that the first row of images above the central system slide at a velocity \(+U_y\) relatative to \(CS\), a second row above would slide at a velocity \(+2U_y\) relative to \(CS\), and so on. Likewise, successive rows below slide at \(-U_y,-2U_y\) and so on relative to the central system.

Schematic illustration of a shearing periodic system

We assume that at time \(t = 0\) the system is aligned in the standard periodic configuration where there is no relative displacement of the rows. In general, a particle leaving the top of the computational system \(CS\) with coordinate \(y\) would enter the adjacent image at a corresponding position \(y' = y - U_yt\). The corresponding coordinate in the true computational system is then \(y'' = y - 2U_yt\). Displaced positions \(y + U_y t\) may always be reduced by modulo \(L_y\) to return the relevant particle coordinate to the true computational system.

The \(y\)-component of the velocity of a particle leaving the top of the central system is transformed via \(v_y' = v_y - U_y\). The corresponding component of velocity of a particle at the lower boundary of the true system would be \(v_y'' = v_y - 2U_y\).

For a particle leaving the bottom of the central system the sign of \(U_y\) in the corresponding transformations is reversed.

5.5.2. Background: Shear for lattice Boltzmann

Wagner and Pagonabarraga [WagnerPagonabarraga2002] were the first to introduce the idea of Lees Edwards sliding periodic boundaries in the lattice Boltzmann picture. The idea here is slightly different to the case of molecular dynamics. We consider a single computational system which is sub-divided into one or more blocks by the introduction of internal Lees Edwards boundaries. (The boundaries we often refer to as Lees Edwards “planes”.) In the the illustration below there are four boundaries, or planes, with positions \(x_1, \ldots, x_4\).

Schematic illustration for the multi-block picture of Wagner and Paganobarraga

At each boundary, a notional particle crossing in the positive \(x\)-direction undergoes the Galilean transformation \(v_y' = v_y + U_y\), with corresponding transformation \(v_y = v_y' - U_y\) if travelling in the negative \(x\)-direction. There is an analogous position transformation as a function of time of \(y' = y \pm U_yt\). So the picture is one in which we have discrete blocks which are - conceptually - in relative motion, with the central block considered a rest frame.

The net result is a single system in which we can impose an overall shear rate \(\dot{\gamma} = NU_y/L_x\), where \(N\) is the number of planes. This allows a larger shear rate than would be possible with a single set of sliding periodic boundaries (as in the molecular dynamics case) as the largest speed relative to the lattice in any block is restricted to \(v_y \sim U_y\). This respects the low Mach number constraint of lattice Boltzmann.

There are a number of points to note:

  1. We retain the picture in which \(x\) is the velocity gradient, or shear, direction, and \(y\) is the flow direction. While not shown in the illustrations, the \(z\)-direction is the voriticity direction. As there are no additional special considerations for the vorticity direction, it is easiest to use two-dimensional examples.

  2. There is a “standard” periodic boundary in the \(x\)-direction, and in the other two dimensions.

The picture for an odd number of planes is slightly different, and is not discussed here; we recommend sticking to an even number of planes. (There is no particular reason to use an odd number of planes, and it is easier to arrange the relevant spacing and parallel decomposition with an even number.)

5.5.3. Practicalities: generating and analysing results

The input for the example discussed below is found in the repository in the docs/tutorial/lees-edwards in the repository along with some sample output. The code has been compiled with -D_D2Q9_.

See also the Lees Edwards SPBCs page for details of input key value pairs for Lees Edwards boundaries.

5.5.3.1. Raw results

The picture above is one of a single system with a number of blocks sliding relative to one another in the flow direction. In practical implementation, of course, data do not move in memory, which is fixed. We can see this if we examine the raw output, where we expect to see a “displaced” picture as a function of \(U_y t\).

For the purposes of illustration we consider the simple example of a binary droplet in two dimensions. The droplet is initialised in the centre of a system with two Lees Edwards planes having \(U_y = 0.005\) (the shear rate \(\dot{\gamma} = 2U_y/L_x\)). We show below the unprocessed output for times \(t = 0\), and \(t = 3200\), the latter corresponding to \(U_yt = 16.0\). We’ve adopted the same coordinate orientation as in the schematics above.

Unprocessed results from a simulation with two planes

It can be seen that the planes (two horizontal lines) separarate different blocks in the computational system. At first sight, the displacement of the structure may not be what is expected; one must remember that, conceptually, the whole of the upper block is translating to the right, and the whole of the lower block to the left.

If the displacement were a whole number of system lengths \(L_y\), the picture would actually be the correct one. However, to be useful in the general case, we need to adjust the raw output to allow for the relevant block displacements. This process is referred to as “unrolling”, and can be performed using a utility provided.

5.5.3.2. “Unrolling” the output

Information about the system in the meta data file includes the details of the Lees Edwards planes. For the binary droplet case, we would look for phi-metadata.001-001 which is produced automatically at run time. In out illustrative case, we should see the JSON section:

"lees_edwards": {
        "Number of planes":     2,
        "Shear type":   "STEADY",
        "Reference time":       0,
        "Plane speed":  0.005
}

This, along with the time step (which is encoded in individual file names) allows the appropriate displacement \(U_y t\) to be made by the extraction utility.

To perform the unrolling for the raw output phi-000003200.001-001 one would run, e.g.,

$ ./extract phi-000003200.001-001
...
Writing result to phi-000003200
Completed processing for phi-000003200.001-001

The utility will identify this file as the scalar order parameter field for composition from the phi- stub, and look for the corresponding metadata file phi-metadata.001-001. Files for different time steps (also identified from the file name) share the same metadata file. The unrolled file will be in the same format as the original, and can the viewed in the same way.

Unrolled Lees-Edwards output

There are a number of points to note about this process:

  1. If \(U_y t\) is not an integer (the general case) data are interpolated using cubic interpolation to bring values back to the lattice points.

  2. The unrolled data should be used for analysis only; the original raw data is always used for restart purposes.

5.5.3.3. The velocity field in the presence of planes

A similar discussion is relevant for the velocity field, with the additional consideration that the component of the velocity in the flow direction must be corrected for the relative motion of the blocks. This correction is made automatically by the extraction utility if it detects that the input is the velocity field.

$ ./extract vel-0000032000.001-001
...
Data is velocity field
Writing result to vel-000003200
Complete processing for vel-000003200.001-001

To illustrate, we plot a single profile of the fluid velocity component \(u_y(x,L_y)\) using the same binary droplet example, but use a position \(y = L_y\) away from the droplet. This is shown below for both (left) the uncorrected values and (right) the corrected values. The uncorrected values are the components relative to the lattice, whereas the corrected values are the view from the central rest frame.

Flow at the lattice level and corrected for two Lees Edwards boundaries.

In a fully-developed flow in a simple fluid, this would be a linear shear profile.

5.5.3.4. Choice of rest frame

In the example above, the data are unrolled with respect to the central block. This gives rise to a picture in which the viewer sees the flow moving to the right in the upper half of the system, and to the left in the lower half of the system.

The choice of the rest frame is actually arbitrary. One could chose a different block in the computational system to be the rest frame, e.g., the bottommost block at \(x = 0\). This gives rise to a picture in with the viewer sees the flow stationary at the bottom and moving with \(u_y \simeq NU_y\) at the top. Experience suggests this can be slightly “unsettling” for the viewer. We therefore always use the central block as the rest frame.

[LeesEdwards1972]

A.W. Lees and S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C: Solid State Phys., 5 1921-1929 (1972).

[WagnerPagonabarraga2002]

A.J. Wagner and I. Pagaonabarra Lees-Edwards boundary conditions for lattice Boltzmann, J. Stat. Phys., 107 521-537 (2002).