Rheological response and dynamics of the amphiphilic diamond phase from kinetic lattice–Boltzmann simulations

R.S. Saksena, P.V. Coveney


The purpose of the present paper is to report on the first computational study of the dynamical and rheological response of a self-assembled diamond mesophase under Couette flow in a ternary mixture composed of oil, water and an amphiphilic species. The amphiphilic diamond mesophase arises in a wide range of chemical and biological systems, and a knowledge of its rheological response has important implications in materials science and biotechnological applications. The simulations reported here are performed using a kinetic lattice–Boltzmann method. Lyotropic liquid crystals exhibit characteristic rheological responses in experiments that include shear-banding and a non-Newtonian flow curve as well as viscoelasticity under oscillatory shear. Their behaviour under steady and oscillatory shear is correctly reproduced in our simulations. On cessation of shear, as the morphology returns to the diamond phase, the relaxation of the stress response follows a stretched-exponential form for low initial strain rates.

1. Introduction

Complex fluids composed of oil, water and amphiphilic species can undergo self-assembly to form triply periodic, bicontinuous structures at specific concentrations of the components and/or under particular temperature and pressure conditions (González-Segredo & Coveney 2004; Conn et al. 2006). The most commonly studied cubic phases are the bicontinuous gyroid (G, spacegroup Ia3d), diamond (D, spacegroup Pn3m) and primitive (P, spacegroup Im3m) phases. Here, we will consider ternary amphiphilic mixtures symmetric in the oil and water composition. The periodic bicontinuous structures consist of a three-dimensional network of interpenetrating yet non-intersecting channels filled with oil and water separated by an amphiphilic monolayer (Chin & Coveney 2006) in which the hydrophobic and hydrophilic ends of the amphiphile particles are oriented towards oil and water, respectively. The stability of these structures arises due to low bending energy as a result of the characteristic saddle-shaped surfaces of zero spontaneous curvature and the reduction of surface tension owing to the presence of amphiphiles at the interfaces (Lynch & Spicer 2005). The first bicontinuous structures were observed experimentally in soap–water systems and later in other systems including biological lipids. These phases are ubiquitous in nature; for example, they are observed in biomineralization processes and in lung tissue, being formed by both low molecular weight surfactants and amphiphilic polymers (Schwarz & Gompper 2002; Lynch & Spicer 2005). More recently, triply periodic minimal surface-containing cubic phases have been finding applications in consumer products, food technology (Mezzenga et al. 2005b), in materials science for the synthesis of mesoporous nanomaterials with interesting structural and electronic properties (Schwarz & Gompper 2002), in biosensing and ultrafiltration, protein crystallization (Grabe et al. 2003) as well as drug delivery applications in biotechnology (Shearman et al. 2006). Furthermore, elucidation of the mechanism and kinetics of self-assembly and phase transitions in cubic phases should provide insight into cellular processes, an active area of experimental investigation recently (Squires et al. 2005; Angius et al. 2006; Conn et al. 2006; Seddon et al. 2006).

The non-equilibrium behaviour of the amphiphilic cubic phases and their response to external fields such as an applied shear flow are commonly observed in experiments but not well understood theoretically (Chen et al. 1997). The rheological response of these fluids arises as a result of coupling between the microscopic dynamical processes governed by local interactions and the macroscopic hydrodynamic flow field. Applied shear may induce macroscopic transitions to morphologically and/or orientationally distinct self-assembled states; it can induce order–disorder transitions as well as order–order transitions from one phase with long-range order to a distinct phase also with long-range order (Chen et al. 1997; Mezzenga et al. 2005a,b; González-Segredo et al. 2006; You et al. 2007). In addition to oscillatory and steady shear, other flow fields, such as roll casting and extensional flow, are used in the processing of various self-assembled systems (Chen et al. 1997). Applied fields can be used to process materials, for example leading to defect-free assembly of mono-domains and mono-crystalline samples (Williams & Mackintosh 1994) and in the selection of particular orientations of the amphiphilic liquid crystal relative to the applied shear flow. A knowledge of the variational response of liquid crystals to a range of processing conditions opens up possibilities for the rational design and fabrication of functional materials for switchable optical, transport and electronic properties and of ordered nanostructures such as nanopores, nanosheets and nanofibres with desired directional properties (Fahmi et al. 2006). However, in order to exploit the characteristic response of amphiphilic liquid crystals to various types of external fields, one has to understand the physical basis of the effects of processing conditions such as the shear frequency, strain, temperature and pressure. A universal description applicable to a wide range of amphiphilic liquid crystals may be independent of the microscopic characteristics of the amphiphilic system and instead depend on the mesoscopic response of the system to the applied shear field.

The self-assembly of ternary amphiphilic mixtures into the gyroid, diamond, primitive and sponge phases has been studied previously using our lattice–Boltzmann (LB) approach (González-Segredo & Coveney 2004; Chin & Coveney 2006; Giupponi et al. 2006; González-Segredo et al. 2006; Saksena & Coveney 2008). Rheological simulations of the gyroid phase under Couette flow for 643 and 1283 lattice-site models exhibit non-Newtonian viscoelastic behaviour (Giupponi & Coveney 2006; Giupponi et al. 2006; González-Segredo et al. 2006), as is widely observed in rheological experiments of self-assembled amphiphilic systems (Bonn et al. 1998; Pitzalis et al. 2000; Porcar et al. 2003; Mezzenga et al. 2005a; Pouzot et al. 2007); in these LB simulations, the gyroid mesophase is found to possess enhanced viscosity as compared with the sponge and lamellar mesophases. The cubic diamond phase is observed to form quite commonly in amphiphilic fluid mixtures (Squires et al. 2005; Conn et al. 2006; Seddon et al. 2006; Shearman et al. 2006) and some experimental studies of its rheological properties have also been reported in the literature (Pitzalis et al. 2000; Mezzenga et al. 2005a). However, the self-assembly and rheology of the diamond phase have not received as much attention in simulation studies (Marrink & Tieleman 2001, 2002; Saksena & Coveney 2008) as the gyroid phase. We note here that, although treatment of spontaneous self-assembly of amphiphilic species into mesoscale structures is currently beyond all-atom molecular dynamics simulations because of the large length and time scales involved, considerable progress has been made in the application of coarse-grained molecular dynamics (CGMD) simulations to such problems (Michel & Cleaver 2007; Klein & Shinoda 2008). Ellison et al. (2006) have simulated the spontaneous self-assembly of the gyroid phase for hard- and soft-repulsive pear-shaped particles using coarse-grained Monte Carlo and molecular dynamics methods, respectively. CGMD has also been used to treat interactions of self-assembled lipid membranes with embedded proteins (e.g. Ayton et al. 2007; Treptow et al. 2008).

