Oscillatory Shear Response of the Rigid Rod Model: Microstructural

Jun 18, 2019 - In this study, colloidal liquid crystals are described by the Doi–Hess model and their response to an external transient flow (oscill...
0 downloads 0 Views 1MB Size
Download PDF
Article Cite This: Macromolecules 2019, 52, 4907−4915

pubs.acs.org/Macromolecules

Oscillatory Shear Response of the Rigid Rod Model: Microstructural Evolution Marco De Corato*,† and Giovanniantonio Natale*,‡ †

Department of Chemical Engineering, Imperial College London, London W72AZ, United Kingdom Chemical and Petroleum Engineering, University of Calgary, 2500 University Drive NW, Calgary T2N 1N4, Canada

Downloaded via STOCKHOLM UNIV on August 22, 2019 at 07:24:24 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: In this study, colloidal liquid crystals are described by the Doi−Hess model and their response to an external transient flow (oscillatory shear) is investigated. The nonlinear partial differential equation governing the probability distribution of rods is solved numerically via an expansion in spherical harmonics and by solving a set of equivalent stochastic differential equations. These approaches provide microstructural insights into the behavior of systems with broken symmetries far from equilibrium. Thanks to this level of microstructural details, we here propose a new methodology to switch between the nematic and isotropic orientation state thanks to the transient nature of the externally imposed flow. Moreover, we show that the oscillatory shear flow is more efficient than the simple shear flow to capture the full microstructural dynamics that these systems can exhibit. Specifically, when a sinusoidal shear rate with large amplitude and low frequency (compared to the relaxation time of the system) is applied, a single oscillation period is sufficient to obtain a succession of microstructural dynamics (e.g., tumbling and wagging). These states would be obtained in a simple shear flow only by applying multiple shear rates to the system in multiple experiments. Hence, this methodology can provide a new strategy for experimental characterization.

■

3D.13 Later, it has been shown that the model also predicted a quite complex dynamical response in shear flow, and depending on the interplay between the intensity of the excluded volume potential and the applied shear flow, different dynamics can be obtained (log-rolling, wagging, and tumbling regimes).14−16 These theoretical predictions were also reported experimentally for the case of nematic solutions of poly(benzyl-glutamate). Burghardt and Fuller17,18 showed the transition from tumbling to shear alignment via rheo-optical techniques. Successively, log-rolling and wagging dynamics were captured experimentally19 and sustained oscillation of the shear stress component was also measured.20 However, only few experimental works were able to verify the predictions of the DH theory. This is mainly due to the presence in the experiments of spatial dependency of the director and textural aspects. These two phenomena contribute directly to the overall stress tensor via Frank elasticity contributions and viscous interactions between nematic domains. The theoretical development to include such issues in the original theory has also been proposed.21,22 More recently, the prediction of the model was found to quantitatively agree with the behavior of a roughly monodisperse suspension of fd-viruses in the nematic state.23−25 However, the authors also observed via optical microscopy texture formation during flow. Hence, even in the

INTRODUCTION Colloidal liquid crystals have been studied for their applications to gas-permeable membranes for molecular separations and large-area flexible display materials.1−3 Processing of liquid crystals presents rheological complexities given the intrinsic anisotropy of their constitutive elements (particles or molecules) and spatial variation of an average molecular orientation (director) in the bulk. Theoretically, the Doi−Hess (DH) model4,5 is successful in predicting correctly the microstructural dynamics and the macroscopic rheological response in the simple shear flow of a nematic system. This theory describes the evolution of a suspension of thin rigid rods dispersed in a Newtonian fluid and interacting via excluded volume potential (hard rods) under the effect of an external flow field. In static conditions, the DH model recovers the Onsager theory,6 and it can be used to describe the phase transition from isotropic to the nematic (I−N) state in the case of hard rods.6,7 In the presence of the flow, the DH model predicts that the rheological properties are highly nonlinear as a function of the Péclet number (Pe) defined as the ratio between the Brownian (τ = 1/Dr) and flow (1/γ̇) time scale, where γ̇ is the shear rate and Dr the rotary Brownian motion coefficient of the rods. In simple shear flow, the flow of nematic systems presents a characteristic change in sign from positive to negative of the first normal stress difference N1 increasing from low to intermediate shear rates and then back to positive at high shear rates.8−11 This rheological signature was also predicted by the DH model first in 2D12 and then in © 2019 American Chemical Society

Received: February 27, 2019 Revised: June 4, 2019 Published: June 18, 2019 4907

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules

can be used as a technique to switch between phases and control microstructures.

optimal conditions for the experimental system (monodispersed, no electrostatic, rigid particles, etc.), the matching between theoretical and experimental results is challenging. In oscillatory shear flow, the only work (to our knowledge) that explored the response of the DH model is by Dhont and Briels.7 The authors derived the DH theory from microscopic considerations, and they found a Maxwellian behavior of the viscoelastic response function in small amplitude oscillatory shear. At higher shear rates, they found that the viscoelastic response function for the nematic state does not have stationary solutions. However, they did not explore in detail the microstructural evolution predicted by the model in shear flow. Oscillatory shear flow can be considered as a model transient flow field which introduces a transient and periodic perturbation to the system.26 Specifically, a sinusoidal strain, γ(t) = γ0 sin(ωt), is imposed to the material, where γ0 is the strain amplitude and ω is the angular frequency. The amplitude of the associated shear rate is given by γ̇0 = γ0ω. Here, in addition to the flow time scale (which is equal to 1/ γ̇0), a second time scale (which is not present in simple shear flow) is introduced by the imposed frequency called the oscillation period, 2π/ω. For the case of a monodisperse Brownian suspension of elongated particles, the above two time scales (1/γ̇0 and 2π/ω) can be made dimensionless by combining them with τ. The ratio between the relaxation and the flow time scales is called the Weissember number, Wi = τγ̇0, whereas the Deborah number is the ratio of the relaxation time to the oscillation period, De = τω. If Wi and Wi/De are larger than 1, the system is driven out of equilibrium and the simple linear proportionality of the shear stress with the amplitude is lost.27−29 This regime is called large amplitude oscillatory shear (LAOS) and has recently attracted major interest as a characterization methodology to fully explore material properties. In the case of colloidal liquid crystals, LAOS coupled with small-angle X-ray scattering was used to study the time evolution of the microstructure of 4-cyano-40-octylbiphenyl (8CB) in the lamellar smectic A phase.30 The authors observed that new stable states of liquid crystal 8CB can be induced by the nonlinear shear conditions. A similar technique was also used to investigate the behavior of nanoclay suspensions in the nematic state.31 For these platelet suspensions, they identified a yielding transition between elastic and plastic deformation accompanied by a flipping of the director at large amplitude of the applied deformation. However, the interpretation at the microstructural level of the LAOS response is in general still challenging. Here, we investigate the response of the DH model in oscillatory shear flow, which can be a useful tool for colloidal suspensions. We perform numerical simulations of the DH equation for the molecular orientational distribution function using Brownian dynamics and an expansion in spherical harmonics. We consider the case of a single crystalline domain and we neglect the presence of textures and distortions, which could be relevant when the periodic flow attains small shear rates.32,33 Including spatial variations of the molecular orientation within the present model is the objective of future work. We show for the first time that oscillatory shear can be a more efficient characterization technique to capture the complex dynamical response of the nematic system. Moreover, we demonstrate that in specific regions of the phase diagram, oscillatory flow

■

RIGID ROD MODEL The nematic transition of a suspension of rigid rods is well described by the nonlinear DH model.4,5 This microstructural model provides an equation for the evolution of the orientation distribution function ψ(pα,t) in the presence of flow and accounting for interaction between rods via excluded volume potential Vex(pα)34

ij ∂ψ ∂ψ ∂ ∂ ∂ V yz = − α ·(ṗα ψ ) + α ·Dr jjj α + ψ α ex zzz j ∂t ∂p ∂p ∂p kBT z{ (1) k ∂p where Dr is the rotational Brownian motion coefficient and kB is the Boltzmann constant. The first term on the right represents the convective flux of the particle probability in the orientation space because of the torque induced by the external flow field. The second and third terms represent the effects of Brownian motion and of the excluded volume potential, respectively. The rotational velocity of a single particle ṗα follows Jeffery’s equation35 ṗα = Ω·pα + λ(E·pα − E: pα pα pα )

(2)

where E is the rate of strain tensor, Ω is the vorticity tensor, and λ is the particle form factor,36 which is equal to 1 for very long thin rods considered here. Without any loss of generality, we assume a velocity field v = γ̇(t)yx̂, which fixes the vorticity axis to ẑ, the flow direction to x̂, and the flow gradient to ŷ. The shear plane is thus given by the x−y plane. The Maier−Saupe formulation of the excluded volume potential15,16 is employed here, which can be regarded as the first term of a spherical harmonic expansion of the Onsager potential6 Vex(pα ) =

∫ C(pα , p β )ψ(p β ) dp β

C(pα , p β ) = const −

3 UkBT (pα ·p β )2 2

(3)

where U is an adimensional rod concentration (U = cd2L with c concentration, d rod diameter, and L rod length) and it represents the intensity of the excluded volume potential.34 It is important to underline here that the Maier−Saupe formulation of the excluded volume potential is used here with a constant rotational diffusion coefficient, while in the original rigid rod theory,34 the rotational diffusion coefficient is a function of the orientation state of the system. However, such approximation does not change the character and dynamical response of the system as previously shown by Larson and Ottinger.14 From the orientation distribution, it is possible to calculate the second-order orientation tensor, which is the second moment of the orientation distribution defined as ⟨pα pα ⟩ =

∫ pα pα ψ dpα

(4)

Any macroscopic properties (e.g., stress tensor, optical properties, and thermal and electrical conductivities) of the system depend on ⟨pαpα⟩.26,34,37 From this, the instantaneous order parameter tensor can be defined as S = ⟨pαpα⟩ − δ/3, where δ is the unit tensor. A useful scalar that is used here to characterize the state of the system is the scalar order 4908

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules parameter S = [(3/2)S:S]1/2, which is an invariant for the order parameter tensor. S tends to 1 when high level of rod alignment is present in the system and it is zero for an isotropic suspension. At equilibrium, if the scalar order parameter is different from zero, the system is in the nematic state.12−14 In the literature,7,38,39 the scalar order parameter is often defined as the largest eigenvalue of the order parameter tensor. The definition used here is a sufficient approximation given that we are working in LAOS. To solve the nonlinear problem numerically defined by eqs 1−3, two methods are employed here following the approaches used previously.14 The first method consists in a spherical harmonic expansion of the orientation distribution15 while the second method integrates directly the stochastic equation of motion of each rod y D ∂Vev ji αz zdt + dpα = jjjṗα − r zz α − 2Dr p z j ∂ k T p B k {

Figure 1. Equilibrium values of the scalar order parameter S as a function of U. Solid and dashed lines represent stable and unstable solutions, respectively.

Dr (I − pα pα )· dW (t )

the stable roots increases with U. These solutions represent the nematic state that forms at high dimensionless concentration U. With increasing intensity of the excluded volume potential, the suspensions have a lower and lower spreading around the nematic director, which represents the average molecular orientation. To characterize the relevant time scales present in this nonlinear model, one possibility is to investigate the rod orientation dynamics in the absence of external fields. This dynamics can be obtained by calculating the mean-square angular displacement (MSAD), defined as MSAD(t) = ⟨(p(t) − p(0))2⟩. Because the particles are rigid rods, the unit vector p(t) is confined on the surface of a unit sphere. On the singleparticle perspective, the rotational motion of a rod is governed by a nonlinear Langevin equation. In Figure 2, we report the MSADs obtained at different intensities of the nematic potential (U).

(5)

In the right-hand side of eq 5, the first bracket represents the deterministic contribution to the displacement of the orientation vector due to the shear flow and the excluded volume potential. The last term in the right-hand side represents the random displacements due to orientational diffusion, which is proportional to the increments of a vectorial wiener process dW(t).40 By computing the evolution of the orientation of a large number of rods, it is possible to approximate the dynamics of the entire population, ψ, and of its moments. More details about the stochastic differential equation (SDE) approach and the spherical harmonics expansion are reported in Appendix A and Appendix B, respectively. The two numerical methods have different advantages with one being more convenient than the other depending on the initial conditions and the dimensionless parameters. The results presented here can be obtained with either of the two methods, but care needs to be taken in truncating the number of spherical harmonics or the number of rods. Hence, we verified the results presented below with both strategies.

■

RESULTS In this section, we present the microstructural prediction of the rigid rod model in oscillatory shear flow. We first present the equilibrium properties that the model predicts and we analyze the Brownian orientation dynamics of a test rod in a nematic suspension at equilibrium. Then, we proceed in moving away from equilibrium conditions and we investigate the effect of the intensity of the shear flow (Wi) in LAOS. Finally, we conclude demonstrating a possibility to switch between microstructural states at low De numbers and explain the advantages of LAOS to experimentally probe nematic colloidal systems. Predictions at Equilibrium. In the absence of external flow, the DH model is controlled by the Brownian term and the excluded volume potential.34 Numerically, it is possible to extract the equilibrium solutions predicted by the model as a function of the intensity of the potential U. The solutions are summarized in Figure 1. Solid and dashed lines represent stable and unstable solutions, respectively. At equilibrium, the model predicts a single stable solution for U < 4.7, two stable solutions and one unstable solution for 4.7 < U < 5, and one stable solution for U > 5. For U > 5, the value of S identified by

Figure 2. Rotational mean-square displacement for different values of U.

In the absence of the excluded volume potential, (U = 0), the rods do not interact with each other and the Brownian motions of their orientation vectors are given by a random walk on the unit sphere.34,41 In this case, the MSAD of the rods has a closed form.41 In Figure 2, we show that the numerical simulations performed by setting U = 0 recover the theoretical solution. The MSAD reaches a plateau at long time because the motion of the orientation vector is confined to a sphere. Figure 2 shows that, in the case of a system in a nematic state (U > 5), the orientational dynamics is more complex. The MSADs present an initial linear behavior at low time (ballistic 4909

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules motion is absent because the particles are inertialess), an intermediate plateau, and a longtime growth plateau because of the orientational confinement on the surface of the unit sphere. At short times, the MSAD is the same as the isotropic case, regardless of the value of the excluded volume potential. For small time intervals, a rod does not interact with other neighboring rods and diffuses as if it was unconfined. At intermediate times, the excluded volume potential introduces a first confinement because the rod is “trapped” by the topological cage developed by the other particles. Interestingly, the intensity of the potential controls both the strength of the cage and the time scale associated with it. At longer time, a diffusive behavior reappears. This is representative of the coordinated Brownian motion of the entire population. The MSAD of a system at U = 6.66 is well described by the sum of two exponential function, MSAD(t) ≈ a(1 − e−4Drt/a) + (2 − a)(1 − e−bDrt), with the parameter a given by the value of the first plateau. The parameter b sets the time scale over which the collective diffusion of the rods allows the individual rods to explore the entire unit sphere. To draw a parallel with optical tweezers, this MSAD evolution is equivalent to that of a particle trapped in a potential well whose center is diffusing in time as well.42 In the cases of strong excluded volume potentials U = 10 and U = 13.32, the longtime behavior of the MSAD is no longer captured by an exponential growth with a single parameter, thus suggesting an anomalous collective diffusion.43 From Figure 2, it can also be concluded that in the nematic state, there are two effective time scales. The first corresponds to the rotary Brownian motion of a rod within its cage and it is given by D−1 r . The second is given by the collective diffusion of the rods and can be orders of magnitude larger than D−1 r in cases where the excluded volume, U, is much larger than 6.66. Predictions at Finite Wi. In this section, we investigate the microstructural perturbation to the system induced by the presence of oscillatory shear flow to shed some light in the microstructural evolution linked with viscoelastic characterization. From the equilibrium results, slow dynamics associated with the diffusion of the director is expected. Hence, numerical simulations are run for time of the order of 103Dr−1. In the absence of excluded volume potential, the microstructural dynamics was previously investigated.44−46 Wagging dynamics was observed at finite Wi numbers and at De number equal to or above Dr.45,46 In the presence of the excluded volume potential, the situation changes drastically as soon as we reach nematic state. In this section, all the results are obtained in nematic state U = 6.66 and at the onset of LAOS, Wi = De = 1. The results are qualitatively similar for systems in nematic state as long as Wi < U. Because the system is nematic, it is important to investigate the dynamics of the system for an average molecular initial orientation both in the plane of shear and outside because this can happen experimentally.14,15 In Figure 3, the initial orientation of the director is placed at θ = π/2 and ϕ = π/4 on the shear plane where θ and ϕ are the polar and azimuthal angle, respectively. Under the effect of oscillatory shear flow, the system presents large oscillations. The rods stay on the plane of shear and continue to oscillate between two positions around the flow direction. This dynamic is known as wagging and it is not expected in simple shear flow at the same flow conditions.15 The time periodic flow field traps the rods (and consequently the director) in this continuous periodic motion.

Figure 3. Initial condition θ = π/2 and ϕ = π/4.

The dynamics of the system changes drastically once the initial orientation is placed outside the shear plane. In Figure 4, the initial orientation is fixed at θ = π/4 and ϕ = π/4. The flow conditions and the intensity of the potential are here unchanged with respect to the one used in Figure 3. At short time (Figure 4 left), the system oscillates with the same period shown in Figure 3 dictated by the intensity of the flow field Wi. However, the amplitude of the oscillations of the components of ⟨pp⟩ is lower compared to that shown in Figure 3. This behavior represents the periodic orbits of the director outside of the shear plane. If the system is followed at very long time (Figure 4, right), an interesting dynamic is observed. The rods slowly migrate toward the vorticity direction and once they reach there, their oscillations are damped because the hydrodynamic torque becomes null. Hence, the vorticity direction represents an attractor point for the rods. In the absence of orientational diffusion, once a rod reaches that position, it will be there indefinitely: log-rolling regime. The excluded volume potential helps in attracting more and more particles in the vorticity eventually aligning the majority of the rods in that orientation. Not all the rods are aligned along the vorticity because there the hydrodynamic torque is small and thermal noise spreads the orientation of the rods around it. Such migration in the vorticity direction was also previously observed in simple shear flow but at different values of U.14,15 The difference here is induced by the fact that the intensity of the flow is dynamically changing with time and hence on average, the potential can become dominant. The migration toward the vorticity axis becomes faster with the increase of Wi until Wi ≈ U (data not shown). For values of Wi > U, the flow dominates over the interaction potential and the particles are always pulled in the shear plane. Predictions at De ≪ 1. It is interesting now to investigate the predictions of the model at low De numbers where elastic and viscous effects compete. The analysis at De ≫ 1 would provide only the elastic response of the model. Faraoni et al.15 predicted the presence of a biphasic coexisting region of a nematic and isotropic state at intensity of the potential U = 4.45 and for Pe values between 0.05 and 0.12. The authors find that a steady shear flow drives the system, which is isotropic at equilibrium, into a nematic state (para-nematic state). We probe the system response to an oscillatory shear flow around this biphasic state. In Figure 5, we report the model predictions at U = 4.45, Wi = 0.25, and De = 0.01. The components of the second-order tensor show oscillations between two values. Following ⟨pxpx⟩, it oscillates 4910

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules

Figure 4. Initial condition θ = π/4 and ϕ = π/4. Left figure shows the short time dynamics. Right figure shows the average location of the rods during the oscillations over a long period of time 3·103Dr.

De, the system oscillates around the nematic solution and is never able to relax to the isotropic state. From the work of Faraoni et al.,15 it is also expected that at U = 6.66 and in simple shear De = 0, a series of different dynamics are obtained by varying Wi. Moving from Wi = 0 to ∼9, an unstable tumbling solution is predicted coexisting with a stable log-rolling behavior that persists up to values of Wi = 15.82. In addition, a stable wagging branch is present for 9.50 < Wi < 16.15. This branch dies out to a flow-aligned stable solution for Wi > 16.15. The predictions of the DH model were verified experimentally in a more recent work.20 The persisted oscillations in shear stress were predicted by the theory; however, it was very challenging to capture this behavior experimentally as mentioned in the Introduction section. Inspired by the previous results in the biphasic state, we tested if the oscillatory shear flow at low De and at large Wi can capture in a single period all the dynamics predicted by the model. We fixed the intensity of the potential at U = 6.66 and De = 0.01 while the maximum intensity of the flow is Wi = 18 above the transition from wagging to flow aligned. The results for the components of the second-order tensor and the scalar order parameter are reported over one period in Figures 7 and 8, respectively. The dynamics reported in Figure 7 is predominately in the plane of shear. The components of ⟨pipj⟩ are aligned to the flow direction during most of the period but display oscillations with varying frequencies and amplitudes as the absolute value of the applied shear rate changes. Although the absolute value of the shear rate is decreasing during the oscillations, the system moves through wagging and full tumbling. These oscillations dampen once the shear rate has changed sign. The dynamics appears to be symmetric with respect to half of the period around zero shear rate. In Figure 8, we report the scalar order parameter over the period. The variation in S shows that the orientational spreading is modified because of the applied shear rate during the period. As expected, the highest value of S is obtained when the value of shear rate moves around zero, allowing the excluded volume potential to dominate. Once the shear rate approaches its maximum or minimum magnitude during the oscillation period, S reaches values around 0.8, again showing the competition between the potential and flow term. The possibility to capture the dynamics predicted by the model in a single cycle represents a unique opportunity for experimental characterization. Indeed, it should be possible to capture these oscillations by rheometry in a single experiment by simply working at low De numbers and high Wi.

Figure 5. Components of the orientation tensor showing oscillation between two states at U = 4.45.

between values of roughly 0.58 (partially aligned in the flow direction) to 1/3 (isotropic). An overshoot is clearly observable in the mixed component ⟨pxpy⟩ and it is located at the maximum value of the shear rate during the oscillations always before the aligned state. The overshoot is a clear contribution of the elastic terms (Brownian and the excluded volume potential) introducing an elastic recoil.12,21 To verify effectively if we are changing the orientational spreading of the particles and not only their average orientation, it is necessary to calculate the scalar order parameter S. In Figure 6, the oscillations between the two

Figure 6. Scalar order parameter showing oscillation between the nematic and the isotropic solution at U = 4.45.

states are clearly shown by S. Here, the oriented state S ≈ 0.5 is switched on and off by the externally applied oscillatory shear flow. This result is very relevant because it demonstrates a new methodology to switch the microstructure between phases with oscillatory shear flow. It is worth to underline that the switching is only possible at low De numbers. Low frequencies of the oscillations are necessary to allow the system sufficient time to adapt to the applied (transient) shear rate. At higher 4911

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules

Figure 7. Components of the orientation tensor over a full (a) and half (b) oscillation period showing the possibility to capture multiple dynamics via LAOS at U = 6.66, Wi = 18, and De = 0.01.

■

APPENDIX A

SDE Approach

Instead of directly solving the Smoluchowski equation, given by eq 1, one can solve the associated SDE.14,40 The solution of SDE yields a trajectory of the stochastic process p(t) belonging to the probability distribution ψ(p,t). By ensemble averaging over many trajectories, we approximate the evolution of the population ψ(p,t). The Ito version of SDE associated with the diffusion equation reads40 dp = A(p)dt + B(p) ·dW (t )

(6)

In eq 6, dW(t) represents a vectorial Wiener process whose amplitude is determined by the matrix B(p) that is the Cholesky decomposition of the diffusion tensor appearing in the Smoluchowski equation. The diffusion tensor is related to the SDE as

Figure 8. Scalar order parameter shown during the different dynamics at U = 6.66.

■

CONCLUSIONS We analyzed the microstructural response of the DH model to an externally imposed oscillatory shear flow. Numerical solutions of the nonlinear DH model were obtained by solving a set of SDEs and by using a spherical harmonic expansion. We demonstrated that oscillatory shear flow can be a powerful technique to probe the microstructural evolution of a nematic colloidal suspension. In addition, we showed that at low De, the system undergoes switchable phase transitions. This finding can be applied to develop new processing schemes to control the final microstructure in colloidal suspension-based materials. Finally, we demonstrated that the complex dynamics predicted by the model increasing the flow intensity Wi could be captured in a single period, provided that we worked at relatively low De numbers in oscillatory shear flow thanks to the transient nature of this flow field. These findings also support that LAOS can be used to facilitate the experimental work on the characterization of nematic colloidal suspensions. Indeed, only a single experiment would enable a complete characterization of the predicted flow instabilities. This is a major advantage because experimental conditions are often affected by the presence of poly-domains, making it challenging to capture the model-predicted behaviors. In addition, the model showed migration of the director for Wi ≈ U. This phenomenon can significantly impact the measurements of the rheological properties in LAOS because the system can be changing during the measurement. These results also motivate future works in investigating the link between these microstructural responses and the rheological response of nematic colloidal systems in LAOS. Including spatial variations of the molecular orientation, which can be relevant for time periodic flows, is also part of future work.

D(p) = Dr (I − pp)

(7)

The matrix B(p) appearing in eq 6 is defined as the Cholesky decomposition of the tensor D(p) D(p) = B(p) ·BT(p)

(8)

We remark that the tensor I − pp is singular; therefore, the decomposition eq 8 is not unique. Without any loss of generality, we choose the tensor B(p) as follows B(p) =

D(p) (I − pp)

(9)

It is straightforward to show that the above definition satisfies eq 8. The vector function A(p) contains the deterministic contribution to the orientation vector displacement which is given by D ∂Vev 1 1 1 ∂ A(p) = − ω· p + (γ ·̇ p + γ :̇ ppp) − r + · D(p) 2 2 kBT ∂p 2 ∂p (10)

The first two terms in eq 10 represent the orientational velocity induced by the time-dependent shear flow. The third term represents the effect of the excluded volume interaction between rods. The last term in the right-hand side is required for the Ito SDE eq 6 to be consistent with the Smoluchowski equation, given by eq 1.47,48 This term can be explicitly calculated using the definition of the diffusion tensor eq 7, which yields 1 ∂ ·D(p) = −2Dr p 2 ∂p 4912

(11) DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules The average of any quantity X(p) over the orientation distribution ψ(p,t) can be calculated from the population of rods as

To evaluate the convection inner product in eq 16 for simple shear with flow in the x direction and gradient in the y direction, we write the convection term as ∂ ·(ṗα ψ ) = γ ̇Γψ ∂pα

Nf

∫ X(p)ψ(p , t )dp = N1 ∑ X(pα ) f α=1

(12)

where Γ is the following operator

The initial orientation distribution, ψ(p,0), is approximated considering Nf rods. Equation 6 is discretized using a stochastic first-order explicit scheme (Euler−Maruyama)40 with time step Δt. The discretization introduces an error of 6(Δt 2) in the norm of pα, which will be accumulated during the time integration. We renormalize the orientation vector of each rod with every time step.47,48 To obtain a convergent approximation of ψ(p,t), we used Nf = 8 × 105 rods. The Brownian dynamics simulations were carried out in parallel over 10 or more cores.

■

λ ∂ Γ = −3λ sin(θ )2 sin(ϕ)cos(ϕ) + sin(2θ )sin(2ϕ) ∂θ 4 ÄÅ ÉÑ ÅÅi λ + 1 y Ñ ∂ 2Ñ j z Å j z + ÅÅj − sin(ϕ) ÑÑÑ z ÅÅÇk 2 { ÑÑÖ ∂ϕ (20)

Γ can be written in terms of spherical harmonic functions and the angular momentum operators. The angular momentum operators Lj are defined as ∂ zy ji Lj ≡ −ijjjpα × α zzz ; j ∂p z{ j k

APPENDIX B

Spherical Harmonic Solution

∞

Γ = −i −i

(13)

where bl,m is the complex function of time. The restriction to even values of l is a reflection of the even parity of ψ, namely, ψ(pα,t) = ψ(−pα,t). The results presented in this work are obtained with lmax = 16 truncating the summation in eq 13. Furthermore, because ψ is real, the following condition needs to hold bl , −m = ( −1)m bl*, m

L− ≡ Lx − iLy

(22)

Γ is then given by

bl , m(t )Ylm(θ , ϕ)

l = 0,even m =−l

(21)

L+ ≡ Lx + iLy ;

m=l

∑ ∑

j = x, y, z

Here, i is the imaginary unit. Furthermore, we define

In the present section, the Fokker−Planck equation (eq 1 in the main text) is reduced to a set of first-order ordinary differential equations by expanding the orientational distribution function in spherical harmonics and then applying a Galerkin procedure.13−15 The orientation distribution function can be written as ψ (pα , t ) =

(19)

+i

π π 6π λL−Y2−1 + i (Y2 2 − Y2−2) + i L+Y21 5 30 30

π π 2π λL+Y21 − i λLz(Y2−2 + Y2 2) L+Y2−1 − i 30 30 15 2 π 2 LzY 00 − i 3 3

π LzY 20 5 (23)

A shorthand notation is defined for the following operators for conciseness l′

[l , m , p , q , l′]− ≡

(14)

∑

⟨l , m|Y pq(L+ − L−)|l′, m′⟩bl ′ , m′

m ′= 0

The asterisk denotes a complex conjugate. We define the inner product as ⟨l , m|X ⟩ =

∫

Ylm *X

(24) l′

[l , m , p , q , l′]+ ≡

α

dp

(15)

(25) l′

ÄÅ ÉÑ ∂ ÅÅ ∂ψ ÑÑ ∂ D ⟨l , m|ψ ⟩ = −⟨l , m| α · (ṗα ψ )⟩ + Dr⟨l , m| α · ÅÅÅÅ α ÑÑÑÑ⟩ ∂p ∂p ÅÅÇ ∂p ÑÑÖ Dt ÉÑ ÅÄÅ Dr ÅÅ ∂Vex ÑÑÑ ∂ + ⟨l , m| α · ψ ÅÅÅ α ÑÑÑ⟩ ÅÇÅ ∂p ÑÖÑ ∂p kBT (16)

[l , m , p , q , l′]z ≡

∂ ∂ψ ⟨l , m| α · α ⟩ = −l(l + 1)bl , m ∂p ∂p

∑

⟨l , m|Y pqLz|l′, m′⟩bl ′ , m′

m ′= 0

(26)

l′

[l , m , p , q , l′]0 ≡

∑

⟨l , m|Y pq|l′, m′⟩bl ′ , m′

m ′= 0

(27)

where the angular momentum operators act on |l′,m′⟩, which denotes Ylm′ ′ . These operators can be simplified using the following identities

By developing the inner products in eq 16, the set of equations for the bl,m is derived. Some of the inner products are easily calculated by exploiting the orthogonality of the spherical harmonic functions over the unit sphere

∫ Ylm*ψ dpα = bl ,m

⟨l , m|Y pq(L+ + L−)|l′, m′⟩bl ′ , m′

m ′= 0

where X is some expression involving ψ. Multiplying eq 1 by Yml * and integrating over the unit sphere, we obtain

⟨l , m|ψ ⟩ =

∑

Lz|l′, m′⟩ = m′|l′, m′⟩

(28)

L+|l′, m′⟩ = [l′(l′ + 1) − m′(m′ + 1)]1/2 |l′, m′ + 1⟩

(17)

(29)

L−|l′, m′⟩ = [l′(l′ + 1) − m′(m′ − 1)]1/2 |l′, m′ − 1⟩

(18)

(30) 4913

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules At this point, we can finally write the convective term in eq 16 as

Now, we can expand the integrals involved in eq 35, and as an example, we focus on the first integral of the right side of eq 35

1 ∂ − ⟨l , m| α ·(ṗα ψ )⟩ = −⟨l , m|Γψ ⟩ γ̇ ∂p ∞ π = ∑ −i [l , m , 2, −1, l′]− 30 l ′= 0,even π 2π −i [l , m , 2,1, l′]− − iλ [l , m , 2, −2, l′]z 30 15 2π 2 π [l , m , 0, 0, l′]z − iλ [l , m , 2, 2, l′]z + i 15 3 2 2π π [l , m , 2,0, l′]z + 3iλ −i [l , m , 2, −2, l′]0 3 15 2π − 3i [l , m , 2,2, l′]0 (31) 15

ÅÄÅ 2 ÅÅ ÅÅ ∑ ÅÅ ÅÅ ÅÅÇ l ′= 0,even

=

=

∫

∫ (33)

■

(−)

Lz[Ylm *]ψLz[Vex ]dpα =

p

l′

2

∑ ∑ ∑

∑

(− mm′)c(l′)

(40)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (M.D.C.). *E-mail: [email protected] (G.N.). ORCID

Marco De Corato: 0000-0002-9361-4794 Giovanniantonio Natale: 0000-0001-5331-0932 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Pier Luca Maffettone for very useful discussions. G.N. acknowledges funding from NSERC grant no. RGPIN-2017-03783. The authors declare no competing financial interest.

■

c(l′)Ylm′ ′*(pα )Ylm′ ′(p β )c(l′)

REFERENCES

(1) Kikuchi, H.; Nishiwaki, J.-i.; Kajiyama, T. Mechanism of electrooptical switching hysteresis for (polymer/liquid crystal) composite films. Polym. J. 1995, 27, 1246. (2) Kajiyama, T.; Miyamoto, A.; Kikuchi, H.; Morimura, Y. Aggregation states and electro-optical properties based on light scattering of polymer/(liquid crystal) composite films. Chem. Lett. 1989, 18, 813−816. (3) Miyamoto, A.; Kikuchi, H.; Kobayashi, S.; Morimura, Y.; Kajiyama, T. Dielectric property-electrooptical effect relationships of polymer/liquid-crystal composite films. Macromolecules 1991, 24, 3915−3920.

(36)

By exploiting the orthonormal properties of the spherical harmonics, the excluded volume potential takes the following form m=l

l = 0,even m =−l

∑

The sum of terms in eqs 38−40 (taking into account eq 35) is the contribution due to the nematic term defined in eq 32.

In the Maier−Saupe formulation of the excluded potential, the term |pα × pβ|2 can be rewritten as

2

l′

2

bp , qbl ′ , m′⟨l , m|p , q|l′, m′⟩

(35)

∑ ∑

p

∑ ∑ ∑

p = 0,even q =−p l ′= 0,even m ′=−l′

∫ Lx[Ylm*]ψLx[Vex ]dpα + ∫ Ly[Ylm*]ψLy[Vex ]dpα 1 1 = ∫ L+[Ylm *]ψL−[Vex ]dpα + ∫ L−[Ylm *]ψL+[Vex ]dpα 2 2

Vex =

L−[Ylm *]ψL+[Vex ]dpα =

∞

By using the identities eq 22, we have

ÑÑ ÑÑ ÑÑ ÑÑÖ

c(l′)bp , qbl ′ , m′⟨l , m − 1|p , q|l′, m′ − 1⟩

The term related to z-component in eq 34 is

(34)

yzÅÅÅ 2( −1) zzÅ 1 {ÅÅÅÅÇ 3

(− )[l(l + 1) − m(m − 1)]1/2

∑

[l(l + 1) − m(m − 1)] [l′(l′ + 1) − m′(m′ − 1)]1/2 c(l′)bp , qbl ′ , m′⟨l , m + 1|p , q|l′, m′ + 1⟩ (39)

∫ Ly[Ylm*]ψLy[Vex ]dpα + ∫ Lz[Ylm*]ψLz[Vex ]dpα

i 4π = jjj k 2l′ +

l′

2

∑ ∑ ∑

p = 0,even q =−p l ′= 0,even m ′=−l′ 1/2

∫

l ′= 0,even m =−l′ É ÅÄÅ l ′ /2 Ñ Ñ

c(l′)bl ′ , m′[l′(l′ + 1) − m′(m′ − 1)]1/2 Ylm′ ′− 1 dpα

∞

∫

∑

l′

∑

Analogously, we can calculate the second integral in eq 35

ÄÅ ÉÑ ij ∂Y m * yz ji ∂V zy ∂ ÅÅÅ ∂Vex ÑÑÑ Å ⟨l , m| α ·ÅÅψ α ÑÑÑ⟩ = − jjjjpα × l α zzzz·ψ jjjpα × exα zzzd j Å Ñ ∂p ÅÇ ∂p ÑÖ ∂p { k ∂p z{ k pα = Lx[Ylm *]ψLx [Vex ]dpα

m = l′

bp , qY pq

(38)

we obtain

2

2

∑

[l′(l′ + 1) − m′(m′ − 1)]

By invoking the identity

∑

p

∞

∫ (−)m [l(l + 1) + m(−m + 1)]1/2 Yl−m+1 ∑ ∑

p = 0,even q =−p l ′= 0,even m ′=−l 1/2

∫

|pα × p β |2 =

l′

∑

l ′= 0,even m ′=−l′ ∞ p

∫

+

ÉÑ ÑÑ Ñ c(l′)bl ′ , m′(t )Ylm′ ′(pα )ÑÑÑÑdpα ÑÑ m ′=−l′ ÑÑÖ

bp , q(t )Y pq(pα )L−

p = 0,even q =−p

p = 0,even q =−p

We next evaluate the term containing the excluded volume potential in eq 16 ÄÅ ÑÉ ÅÄÅ ÑÉ ∂Vex ÑÑÑ α ∂ ÅÅÅ ∂Vex ÑÑÑ m* ∂ Å Å Å Ñ Å ÑÑdp ⟨l , m| α ·ÅÅψ α ÑÑ⟩ = Yl ·Åψ ∂p ÅÅÇ ∂p ÑÑÖ ∂pα ÅÅÅÇ ∂pα ÑÑÑÖ ÄÅ ÉÑ ∂Ylm * ÅÅÅ ∂Vex ÑÑÑ α Å ÑÑdp =− ·Åψ ∂pα ÅÅÅÇ ∂pα ÑÑÑÖ (32) ij α y yi ∂f ∂g jp × ∂f zzzjjjpα × ∂g zzz j z j α α = j α α j ∂p ∂p ∂p z{jk ∂p zz{ k

p

∞

∫ L+[Ylm*]ψL−[Vex ]dpα = ∫ (−)m L+[Yl−m] ∑ ∑

c(l)blmYlm(pα ) (37) 4914

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915

Article

Macromolecules (4) Hess, S. Fokker-Planck-equation approach to flow alignment in liquid crystals. Z. Naturforsch., A: Phys. Sci. 1976, 31, 1034−1037. (5) Doi, M. Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid crystalline phases. J. Polym. Sci., Polym. Phys. Ed. 1981, 19, 229−243. (6) Onsager, L. The effects of shape on the interaction of colloidal particles. Ann. N.Y. Acad. Sci. 1949, 51, 627−659. (7) Dhont, J. K. G.; Briels, W. J. Viscoelasticity of suspensions of long, rigid rods. Colloids Surf., A 2003, 213, 131−156. (8) Kiss, G., Porter, R. S. Rheology of concentrated solutions of poly (benzylglutamate). Journal of Polymer Science: Polymer Symposia; Wiley Subscription Services, Inc., A Wiley Company: New York, 1978; Vol. 65, 1, pp. 193−211. (9) Kiss, G.; Porter, R. S. Rheology of concentrated solutions of helical polypeptides. J. Polym. Sci., Polym. Phys. Ed. 1980, 18, 361− 388. (10) Navard, P. Formation of band textures in hydroxypropylcellulose liquid crystals. J. Polym. Sci., Part B: Polym. Phys. 1986, 24, 435− 442. (11) Moldenaers, P.; Mewis, J. Transient Behavior of Liquid Crystalline Solutions of Poly(benzylglutamate). J. Rheol. 1986, 30, 567−584. (12) Marrucci, G.; Maffettone, P. L. A description of the liquidcrystalline phase of rodlike polymers at high shear rates. Macromolecules 1989, 22, 4076−4082. (13) Larson, R. G. Arrested tumbling in shearing flows of liquidcrystal polymers. Macromolecules 1990, 23, 3983−3992. (14) Larson, R. G.; Ottinger, H. C. Effect of molecular elasticity on out-of-plane orientations in shearing flows of liquid-crystalline polymers. Macromolecules 1991, 24, 6270−6282. (15) Faraoni, V.; Grosso, M.; Crescitelli, S.; Maffettone, P. L. The rigid-rod model for nematic polymers: an analysis of the shear flow problem. J. Rheol. 1999, 43, 829−843. (16) Grosso, M.; Keunings, R.; Crescitelli, S.; Maffettone, P. L. Prediction of chaotic dynamics in sheared liquid crystalline polymers. Phys. Rev. Lett. 2001, 86, 3184. (17) Burghardt, W. R.; Fuller, G. G. Transient shear flow of nematic liquid crystals: Manifestations of director tumbling. J. Rheol. 1990, 34, 959−992. (18) Burghardt, W. R.; Fuller, G. G. Role of director tumbling in the rheology of polymer liquid crystal solutions. Macromolecules 1991, 24, 2546−2555. (19) Mewis, J.; Mortier, M.; Vermant, J.; Moldenaers, P. Experimental evidence for the existence of a wagging regime in polymeric liquid crystals. Macromolecules 1997, 30, 1323−1328. (20) Grosso, M.; Crescitelli, S.; Somma, E.; Vermant, J.; Moldenaers, P.; Maffettone, P. L. Prediction and observation of sustained oscillations in a sheared liquid crystalline polymer. Phys. Rev. Lett. 2003, 90, 098304. (21) Marrucci, G., Greco, F. Flow behavior of liquid crystalline polymers. Advances in Chemical Physics; Wiley, 1993; Vol. 86, pp 331−404. (22) Tao, Y.-G.; den Otter, W. K.; Briels, W. J. Kayaking and wagging of rods in shear flow. Phys. Rev. Lett. 2005, 95, 237802. (23) Lettinga, M. P.; Dhont, J. K. G. Non-equilibrium phase behaviour of rod-like viruses under shear flow. J. Phys.: Condens. Matter 2004, 16, S3929. (24) Lettinga, M. P.; Dogic, Z.; Wang, H.; Vermant, J. Flow behavior of colloidal rodlike viruses in the nematic phase. Langmuir 2005, 21, 8048−8057. (25) Lagerwall, J. P. F.; Scalia, G. A new era for liquid crystal research: Applications of liquid crystals in soft matter nano-, bio- and microtechnology. Curr. Appl. Phys. 2012, 12, 1387−1412. (26) Bird, R. B.; Armstrong, R. C.; Hassager, O.; Curtiss, C. F.; Middleman, S. Dynamics of Polymeric Liquids, vol. 2. Phys. Today 1978, 31, 54. (27) Khair, A. S. Large amplitude oscillatory shear of the Giesekus model. J. Rheol. 2016, 60, 257−266.

(28) Khair, A. S. On a suspension of nearly spherical colloidal particles under large-amplitude oscillatory shear flow. J. Fluid Mech. 2016, 791, R5. (29) Hyun, K.; Wilhelm, M.; Klein, C. O.; Cho, K. S.; Nam, J. G.; Ahn, K. H.; Lee, S. J.; Ewoldt, R. H.; McKinley, G. H. A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS). Prog. Polym. Sci. 2011, 36, 1697−1753. (30) Struth, B.; Hyun, K.; Kats, E.; Meins, T.; Walther, M.; Wilhelm, M.; Grübel, G. Observation of new states of liquid crystal 8CB under nonlinear shear conditions as observed via a novel and unique rheology/small-angle x-ray scattering combination. Langmuir 2011, 27, 2880−2887. (31) Lettinga, M. P.; Holmqvist, P.; Ballesta, P.; Rogers, S.; Kleshchanok, D.; Struth, B. Nonlinear behavior of nematic platelet dispersions in shear flow. Phys. Rev. Lett. 2012, 109, 246001. (32) Marrucci, G.; Greco, F. A molecular approach to the polydomain structure of LCPs in weak shear flows. J. Non-Newtonian Fluid Mech. 1992, 44, 1−13. (33) Greco, F.; Marrucci, G. Molecular structure of the hedgehog point defect in nematics. Mol. Cryst. Liq. Cryst. 1992, 210, 129−141. (34) Doi, M., Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press, 1988; Vol. 73. (35) Jeffery, G. B. The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. London, Ser. A 1922, 102, 161−179. (36) Bretherton, F. P. The motion of rigid particles in a shear flow at low Reynolds number. J. Fluid Mech. 1962, 14, 284−304. (37) Fuller, G. G. Optical Rheometry of Complex Fluids; Oxford University Press on Demand, 1995. (38) Dhont, J. K. G.; Briels, W. J. Isotropic-nematic spinodal decomposition kinetics. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2005, 72, 031404. (39) Sonnet, A. M.; Virga, E. G.; Durand, G. E. Dielectric shape dispersion and biaxial transitions in nematic liquid crystals. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2003, 67, 061701. (40) Ö ttinger, H. C. Stochastic Processes in Polymeric Fluids: Tools and Examples for Developing Simulation Algorithms; Springer Science Business Media, 2012. (41) Dhont, J. K. An Introduction to Dynamics of Colloids; Elsevier, 1996; Vol. 2. (42) Khan, M.; Mason, T. G. Trajectories of probe spheres in generalized linear viscoelastic complex fluids. Soft Matter 2014, 10, 9073−9081. (43) Metzler, R.; Klafter, J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 2000, 339, 1−77. (44) Férec, J.; Heuzey, M. C.; Ausias, G.; Carreau, P. J. Rheological behavior of fiber-filled polymers under large amplitude oscillatory shear flow. J. Non-Newtonian Fluid Mech. 2008, 151, 89−100. (45) Natale, G.; Reddy, N. K.; Prud’homme, R. K.; Vermant, J. Orientation dynamics of dilute functionalized graphene suspensions in oscillatory flow. Phys. Rev. Fluid. 2018, 3, 063303. (46) Bechtel, T. M.; Khair, A. S. Nonlinear viscoelasticity of a dilute suspension of Brownian spheroids in oscillatory shear flow. J. Rheol. 2018, 62, 1457−1483. (47) Cobb, P. D.; Butler, J. E. Simulations of concentrated suspensions of rigid fibers: Relationship between short-time diffusivities and the long-time rotational diffusion. J. Chem. Phys. 2005, 123, 054908. (48) De Corato, M.; Greco, F.; D’Avino, G.; Maffettone, P. L. Hydrodynamics and Brownian motions of a spheroid near a rigid wall. J. Chem. Phys. 2015, 142, 194901.

4915

DOI: 10.1021/acs.macromol.9b00407 Macromolecules 2019, 52, 4907−4915