Distortion and flow of nematics simulated by dissipative particle dynamics Tongyang Zhao and Xiaogong Wang Citation: The Journal of Chemical Physics 140, 184902 (2014); doi: 10.1063/1.4873699 View online: http://dx.doi.org/10.1063/1.4873699 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/18?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Static and dynamic properties of magnetic nanowires in nematic fluids (invited) J. Appl. Phys. 97, 10Q304 (2005); 10.1063/1.1852171 Molecular dynamics simulations of a ferroelectric nematic liquid under shear flow J. Chem. Phys. 117, 8551 (2002); 10.1063/1.1512275 Erratum: “A theory for flowing nematic polymers with orientational distortion” [J. Rheol. 44, 1085–1101 (2000)] J. Rheol. 44, 1435 (2000); 10.1122/1.1323466 A theory for flowing nematic polymers with orientational distortion J. Rheol. 44, 1085 (2000); 10.1122/1.1289278 Shear flow simulations of biaxial nematic liquid crystals J. Chem. Phys. 107, 3144 (1997); 10.1063/1.474666

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

THE JOURNAL OF CHEMICAL PHYSICS 140, 184902 (2014)

Distortion and flow of nematics simulated by dissipative particle dynamics Tongyang Zhao and Xiaogong Wanga) Laboratory of Advanced Materials (MOE), Department of Chemical Engineering, Tsinghua University, Beijing 100084, People’s Republic of China

(Received 3 March 2014; accepted 17 April 2014; published online 8 May 2014) In this study, we simulated distortion and flow of nematics by dissipative particle dynamics (DPD). The nematics were modeled by a binary mixture that contained rigid rods composed of DPD particles as mesogenic units and normal DPD particles as solvent. Elastic distortions were investigated by monitoring director orientation in space under influences of boundary anchoring and external fields. Static distortion demonstrated by the simulation is consistent with the prediction of Frank elastic theory. Spatial distortion profile of the director was examined to obtain static elastic constants. Rotational motions of the director under influence of the external field were simulated to understand the dynamic process. The rules revealed by the simulation are in a good agreement with those obtained from dynamical experiments and classical theories for nematics. Three Miesowicz viscosities were obtained by using external fields to hold the orientation of the rods in shear flows. The simulation showed that the Miesowicz viscosities have the order of ηc > ηa > ηb and the rotational viscosity γ 1 is about two orders larger than the Miesowicz viscosity ηb . The DPD simulation correctly reproduced the non-monotonic concentration dependence of viscosity, which is a unique property of lyotropic nematic fluids. By comparing simulation results with classical theories for nematics and experiments, the DPD nematic fluids are proved to be a valid model to investigate the distortion and flow of lyotropic nematics. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4873699] I. INTRODUCTION

Dissipative particle dynamics (DPD) is a mesoscopic particle-based simulation method originally proposed to study the hydrodynamics of complex fluids.1, 2 The main advantage of DPD method is that this particle-based simulation method is able to model both thermal fluctuations and hydrodynamic interactions. It has been proven that mimetic fluids composed of DPD particles are able to reproduce the fluid mechanics of simple fluids, which are described by the Navier-Stokes equations in the large scale limit.3–6 By introducing different kinds of inter-particle links, DPD method has been extended to model various kinds of complex fluids, such as colloidal suspensions,7–9 polymer solutions,10–14 and polymer melts.15–18 In addition to the isotropic fluids, DPD method can also be used to simulate anisotropic systems. DPD fluids containing rodlike chains are able to form anisotropic phases such as nematic and smectic liquid crystals (LCs).19–22 We have used binary mixtures composed of rodlike chains and normal DPD particles to model lyotropic LCs, which are in nematic phase under certain concentration and temperature range.23–26 Theoretically, DPD method exhibits the potential to model nematic fluids and can be further explored to simulate other important characteristics of the anisotropic fluids. Two of the most important responses of nematics to the external forces are distortion and flow, which have been theoretically studied for years.27 In the continuum theory, a field of the director n(r) is used to represent the direction of preferred orientation of molecules in the neighborhood of a point a) Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2014/140(18)/184902/10/$30.00

r. A state of static distortion is described by a spatially varied vector field of n(r). According to the Frank elastic theory, the free energy density (per cm3 of a nematic material) of the spatial distortions of n(r) can be given by the following equation: 1 1 1 K1 (∇ · n)2 + K2 (n · ∇ × n)2 + K3 (n × ∇ × n)2 , 2 2 2 (1) where K1 , K2 , and K3 are the splay, twist, and bend elastic constants, respectively.27 In this equation and the following discussions, no distinction will be made between n(r) and n, where n will be understood to denote the vector field n(r) at representative point r. The equilibrium state corresponds to the minimum of the total distortion energy considering the effects of boundary anchoring and external fields. Rotational motion of the director under the influence of external fields has been theoretically treated as a dynamic process with energy dissipation.27 In addition to the distortion and rotation, the translational motion of nematics is another important dynamic process. Flow of a nematic fluid is a more complicated issue compared with flow of isotropic liquids. One of the reasons is that the translational motions are coupled with orientational motions of the molecules in nematic fluids. For the nematic fluids, the flow behavior is usually described by a series of phenomenological equations.27 The relationship between viscous stress tensor σ  αβ and flow is related to n by Leslie-Ericksen equation:28–30 F =

 σαβ = α1 nα nβ nμ nρ Aμρ + α4 Aαβ + α5 nα nμ Aμβ

+α6 nβ nμ Aμα + α2 nα Nβ + α3 nβ Nα .

(2)

In this equation, A is the symmetric part of the velocity gradient tensor and N represents the rate of change of n with