In this work, we focus on the rheology of the self-assembled diamond mesophase simulated using the kinetic LB method under applied Couette flow. Briefly, the kinetic LB method takes into account the hydrodynamic flow under shear as well as quasi-molecular interactions that describe the coupling between the applied shear field and microscopic dynamical processes within the system. Thus, the method provides access to chemical detail as well as longer (continuum) time scales compared with microscopic simulation methods such as molecular dynamics. The kinetic LB method differs from other LB methods based on the relaxation of a suitably chosen free energy functional (Theissen et al. 1998; Lamura et al. 1999; Denniston et al. 2001, 2004; Cates et al. 2005; Yeomans 2006). In the kinetic LB method, the system dynamics and its long-time behaviour (including the approach to equilibrium) are emergent from the quasi-molecular interactions specified in the particular model: we do not a priori impose a form for the free energy functional that would govern the self-assembly dynamics or the rheological response.

The paper is organized as follows. In §2 we describe our LB method for simulation of ternary amphiphilic fluids, as well as the visualization and analysis tools and computational infrastructure used in this work. In §3, we describe results from our simulations of diamond phase self-assembly, while, in §4, we present rheological simulations under both steady and oscillatory shear. In §5, we present our conclusions from this work.

2. Modelling and simulation

In this section we describe the ternary amphiphilic LB fluid model and simulation methodology. We also discuss visualization and analysis tools and the computing infrastructure used in this work.

(a) LB simulation method

We use a kinetic, mesoscale LB method that, in the limit of low Mach number, corresponds to the Navier–Stokes equations for nearly incompressible single-phase fluids. The amphiphilic species is modelled as a point dipole with one end having affinity for the oil species and the other end having an affinity for the water species. The strength of the amphiphile–amphiphile coupling, the repulsion strength between the oil and water species and the coupling of amphiphile to the oil and water species are governed by three independent interaction parameters. The kinetic LB model has been implemented in a Fortran 90 code called LB3D, which has been parallelized with a message passing interface using a domain decomposition strategy. The code has been described in detail previously (Love et al. 2003; Giupponi et al. 2006) and here we briefly mention salient points of the model. The evolution of the LB single-particle velocity distribution functions, Embedded Image and Embedded Image belonging to component σ (oil or water) and amphiphile (s), respectively, and with discrete velocity ci, is given by the following set of equations:Embedded Image(2.1)Embedded Image(2.2)where τα is the relaxation time for the species α=(σ, s). The terms Embedded Image, Embedded Image, Embedded Image and Embedded Image are matrix elements of the collision operators that result from mean-field interactions among different fluid components (Nekovee et al. 2000). The local equilibrium distribution functions, Embedded Image, are calculated by shifting the average velocity by an increment proportional to the net interaction force on the corresponding fluid species,Embedded Image(2.3)The time evolution of amphiphilic dipoles d(x, t) is given byEmbedded Image(2.4)Embedded Image(2.5)where d is the average dipole vector per site, whose relaxation is governed by the Bhatnagar–Gross–Krook (BGK) equation (2.5). Embedded Image is the average dipole vector at site x prior to collision and is calculated according to equation (2.4). This simple model of a ternary amphiphilic mixture gives rise to self-assembly into a range of periodic structures and complex transition behaviour under steady shear, which is consistent with experiments and simulations (Chin & Coveney 2006; Giupponi et al. 2006; González-Segredo et al. 2006; Saksena & Coveney 2008). These simulations allow us to directly observe the evolution of the initial mesophase as it undergoes morphological transitions under shear and to determine how the interactions between the component species affect the dynamics. In Couette flow simulations, shear is applied in the xz-plane using the Lees–Edwards boundary condition (Giupponi et al. 2006) with the shear velocity in the z-direction and the shear gradient in the x-direction. For a particle subject to Couette flow under Lees–Edwards boundary conditions, the new position and velocity along the flow direction are given, respectively, by the following:Embedded Image(2.6)Embedded Image(2.7)where ΔzUΔt; U is the shear velocity; uz is the z-component of the average fluid velocity u; and Nx is the system length in the x-direction. Additionally, an interpolation scheme suggested by Wagner & Pagonabarraga (2002) is used as Δz is not generally a multiple of the lattice site. For the cubic L3 lattice site systems we use in these simulations, the relationship between shear rate and applied shear velocity is given by Embedded Image. The minimum shear velocity we have used in our simulations is U=0.05. The maximum applied shear rate is determined by the requirement that the shear velocity should be much smaller than the speed of sound Embedded Image in the D3Q25 lattice used here (González-Segredo et al. 2006) for the stability of the LB algorithm.

The components of the pressure tensor are computed as the sum of a kinetic term and a virial mean-field term (González-Segredo et al. 2006)Embedded Image(2.8)The summation in the second term on the right-hand side runs over all combinations of the three fluid components, oil (r), water (b) and amphiphile (s). Since the interaction matrix Embedded Image is symmetric with only the nearest neighbour interactions considered, the virial term reduces toEmbedded Image(2.9)Under applied Couette flow, the negative of the off-diagonal component of the pressure tensor in the shear plane, −Pxz above, is equal to the stress response, σxz, of the sheared mesophase to the applied strain.

(b) Mapping of lattice units to physical units

