Stress response and structural transitions in sheared gyroidal and lamellar amphiphilic mesophases: lattice-Boltzmann simulations

0.54 MB·

by  Nelido Gonzalez-Segredo

Stress response and structural transitions in sheared gyroidal and lamellar amphiphilic
mesophases: lattice-Boltzmann simulations
N´elido Gonz´alez-Segredo∗
FOM Institute for Atomic and Molecular Physics (AMOLF),
P. O. Box 41883, 1009 DB Amsterdam, The Netherlands
6 Jens Harting†
0
Institut fu¨r Computerphysik, Pfaffenwaldring 27, D-70569 Stuttgart, Germany
0
2
Giovanni Giupponi and Peter V. Coveney‡
n
Centre for Computational Science, Department of Chemistry,
a
J University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom
0 (Dated: February 6, 2008)
2
Wereportonthestressresponseofgyroidalandlamellaramphiphilicmesophasestosteadyshear

] simulated using a bottom-up lattice-Boltzmann model for amphiphilic fluids and sliding periodic
t (Lees-Edwards) boundary conditions. Westudy the gyroid per se (above the sponge-gyroid transi-
f
o tion,of high crystallinity) andthemolten gyroid (within such atransition, of shorter-range order).
s
We find that both mesophases exhibit shear-thinning, more pronounced and at lower strain rates
.
t for the molten gyroid. At late times after the onset of shear, the skeleton of the crystalline gyroid
a
m becomes a structureof interconnected irregular tubesand toroidal rings, mostly oriented along the
velocityrampimposedbytheshear,incontradistinctionwithfree-energyLangevin-diffusionstudies
-
d which yield a much simpler structure of disentangled tubes. We also compare the shear stress and
n deformationoflamellarmesophaseswithandwithoutamphiphilewhensubjectedtothesameshear
o flow applied normal to the lamellae. We find that the presence of amphiphile allows (a) the shear
c
stressat late timestobehigherthan inthecase without amphiphile,and (b)theformation ofrich
[
patternson thesheared interface, characterised byalternating regions of high and low curvature.

1
v
6 I. INTRODUCTION ous and hydrocarbon species. These sheets or labyrinths
7
form a triply periodic minimal surface (TPMS) whose
4
unit cell is of cubic symmetry, has zero mean curvature,
1 The study of the response to shear in amphiphilic
0 no two points on it are connected by a straight segment,
mesophases has been the subject of attention for numer-
6 and has no reflexion symmetries. Their skeletons, i.e.
0 ical modellers only in recent years. The interest in the
the locus bounded by the TPMS, for each immiscible
/ subjectissustainednotonlybythewiderangeofapplica-
t phase, form double (inter-weaving), chirally symmetric
a tions in materials’ science and chemical engineering, but
m 3-fold coordinated lattices. There are lyotropic [14, 15]
also by the need to gain a fundamental understanding of
- the universal laws governing the self-assembly processes andthermotropictransitionsbetweenthe gyroidandthe
d microemulsion mesophase, the latter being a bicontinu-
n and competing mechanisms present.
ous mesophase of short-range order. The morphologies
o
Hitherto, studies have focused mainly on the struc-
c inthecrossoverregionsofthesetransitionsshowshorter-
: tural changes induced by steady and oscillatory shear, rangeorderthanthegyroid’sandlonger-rangeorderthan
v
near and far from critical points, in polymer sys-
i the microemulsion’s, reasons for which they are termed
X tems [1–6]. The morphologies studied include cubic-
‘molten gyroids.’
r andwormlike-micellar,lamellarandhexagonally-packed-
a Bicontinuous cubic mesophases of monoglycerides and
tubular mesophases;more complexstructures arethe so-
the lipid extract from archaebacterium Sulfolobus solfa-
called bicontinuous mesophases, of which those liquid-
taricushavebeenfoundatphysiologicalconditionsincell
crystalline of cubic symmetry have thus far been consid-
organelles and physiological transient processes such as
ered in far less detail.
membrane budding, cell permeation and the digestion of
Theamphiphilic gyroid[14,15]isabicontinuouscubic fats [7]. They can also be synthesised for important ap-
liquid crystalconsistingof multi- or mono-layersheets of plicationsinmembraneproteincrystallisation,controlled
self-assembledamphiphiledividingtworegions,eachcon- drug release and biosensors [8, 9].
tainingphaseswhicharemutuallyimmisicible,e.g. aque- The purpose of this paper is to report on the re-
sponse to shear of gyroid (G), molten-gyroid (MG) and
lamellar (L ) amphiphilic mesophases simulated using a
α
bottom-up kinetic-theoretic model for fluid flow. The
∗Email: [email protected]
†Email: [email protected] model is based on a lattice-Boltzmann (LB) method,
‡Email: [email protected] which has proved to be a modelling tool alternative to
2
and more efficient and robust than sofisticated methods whichallowsustosimulate,forlargeenoughlattices[19],
based on continuum equations. This LB method ad- the Navier-Stokes (NS) momentum-balance equation in
heres to a bottom-up complexity paradigm [13] in the thebulkofeachimmisciblefluidspecies,namely“oil”(or
sense that it is simple, fully particulate and no hypothe- “red”,‘r’)and“water”(or“blue”,‘b’). Themodelallows
ses of desirable macroscopic behaviour are imposed on thesimulationofcorrectphase-segregationkineticsinthe
the microdynamics—yet, we have shown in the past its absence [20] and presence [15] of a third, amphiphilic
abilitytosimulatecorrectsegregationkineticsforimmis- (surfactant-like,“s”)dipolarspecies. Themodelcontrols
cible fluids [20] and non-equilibrium self-assembly into the inter-particle forces between r, b, and s species via
amphiphilic mesophases [14, 15]. Knowing that such a couplingparameters(g ,g ,g ),andtransientsarecon-
br bs ss
simple model is capable of simulating these kinetic pro- trolled via relaxation times for densities (τb,τr,τs) with
cesses from a purely bottom-up dynamics, in this pa- an additional relaxation time for the orientation of the
per we investigate how hydrodynamic interactions cou- amphiphiledipoles(τd). Inaddition,themodelsimulates
ple with self-assembly and modify the stability and mor- the nonequilibrium self-assembly and relaxation dynam-
phologyof the mesophases. The noveltyof ourworkalso ics of sponge (L ) and gyroid mesophases [14, 15]. The
3
restsinthemodel’scapabilitytoreproducemorphological gyroidsthat it simulates show rigidity, arisingfrom their
transitionswithouthavingtoassumeamacroscopic,free- crystallineordering,whichdecreasesastheconcentration
energy model, used in other LB methods [3, 10] to com- of amphiphile is reduced; indeed, a lyotropic transition
pute the diffusive currents substantiating self-assembly. causes the correlation length to decrease towards that
In addition, since our method models amphiphilic of a sponge mesophase through a molten-gyroid state.
molecules as point dipoles—the simplest possible partic- Thisidea iscentraltothe workwepresenthere: weshall
ulate model for an amphiphile—the rheological features seethatthemesophase’scrystallineorderingenhancesits
emergentfromitareexpectedtobeuniversalforabroad stress response; indeed, we find shear-thinning to occur
range of amphiphilic systems. Finally, most of the nu- at higher strain rates for gyroids than for sponges.
mericalstudies measuringthe stressresponse of complex The Lees-Edwards boundary conditions (LEBC) were
fluids to shearreportedinthe literaturedealwithphase- originallyproposedbyLeesandEdwardsinthecontextof
segregating fluids on one side [10], and the more compli- molecular dynamics simulations [16]. They showed that
catedpolymeric[11]andglassysystemsontheother[12]. theseboundaryconditionswouldgiverisetoadesiredlin-
In this respect, the present article stands somewhere in ear, wedged velocity profile whilst avoiding the trouble-
between these two. some spatial inhomogeneities appearing when solid walls
Our paper is structured as follows. In the next section are used to induce the shear flow [21]. A particular
we briefly introduce the model and describe the bound- realisation of the LEBC on a cartesian simulation box
ary conditions for the imposition of shear. In Section III [0,Nx] [0,Ny] [0,Nz] is established by letting the pe-
× ×
we report on simulation data and conclude that shear- riodic images, at Nx < x 2Nx and Nx x 0,
≤ − ≤ ≤
thinning occurs for both G and MG mesophases, lead- move parallel to unit vectors ˆz, respectively, both with
±
ing to a transition to a mesophase consisting of tubu- speedU. TheLEBC,intheiroriginal,moleculardynam-
lar and ring-like structures as the strain increases. In ics form, are expressed as a Galilean transformation on
section IV we reveal how the presence of amphiphile in the position (x,y,z) and velocity (ξx,ξy,ξz) co-ordinates
lamellar mesophases induces the formation of rich inter- of a molecule, as follows
facialpatternssurvivingshearandallowshighervaluesof
x′ xmodN
x
stress than in lamellar mesophases without amphiphile. ≡
y′ ymodN
Finally, we provide our conclusions in Section V. y
≡
(z+∆ )modN , x>N ,
z z x
z′  zmodN , 0 x N , (1)
≡  z ≤ ≤ x
(z ∆ )modN , x<0,
II. THE MODEL AND THE LEES-EDWARDS z z
 −
BOUNDARY CONDITIONS ξ′ ξ
x ≡ x
ξ′ ξ
y ≡ y
We utilised an existing bottom-up lattice-Boltzmann
ξ +U , x>N ,
z x
(LB) model for amphiphilic fluids [14, 15], extended to ξ′  ξ , 0 x N , (2)
simulate shear flow by means of Lees-Edwards bound- z ≡  z ≤ ≤ x
ξ U , x<0,
z
ary conditions [16]. The model is in turn based on an  −
extension made to the Shan-Chen bottom-up LB model where∆ U∆tisthe image’sshiftattime ∆t afterthe
z
≡
forimmisciblefluidstomodelamphiphilic-fluidflow,and onset of shear.
employs25 microscopicvelocities,ofspeeds 0, 1and √2, An implementation of the LEBC on our LB dynamics
in three dimensions (D3Q25 lattice) [17, 18]. The model (LB-LEBC)differsfromthatusedinmoleculardynamics
usesaBGK(relaxation-time)approximationtothecolli- (MD-LEBC)in that the shift ∆ is notin generala mul-
z
sion term of the Boltzmann equation for fluid transport, tiple of the lattice unit, as Wagner and Pagonabarraga
3
havepointedout[21], andhenceaninterpolationscheme have previously found that the lattice size guaranteeing
is needed. This interpolation scheme streams the am- condition (b) is 1283 for the parameters generating the
phiphile dipoles d(x) and mass densities nα(x) located mesophases investigated here [14, 15].
k
at position x on the shearing wall, where c is the rele-
k
vant discrete molecular velocity, k = 1,...,25, for each
(fluid and amphiphilic) species α. III. SHEARING GYROIDAL MESOPHASES
In our LB-LEBC, while the spatial displacement fol-
lowsEqs.(1),thevelocityshiftcannotbeenforcedbyre- We sheared two gyroidal mesophases differing in the
placing the continuum velocity component ξz in Eqs. (2) amountof amphiphile presentand the value of the inter-
withthediscretemicroscopicspeedsck ˆz,sincetheveloc- amphiphileinteractioncouplingparameter. Eachofthese
·
ities ck areconstantvectors. Instead, this accelerationis structures was allowed to self-assemble from homoge-
enforced on the macroscopic fluid velocity around which neous mixtures of oil, water and amphiphile using pe-
thelocal-Maxwelliandistributesmoleculesatequilibrium riodic boundary conditions. They have been appropri-
and towardswhich the BGK scheme relaxes,similarly to ately characterised by probing direct and Fourier-space
how immiscibility forces are implemented [15, 20]. This late-time snapshots of the density order parameter φ
≡
procedure guarantees that all accelerations in the fluid ρoil ρwater;moreprecisely,theycorrespondtogyroid(cf.
−
are ruled by the same BGK process, controllable via the Fig. 5(a)) and molten gyroid mesophases, as previously
shapeofthedistributionfunctionandtherelaxation-time reported by us [14, 15].
parameter,includingtheaccelerationduetotheshearing The common parameters used for both gyroids were
walls. oil and water densities flatly distributed in the range
The MD-LEBC give rise, at steady state (late times), 0 < n(0)b = n(0)r < 0.7, coupling strengths g =0.08,
br
toashearstatewhichisGalilean-invariant,i.e. nopartic- g = 0.006, relaxation times τb = τr = τs = τd = 1,
bs
−
ularplaneinthesystemisfavouredoveranother. Thisis and, for the amphiphile’s dipoles, β =10 and d =1.
0
asine qua nonforanyshearingmethod, andourmethod Their differing parameters were surfactant densities,
satisfies it too. As regards the unsteady, transient initial flatly distributed in the initial homogeneous mixture, in
states, the MD-LEBC are unphysical since they cannot the ranges 0<n(0)s <0.9 for the gyroidand 0<n(0)s <
provide the molecular specificity (e.g. wall roughness) 0.6 for the molten gyroid, with coupling strengths g =
ss
required in an atomistic approach to boundary effects, 0.0045 for the gyroid and g = 0.003 for the molten
ss
− −
suchasdensitylayeringandslipatwall. However,meso- gyroid. These values for the gyroid are 50% higher than
scopic methods—LB is one of them—in general only de- those for the molten gyroid.
scribe low wavenumbers and frequencies, which means While the gyroid relaxes to a highly crystalline struc-
that, with respect to MD, (a) the atomistic detail of the ture [23], the molten gyroid shows both shorter-range
shearing walls is largely coarse-grained and (b) the fluid order and stronger temporal fluctuations than the for-
structure and dynamics are much less sensitive to the mer [15]. In order to obtain a sufficiently relaxed molten
atomistic detail of the walls. Since most boundary ef- gyroid as an initial condition for the shear, we took the
fects present in MD are absent in LB, the fact that the structure as evolved up to time step 32500; regarding
LE boundary conditions eliminate them does not pose the gyroid, the time slice chosen was time step 15000.
a problem. This should be taken with a caveat: our gy- For practical reasons, instead of letting the molten gy-
roidalmesophasesmeltwhenplacedinasolidbox,which roid self-assemble starting from a homogeneous initial
means that the approach to equilibrium is sensitive to mixture,weupscaledasmallermoltengyroid,previously
momentum transfer with the walls,and therefore the LE self-assembled using the same parameters on a 643 lat-
boundary conditions do not mimic shearinga mesophase tice[15],toa1283lattice. Upscalingconsistedinreplicat-
in confinement. (To our knowledge, no bottom-up sim- ing identicalcopies ofthe system: the periodic boundary
ulations have ever reported mesophase self-assembly in conditionsusedtogeneratethe 643 system(a)guarantee
confinement.) Rather, the LE boundary conditions in a thatthe densityfieldissmoothacrossthe replicabound-
LBmodelmimicshearingwithwallswhicharefarenough aries, yet, for this same reason, (b) produce a molten
from the locations in the system where observables are gyroid with an additional, undesirable long-wavelength
probed such that microscopic boundary effects are ab- fluctuation whose periodicity is half the lattice size. The
sent. amplitude of this undesired long-wavelength fluctuation
OurLEBCimplementationisembeddedwithinaneffi- relaxes in time to a vanishingly small value, fact which
cientparallelLBalgorithm[22]whichallowsustoemploy provides us with the 1283 mesophase we seek. We ob-
large lattices and hence reach the small Knudsen num- served, however, that this relaxation takes place in less
ber limit where (a) regions away from interfaces satisfy than1000timesteps[23],i.e. itisafasttransientwhich,
the incompressible NS equation in the limit of low Mach therefore, does not affect the shear response at the late
numbers(Ma)[20],and(b)observablesvarybylessthan times that we are interested in. In other words,the late-
10% when the lateral lattice dimension is doubled. We time shearresponseisinsensitive toasmallperturbation
4
in the initial condition. This allowed us to take the up- onlynearestneighbourinteractionsarebeingconsidered,
scaled,unrelaxedstructureastheinitialconditionforthe the virial term reduces to
molten gyroid.
1
Itis worthnotingthatwedidnotrequireanelongated 2 Xgαα¯Xψα(x)ψα¯(x+ck)ckck. (5)
aspect ratio for the lattices along the direction parallel α6=α¯ k
to the translationof the shearing walls since spatialden-
In the incompressible, low Mach number limit, our LB
sity fluctuations were much smaller than the lattice size.
model reproduces the NS equation away from inter-
This is not the case whenshearingphase-segregatingflu-
faces [17, 19], which describes a Newtonian fluid with
ids without an amphiphilic, growth-arresting species, as
a viscosity being a well known function of the relaxation
has been previously reported using LB lattices of up to
time. The presence of an interface, characterised by an
128:128:512 sizes and aspect ratio [24].
interfacialtensionandabending rigidity,however,intro-
ducesanisotropiesinthefluid’sstresstensorwhichcanbe
accounted for by a tensorial effective viscosity. Since the
A. Stress response and transients
interface may move, at a speed growing with the strain
rate,theseanisotropiescanbecomeunsteady. Ouraimis
Shearthinningissaidtooccurwhentheshearviscosity thentomeasurehowthisviscosityevolveswiththestrain
drops as the strain rate increases. For structured fluids and the strain rate.
such as those we study in this paper, the dynamic shear In order to probe the function η = η(γ˙) for both gy-
viscosity,η,isnotexpectedtobeaconstantofthe strain roidal mesophases, we measured P for a number of
xz
rate γ˙ ≡ 12(∂xuz+∂zux) as is true of Newtonian fluids, different applied shear rates. The chosen values for U
for which, were such that they remained within the incompressibil-
ity limit, i.e. small compared to the speed of sound on
Pxz =±2ηγ˙ , η 6=η(γ˙). (3) theD3Q25lattice,cs =3−1/2 ≈0.58. Valueschosenwere
U = 0.05,0.10,0.15,0.20, corresponding to Mach num-
Here P is one off-diagonal component of the pressure
xz bersMa U/c =0.086,0.17,0.26,0.34,respectively. All
s
(or stress) tensor, and the sign, by convention, indicates ≡
observables we report in this paper are spatial averages,
that the pressure is exerted by the fluid element on the
at least on x = const. planes where a simple fluid un-
surroundings(‘+’)orfromthelatterontheformer(‘ ’),
− der the same shear would show translational symmetry
respectively. We adhere in this paper to the secondcase.
for the velocity field, i.e., perpendicular to the velocity
In our simulations, we apply the steady shear described
gradient. Since, for reasonsof computationalcost, we do
inSectionII,i.e.,theshearisgeneratedbythetwoimage
not perform averages over the seed used to generate the
cells of the LB lattice located along the x-axis moving in
pseudo-random initial configuration mimicking a homo-
opposite directions. As a consequence,∂ u becomes the
x z geneous ternary mixture, we do not provide error bars
only non-vanishing component of the velocity gradient,
around averages.
which is also true for the P component of the stress
xz Figure1showsthe profileofthe stress,forthe sheared
tensor (and P , since the physical requirement that the
zx gyroid,alongtheappliedvelocitygradientdirection. Sev-
vorticity, W 1(∂ u ∂ u ), remains upper bounded
≡ 2 x z − z x eralcurvesthereindepictthetransportofmomentumto-
requires the stress tensor to be symmetric).
wardsthecore(i.e.,theplanex=64)asthestraingrows
As we have likewise done previously while computing
asafunctionoftime. Distinctively,theprofileshavespa-
diagonal components of the pressure tensor [14, 15, 20],
tialfluctuations,aconsequenceofthegyroid’sconvoluted
here we measured P from its definition as the sum of a
xz structure whose interfacial tension locally modifies the
kinetic term plus a virial mean-field term accounting for
viscosity expected for a simple fluid. The u component
z
interactions and giving rise to non-ideal gas behaviour,
of the velocity field, shown in Fig. 2 and averagedin the
namely,
samewayas statedfor P in the captionofFig.1,is
xz
h− i
however not inhomogeneous but follows a transient sim-
P(x) ρα(x)(c u(x))(c u(x))
≡ XX k k− k− ilar to that expected for a simple fluid: we observe the
α k setting up of a steady, smooth and wedge-shaped pro-
1
+ g ψα(x)ψα¯(x′)+ψα¯(x)ψα(x′) file, except at the borders. Figure 2 also includes the
4X αα¯Xh i×
α,α¯ x′ behaviour of the averaged velocity profile for the molten
(x x′)(x x′), (4) gyroidMGatlatetimes,andisseentomatchthatofthe
− −
gyroid G.
whereψhastheformψ 1 exp[ n(x)],whichsaturates RemainingwiththeGmesophase,weshowinFig.3the
≡ − −
athighdensityvaluesinordertoavoidunboundedinter- temporal evolution of the stress displayed in Fig. 1; the
particle forces whilst reproducing a meaningful equation values plotted are averages of the latter on the 8 x
≤ ≤
of state [15]. Since the interaction matrix g is sym- N 8 = 120 interval, which amounts to averaging over
αα¯ x
{ } −
metric with all diagonal elements identically zero, and thewholelatticeexceptthinslabsadjacenttothebound-
5
aries. In addition to Fig. 1, we include higher and lower lost any resemblance with the initial gyroid, except for
shear velocities, namely U =0.05, 0.15, 0.20. Note that the persistence of the toroidal rings, see Fig. 5(c), which
thetimeevolutionoftheaveragedstressisasuccessionof are defects in G. Also, the structure at ∆t = 21000 is
peaks and troughs, denoting successive intervals of yield essentially the same as that at time step ∆t = 5000—it
andrecoil,whichisacanonicalfeatureofviscoelasticbe- is a nonequilibrium steady state for at least the previous
haviour. Werethestrainrateatwhichthegyroiddeforms 16000timesteps,atimelongerthanthatrequiredforthe
coincidentwiththeappliedshearrate,thesecurveswould initialconfigurationtoself-assemblefromahomogeneous
implyshearthinning. Infact,whiletheincrementsinap- mixture of oil, water and amphiphile. The structure at
plied shear rate between these curves are kept constant, ∆t = 21000 consists of an irregular network of mainly
the increments in the (absolute) values of the stress at the same structuralelements characterisingthe defective
late times do not remain so but decrease. In Fig. 4 we regions before the onset of shear, namely, (a) elongated
showthestressaveragedovertimesteps24000to28000, tubules, withatendency to alignalonga directionwhich
plotted againstthe true strain rate, where the latter was is a linear combinationofdirections (1,0,0)and(0,0,1),
measured from the linear velocity profile generated at and(b)toroidal,ring-likestructures. Thisdescriptionis,
∆t 9000 (t 24000), as displayed in Fig. 2. Figure 4 by visual inspection, similar for every subvolume of the
≥ ≥
clearly shows shear thinning: the slope, i.e., the effective lattice visualised.
viscosity ηeff ∂Pxz/∂γ˙, decreases with the strain rate. Wealsolookedintothestructureoftheshearedmolten
≡
Figure4alsocontainstheanalogouscurveforthemolten gyroid at late times. In contradistinction to the gyroid’s
gyroid,whichshowsshearthinningforthelatteratlower state at high strain, showing tubules of shape similar to
strain rates than those at which the gyroid does, and at thatdepictedinFig.6andatananglewiththex=const.
higher intensity, i.e. planes, the highly strained molten gyroid displays tubes
∂ηeff ∂ηeff which are more stretched and aligned along the ˆz direc-
∂γ˙ (cid:12)(cid:12)molten < ∂γ˙ (cid:12)(cid:12)gyroid <0. (6) tion. The toroidal rings, also present for the molten gy-
(cid:12) (cid:12) roid before shear, represent a much smaller volume frac-
This is the first indication of shear thinning reported by
tion for the sheared molten gyroid than for the sheared
means of a bottom-up kinetic-theoretic model for fluid
gyroid.
flow.
Figure 7 shows the summed structure function
S(k,t), or scattering pattern, of the sheared gy-
Pky
B. Morphological transitions roid mesophase, showing stages of its plastic deforma-
tion. Here, S(k,t) is the structure function, computed
according to [15, 20]
Figure 5 shows the configuration of the gyroid in the
40 y 52slabof the 128:128:128lattice, before and
≤ ≤ ς 2
at late times after applying a shear of U = 0.20. The S(k,t) φ′(t) . (7)
volume-rendering graphicalrepresentationemployed [25] ≡ V (cid:12)(cid:12) k (cid:12)(cid:12)
(cid:12) (cid:12)
makes regions where φ 0.37 opaque to the lighting Here,kisthediscretewavevector,V isthelatticevolume,
≥
rays, assumed to shine normal to the plane of the text ς is the unit cellvolume forthe D3Q25lattice, andφ′(t)
k
and inwards; since 0.79 φ 0.79 over the entire is the Fourier transform of the fluctuations of φ. S(k,t)
− ≤ ≤
system, these regionsare the high-density locus ofone of is the Fourier transform of the autocorrelation function
thespecies(say,oil). Beforeshear,thestructurecontains for the order parameter,
highlyorderedsubvolumesofgyroidsymmetryanddiago-
nallength fromabout32 to 64 lattice sites, cf. Fig. 5(a). C (r,t) φ((x,t))φ(x+r,t) (8)
φφ
≡h i
This gyroid is hence a collection of subvolumes with a
regular tubular structure making up two three-fold co- where r is a vector lag and the brackets indicate an av-
ordinated,interweavingchirallattices ofwhichwedepict erage over the spatial coordinate x. Figures 7(a), (b)
onlyone. SincethesizeoftheGunitcellisapproximately and (d) are the xz ‘scattering patterns’ of the structures
5 to 6 lattice units, the depth (y-dimension) of the slabs in Fig. 5, produced by summing up the structure func-
shown in Fig. 5 is of about two gyroid unit cells. As can tion along the x direction. At ∆t = 1000 (not shown),
be seen in Fig. 5(a), the interfaces between these gyroid the maximum intensity is reduced to 29% of its value
subvolumes are defective regions where long-range order at ∆t = 0, while there appear horizontal ‘smeared out
andsymmetry appearsto be drasticallyreduced[14, 15]. filaments’ of very weak intensity, intrinsically related to
Two features characterising them are the spatial varia- the shearing process, as we shall conclude from Fig. 8.
tion in coordination number and chirality, seen by the At ∆t = 5000 a clear cardioid shape has developed; the
presence of elongated tubules and toroidal rings, cf. fig- factthatitpersistsfortherestofthesimulationconfirms
ure 6. ourobservationthatthesystemreachesasteadystateat
At ∆t = 21000, which is a late time after the onset time step ∆t = 5000. In addition, there is no trace of
of shear and we take as steady state, the structure has gyroidalpatterns along the x-direction.
6
Inordertoinvestigatethe originofthecardioidshape, The initial configuration employed was a cubic 1283
we computed the scattering pattern for a ‘synthetic gy- lattice with 16 lamellae, stacked perpendicularly to unit
roid’, vectorzˆ. Thelamellaewereofalternating,oil-watercom-
positions, separated by a thin monolayer of amphiphile;
G(x) sinqxcosqy+sinqycos(qz δ(x))+
the thickness of the immiscible and amphiphilic lamel-
≡ −
sin(qz δ(x))cosqx. (9) lae were 7 and 1 lattice sites, respectively. We popu-
−
lated each lattice site with a value of density kept con-
where δ(x) = (x N /2)δ is a spatially-varying de-
x max stant over the region corresponding to a given species;
−
phase used to obtain a linear strain on the morphology
each microscopic velocity is assigned the same fraction
(itsmaximumvalue,δ ,isreachedatthelatticebound-
max of this value. We gave amphiphilic regions the densities
aries), and q = const. is a wavenumber controlling the n(0)s =0, 0.80, 0.95,andoilandwaterregionsthe densi-
size of the surface’s unit cell. It is known that G(x) = 0 ties n(0)r =n(0)b =0.7. Shearwas appliedperpendicular
for δ 0 is a good approximationto the Schoen “G”
max tothelamellaewiththesameLEBCemployedinthelast
≡
triply periodic minimal surface of Ia3d cubic symmetry,
section, with speed U =0.10.
referred to as ‘the ideal gyroid’ hereafter [26]. Figure 8
Beforethe onset ofshear,the case without amphiphile
showsthescatteringpatternsfortheunstrainedmorphol-
for the lamellar initial condition just described is, a pri-
ogy and for dephases δ =8, 16.
max ori, a metastable state in our LB model. In fact, the
ComparingstructurefunctionmapsinFigs.7and8,at
structure has a stationary morphology since short-range
the same value of the strain rate, proves useful. For the
oil-water forces and the absence of fluctuations main-
synthetic gyroid, the strain is controlled by the number
tainimmiscibility, i.e. avalue forthe interfacesteepness,
ofunitcellsthatthedephasecausesthestructuretoshift
φ; however, a large enough perturbation in φ may al-
at the lattice boundary, following a linear profile as we |∇ |
lowafluctuationinsurfacetensionwhichdrivestheentire
approach the other boundary going through zero strain
interface to a radically different shape. Another factor
atthelatticecore. Foroursimulatedamphiphilicgyroid,
disrupting this lamellar morphology is shear, which may
however,thestraindoesnotfollowalinearprofileatearly
workagainsttheinterfacialtensionbyreducing φ;this
times; instead, the strain at time t would need to be |∇ |
can lead to miscibility (φ 0) for high enough strain
computed from the integral 1 t Nxdt′dx∂ u (x,t′), ≡
Nx R0 R0 x z rates. Despite these arguments,we observedstability for
where t′ is the time parameter. For the purposes of this
the sheared lamellar mesophase without amphiphile, as
paper, however, such an analysis would be superfluous;
we report next.
in fact, Fig. 8 already provides us with enough informa-
Figure 9 shows the stress as measured in the same
tion to understand the origin of the cardioid shape. For
fashion performed on the data plotted in Fig. 3, for sev-
all panels, (a), (b) and (c) therein, the position of the
eral amphiphile densities. The behaviour observed is di-
peaks at k =0 (k /(2π/N) 14, 15, where N =128)
z x ≈− verse. For zero amphiphile concentration (solid curve),
are invariant under the strain (dephase); not so with the
the stress reaches a peak at early times before it pro-
peaks at k = 0, which shift leftwards. (The shift would
z 6 ceeds to a second, lower maximum at late times, going
be rightwards were ∂ u < 0 or δ .) The shape of
x z max through a trough at intermediate times due to the fact
the maps in Figs. 7(c) and 7(d) is that of a transformed
that φ experiences a transient decrease.
scatteringpatternshiftedleftwards. This transformation |∇ |
The high-density regions of one of the immiscible
occursearly,between∆t=0and∆t=3000,andischar-
species(say,oil)isshowninFig.10(a)atlatetimes,∆t=
acterisedbytwo(strong,S 700)peakssimilarto those
≥ 8000; these are representative of the shape of the oil-
ofthegyroidatk =0,andtwo(weaker,200 S <700)
x
≤ water interface. Away from the boundaries (x=0, 128),
peaks at k =0.
z
thereisalargeinterfacialareawithzerocurvature,where
we define the curvature as H ∂2 x (z), x (z) being
≡ zz φ φ
IV. A SIMPLER CASE: SHEARING THE the curve resulting from projecting the φ = 0.18 surface
LAMELLAR MESOPHASE onto the xz plane. Curiously, we observe three changes
of curvature as we follow the curve x (z) for y =const.,
φ
Inthelastsectionwereportedonthegyroiddisplaying namely, H < 0,H > 0,H < 0,H > 0; instead, we would
higher shear stress than the molten gyroid. Since the have expected the steady, late-time configuration for the
structural transition between these two mesophases can sheared lamellar mesophase to rather minimise the in-
be driven by both the amphiphile density and the inter- terfacial area, leaving only one inflexion point. For fluid
amphiphile coupling parameter, as we have reported in regimes under conditions of local thermodynamic equi-
the past [15], our aim in this section is to elucidate the librium, we can think of H2 as an interfacial free energy
roleoftheamphiphiledensityaloneonthestressresponse densityassociatedwiththecurvature[27];inthiscase,we
toshear;wechoosethelamellarmesophaseasthesubject would have expected the steady, late-time configuration
of study, since this is the mesophase with the simplest to also minimise the interfacial free energy.
possible internal interface. The stress curve corresponding to n(0)s = 0.80
7
(cf. Fig. 9) shows the absence of large troughs, as occurs gyroidal morphology (G and MG) to a mesophase char-
for the n(0)s = 0 case, despite the fact that interfacial acterised by coexisting elongated tubules and toroidal,
tension is drastically reduced by the presence of the am- ring-like structures. The features of this mesophase is in
phiphile. In addition, the stress grows at late times to contrast to those of the mesophase reported by Zvelin-
higher values than those achieved by the n(0)s = 0 case. dovsky et al. using free-energy Langevin-diffusion meth-
The late-time order-parameterconfigurationis displayed ods by shearing a bicontinuous structure reminiscent of
in Fig. 10(b), showing a rich interfacial pattern. Using a molten gyroid [3]. The structure they found is of a
the same arguments as those of the last paragraph, this shorter-range ordering than that of the MG mesophase
structurecouldbecharacterisedbyahighercurvatureen- described here, and the high-strain structure consists of
ergy, d2xH2, where d2x is a measure on the oil-water coexisting lamellae and hexagonally packed tubes elon-
R
interface, and H is now defined as the inverse radius of gated along the direction of the imposed shear veloc-
curvature,parametrisedonthearclength,s. Figure10(b) ity. Our sheared mesophases also show enlongated tubes
shows similar regions of high curvature at an equal dis- along this direction, but the structure is far more com-
tance from the shearing walls, where u = const., which plicated than that found by Zvelindovsky et al. in that
z
weshallcallnodalplanes. Alsonotethattheinterface,as it exhibits remnant toroidal rings and ‘hard shoulders’
approximately depicted by the boundary of the φ 0.22 reminiscent of gyroidalskeletons; hexagonalpacking and
≥
volume,joinsthelatticeboundaryatananglecloseto90 coexisting lamellae are, on the other hand, absent.
degrees. The shear performs a plastic deformation which ef-
The stress curve for the n(0)s =0.95 case shows a dra- fectively breaks the links of the gyroidal skeleton; this
maticallydifferentsituationforthefirst5000timesteps: happens as these links interpose an (oil-water) interface
the presence of a trough, deeper than that present for whose normal, n, is parallel or anti-parallel to the flow
the n(0)s = 0 density. After that, there appears a shoot- streamlines, u. In other words, shear effectively applies
off whereby the stress rapidly grows and equals the late a ‘mixing’ force which is in competition with the inter-
time value achieved in the n(0)s = 0.80 case, while the particle forces keeping the mesophase in place, namely,
order-parameterdisplaysaconfigurationanalogoustothe those controlled by coupling parameters g , g and g .
br bs ss
n(0)s = 0.80 case, cf. Fig. 10(c). By looking at the am- Our hypothesis is that adsorbed dipoles sitting on in-
phiphile density field, ρs(x), for the case n(0)s =0.95, we terfacial regions at an angle θ ∠(u,n) other than
≡
observed that the high curvature regions arise close to θ = 0,π require more work from the shear forces to
the boundaries first (∆t<1000), and then rapidly move be drawn away from the interface than those regions on
away from them as the strain progresses. whichthestreamlineimpingesnormally,sincethemixing
force goes as cosθ. In particular, since the mixing force
vanishes for θ = π/2, considerably longer interfaces can
V. CONCLUSIONS survive the flow—shear induces a preferential direction
along which the long-range order present before the on-
In this paper we have reported on the shear stress set of shear is not reduced. This explains not only the
response of two gyroidal cubic amphiphilic mesophases formation of the elongated tubules but also their recon-
previously self-assembled using the same bottom-up LB nection (increase in coordination number). In fact, the
model we employ here, namely, the gyroid per se, G, toroidal,ring-like structuresare not only vestigialgyroid
whichshowshighcrystallinityatlateself-assemblytimes, defects which have survived the gradient u, but are
∇
and the molten gyroid,MG, endowedwith shorter-range also born anew as a result from reconnections. It is rele-
order and located within the sponge-gyroid lyotropic vanttopointoutthatPaddingandBoek,usingacoarse-
structural transition [15]. Shear was imposed via sliding grainedmoleculardynamicsmodelforwormlikemicelles,
periodic(Lees-Edwards)boundaryconditions,andwein- reportedon the formation of rings when applying steady
vestigatedtheresponsetoseveralvaluesofthestrainrate. shear to a wormlike micellar mesophase [6]—this is an
Inaddition,inordertoinvestigatethe dependence ofthe ‘amphiphile-in-water’ binary mesophase, in contrast to
shear stress on the amphiphile density, we also sheared the ‘oil-amphiphile-water’ ternary mesophases we study
a lamellar mesophase, of much simpler morphology than in this paper.
the gyroidal mesophases. By applying shear to a lamellar mesophase we found
We found that the gyroidal mesophases exhibit shear thatthepresenceofamphiphileontheoil-waterinterface
thinning, more pronounced and at lower strain rates for ofthemesophasecausestheinterfacetofoldintoawealth
the MG mesophase than the G mesophase. In other of structures with a (discrete) translational symmetry
words, momentum introduced into the system due to on planes equidistant to the shearing walls and along
shear is transported more easily for the mesophase con- the direction of the shear velocity. In other words, the
taining more amphiphile, with longer-rangeordering,i.e. inter-amphiphile force couples the adsorbed amphiphilic
the effective viscosity is higher for the G mesophase. dipolessothattheinterfacelocallyincreasesitscurvature
Wealsofoundashear-inducedtransitionfromaninitial energy density. It is worth investigating whether this
8
local increase is due to the amphiphile being incapable ening discussions. We acknowledge Iain Murray, Elena
of sustaining interfacial regions of low curvature under Breitmoser and Jonathan Chin for their involvement in
shear, i.e. whether the ‘breaking’ mechanism induced algorithm implementation and optimisation within the
by shearis counteractedby regionsof high curvature en- RealityGrid and TeraGyroid projects. This work was
ergydensity. Regardingtheshearstress,ouramphiphile- supported by the UK EPSRC under grant RealityGrid
containing lamellae responded with higher stress at late GR/R67699whichalsoprovidedaccesstoa512-processor
timesthanthosewithoutamphiphile. Thiscontrastswith SGI Origin3800 platform at Computer Services for Aca-
theresultsfoundforthegyroidalmesophases,andletsus demic Research (CSAR), Manchester Computing, UK.
conclude that it is the gyroid’s cubic morphology that WealsothanktheHigherEducationFundingCouncilfor
allows this structure to be stiffer. Understanding the be- England (HEFCE) for our on-site 16-node SGI Onyx2
haviour of the lamellar mesophase under shear requires graphical supercomputer.
the study of amphiphile self-assembly under shear, in-
cluding in and out of plane amphiphilic and associated
Marangoni currents, and their coupling to the imposed
flow.
VI. ACKNOWLEDGMENTS
We thank Dr. Rafael Delgado-Buscalioni, Prof. Anto-
nio Coniglio and Prof. Francesco Sciortino for enlight-
[1] A.V.M.Zvelindovsky,B.A.C.vanVlimmeren,G.J.A. [16] A.W.LeesandS.F.Edwards,J.Phys.C5,1921(1972).
Sevink,N.M.Maurits,andJ.G.E.M.Fraaije,J.Chem. [17] H. Chen, B. M. Boghosian, P. V. Coveney, and M.
Phys.109, 8751 (1998). Nekovee,Proc. R. Soc. London A 456, 2043 (2000).
[2] A. V. M. Zvelindovsky, G. J. A. Sevink, B. A. C. van [18] An example of the use of the D3Q25 lattice in lattice-
Vlimmeren,N.M.Maurits,andJ.G.E.M.Fraaije,Phys. gas methods is: P. J. Love, P. V. Coveney, and B. M.
Rev.E 57, R4879 (1998). Boghosian, Phys. Rev.E 64, 021503 (2001).
[3] A.V.M. Zvelindovsky,G.J. A. Sevink,and J. G. E. M. [19] S. Succi, “The lattice-Boltzmann equation–for fluid dy-
Fraaije, Phys. Rev.E 62, R3063 (2000). namics and beyond.” (Oxford University Press, Oxford,
[4] A. N. Morozov, A. V. M. Zvelindovsky, and J. G. E. M. 2001.)
Fraaije, Phys. Rev.E 64, 051803 (2001). [20] N. Gonz´alez-Segredo, M. Nekovee, and P. V. Coveney,
[5] A.N.MorozovandJ.G.E.M.Fraaije,Phys.Rev.E65, Phys. Rev.E 67, 046304 (2003).
031803 (2002). [21] A. J. Wagner and I. Pagonabarraga, J. Stat. Phys. 107,
[6] J. T. Padding and E. S. Boek, Phys. Rev. E 70, 031502 521 (2002).
(2004). [22] P. J. Love, M. Nekovee, P. V. Coveney, J. Chin, N.
[7] P. Mariani, V. Luzzati and H. Delacroix, J. Mol. Biol. Gonz´alez-Segredo, and J. M. R. Martin, Comp. Phys.
204, 165 (1988). Commun. 153(3), 340 (2003).
[8] S.-J. Marrink and D. P. Tieleman, J. Am. Chem. Soc. [23] J. Harting, M. J. Harvey, J. Chin and P. V. Coveney,
123, 12383 (2001). Comp. Phys. Commun. 165, 97 (2005); J. Harting, J.
[9] V. Luzzati, R. Vargas, P. Mariani, A. Gulik, and H. Chin, M. Venturoli, P. V. Coveney, “Large-scale lat-
Delacroix, J. Mol. Biol 229, 540 (1993). tice Boltzmann simulations of complex fluids: advances
[10] A. Xu, G. Gonella, and A. Lamura, Phys. Rev. E 67, throughtheadventofcomputationalgrids.” Phil.Trans.
056105 (2003). R. Soc. London Series A,in press (2005).
[11] R.D. Groot and W. G. M. Agterof, Macromolecules 28, [24] J. Harting,M.VenturoliandP.V.Coveney,Phil.Trans.
6284 (1995). R. Soc. London Series A 362, 1703 (2004).
[12] P.Sollich, Phys. Rev.E 58, 738 (1998). [25] For an explanation of volume rendering, see, e.g.,
[13] P. V. Coveney, Phil. Trans. R. Soc. Lond. A 361, 1057 W. Schroeder, K. Martin, B. Lorensen, “The Visualiza-
(2003). tionToolkit: AnObjectOrientedApproachto3DGraph-
[14] N.Gonz´alez-SegredoandP.V.Coveney,Europhys.Lett. ics.” Kitware Inc., ISBN 1930934076, 3rd Edition, 2003;
65, 795 (2004). URL:http://public.kitware.com/VTK/ .
[15] N. Gonz´alez-Segredo and P. V. Coveney, Phys. Rev. E [26] M.Wohlgemuth,N.Yufa,J.HoffmanandE.L.Thomas,
69,061501(2004);Virt.J.NanoscaleSci.&Tech.9,(Is- Macromolecules 34, 6083 (2001).
sue 23, June 14 2004), Supramolecular and Biochemical [27] R. Lipowsky, Nature 349, 475 (1991).
Assembly.
9
FIG. 1: Shear stress response of a gyroid mesophase along the
direction of the velocity gradient. As initial condition, we have
taken a gyroid on a NxNyNz = 1283 cubic lattice at time step
t=15000of self-assembly[14, 15]. The Lees-Edwardswallsmove
with speed U = 0.10 (Ma = 0.17). For each x coordinate, the
originalfieldhasbeenaveragedontheplane[1,Ny]×[16,Nz−16],
where the excluded interval on the z-axis accounts for wrapped-
rounddensities. Standarderrorsoftheaveragesareabout6×10−8
throughout, and arenot shown. Each linerepresents the response
at ∆t time steps after the start of steady shear: ∆t = 0 (dotted
line), ∆t=100 (dash-dotted), ∆t=800 (dashed) and ∆t=9000
(solid), where the last is ca. the time at which the core (i.e., the
plane x = 64) fully responds. From the figure we can see that
momentumtransferdecreasesasitreachesthecorefromthewalls.
Also,note that the stress inverts its signatlate times adjacent to
theboundaries,|x−x0|≤2(x0=0,128). Allquantitiesreported
areinlatticeunits.
10
FIG.2: Spatiallyaveragedvelocitycomponentuz forthemolten
gyroid and the gyroid mesophases sheared with U = 0.10, at late
times and over the x ≥ 64 half of the system. The dashed thin
and thick curves correspond to the molten gyroid at time steps
∆t=9000and13000,respectively. Thesolidthinandthickcurves
correspond to the gyroid at time steps ∆t = 9000 and 13000,
respectively. Theaverageisoverthesametwo-dimensionaldomain
as described in Fig. 1, for each x, and its standard error is shown
as negligibleerror bars. Note that the velocity shows amaximum
located from2to 4sites away fromthe boundary, unlikea simple
fluidwhichwoulddisplayitexactlyattheboundary. Thevalueof
thismaximumcoincideswiththeactualvelocityatwhichtheBGK
relaxation process of our LB model is forcing the fluid to move,
whichneedsnotcoincidewiththeinputparameterU =0.10. Note
that the inversion in the sign of the stress that we reported in
Fig. 1 occurs precisely for |x−x0|≤ 2, x0 = 0,128 and at (late)
times close to and after ∆t = 9000. The behaviour at the other
boundary region is similar and symmetric to that displayed here.
Allquantitiesreportedareinlatticeunits.