140, 184902-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-2

T. Zhao and X. Wang

J. Chem. Phys. 140, 184902 (2014)

respect to the background fluid: N=

Dn − ω × n, Dt

Newton’s equations of motion: (3)

where ω is the anti-symmetric part of the velocity gradient tensor. The operator D/Dt stands for the material derivative. Together with Eq. (3), time evolution of the director n is controlled by h(r) = γ1 N + γ2 An.

(4)

The parameters α1 ∼ α6 and γ 1 , γ 2 are the Leslie coefficients with the dimension of viscosity. h(r) is the “molecular field” that mediates the effects of the elastic energy and other external applied fields.27 Although above theories give a general understanding of nematic fluids, based on the concepts such as the Frank elasticity, dissipative rotational motion of the director, and anisotropic viscosities, quantitative verification by experiments can never be a simple task. The Frank elastic theory has been experimentally verified mainly for thermotropic LCs of small molecules. The experimental investigations are rare for lyotropic liquid crystals, which are formed by dissolving rodlike macromolecules in a solvent. Due to the experimental limitation, the Leslie coefficients are difficult to be obtained to solve the flow problems of nematic fluids. As a powerful supplement to the theoretical and experimental approaches, DPD method has been purposely designed to investigate fluid mechanics.1, 6 However, up until now, most DPD studies involving nematics are concentrated on the topics of static equilibrium properties such as phase behavior and selfassembly. Few efforts have been devoted to the study of DPD nematics under the influence of external perturbations.31, 32 To our knowledge, no systematical investigation on the distortion and flow of DPD nematics has been reported in the literature. Such studies are of particular significance for understanding rheological properties of nematic fluids by DPD modeling. In this study, we use DPD simulations to investigate the distortion and flow of DPD nematics. The paper is organized as follows: In Sec. II, we present the DPD simulation methods, the models for nematic fluids, and the simulation techniques. In Sec. III, simulation results are given in three subsections. The static distortion investigated by the simulation is presented in Sec. III A and the results are discussed by comparing with the Frank elastic theory. The results of the rotational motion of n are presented in Sec. III B, which are compared with dynamical experiments and classical theories for nematics. In Sec. III C, simulation results about the flow are given and the anisotropic shear flow is characterized by Miesowicz viscosities obtained from the simulation. A conclusion of the studies is presented in Sec. IV after the discussions.

dvi dr = vi , mi = fi , (5) dt dt where ri , mi , and vi are the position, mass, and velocity of the ith particle; fi is the total force experienced by the ith particle. The total force exerted on the ith particle by all other particles j is given by   R FCij + FD (6) fi = ij + Fij , j =1 R where FCij , FD ij , and Fij represent the conservative, dissipative, and random forces, respectively. These pairwise additive forces are given as

FCij = aij w(rij )ˆrij ,

(7)