We determine the mapping of simulation units to physical units by comparing the simulated diamond phase self-assembly with the lyotropic diamond phase self-assembly observed in experiments (Squires et al. 2005; Rappolt et al. 2006). Based on the experimental lattice parameter size of the order of 10 nm, the simulation resolution is approximately equal to 2 nm per lattice unit; this mapping is also in agreement with the previous simulations of the gyroid-phase assembly (González-Segredo & Coveney 2004). In experiments on diamond-phase self-assembly, for example, as induced by changing the concentration of water or on application of a thermodynamic jump, the self-assembly times are found to be of the order of minutes; comparing this with simulations, we estimate that each simulation time step corresponds to approximately 4–40 ms of physical time. Similarly, the mass in lattice units corresponds to approximately six surfactant molecules (of the order of 10−24 kg), based on observations of particle distribution functions around a single lattice site, which approximate the physical liquid crystal lattice site. According to the above mapping between lattice units and physical units, the applied strain rates in the Couette flow simulations are in the range of 0.01–0.1 s−1 in physical units with the frequency of applied oscillatory shear between 0.001 and 0.2 rad s−1. All simulation results in §§3 and 4 are reported in lattice units.

The mapping of lattice to physical units using the observed behaviour of the system, as described here, is likely to lead to optimistic estimates of the length and time scale covered by the simulation. One might attempt a more precise mapping by comparing dimensionless groups such as the Reynolds and capillary numbers. However, for complex fluids such as those in the present work, full resolution of all the length, time and energy scales is not generally possible with the LB method and, at best, one expects to be able to correctly treat the hierarchy of scales in order to reproduce the physics of interest (Cates et al. 2005).

(c) Isosurface visualization and topological characteristic

The interface between the oil and water phases is defined by the implicit surface of zero colour order parameter, ϕ(r)=ρr(r)−ρb(r)=0, for rEmbedded Image3 (here, we use superscripts σ=b for water species and σ=r for the oil species). As described by Chin & Coveney (2006), the continuum-order parameter field on the discrete simulation lattice is estimated by trilinear interpolation. In order to graphically visualize the interface, a polygon mesh approximation is generated using the marching cubes algorithm implemented using our own isosurfacing software. The marching cubes algorithm detects changes in the sign of ϕ(r) corresponding to a transition from an oil-rich site to a water-rich one, indicating the point where the link between the sites intersects the isosurface. The marching cubes algorithm proceeds by considering the sign of ϕ(r) at the vertices of each cube on the lattice to identify a set of polygons approximating the oil–water interface within the cube. The resulting polygon mesh representation of the isosurface of zero order parameter provides direct visualization of the interface between the oil and water phases. Given the polygonal mesh representation of a surface, we can calculate its topological Euler characteristic, Χsurf=VEF, where the surface has V vertices, E edges and F polygon faces. The Euler characteristic depends only on the topology of the surface and not on its shape or the particular polygon mesh representation (Kinsey 1993). Hence, the topological Euler characteristic of the surface or, equivalently, the genus g, related to Χ by Χ=2(1−g), is a useful metric to identify and compare cubic surfaces and has been used in this work. The Euler characteristic per unit cell is –8 for the gyroid surface, –16 for the diamond surface and –4 for the primitive surface (Schwarz & Gompper 1999). Our isosurfacing software runs on a single processor and has been used to analyse lattice sizes of up to 2563 lattice sites.

The morphology of the self-assembled cubic phases is verified with a numerically exact, ideal structure generated using the periodic nodal approximation (PNA; Schwarz & Gompper 1999; Squires et al. 2005). The PNAs for the three simulated cubic surfaces are

  1. gyroid surface is Embedded Image,

  2. diamond surface is Embedded Image, and

  3. primitive surface is Embedded Image.

The topological characteristics of the simulated self-assembled mesophase can be compared with the numerically generated surface, using the corresponding PNA on a lattice of the same size, in order to compare and contrast the simulated morphologies.

(d) Grid infrastructure