D rij · vij )ˆrij , FD ij = −γ w (rij )(ˆ

(8)

FRij = σ w R (rij )θij ( t)−1/2 rˆ ij ,

(9)

where rij = ri − rj , rˆ ij = rij /|rij |, vij = vi − vj , t is the time step, and aij is the repulsion parameter. In this study, the same value was adopted for aij between all particles (aij = a). In Eqs. (7)–(9), γ is the friction coefficient, σ is the noise amplitude, wD (r) and w R (r) are dimensionless weight functions for the dissipative and random processes, and θ ij is a delta-correlated Gaussian random variable (θ ij (t) = 0 and θ ij (t)θ kl (t ) = (δ ik δ jl + δ il δ jk ) δ(t − t )). To be consistent with the fluctuation-dissipation theorem, w D (r), w R (r), γ , and σ should meet the conditions, w D (r) = [wR (r)]2 and σ 2 = 2γ kB T (kB is the Boltzmann constant and T is the temperature).33 Under these conditions, the equilibrium Gibbs-Boltzmann distribution is obtained when the dissipative and noise terms are switched on. As a common practice in DPD, the soft conservative interaction meets the following condition:  1 − r/rc (r ≤ rc ) , (10) w(r) = 0 r > rc where rc is the maximum interaction radius. Typically, the weight function for the random force w R (r) is the same function as for w(r). The particle mass m, maximum interaction radius rc , and kB T are, respectively, taken as the units of mass, length, and energy, i.e., m = rc = kB T = 1. Thus, the unit of time is given by the dimension (mrc 2 /kB T)1/2 = 1. In this study, a velocity-Verlet integration was employed to update the position and velocity of DPD particles, and the time-step for integration was t = 0.04. The repulsive parameter a was set to be 25 and the friction coefficient γ was set to be 4.5, which are common choices used in DPD simulations.33 A free source code LAMMPS was employed for the DPD simulations,34 with some modifications for implementing boundary conditions.

II. SIMULATION METHOD AND MODELS A. DPD method

B. Model of nematic fluid

In DPD simulations, a group of atoms and molecules are coarse-grained into a DPD particle.1 DPD particles obey

In this model, the nematic fluids are composed of rigid rods and normal DPD particles (Fig. 1(a)). Each rigid rod is

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-3

T. Zhao and X. Wang

J. Chem. Phys. 140, 184902 (2014)

where u is the unit vector along the long axis of the rods, I is the unit tensor, and brackets denote the average over rigid rods in the body. If there is no distortion, n represents the preferential orientation of the simulation volume. Otherwise, it only indicates the preferential orientation at the representative points in selected regions. The scalar order parameter S2 , defined as 3/2 of the largest eigenvalue of Q, is used to characterize the degree of orientation of the rods, which equals to 0 for random distribution and 1 for a perfect alignment. C. External fields

In the simulations, external fields were applied through a pair of forces acting on two terminal particles of each rigid rod. The forces were set to have the same magnitude fext but in the opposite direction. The translational motion was not influenced by the forces, but the forces could exert a torque on the rigid rod, which is equivalent to a permanent dipole (along the long axis u) interacting with an external field. The torque (Tu ) can be represented by the following equation: Tu = fext Lu × Hext ,

FIG. 1. (a) Illustration of the DPD liquid model. (b) The setup to study the twist distortion and flow.

formed by several fused DPD particles and the distance between consecutive particles is fixed at 2/3rc . The rods represent the mesogenic units in lyotropic systems and the normal DPD particles are used to represent coarse-grained solvent.23 The AlSunaidi’s framework for mesogens is adopted in this model.19 After each time-step, the total force and torque on each rigid rod are computed as the sum of the forces on its constituent particles; therefore, the particles in a rigid rod move as a single rigid-body. The ratio of the number of particles in rigid-rods to the total number of particles is referred to as the solution concentration (c). Phase behavior of this binary DPD fluids has been studied in our previous work.23 Nematic phase is formed in an appropriate concentration and temperature range. In this study, rigid rods containing 9 DPD particles are used for the simulations typically under proper concentration (c = 0.9) and temperature (T = 1.0), which has been proven to exhibit the nematic phase.23 The rodsolvent binary system at this state is used as a representative case to investigate static and dynamic behavior of the DPD nematic fluids. In some cases, the solution concentration is also adjusted to study the dependence of shear viscosity on concentration. Common orientation axis of the rigid rods in a body of the DPD nematic fluid is represented by the director n calculated as the principal axis of order tensor:   1 (11) Q(r) = uu − I , 3

(12)

where L is length of the rigid rod, Hext is an unit vector representing direction of the applied external field, fext represents magnitude of the field in the dimension of force kB T/rc . The external fields with the strength in the range from 0.01 to 0.05 (fext = 0.01 to 0.05) are used for the simulations. As it can be seen in supplementary material,35 the influence of the weak external fields on the order parameter S2 of the nematics can be ignored. The applied external fields play an important role in this study. In Sec. III A, distortions are induced through the external fields applied on the rods to investigate elastic properties of the DPD nematics. In Sec. III B, the orientation direction of rods is purposely changed by the external fields. Through this process, the dynamic of the rotational distortion is investigated. In Sec. III C, the external forces are used to hold the n orientation so that the anisotropy in viscous properties can be studied in the presence of shear flow. D. Distortion simulation

Twist distortions are employed to study the static responses of the DPD nematics, which is illustrated by Fig. 1(b). A procedure similar to the experiment of Frederiks transition is used to generate the distortions. Initially, the whole system (including rigid rods in anchor zones and simulation volume) forms a homogenous nematic single crystal with n along the z direction. Twist distortions are produced in the simulation volume by applying an external field that tends to orientate rigid rods towards y direction, where the boundary constraint tries to keep the rigid rods in z direction at x = 0 and x = d. The boundary constraint is achieved by the interaction between rods in anchor zones and simulation volume. The anchor zone is composed of two regions at the two opposite sides of the simulation box, which are perpendicular to the x axis. A much stronger specific external field is applied to the rods in anchor zones to fix their orientation direction in the z direction. In the nematic phase, the rigid rods tend to be

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-4

T. Zhao and X. Wang

J. Chem. Phys. 140, 184902 (2014)

parallel with the neighboring ones, orientation of rods at the boundary area between anchor and simulation volume is kept along z direction. In order to lower the potential energy in the applied external field, rods in simulation volume tend to orientate along y direction. At the same time, the rigid rods try parallel to each other to maintain the orientation order and favor the orientation in z direction. The competition occurs between the alignment favored by the anchor zone and that favored by the external field in the simulation volume. A curved spatial profile of n along the x direction is produced by the simulations as the result of final balance, which corresponds to the twist distortion described by the Frank elasticity theory. By examining the spatial profile of the director as function of x, the elastic constants of the DPD nematics are obtained from the simulation results.

Two layers of wall particles outside of the simulation domain are used to reduce the density fluctuations near the boundary.36, 37 The density distortions near the boundary are minimized by adjusting the distance of the wall layers from the boundary and the repulsive parameter between fluid particles and wall particles. We choose α 1 = 0.125rc and α 2 = 0.2rc , and the repulsive parameter between wall particles and fluid particles is set to be aw = 10/25a = 10, where a is the repulsive parameter between fluid particles. Particles in the boundary layer at the side x = d are set to move with velocity vwall in Eq. (15) to produce the correct boundary condition of the velocity.

III. RESULTS AND DISCUSSIONS A. Static distortion and elastic properties

E. Velocity gradient and boundary condition

When studying shear viscosities, a bounce-back boundary condition is used to implement velocity gradient.36 The rectangular Cartesian coordinates are defined by the same way as given in Fig. 1(b). Fixed boundaries are used for the walls with the normal in the x direction, while periodic boundaries are used for the walls with the normals in the y and z directions. If a particle at position r = (x, y, z) travels out of the simulation domain in the x direction (x > d, or x < 0), it will be bounced back into the simulation domain with a new position given as rnew = rold + 2dx nwall .

(13)

Here, nwall is the vector normal to the wall and pointing into the simulation domain, and dx is the distance between the particle and the corresponding boundary. Due to the fact that mesogenic units are modeled by rigid rods, the bounce-back rule for the particles in rigid rods needs to be modified by the following protocol. If the first of the particles in a rigid rod travels out of fluid domain, the whole rod will be bounced back, which means that all particles in that rod will have a new position, rnew,i = rold,i + 2dnwall,1 , i = 1...m,

(14)

m is the number of particles in a rigid rod. In order to get the desired velocity profile at the boundaries, the velocity of the particle bounced back is also changed at the same time, vnew = 2vwall − vold ,

(15)

which is to implement a no-slip boundary condition for both still and moving walls. Couette flow is produced by setting different wall velocity at the two boundaries in x axis. At one side where x = 0, vwall = (0, 0, 0), which makes the equilibrium average velocity at the boundary is zero. At the other side where x = d, the wall velocity is set to be vwall = (0, 0, v) or vwall = (0, v, 0), so that particles in vicinity of the boundary have an average velocity with non-zero z component or y component. The velocity gradient (shear rate) can be calculated as v (16) γ˙ = . d

Experimentally, static distortion of a nematic sample is studied by examining director alignment under competing influences of boundary anchoring and external fields. A typical setup for these investigations was developed by Frederiks and has been widely used for a long time.27 In this study, as described in Sec. II D, we developed a similar computational model to investigate this problem (Fig. 1(b)). In this typical case, the easy axis of the rigid rods is parallel to the z direction at x = 0 and x = d. Initially, the whole system (including rigid rods in anchor zones and simulation volume) forms a homogenous nematic single crystal domain with the director n along the z direction. An external field, mimicking magnetic or electric field, is applied in the y direction normal to the x-z plane. This case corresponds to the twist distortion of the Frederiks setup.27 For a nematic, a single parameter θ can be used to characterize the director alignment, which has the relation n = (0, sinθ , cosθ ). In the simulation process, the director alignment in each thin layer perpendicular to the x direction is examined to investigate director distortion in the simulation box. For a typical case, simulations were performed in a (26rc × 40rc × 40rc ) simulation box, with one (20rc × 40rc × 40rc ) simulation volume and two (3rc × 40rc × 40rc ) anchor zones. The enlarged box dimensions in y and z direction were used to get more rods in the system and obtain better statistics. The simulation volume was divided into 100 slab-like layers (with an area of 40rc × 40rc ) normal to the x direction. As the DPD particle density was taken to be 3rc −3 , the rod concentration was 0.9 when each rod contained 9 particles, there were on average 96 rods in each layer. After reaching equilibrium (when θ did not change obvious with time), data of 20 time-steps were used to obtain each data point of θ , which represented the average orientation of about 1920 rods. For each layer, the director n and value of θ were calculated using the method described in Sec. II B. The typical result obtained from the simulations is given in Fig. 2, which was obtained for the system with concentration c = 0.9 and temperature T = 1.0. It shows the variation of θ along x direction under the influence of external fields of different magnitudes. The x coordinate in the figure is normalized by the length of simulation box in x direction (d = 20 rc ). It can be seen that n(r) of the DPD nematic

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-5

T. Zhao and X. Wang

J. Chem. Phys. 140, 184902 (2014)

FIG. 2. Simulation results of the twist distortions together with fitted curves obtained from the Frank elastic continuum theory with K2 as the fitting parameter.

varies continuously along the x direction. It means that n(r) in the bulk tends to take the preferential direction of the external field. On the other hand, n(r) near both surfaces takes the easy direction imposed by the strong anchoring. The general characteristic revealed by the simulation is consistent with the experimental result obtained from the Frederiks setup. It indicates that the DPD nematics show long-range elasticity and can transmit torques, which can be described by the Frank elastic theory. According to the elastic continuum theory, a certain state of a nematic system is controlled by the distortion free energy, which can be given by an equation of n with three elastic constants (Eq. (1)). The equilibrium state is reached when the distortion free energy is minimized. For the twist distortion, the director n(r) rotates within y-z plan and can be represented by θ through n(r) = (0, sinθ , cosθ ). By minimizing the distortion free energy, the equilibrium equation can be obtained:27 K2

∂ 2θ + ρfext L cos θ = 0, ∂x 2

orientate the director of DPD nematics, but are far below the value to influence the orientation of an individual rod in the nematics and isotropic phase. As shown in Fig. 2, the curves obtained from Eq. (17) can well fit the simulation results of the twist distortions using K2 as the fitting parameters, where the open symbols are simulation results and solid lines are fitted curves. The open symbols are very close to the fitted curves which are the solutions of Eq. (17). It means that the twist distortion behavior of DPD nematics is similar to those described by the Frank elasticity theory. By examining the fitting results, we can obtain the value of the twist elastic constants K2 . As the number density of DPD particles is 3 rc −3 , the number of particles in each rigid rod is 9, and the concentration of the rigid rods is 0.9, the number density of rods can be calculated to be ρ = 3/9*0.9 = 0.3 rc −3 . Taking length of the rigid rod to be L ≈ (91)*2/3+1 = 6.33, the fitting results of K2 for different external fields are shown in Fig. 3. The value of K2 was obtained from a group of parallel simulations and the error bars represent the standard variation. By comparing the result given in Fig. 3 with the Frank elastic theory, one important difference can be noticed. The K2 obtained from the simulation is no longer a constant but dependent on the strength of the applied external fields. This point can be understood by considering the following factors. Equation (1) is deduced under the condition that the spatial gradient of the director field is small in the length scale of rigid rods (L∇n 1). In our simulations, the distortion of n occurs in a length scale (d = 20rc ) and is actually comparable with the length of the rigid rods (L = 6.33 rc ). As the external field increases, the variation of n becomes larger and the assumption of slow variation is no longer suitable for the DPD nematic. Consequently, the apparent value of K2 obtained by fitting Eq. (17) will become a function of the magnitude of external fields. As shown in Fig. 3, K2 decreases with the fext increase following the linear relationship. The value of K2 in the limit of weak distortion can be obtained by extrapolating the K2 ∼ fext curve to fext = 0, which is 3.38 kb T/rc in this case. This rule revealed by the DPD simulation might be able

(17)

where the first term is the torque caused by elastic distortion and K2 is the twist elastic constant, the second term is the toque applied by external field. In a nematic phase with the high orientational order, the long axis u of most rigid rods orientates within a small solid angle around n. Therefore, the torques on n can be approximately estimated by summation of torques on individual rods in a unit volume: Tn = ρfext Ln × Hext = ρfext L cos θ,

(18)

where ρ is the number density of the rigid rods. The variation profile of director n along the x direction θ (x) can be obtained by solving the equilibrium equation. The results obtained from the simulation for the twist distortions are compared with the curves calculated by using Eq. (17) to obtain K2 . Weak external fields (fext ranges from 0.01 to 0.05 kB T/rc ) are used in this simulation to meet the conditions of the Frederiks setup. As it is shown in Sec. II C and will be further discussed in Sec. III B, the external fields can only

FIG. 3. The fitting results of K2 for different external fields, the number density, and length of the rods ρ = 0.3rc −3 and L = 6.33.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-6

T. Zhao and X. Wang

J. Chem. Phys. 140, 184902 (2014)

to provide a way to understand the large scale distortion in the space. B. Dynamics of orientational motions

An obvious advantage of the DPD simulation is that it can provide a feasible way to study the dynamical properties of the orientational motions as discussed below. If the rigid rods are initially orientated along the z direction and a field in the y direction is applied, the rods will gradually rotate from the z direction to the y direction. The director n, representing the average orientation direction of the rods, will consequently rotate from the z direction to the y direction. Orientation of n can be represented as a function of time t, n = (0, sinθ (t), cosθ (t)). During the orientational motion, n will be driven by the torque caused by external fields. For DPD nematics, the rigid rods tend to rotate collectively due to the nematic potential. Therefore, the torques on n can be approximately estimated by summation of torques on individual rods in a unit volume, Tn = ρfext Lcos θ . On the other hand, the orientational motion of n is resisted by the background fluid. According to the Leslie-Ericksen theory (Eq. (4)), this process corresponds to the energy dissipation by rotation of n. This resistance is proportional to the rate of the change of the director with respect to the background fluid and the proportional coefficient is the rotational friction coefficient γ 1 . Ignoring inertia, the energy input by external field and the energy dissipation by friction should be equal to each other. Therefore, time evolution of n in this process is represented by the equation ∂θ , (19) ∂t where the left side is the torque of the external field and the right side is the friction torque. With initial condition θ (0) = 0, letting  π θ γ1 y = tan − ,τ = , (20) 4 2 ρfext L ρfext L cos θ = γ1

Eq. (19) can be solved analytically as  t y = exp − . τ

FIG. 4. The simulation results of y varying with time and y ∼ t fitting curves for the different external fields, where open symbols are simulation results and solid lines represent the numerical solutions of Eq. (21).

friction coefficient γ 1 can be calculated. Taking the number density of rods ρ = 0.3rc −3 and L = 6.33, the rotational friction coefficient of the DPD nematic fluid is γ 1 = 203 (mkB T)1/2 /rc 2 . In this section as well as in Sec. III A, the external fields for the simulations are relatively weak. The largest potential energy due to external fields can be estimated as fext L, which is about 0.3 for fext = 0.05. Compared with the kinetic energy of each rigid rod (mkb T = 9), the external fields are not strong enough to completely suppress the thermal orientational motions of individual rigid rods. As shown in Sec. II C, for the DPD nematics, the influence of the weak external fields on the order parameter S2 can be ignored. This point can further be proved by applying the same external fields on the DPD fluids in an isotropic phase (with high temperatures or low concentrations). The simulation result shows that no orientational motion of the rods can be induced by the external fields. It indicates that the DPD nematics possess the unique characteristic of liquid crystals, for which the orientation is

(21)

On the basis of this understanding, the orientational motions of nematics were investigated by the DPD simulation. The DPD nematics used for the simulation were composed of rigid rods with 9 consecutive particles and normal DPD particles. Concentration of the system was set to be c = 0.9 and temperature T = 1.0, which has been proved to possess the nematic phase.23 From simulation results, we first calculate n for all the rods in the simulation box and calculate quantity y according to Eq. (20). Results under different strengths of the external fields (fext ranges from 0.01 to 0.05) are presented in Fig. 4. It can be seen that the time evolution of y can be well represented by Eq. (21). By fitting the simulation data with Eq. (21), we get a series of relaxation time τ , which is reciprocally proportional to the strength of external fields. The relationship between 1/τ and fext is shown in Fig. 5, which shows the linear 1/τ ∼ fext relationship passing the origin point. From the slope, the value of rotational

FIG. 5. The relationship between 1/τ and fext , which shows a linear correlation passing the origin point.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-7

T. Zhao and X. Wang

sensitive to weak external perturbations.27 As a result of the formation of the nematic phase, rigid rods tend to rotate collectively. The collective rotational motion of the rigid rods shown by the DPD simulation can well represent the orientational response of nematics to external fields. According to the elastic continuum theory, the static twist distortion in nematic can be described by Eq. (17). If considering friction and energy dissipation, the time-dependent rotational motions of the director should be described by the following equation: ∂ 2θ ∂θ (22) + ρfext L cos θ = γ1 , 2 ∂x ∂t where the friction torque is added in the right side of the equation.27 The friction coefficient γ 1 in above equation has been treated as the rotational viscosity in the field theory. In order to compare the DPD results with the field theory, we use the same symbol γ 1 without further distinction. By taking K2 values obtained at corresponding fext shown in Fig. 3 and γ 1 obtained from the DPD simulation reported in this section, Eq. (22) can be solved numerically. As indicated in Sec. III A, K2 obtained from the DPD simulation possesses different values at the corresponding fext . In Fig. 6, we show the variation of θ at x = d/2 as function of time. Open symbols are simulation results obtained by averaging the principal orientation direction of rods in vicinity of x = d/2 and solid lines represent the numerical solutions of Eq. (22). Simulation results and theoretical solutions show a satisfactory agreement. As the values of γ 1 and K2 are obtained independently from the simulations, it further validates the DPD model used to simulate both the static distortions and the orientational motions of nematics. K2

C. Flow of DPD nematics

Generally, the flow regime of nematics is complicated due to the coupling of the translational motions to the orientational motions.27 In Sec. III B, we present the results obtained from the investigation on the orientational motions by

FIG. 6. The variation of θ at x = d/2 as a function of time, open symbols are simulation results obtained by averaging the principal orientation direction of rods in vicinity of x = d/2 and solid lines represent the numerical solutions of Eq. (22).

J. Chem. Phys. 140, 184902 (2014)

DPD method, where the motions are driven by applying external field perpendicular to the initial orientational direction. In this section, we will concentrate on flow behavior of the DPD nematics for the idealized cases used in the Miesowicz experiment, where the director n is kept in some fixed orientations. Miesowicz experiment is a typical way to study the flow behavior of nematics under well-controlled conditions.27 Three Miesowicz viscosities are used to characterize the viscous response under three special conditions with a fixed orientation of n. Miesowicz viscosities ηa , ηb , and ηc are defined as the apparent shear viscosities when the director is steadily aligned in the vorticity, velocity, and velocity gradient direction, respectively.27 The Miesowicz viscosities have been experimentally measured to characterize the anisotropic viscous properties of the nematic fluids. In this simulation, similar conditions were adopted and techniques implementing velocity gradient and external fields are described in Sec II. To obtain the shear viscosities of different conditions, stress tensor was calculated by Eq. (23) and then averaged over time after equilibrium state is reached:

 1  mk vki vkj + rki fkj . (23) σij = − V k k In above equation, V is the volume of the simulation box, mk is the mass of the kth particle; vki , rki , and fki are the ith component of velocity, position, and total force applied on the kth particle, respectively. As physical quantities in the simulation are all represented in reduced units, stress is given in the dimension of energy/length3 (kB T/rc 3 ). Shear rate is in dimension of time−1 ((kB T/ mrc 2 ) 1/2 ). Viscosity equals to stress divided by shear rate so that it is in the dimension of (mkB T)1/2 /rc 2 . Simulations were performed in a (30rc × 30rc × 30rc ) cubic simulation box. The velocity gradients were generated in the x direction, using boundary conditions described in Sec. II D. Periodic boundary conditions were used for the walls with the normal in the y and z directions. External fields were employed to hold n unchanged in the shear field to obtain the Miesowicz viscosities. By the definition, if the director n is fixed by external fields in the z direction and velocity direction is in the y direction, apparent viscosity calculated from the ratio between σ yx and γ˙ will be the Miesowicz viscosity ηa . Accordingly, ηb is obtained when n is kept in the z direction and the fluids flow is along the z direction, and ηc is obtained by holding n in the x direction and letting the fluids flow along the y direction. The DPD nematic fluids used in this section are composed of the rigid rods with 9 consecutive particles and normal DPD particles. Concentration of the system is c = 0.9 and temperature takes T = 1.0, which can keep the system in the nematic phase.23 Figure 7 shows the σ − γ˙ relation for obtaining ηa under the Miesowicz condition. The external field is used to hold the director from being rotated by the shear flow, which is set to be fext = 0.1 as a typical condition. The result shows that the shear stress σ grows linearly with the shear rate γ˙ . From the slope of the fitted line, the Miesowicz viscosity ηa is obtained. In Fig. 7, we also give the result when

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-8

T. Zhao and X. Wang

J. Chem. Phys. 140, 184902 (2014)

FIG. 7. The σ − γ˙ relation obtained from the simulation result for calculating Miesowicz viscosity ηa , c = 0.9, and T = 1.0.

FIG. 9. The σ − γ˙ relation obtained from the simulation result for calculating Miesowicz viscosity ηb , c = 0.9, and T = 1.0.

a stronger external field (fext = 0.5) is applied for comparison. It can be seen that the strength of the external fields varying in the range has a trivial influence on the σ − γ˙ relation, which indicates that the fields used here are strong enough to prevent the director n from being rotationally disturbed, but far from the criteria to influence structure of the fluids. Therefore, the value of ηa obtained by this method can be considered as an intrinsic property of the DPD nematic fluid, which is independent of the external field strength to obtain the result. Similarly, ηc is obtained when n is kept in the velocity gradient direction. The σ − γ˙ relations obtained from the simulation for the calculation of ηc is shown in Fig. 8. In a shear flow, a rod with the long axis along the velocity gradient direction will feel the largest torque compared with other orientations. In the simulation, to avoid the rods from rotation, a stronger external field (fext = 1.0) is needed to hold the director n. In the tested shear rate range, the shear stress σ grows linearly with the shear rate γ˙ . In Fig. 8, we also present the simulation results when weaker external fields are used. It can be seen that the linearity of the σ − γ˙ relations fails at the high shear rate. The reason for the failure is that the external fields are no longer strong enough to hold the director

n unchanged in the velocity gradient direction. As shown in Fig. 8, the value of the shear rate where linearity starts to fail is related to the strength of the external field, but the slope of the linear part almost remains the same. Therefore, ηc can still be estimated from the slope of the linear parts of σ − γ˙ curves. Different with above two cases, ηb can be obtained without applying an external field. The reason is that the director n is initially in the velocity direction and it is not changed within the time period when the shear stresses are calculated. The σ − γ˙ relation obtained under this condition is shown in Fig. 9. It can be seen that the shear stress also grows linearly with shear rate, and the value of ηb is obtained by fitting the line. From above simulations, we obtain three Miesowicz viscosities for a representative DPD nematic fluid and the results are summarized in Table I. The same as shown by the continuum model for nematic liquid crystals, the DPD nematic fluids own anisotropic viscous properties depending on the direction of n relative to the flow field. The results show that the Miesowicz viscosities of the DPD nematic fluid have the order ηc > ηa > ηb , which is a typical result observed by the experiments for nematic fluids. However, considering the relative magnitude of the three viscosities, ηc is much larger than ηa and ηb , which is different from available experimental results obtained for low-molecular-weight nematics.27, 38 Although experimental data about polymeric nematics are rare, a relationship of ηc ηa > ηb has been obtained from Doi’s molecular theory.39 It indicates that the DPD nematic fluid is a model more suitable for lyotropic polymeric nematic fluid. It seems to be a natural consequence as the DPD fluid used TABLE I. Fitting results obtained from the stress-shear rate curves shown Figs. 7–9. Miesowicz viscosities are calculated from the slopes of the lines.

FIG. 8. The σ − γ˙ relation obtained from the simulation result for calculating Miesowicz viscosity ηc , c = 0.9, and T = 1.0.

ηa ηb ηc

Slope

R2

3.77 ± 0.06 1.68 ± 0.05 27.4 ± 0.3

0.997 0.992 0.999

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-9

T. Zhao and X. Wang

FIG. 10. The relationship between the shear viscosity and the solution concentration at the shear rate γ˙ = 3.3 × 10−3 .

here is a binary mixture of the rigid rods and normal DPD particles. One remarkable characteristic of lyotropic liquid crystalline polymer systems is the non-monotonic dependence of shear viscosity on polymer concentration. As polymer concentration increases, the system turns from the isotropic to anisotropic phase and the maximum viscosity appears around the critical concentration of the isotropic-anisotropic phase transition.39 According to the phase diagram that we obtained previously,23 at temperature T = 1.0, as concentration c goes from 0.1 to 1.0, the system changes from the isotropic phase, to the coexistence of isotropic phase and nematic phase, and finally to the homogenous nematic phase. In Fig. 10, we present the relationship between the shear viscosity and concentration at the shear rate γ˙ = 3.3 × 10−3 . It can be seen that the viscosity of the system first increases and then falls with the concentration increase. The maximum viscosity value appears at the critical concentration, at which the nematic phase emerges. This simulation result well reflects the general tendency already revealed by the experiments.40, 41 In this section, we present the viscous response of the DPD nematic fluids to the shear flow fields. As the same DPD nematic fluids are used for the investigations on rotational viscosity and Miesowicz viscosities, a comparison between the relative magnitude of rotational viscosity and Miesowicz viscosities can be made. The result indicates that the rotational viscosity γ 1 is about two orders larger than the Miesowicz viscosity ηb . This unique feature is consistent with the results obtained for polymeric nematic fluids from the experiments42 and theoretical investigation.43 On the other hand, this magnitude difference between γ 1 and ηb cannot be seen for lowmolecular-weight nematics.27 It also indicates the DPD model used in this study is more suitable for lyotropic polymer systems. IV. CONCLUSIONS

The distortion and flow of nematics are investigated by the DPD simulation in this study. The DPD models are composed of rigid rods and normal DPD particles, which are in the

J. Chem. Phys. 140, 184902 (2014)

nematic phase in the studied concentration and temperature range. The simulations show that the DPD nematics are able to reveal some fundamental distortion and flow characteristics of nematics. Parameters used to characterize the deformation and flow, such as Frank elastic constants, rotational viscosity, and Miesowicz viscosities, are obtained from the simulations. Spatial distortion profile of the director n demonstrated by the simulation is consistent with the prediction of the Frank elastic theory. The elastic constant K2 is obtained by fitting the spatial profile of the director field when the elastic distortions are generated. Orientational motions under the influence of external fields show the characteristics of the friction rotations of n, which are in a good agreement with the dynamical experiments and the classical theories for nematics. For the shear flow, the Miesowicz viscosities obtained from the simulation show a relationship of ηc ηa > ηb . The rotational viscosity γ 1 is about two orders larger than the shear viscosities ηb . These results indicate that the DPD nematic fluids are more suitable to simulate polymeric nematics rather than thermotropic low-molecular-weight nematics. Using these constants and parameters, the relative strength of the viscous forces and elastic forces in a certain flow filed can be estimated by Ericksen number, which is defined by Er = ηVh/K, where η is a characteristic viscosity, K is a Frank elastic coefficient, V is a flow velocity, and h is a length scale of the flow geometry.44 The results indicate that the DPD method, with the links between the particles to form the rigid rods, can be used as a valid model to study distortion and flow of nematics. In this paper, we verify the validity of this method as a computational approach to study nematic fluids. On the other hand, the distortion and flow behavior of DPD nematics is complicated and dependent on other factors such as rod length, concentration, temperature, as well as the DPD parameters. To uncover possible influences of these factors on elastic and viscous properties of the DPD nematics, further comprehensive investigations will be needed. ACKNOWLEDGMENTS

The financial supports from the National Basic Research Program of China (NBRPC) (973 Program) under Project No. 2011CB606102 and the National Science Foundation of China (NSFC) under Project No. 51233002 are gratefully acknowledged. 1 P.

Hoogerbrugge and J. Koelman, Europhys. Lett. 19(3), 155 (1992). Koelman and P. Hoogerbrugge, Europhys. Lett. 21(3), 363 (1993). 3 P. Espanol and P. Warren, Europhys. Lett. 30(4), 191 (1995). 4 C. A. Marsh, G. Backx, and M. H. Ernst, Phys. Rev. E 56(2), 1676–1691 (1997). 5 E. G. Flekkøy, P. V. Coveney, and G. De Fabritiis, Phys. Rev. E 62(2), 2140 (2000). 6 P. Espanol, Phys. Rev. E 52, 1734–1742 (1995). 7 E. Boek, P. Coveney, and H. Lekkerkerker, J. Phys.: Condens. Matter 8(47), 9509 (1996). 8 E. Boek, P. Coveney, H. Lekkerkerker, and P. van der Schoot, Phys. Rev. E 55(3), 3124 (1997). 9 J. M. Kim and R. J. Phillips, Chem. Eng. Sci. 59(20), 4155–4168 (2004). 10 Y. Kong, C. W. Manke, W. G. Madden, and A. G. Schlijper, J. Chem. Phys. 107(2), 592–602 (1997). 11 N. Spenley, Europhys. Lett. 49(4), 534 (2000). 2 J.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

184902-10 12 W.

T. Zhao and X. Wang

Jiang, J. Huang, Y. Wang, and M. Laradji, J. Chem. Phys. 126(4), 044901 (2007). 13 V. Symeonidis, G. E. Karniadakis, and B. Caswell, Phys. Rev. Lett. 95(7), 076001 (2005). 14 E. Sultan, J.-W. van de Meent, E. Somfai, A. N. Morozov, and W. Van Saarloos, Europhys. Lett. 90(6), 64002 (2010). 15 R. D. Groot and T. J. Madden, J. Chem. Phys. 108, 8713 (1998). 16 R. D. Groot, T. J. Madden, and D. J. Tildesley, J. Chem. Phys. 110, 9739 (1999). 17 X. Guerrault, B. Rousseau, and J. Farago, J. Chem. Phys. 121, 6538 (2004). 18 P. Nikunen, I. Vattulainen, and M. Karttunen, Phys. Rev. E 75(3), 036713 (2007). 19 A. AlSunaidi, W. Den Otter, and J. Clarke, Philos. Trans. R. Soc. London, Ser. A 362(1821), 1773 (2004). 20 Y. K. Levine, A. E. Gomes, A. F. Martins, and A. Polimeno, J. Chem. Phys. 122, 144902 (2005). 21 A. AlSunaidi, W. Den Otter, and J. Clarke, J. Chem. Phys. 130, 124910 (2009). 22 Z. Zhang and H. Guo, J. Chem. Phys. 133, 144911 (2010). 23 T. Zhao and X. Wang, J. Chem. Phys. 135(24), 244901 (2011). 24 T. Zhao and X. Wang, J. Chem. Phys. 138(2), 024910 (2013). 25 T. Zhao and X. Wang, Polymer 54(19), 5241–5249 (2013). 26 T. Zhao and X. Wang, J. Chem. Phys. 139(10), 104902 (2013). 27 P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Oxford University Press, 1995). 28 J. Ericksen, J. Rheol. 5, 23 (1961).

J. Chem. Phys. 140, 184902 (2014) 29 J.

Ericksen, Arch. Ration. Mech. Anal. 9(1), 371–378 (1962). Leslie, Arch. Ration. Mech. Anal. 28(4), 265–283 (1968). 31 A. E. Gomes, A. F. Martins, and A. Polimeno, Mol. Cryst. Liq. Cryst. 435(1), 135–152 (2005). 32 Y. Levine and A. Polimeno, Euro. Phys. J. E 23(1), 13–23 (2007). 33 R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997). 34 S. Plimpton, J. Comput. Phys. 117(1), 1–19 (1995). 35 See supplementary material at http://dx.doi.org/10.1063/1.4873699 for the order parameter of nematics. 36 D. Duong-Hong, N. Phan-Thien, and X. Fan, Comput. Mech. 35(1), 24–29 (2004). 37 T. Zhao, X. Wang, L. Jiang, and R. G. Larson, J. Chem. Phys. 139, 084109 (2013). 38 A. M. Jamieson, D. Gu, F. L. Chen, and S. Smith, Prog. Polym. Sci. 21(5), 981–1033 (1996). 39 M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, Oxford, 1986). 40 K. F. Wissbrun, J. Rheol. 25, 619 (1981). 41 M. G. Northolt and D. J. Sikkema, in Separation Techniques Thermodynamics Liquid Crystal Polymers (Springer, Berlin, 1991), Vol. 98, pp. 115–177. 42 V. G. Taratuta, A. J. Hurd, and R. B. Meyer, Phys. Rev. Lett. 55(2), 246– 249 (1985). 43 S. D. Lee, J. Chem. Phys. 88(8), 5196–5201 (1988). 44 R. G. Larson, The Structure and Rheology of Complex Fluids (Oxford University Press, New York, 1999). 30 F.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.28.191.34 On: Mon, 22 Dec 2014 22:50:15

Distortion and flow of nematics simulated by dissipative particle dynamics.

In this study, we simulated distortion and flow of nematics by dissipative particle dynamics (DPD). The nematics were modeled by a binary mixture that...
1MB Sizes 2 Downloads 3 Views