The simulations reported in this work were performed on high-performance transatlantic computing grid infrastructure. The computations were done on the petascale machine Ranger, at the Texas Advanced Computing Centre (http://www.tacc.utexas.edu/resources/hpcsystems/), with a theoretical peak of 0.58 petaflops and the Cray XT4 machine, HECToR, in the UK (http://www.hector.ac.uk), with a theoretical peak of 59 teraflops. The parallel kinetic LB code, LB3D, used in this work demonstrates excellent scaling on both resources, notably performing at a parallel efficiency of close to unity on up to 32 768 cores on Ranger (Saksena et al. in press). The simulation code was hosted in the application-hosting environment that enables distribution of simulations across grid resources (Coveney et al. 2007). Data transfer from the computational resources to local machines at UCL for post-processing and visualization was performed over the JANET-switched lightpath network (Coveney et al. in press; http://www.ja.net/services/lightpath/index.html), which is connected to the Starlight-switched router facility in the USA (http://www.startap.net/starlight/ABOUT/). The disk performance of computational resources is important in determining the overall performance of the LB3D code in our simulations; this is due to the large volume of data produced in this work (approximately equal to 3 TB). We found good scaling of I/O performance on up to 256 cores, the maximum used in this work. Weak-scaling benchmarks on Ranger, however, indicate that parallel scaling of our typical I/O requirements degrades when tens of thousands of cores are used; for example, performance decreases by 10 per cent, instead of remaining constant, on going from 1024 cores to 32 768 cores while keeping a constant load of 643 lattice sites per core.

3. Simulation details

Here we describe LB simulations of the diamond phase in a ternary amphiphilic mixture on two different lattice sizes—643 and 1283—that are chosen to determine the scaling of the rheological response with system size as well as any finite-size effects. The overall computational cost of simulating the entire rheological flow curve increases with lattice size: the computational cost of each simulation scales as L3 where L is the lattice dimension; additionally, the number of simulation time steps required to reach steady state under shear also increases with lattice size (∼L2).

In our LB simulations, the particle distribution functions in the initial random quench, fα(0) (where α can be oil (r), water (b) or amphiphile (s)), are distributed uniformly as 0≤fα(0)≤Fα, with the maximum density, Fα≤1.0, specified as input to the LB3D simulation code. The initial configuration is then allowed to evolve under the collision–advection scheme of the LB method in the presence of interaction forces (see equations (2.1)–(2.5)) for asymptotically long times. In simulations performed at specific values of the amphiphilic fluid composition and interaction parameters of the LB model, the ternary amphiphilic mixture undergoes self-assembly into the corresponding periodic mesophase, irrespective of the lattice sizes used in the simulation. Briefly, the interaction parameter gbr governs the interaction between the non-amphiphilic water (blue, σ=b) and oil (red, σ=r) species. A positive gbr is used to model repulsion between the two species (Chen et al. 2000; Love et al. 2003). A negative value of the water–amphiphile coupling parameter, gbs<0, favours assembly of periodic structures in which the headgroups of the amphiphilic particles are oriented towards water-rich regions and the tail groups are oriented towards oil-rich regions. The amphiphilic coupling parameter gss controls the interaction between amphiphile particles. For gss<0, it describes attraction between two tail groups and repulsion between a head and a tail group, while, for gss>0, the coupling parameter models repulsion between two headgroups, as might arise in ionic surfactants (Bhattacharya & Mahanti 2001). For gss>0, we find that, despite the repulsion between headgroups, the amphiphile particles still align themselves at the oil–water domain interfaces with headgroups pointing towards the water domains; this arises due to their directional affinity to oil and water species and the repulsive interaction (gbr>0) between the oil and water species.

Previously, we have performed systematic parameter sweeps to extend the mesophase diagram explored by the kinetic LB simulation method (Saksena & Coveney 2008). We reported a wide range of mesomorphic behaviours obtained by systematically varying the interaction parameters, gbr and gss, at a fixed species concentration. In these kinetic LB simulations, we were able to observe mesomorphic variations in agreement with experimental studies of lyotropic solutions as (i) the amphiphilic strength gss was varied and (ii) the oil–water repulsive interaction gbr was varied.

(a) Cubic-phase self-assembly

For the same concentration of oil, water and amphiphilic species and different values of gbr, gbs and gss, we observe the self-assembly of the primitive, gyroid and diamond mesophases (Saksena & Coveney 2008). The parameters corresponding to the three phases are listed in table 1. The morphology of the self-assembled cubic domains is verified with an ideal structure generated on a lattice of the same size and with the same number of unit cells using the PNA (Schwarz & Gompper 1999; Squires et al. 2005) by means of (i) visual comparison and (ii) comparison of the topological Euler characteristic (Gandy et al. 1999; Chin & Coveney 2006) with the PNA domain (Squires et al. 2005; Saksena & Coveney 2008). In figure 1, we visualize the three-dimensional isosurfaces of the simulated self-assembled gyroid, diamond and primitive domains and the corresponding PNAs. We have selected sub-domains of the self-assembled cubic mesophases which are closest to the ideal PNA sub-domain. During the simulated self-assembly process, there are small changes in local alignment of the cubic phase with respect to the simulation axes and spontaneous formation of defects, which result in differences between the simulated system morphology and the corresponding ideal PNA structures that are idealized static representations, perfectly aligned with the x, y and z simulation axes.

View this table:
Table 1

Self-assembly of mesophases at oil, water and amphiphile concentrations of Fr=0.7, Fb=0.7, Fs=0.9, respectively. (Variation of the interaction parameters between the three species leads to the self-assembly of the primitive, diamond and gyroid triply periodic cubic phases as well as the disordered sponge phase.)

Figure 1

Comparison of the simulated (a) gyroid, (b) diamond and (c) primitive phases with the numerically generated (d) gyroid PNA, (e) diamond PNA and (f) primitive PNA. Simulation parameters are as per table 1. Isosurface renderings of 163 sub-domains extracted from the simulations of 643 lattice-site systems are shown here. The isosurface represents the interface between the oil and water domains and is generated from points where oil and water densities are equal; the isosurface is coloured red on the oil-majority side and blue on the water-majority side.

In the rest of this paper, we will focus on the diamond phase that undergoes self-assembly in our LB simulations at parameters specified in table 1 (system II). The salient morphologies that are encountered as the system self-assembles into the diamond phase on a 1283 lattice-site system, from a random quench, are shown in figure 2. Systems of both sizes took between 13 000 and 16 000 simulation time steps to self-assemble into the diamond phase from the initial random quench. For the 643 lattice, we observe a self-assembled ideal diamond phase with an Euler characteristic value equal to that of the corresponding PNA surface, while, on the 1283 lattice, the self-assembled diamond phase contains defects and its Euler characteristic differs from that of the corresponding PNA surface.

Figure 2

Salient morphologies that are encountered as the diamond phase (system II of table 1) is self-assembled on a 1283 lattice from a random quench as a function of lattice time step (a) t=0, (b) t=4000, (c) t=8000, (d) t=12 000, (e) t=14 000 and (f) t=16 000. Colour-coding of the isosurface is as per figure 1. Between t=19 000 and 63 000, the system shows further complexity with the formation of intermediate primitive sub-domains that reconvert to the diamond phase (figure 3).

After the self-assembly of the diamond phase, we follow the diamond phase dynamics on the 643 and 1283 lattices over the next 50 000 time steps and find that the system undergoes transitions to a short-lived, intermediate state with primitive domains; this reconverts to a purely diamond morphology within 1500 time steps. Finite-size effects play a major role in the defect dynamics and the subsequent phase transition to primitive domains. In particular, in the simulations on smaller 643 lattices, the primitive sub-domains grow to occupy the entire simulation cell; that is, we see a spontaneous ideal diamond→ideal primitive transition prior to reconversion to the diamond phase. By contrast, in simulations on the larger 1283 lattice, transient primitive-phase sub-domains nucleate around defects and coexist with the diamond and sponge phases before annealing to a lattice with only diamond domains. In figure 3, we plot the dynamical evolution of the Euler characteristic Χ for the 643 and 1283 systems; the ideal line shows the Euler characteristic value of a perfect diamond mesophase with 103 unit cells on a 643 lattice-site system and 203 unit cells on a 1283 lattice-site system. The troughs in the plot of the Euler characteristic correspond to the short-lived intermediate state comprising primitive sub-domains that form and disappear as the simulation progresses. These transitions are hydrodynamically driven (as is confirmed by direct correlations between the fluid velocity field and these non-equilibrium phase transitions; data not shown). Asymptotically, the system is eventually expected to reach an equilibrium state, where such transitions are absent (although it may be noted that interacting LB models neither satisfy detailed balance nor, in turn, support an H-function; Boghosian et al. 2001). However, experimentally, it is known that D and P cubic phases often coexist and thus one may expect the barrier to their interconversion to be small (Czeslik et al. 1995; Winter et al. 1998; Conn et al. 2006).

Figure 3

Plot of the dynamical evolution of the topological Euler characteristic for (a) 643 lattice-site system and (b) 1283 lattice-site system starting from a random quench. After forming the diamond phase, the system undergoes short-lived transitions to primitive domains. The dashed ‘ideal’ line shows the Euler characteristic of the perfect diamond mesophase, while the troughs correspond to short-lived morphologies with primitive sub-domains (solid curve, simulation).

In figure 4, we show isosurface plots of the ternary system in and around the diamond→(primitive+diamond+sponge)→diamond transition for the 1283 lattice-site system. This sequence of isosurfaces illustrates the morphological changes that occur as the diamond phase converts to the intermediate primitive phase. This particular morphological sequence corresponds to the sixth peak in figure 3b. We note that although the genus and Euler characteristic are important metrics for classification and verification of topologies, they do not furnish a sufficient condition to guarantee the existence of a particular cubic phase. However, for systems that are already known to be in a particular cubic morphology, for example the diamond phase in the present case, the topological parameters provide a succinct metric by means of which to follow the morphological dynamics of the periodic mesophase, highlighting points in the morphological landscape which are worthy of more detailed investigation.

Figure 4

Isosurface evolution corresponding to the vicinity of the sixth peak in the genus versus time plot in figure 3b, for the 1283 lattice-site system. Starting from a diamond phase at (a) t=43 300, the following morphologies are observed as the simulation proceeds; (b) t=43 700, primitive+diamond sub-domains; (c) t=44 000, primitive+diamond sub-domains; (d) t=44 600, sponge+primitive+diamond sub-domains; (e) t=44 900, diamond+sponge sub-domains; (f) t=45 200, diamond phase reforms. Colour-coding of the isosurface is as per figure 1.

The relatively high temporal and spatial resolution of the dynamics of the D→P transition that we obtain from our simulations provide an opportunity to analyse the phase transition mechanism. It is notable that these kinetic LB simulations allow us, for the first time, to study cubic→cubic transitions at a high resolution. Preliminary investigation into the origin of the formation of primitive sub-domains indicates that growth of primitive sub-domains is nucleated in the vicinity of defects in the diamond phase.

We turn next to consider the non-equilibrium properties of the diamond phase under an applied Couette flow.

4. Rheology

In this section we report the rheological response of the diamond phase to applied Couette flow. In the LB simulations of Couette flow, a linear shear velocity profile is applied using Lees–Edwards boundary conditions as described in §2. Our simulations are able to correctly simulate the shear-banding effect and non-Newtonian responses of liquid crystals under steady shear as observed in experiments (Bonn et al. 1998; Pitzalis et al. 2000; Rey & Denn 2002; Smith et al. 2002; Porcar et al. 2003; Mezzenga et al. 2005a). The simulated viscoelastic response under an oscillatory shear is also in very good agreement with analogous experiments on the diamond phase (Pitzalis et al. 2000; Mezzenga et al. 2005a).

(a) Approach to steady-state Couette flow and shear banding

The number of simulation time steps required for steady-state Couette flow to be established depends on the lattice size, with a shear velocity profile occurring in approximately 1000 time steps for the 643 lattice and within 4000 time steps for the 1283 lattice. The evolution of the velocity profile from the unsheared system to steady state for the two system sizes is shown in figure 5. The nonlinear velocity profile that is established here at steady state is evidence of a backflow effect (Jewell & Sambles 2002; Smith et al. 2002; Marenduzzo et al. 2003): the movement of the fluid components themselves under the influence of the applied shear velocity as well as local surface reordering causes a modification of the linear velocity profile under Couette flow. This is a mesoscopic phenomenon involving an interplay between the externally applied macroscopic velocity flow field and microscopic dynamical processes governed by the interspecies interactions. The nonlinear velocity profile is accompanied by the appearance of distinct liquid crystal morphology (texture) in regions along the shear velocity gradient within the simulation cell. Such shear banding has been observed in our rheological simulations of the diamond phase for both system sizes. This non-Newtonian effect has also been observed in experiments and other simulations of liquid crystals under shear (Rey & Denn 2002; Marenduzzo et al. 2003).

Figure 5

Evolution of the applied steady-shear velocity profile, U=0.11, starting from the unsheared diamond system (system II in table 1) to steady-state Couette flow for (a) 643 lattice (black curve, t=0; red curve, t=200; green curve, t=500; blue curve, t=700; yellow curve, t=1000; brown curve, t=2000; grey curve, t=3000) and (b) 1283 lattice (red curve, t=0; green curve, t=1000; blue curve, t=2000; yellow curve, t=3000; brown curve, t=4000; grey curve, t=5000; violet curve, t=6000; turquoise curve, t=7000; pink curve, t=8000).

In figure 6, we show the shear-banded texture in the xz (shear)-plane and the xy (transverse)-plane for the 643 and 1283 lattice-site system, at a shear velocity U=0.11 in lattice units. The isosurface visualization of the interfacial structure in the 643 and 1283 lattices is shown here.

Figure 6

Isosurface visualization of the interface between oil (red) and water (blue) using contour plots corresponding to ϕ(r)=ρrρb=0. The view here is of a steady-state shear-banded texture along the xz shear plane for (a) 643 and (b) 1283 lattices and the xy transverse plane for (c) 643 and (d) 1283 lattices at a shear velocity U=0.11. The shear velocity is along the z-direction and the shear velocity gradient is along the x-direction.

In contrast to the morphology at the centre of the simulation cell, the shear-banded texture in the high shear velocity regions comprises a smectic lamellar phase, oriented in the transverse direction relative to the applied shear velocity (figure 6). Such an orientation relative to the applied velocity field is particularly striking as one would expect it to be unstable and rather a parallel or perpendicular lamellar phase to exist at a steady state. However, intermediate smectic domains aligned transverse to the shear flow have indeed been observed in experimental studies of di-block copolymers, during the transition from an isotropic mixture to the parallel/perpendicular lamellar phase under shear (Chen et al. 1997). Similarly, in our simulations, transverse lamellar domains are formed when the diamond channels are subjected to shear flow in the z-direction as the oil and water channels get reoriented under the flow field. We conjecture that this transverse orientation is reinforced due to backflow effects and hence it persists, an effect also observed in twisted nematic devices where a short-lived rotation of the director field occurs in the ‘wrong’ direction upon switching (Smith et al. 2002; Marenduzzo et al. 2003). There are also examples, in numerical studies, of the development of transverse orientational patterns in liquid-crystalline polymer systems under steady shear where the interactions are modelled by the Marrucci–Greco potential or by the simpler Leslie–Ericksen theory (Rey & Denn 2002). Backflow effects will, of course, be significant only when the applied rate of external shear is low enough to be able to couple to intrinsic fluid dynamical processes such as surface reordering; we would expect that, at high shear rates, shear flow will become the dominant factor and break the ordering of lamellae transverse to it. This is indeed what is observed in our simulations beyond a critical value of the strain rate: at such strain rates greater than the critical value, the transverse smectic phase is ruptured at steady state as we will describe in §4b. Although many experimental reports exist in the literature on the non-Newtonian rheology of amphiphilic cubic liquid crystals (Pitzalis et al. 2000; Porcar et al. 2003; Mezzenga et al. 2005a; Pouzot et al. 2007), there are very few papers reporting morphological transitions induced by shear (for example, Porcar et al. (2003) discuss the shear-induced transition from the sponge to the lamellar phase in an amphiphilic solution), although shear-induced morphological transitions in self-assembled block copolymer systems have been studied in more detail (Fredrickson & Bates 1996; Park & Char 2005; Darling 2007). Experimental verification of the formation of such a transverse smectic phase at low strain rates when a self-assembled diamond morphology is subjected to steady shear is thus not currently available; we hope our work will inspire further laboratory studies.

(b) Stress response

Once steady shear had been established, we averaged the stress response over 5000 and 3000 time steps for the 643 and 1283 lattice-site systems, respectively. The average stress-response versus strain-rate profiles and viscosity versus strain rate for the two system sizes are plotted in figures 7 and 8.

Figure 7

Yield curve and viscosity profile characterizing the rheological response of the 643 lattice-site system starting from the diamond phase: (a) stress response versus measured strain rate and (b) viscosity versus measured strain rate. The starting configurations, to which steady shear is applied, are taken from the dynamical evolution of the cubic diamond phase. Error bars indicate the standard deviation of the calculated average stress response and viscosity values at each strain rate shown in the figure (pink curve, t0=30 000; red curve, t0=31 000; green curve, t0=32 000; blue curve, t0=40 000; yellow curve, t0=44 000; brown curve, t0=62 000; black curve, t0=81 000).

Figure 8

Yield curve and viscosity profile characterizing the rheological response of the 1283 lattice-site system starting from the diamond phase: (a) stress response versus measured strain rate and (b) viscosity versus measured strain rate. The starting configurations, to which steady shear is applied, are taken from the dynamical evolution of the cubic diamond phase. Error bars indicate the standard deviation of the calculated average stress response and viscosity values at each strain rate shown in the figure (black curve, t0=21 000; red curve, t0=32 000; green curve, t0=37 000; blue curve, t0=43 000; yellow curve, t0=56 000; brown curve, t0=64 000).

For both lattice sizes, we performed multiple runs to generate the flow curve starting from different initial diamond phase configurations, which have been aged for different lengths of time. At a shear velocity close to U≈0.13, we observe dramatic shear thinning behaviour, i.e. the viscosity of the system reduces sharply as the strain rate is increased. However, at this strain rate, the value of the stress response itself exhibits variations between systems sheared from different initial configurations in the diamond morphology, indicative of hysteresis.

Shear thinning associated with a stress plateau in the flow curve has been widely observed experimentally in a range of self-assembled systems and morphologies such as lyotropic lamellar, multi-lamellar and hexagonal phases, thermotropic liquid crystals and block copolymers. We are unaware of analogous experiments on bicontinuous cubic liquid crystals. Here, we draw comparison with phenomenologically related experiments on other lyotropic self-assembled systems. In experiments involving Couette flow, flow curves have been reported, which are in qualitative agreement with our simulations over the entire range of strain rates simulated. As examples, we mention here experiments comprising fixed strain-rate rheological measurements of (a) the micellar bcc cubic crystal in aqueous triblock copolymers (Eiser et al. 2000) and the lyotropic liquid-crystalline lamellar ‘onion’ phase (Bonn et al. 1998), in which flow curves with a stress plateau together with a sharp drop in viscosity with increasing strain rate have been observed. The existence of the stress plateaux in these experiments has been attributed to the coexistence of high- and low-viscosity phases, which are stable on the lower and higher shear rate branches, respectively. The occurrence of hysteresis in the experiments in the vicinity of the stress-response drop has been attributed to the fact that the system was bistable, with a small energy barrier between the two states, thus resulting in two values of the viscosity for a given strain rate (Bonn et al. 1998). The hysteresis in our simulations would also arise if the system were indeed bistable in the vicinity of this critical value of the strain rate corresponding to a morphological transition in the sheared system. In our diamond phase rheology simulations, the bistable states are (i) the periodic (smectic+diamond) morphology found at lower strain rates and (ii) the disordered morphology, with oil and water channels running parallel to the shear flow, which is found at higher strain rates as described below.

At the transition strain rate (figure 7), the diamond phase undergoes a morphological transition, in which the oil and water channels aligned transverse to the shear flow are disrupted and structure at steady state under shear is no longer composed of the periodic (transverse smectic+diamond) phase. Figure 9 shows the steady-state morphology at strain rates beyond the morphological transition; this can be compared with the steady-state morphology at a lower strain rate corresponding to U=0.11 in figure 6. For the 643 lattice systems and to a lesser extent for the 1283 lattice systems, the disruption in periodicity extends into the centre of the simulation cell.

Figure 9

Isosurface visualization of the interface between oil (red) and water (blue) using contour plots corresponding to ϕ(r)=ρrρb=0. The view here is of a steady-state shear-banded texture along the xz shear plane for (a) 643 and (b) 1283 lattices and the xy transverse plane for (c) 643 and (d) 1283 lattices at a shear velocity U=0.26. The shear velocity is along the z-direction and the shear velocity gradient is along the x-direction.

In the steady-shear experiments of the lamellar ‘onion’ phase (Bonn et al. 1998), the disappearance of the shear-banded morphology was detected at high strain rates by phase contrast microscopy and was attributed to either (i) the actual disappearance of the shear bands or (ii) inability to detect shear bands as they became too small. Our results suggest the former to be the more likely explanation for the disappearance of the shear bands at high shear rates in these experiments (Bonn et al. 1998). We do not observe any strain-rate-dependent variations in the size of the shear bands in simulations of diamond phase rheology.

To the best of our knowledge, there are no experimental results in the literature on the amphiphilic diamond mesophase under steady shear which exhibit the stress plateau and with which we could make a quantitative comparison. However, in the phenomenologically similar steady-strain rheological experiments of the biphasic onion phase (Bonn et al. 1998), analogous extreme shear-thinning behaviour was observed at Embedded Image, while in the steady-strain rheological experiments investigating the conversion of the L3 sponge to the lamellar phase, the transition occurred at Embedded Image with a decrease in viscosity (Porcar et al. 2003); it was also noted that the critical shear rate could be reduced by a factor of 8 by halving the lipid volume fraction.

To summarize, for the two system sizes simulated, we find agreement in the profile of the flow curve comprising a non-Newtonian shear-thinning regime that includes a striking drop in viscosity. For the 643 lattice system, we observe subtle finite-size effects on closely following the shear-induced morphological transition dynamics from the initial unsheared configuration. In contrast to the gradual morphological evolution under shear that is observed in the larger systems, for the 643 lattice-site system, the initial diamond morphology completely breaks up into a random quench as soon as shear is applied; the shear-induced assembly of the transverse smectic+diamond phase is observed subsequently. This indicates that the long-range correlations in the initial diamond liquid crystal present in the smaller 643 lattice-site system are coupled to the periodic boundary conditions: a slight perturbation in the boundary condition causes the liquid-crystalline morphology to break up.

Finally, we find that the steady-state flow curves for the two system sizes have the same shape but differ quantitatively. For the 1283 lattice-site system, the sharp drop in viscosity occurs at a strain rate that is close to half of that for the 643 lattice-site system. There is no theoretical basis for such a scaling of the critical strain rate with system size; however, it is notable that non-overlap of Couette flow curves measured for different Couette cell gap widths and non-trivial correlation of Couette cell gap width with wavelength of vorticity structures have been observed in experiments on amphiphilic fluids (Ganapathy et al. 2008).

(c) Relaxation of the stress response

We also analysed the time dependence of the relaxation of the stress response once the shear flow field is extinguished. In previous LB and lattice-gas simulations, domain growth kinetics during self-assembly of amphiphilic ternary mixtures followed a stretched-exponential functional form (Emerton et al. 1997; Nekovee et al. 2000). The viscoelastic response of the gyroid phase was also found to be consistent with stretched-exponential kinetics (Giupponi et al. 2006). Following earlier work, we fit the relaxation kinetics to the following form:Embedded Image(4.1)We used Matlab (http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml) to perform a least-squares fit of the relaxation curves of the 1283 diamond phase lattice sheared at various initial strain rates to the stretched exponential. The particular sheared system we study here corresponds to the system whose flow curve is plotted in figure 8 with legend t0=43 000. The relaxation curve and corresponding least-squares fit for two initial strain rates are shown in figure 10, while the inset shows the relaxation curves for a range of initial strain rates. The fitted parameters A, B, C and D for the strain rates corresponding to different shear velocities, U, are listed in table 2. The value of the stretched-exponential exponent D, and hence relaxation times, increases with initial applied strain rate and remains close to 1 for the highest strain rates considered here, corresponding to pure exponential relaxation dynamics. The time dependence of the stress relaxation after the cessation of shear flow gives a good indication of the time scales involved in the intrinsic relaxation processes of the self-assembled mesophase. The internal stress response of the system approaches zero and, accompanying this, re-adopts the diamond morphology within 5000 time steps. This is similar to experimental observations of a sheared gyroid phase where the cubic phase recrystallizes after the cessation of shear albeit with a different orientation (Olsson & Mortensen 1995).

Figure 10

Simulated and corresponding stretched-exponential fits for the stress response versus time of a 1283 diamond lattice starting from the time at which the applied shear flow field is terminated. Plots are labelled according to the shear velocities, U, at which the system was initially sheared (black circles, U=0.11; solid curve, fit U=0.11; green circles, U=0.18; dashed curve, fit U=0.18). The inset shows the relaxation profiles for a range of initial shear velocities, U=0.1–0.18 (black curve, U=0.1; red curve, U=0.11; green curve, U=0.12; blue curve, U=0.13; yellow curve, U=0.15; brown curve, U=0.18).

View this table:
Table 2

Parameters A, B, C and D in equation (4.1) are used to fit the simulated relaxation curves of the sheared diamond phase using the least-squares method. (The system studied here is the 1283 lattice-site diamond phase, the flow curve for which is plotted in figure 8 and corresponds to the plot with legend t0=43 000.)

(d) Oscillatory shear and viscoelasticity

Oscillatory shear is simulated by modulating the applied shear velocity according to the equationEmbedded Image(4.2)where ω/2π is the frequency of oscillation and shear velocity, U, is set to 0.1. Here we follow the approach taken by Giupponi et al. (2006) in which such an oscillating shear was applied to the gyroid mesophase to measure its viscoelastic response. The simulations were run for at least three complete oscillations over a range of shear frequencies, ω. We calculated the storage and loss moduli by fitting the stress-response profile to the equationEmbedded Image(4.3)where Embedded Image and G′(ω), G″(ω) are the storage and loss moduli, respectively. The fitted values for G′ and G″ versus ω are plotted in figure 11.

Figure 11

Plots of the dynamic storage G′ (solid curve with circles) and loss G″ (solid curve with crosses) moduli versus applied oscillatory strain frequency ω for the diamond phase modelled on a 1283 lattice-site system (system II in table 1).

At low frequencies, the loss modulus dominates, indicating a viscous response, while at high frequencies the system exhibits a solid-like elastic response. The frequency of oscillation at which the storage and loss moduli curves intersect signals the transition from liquid-like behaviour to solid-like behaviour; it gives an estimate of the characteristic stress relaxation times in the system and provides clear evidence of the viscoelastic behaviour of the simulated liquid-crystalline diamond phase. The form of the dynamic moduli versus strain oscillation frequency is in qualitative agreement with experimental studies of the shear rheology of the diamond phase formed in lyotropic liquid crystals (Pitzalis et al. 2000; Mezzenga et al. 2005a). In figure 11, the crossover between liquid-like and solid-like behaviours occurs at ω≈0.015, which corresponds to a relaxation time τ≈419 in lattice units and 17 s in physical units according to the mapping in §2b. The observed relaxation time is also in good agreement with that obtained for the diamond phase in the experiments of Mezzenga et al. (2005a). It should be noted that, in the experiments of Pitzalis et al. (2000), the diamond phase exhibits a dominant storage modulus and is essentially elastic except for the lowest shear frequencies. We plot the stress-response curve versus simulation time in figure 12 at ω=0.015 where the diamond phase exhibits a viscoelastic response with a phase lag ϕ between 0 (purely elastic) and π/2 (purely viscous).

Figure 12

Plot of stress response versus time for oscillatory shear (U=0.1 and ω=0.015) applied to the diamond phase (system II in table 1). The applied strain (dashed curve) is plotted (scaled) as a guide to the eye (solid curve, 〈σxz〉). Shear moduli are calculated as per equation (4.3). The phase lag is ϕ≈0.25π, intermediate between a 0 and π/2 lag, characteristic of a viscoelastic response.

5. Conclusions

In this paper, we have reported the first simulations of a ternary amphiphilic diamond mesophase characterizing the rheological response under an applied Couette flow. Under shear, the diamond phase exhibits many of the characteristic phenomena of sheared liquid-crystalline systems observed in experiments (Bonn et al. 1998; Pitzalis et al. 2000; Rey & Denn 2002; Smith et al. 2002; Mezzenga et al. 2005a), including a non-Newtonian flow curve and shear banding. Non-Newtonian behaviour arises in complex amphiphilic fluids because of the coupling between the length and time scales of molecular ordering, determined by microscopic interactions and the applied macroscopic flow field. It is notable that such rich phenomenology can be accessed using a minimal amphiphilic model simulated on a three-dimensional lattice with a discrete set of velocity vectors using the LB method. At a critical value of the strain rate, we observe extreme shear thinning corresponding to a steep drop in viscosity accompanied by a morphological transition in which the ordered (smectic+diamond) phase is broken up under shear and channels are formed oriented in the direction of shear. We observe hysteresis in the rheological response of the amphiphilic fluid at the transition strain rate and beyond, particularly in the 643 lattice simulations. Finally, we note that the flow curves for the 1283 and 643 systems exhibit the same shape but do not overlap. The rheological responses of the two system sizes differ in part due to the presence of multi-domain defects in the larger 1283 lattice-site systems.

Finite-size effects are found to be pronounced on following shear-induced morphological transitions in the smaller (643) systems where, under shear, the initial diamond morphology exhibits a striking break-up into a random quench before assembling into the steady-state shear-induced morphology. For larger lattices, the transition to steady-state morphology under shear occurs gradually. This highlights the importance of using a sufficiently large model in simulations of dynamical processes involving local interactions and in which long-range correlations arise due to hydrodynamics, inherent periodicity of the system or are induced by an externally applied field. However, despite the difference in the dynamical approach to steady state, we obtain the same shear-induced steady-state morphology for both smaller and larger system sizes.

On comparing various functional forms (Emerton et al. 1997; Nekovee et al. 2000) for the relaxation of the stress response on the cessation of shear, as the diamond phase reforms, the best fit is obtained by a stretched-exponential function for systems initially sheared at low strain rates. The stretched-exponential form has previously been found to describe domain growth kinetics in simulations of ternary amphiphilic fluids and in experimentally observed micellar self-assembly (Wilcoxon et al. 1995). Stretched-exponential kinetics are characteristic of glassy systems with multiple relaxation times. In our simulations, at high strain rates, as channels develop along the flow direction there is a transition from a high-viscosity glassy steady state to a low-viscosity steady state, for which the stretched exponential gives a poorer fit. The stress-response relaxation provides insight into the probable functional form of the kinetics of intrinsic dynamical processes such as the growth of defects and rearrangement of the domains in the unsheared system.


We are grateful to EPSRC for funding this research through the grant EP/E045111/1, which also provided access to HECToR. This work was also supported by US National Science Foundation MRAC grants TG-DMR070013N and TG-DMR070014N. We are grateful for DECI allocations on the DEISA grid (http://www.deisa.org) provided by the EU (FP6 project 508830 and FP7 project 222919). We would like to thank Davide Marenduzzo for useful discussions.


    • Received November 21, 2008.
    • Accepted March 4, 2009.


View Abstract