- Email: [email protected]
clay nanocomposites under flow
Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
Contents lists available at SciVerse ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: http://www.elsevier.com/locate/jnnfm
Computational model for Brownian dynamics simulation of polymer/clay nanocomposites under flow Takehiro Yamamoto ⇑, Nobuhiro Kanda Department of Mechanical Engineering, Graduate School of Engineering, Osaka University, 2-1, Yamadaoka, Suita, Osaka 565-0871, Japan
a r t i c l e
i n f o
Article history: Received 17 January 2012 Received in revised form 6 June 2012 Accepted 7 June 2012 Available online 28 June 2012 Keywords: Polymer/clay nanocomposites Brownian dynamics simulation Gay–Berne potential Reversible network model Shear flow
a b s t r a c t A Brownian dynamics (BD) simulation model is proposed for simulating the behavior of polymer/clay nanocomposites under flows. In this model, polymer/clay nanocomposites are described by dispersion systems of disk-like particles in polymer melt. The particles are represented by oblate spheroids that interact with each other via the Gay–Berne potential, and the polymer melts are modeled using a reversible network model, in which polymer networks are represented by a set of FENE dumbbells having two states, active and dangling. The particle–bead and the bead–bead interactions are also considered. In the present study, BD simulations for the shear flow of this system was conducted to investigate the orientation behavior of disk-like particles and the effect of polymers on the particle orientation. The numerical results indicate that the disk-like particles are well aligned in flows and that the active/dangling transition of dumbbells affects the orientation behavior of the particles. Furthermore, it is found that the degree of particle orientation is higher when the intensity of entanglement of polymers U 0 is stronger because polymers with larger U 0 tend to be more aligned in the flow direction and affect the particle motion behavior. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Polymer/clay nanocomposites are remarkable materials because of their superior mechanical, thermal, and optical properties. In polymer/clay nanocomposites, nanoscale particles, typically 10–100 nm in size, are dispersed in a polymer matrix as fillers. Because of their high surface-to-volume ratio, nanoscale particles effectively modify the material properties of polymers as compared to usual polymer composites. The effects of nanoscale fillers on material properties, their orientation behavior, the flowinduced structure of the system, and rheological properties have been actively studied, e.g. [1–14]. Okamoto et al. [4] have found a house of cards structure in polypropylene/clay nanocomposites under uniaxial elongational flow by TEM observation, and reported that both strong strain-thickening and rheopexy are originated from the perpendicular alignment of the silicate layers to the stretching direction. Fornes et al. [6] have investigated the dependence of both mechanical properties and rheology on the molecular weight of matrix polymer using organoclay nanocomposites based on three different molecular weight grades of nylon 6. Lele et al. [8] have carried out measurements of both in situ X-ray diffraction and rheology for capillary flow of a system of organically modified hydrophobic montmorillonite clay and syndiotactic ⇑ Corresponding author. Fax: +81 6 6879 7308. E-mail address: [email protected] (T. Yamamoto). 0377-0257/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnnfm.2012.06.005
polypropylene. They investigated the relation between the flow-induced structure and rheology, and reported that yield behavior that might be linked to the clay orientation occurs at high shear rate and that the clay particles are easily oriented by shear to achieve a high degree of orientation after the yielding transition. Lin-Gibson et al. [11] have investigated the shear-induced structure of systems of polyethylene oxide (PEO) and synthetic hectorite clay Laponite using X-ray scattering, light scattering, optical microscopy, and rheometry. This system forms physical gel via the bridging of neighboring clay particles by polymers absorbed on the clay surface. They reported that this system showed the clay particle orientation with their surface normal parallel the vorticity direction in shear flows. Recently, Dykes et al. [14] have studied the shear-induced orientation for systems of montmorillonite or fluorohectrite clays and PDMS and PDPS copolymer matrix using in situ X-ray scattering and rheometry. They introduced an anisotropy factor, which represents the orientation degree of particles and is evaluated from data of small angle X-ray scattering measurements. In a moderately concentrated system, shear rate dependent anisotropy factor and orientation angle, and partial relaxation upon flow cessation were observed. Conversely, a more highly concentrated system showed less rate dependence in orientation state, and no orientation relaxation. Because the properties of nanocomposites depend on both the orientation and distribution of the particles, the analysis of the behavior of particles in flows is an important issue, and numerous
2
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
experimental investigations have been done as mentioned above. On the other hand, numerical simulation is an alternative effective method for the analysis of particle behaviors. Monte-Carlo and self-consistent field simulations have been used to numerically investigate polymer nanocomposites to determine the effects of nanoparticle–polymer interactions in diblock copolymer/nanoparticle mixtures [15–19]. Molecular dynamics simulations have also been used to analyze the effects of nanoparticle–polymer interactions [20–22]. Singh and Rey [23–25] have proposed constitutive equations for discotic nematic liquid crystals and analyzed the flow-induced structure of discotic mesophases using the constitutive equations. Bertevas et al. [26] have proposed an algorithm of simulations for suspension of oblate spheroidal particles in a Newtonian fluid. In their algorithm, the motion of each particle is described by Jeffery’s solution modified by the interactions between the particles. They investigated the relation between rheological properties and the orientation behavior of particles. Yamamoto and coworkers [27,28] have investigated both particle orientation behavior and flow-induced structure in suspensions of disk-like particles using Brownian dynamics (BD) simulation. Yamamoto and Kasama [28] applied a dispersion system of oblate spheroids and elastic dumbbells in Newtonian solvent to a model of polymer/clay nanocomposites. In the present study, we modified this model for the application to a dispersion system of clays in polymer melts. Polymer melts are modeled using a reversible network (RN) model, which was originally developed for associative polymers [29]. In a RN model, polymer network is represented by a set of numerous non-interacting elastic dumbbells, and hence the deformation of network topology is not directly computed. Consequently, this coarse graining representation for polymer network can save much computational cost as compared to a numerical simulation in which the deformation network structure is simulated. We propose this system as a model for dispersion systems of clay particles in polymer melts, and numerically simulated the shear flow of this system to investigate the flow-induced structure and the orientation behavior of both disk-like particles and polymers.
node bead association Pa
We first describe the computational modeling of polymer/clay nanocomposites employed in the present study. As a model of polymer/clay nanocomposites, we consider a dispersion system of disk-like particles in a polymer melt. In the present model, the disk-like particles are approximately represented by oblate spheroids, and the polymer melt is described by a RN model. We modeled polymer melts by modifying the RN model proposed by Hernández Cifre et al. [29], which has been developed for describing the behavior of associative polymers. In this model a polymer network is expressed by a set of active and dangling segments. The former is a chain that has both ends attached to the network, and the latter is a chain that is attached to the network only one end. Active and dangling segments are modeled by two separate sets of non-interacting elastic dumbbells. A dumbbell consists of two beads linked by a finitely extensible non-linear elongational (FENE) spring [30]. In a FENE chain, the spring force between beads i and j; ðF S Þij , is expressed by
ðF Þij ¼
1 ðQ ij =Q 0 Þ2
;
free bead
Fig. 1. Schematic diagram of the association/dissociation process in a reversible network model.
Both ends of an active dumbbell belong to a polymer network and hence the beads of active dumbbell are denoted by ‘node’ beads. A dangling dumbbell has a node bead and a ‘free’ bead, which represents the end of dumbbell that is not linked to a network (Fig. 1). The resistance coefficient of a bead is different between node and free beads. In the present simulation, the resistance coefficient ratio Z defined by Z=fn =ff is larger than unity. Here fn and ff denote the resistance coefficient for node beads and that for free beads, respectively. The resistance coefficient ff is calculated by ff =3pgs rB , where gs is the Newtonian contribution to polymer viscosity and r B is the bead diameter. The transition from an active dumbbell to a dangling dumbbell and the inverse transition are expressed using probability functions Pd and Pa , respectively (Fig. 1). The RN model has been originally developed for associative polymers. However, Sunarso and Yamamoto [31] adopted this model to non-associative polymer melts by modifying the association probability functions. We used this model to adopt the RN model to the numerical simulation for polymer melts. As for the dissociation process, we used a probability function Pd proposed by Hernández Cifre et al. [29] expressed by Eqs. (2)– (4). This function describes a probability that the transition from an active dumbbell to a dangling dumbbell occurs within time interval Dt:
2 Dt ; Pd ¼ 1 exp
ð2Þ
s
2
2.1. Models for polymer melts and clay particles
HQ ij
dangling dumbbell
3
2
H2 d Q 2 7 s ¼ s0 exp 6 4 n o2 5; 4U 0 kB T 1 ðQ =Q 0 Þ2
2. Computational model
S
Pd dissociation
active dumbbell
ð1Þ
where H is the spring constant, Q ij is the connector vector between the beads, Q ij is the length of vector Q ij , and Q 0 is the maximum length of the spring.
s0 ¼ sf exp
ð3Þ
U0 : kB T
ð4Þ
Here s is the lifetime of an associated bead, s0 is the lifetime at U=U 0 , and sf is a fundamental molecular time scale related to the thermal vibration. Tanaka and Edwards [32] estimated that sf is on the order of 1010 s. U 0 =ðkB TÞ indicates the strength of entanglement. d is the depth of entanglement potential well. Considering the numerical results of BD simulations by van den Brule and Hoogerbrugge [33], Hernández Cifre et al. [29] proposed an association probability function Pa , which indicates a probability that the transition from a dangling dumbbell to an active dumbbell occurs in Dt. The probability is an increasing function of the dumbbell length and is expressed by
" ( Pa ¼ 1 exp a0 þ a1
Q
)
1 ðQ =Q 0 Þ2
#
Dt :
ð5Þ
On the other hand, Sunarso and Yamamoto [31] proposed an association probability function for polymer melts:
a2 Q Pa ¼ 1 exp a0 þ a1 Dt ; Q0
ð6Þ
where a0 ; a1 , and a2 are model parameters. This model was proposed to avoid the shear-thickening behavior of a system by
3
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
restricting an abrupt increase in association probability. This model can predict both the strain-thickening at low elongation rates and the strain-thinning at high elongation rates, which are usually observed for polymer melts. We adopted this model to the present numerical simulation. In the present numerical simulations, the transition probability for each dumbbell is calculated and whether the transition occurs in a considering dumbbell is stochastically determined using a computationally created random number. For the dangling transition, it is also determined at random which node bead in a dumbbell changes to a free node. The effect of excluded volume is usually ignored in BD simulations of the RN model, whereas it is considered in our simulation because the excluded volume affects the interaction between polymers and clay particles. The excluded-volume effect is introduced using the repulsive Lennard–Jones (RLJ) potential as explained later. In the present simulation, disk-like particles are approximately represented by oblate spheroids. In the present system, we need to consider the interactions among the particles and polymer chains such as the particle–particle, particle–bead, and bead–bead interactions. These interactions are represented by potential functions. The Gay–Berne (GB) potential [34] is used to represent the particle–particle interaction. The potential function between particles m and n; /GB mn , is expressed by
12 6 r0 r0 ^ ^ ^ ^ ^ ^ ; /GB mn ðum ; un ; r mn Þ ¼ 4eðum ; un ; r mn Þ R R
ð7Þ
^m; u ^ n ; ^rmn Þ þ r0 ; R ¼ r mn rðu
ð8Þ
^ m is a unit vector indicating the orientation of the rotation where u axis of a spheroid particle m; ^ r mn is a unit vector representing the direction from the centroid of particle n to that of particle m, and r mn is the distance between particles m and n. r0 is the major axis of an oblate spheroid particle. In the present paper, we will leave out the explanation of the entire equations for the GB model because they are rather numerous and are available in the literatures, e.g. [34–38]. Those references also give details regarding the strength of the potential e and the inter-particle separation r. We employed the same notation used in previous studies [37,38] for the parameters of the GB potential: e0 ; er ; l, and m. As for oblate particles, the face-to-face configuration, in which a disk surface of a particle faces to that of another particle, is the most stable. A modified Gay–Berne (MGB) potential [39] is adopted for the potential between a disk-like particle m and a bead n. The equations of the MGB potential function /MGB mn are similar in form to those for the GB potential:
(
^ ^ ^ ^ /MGB mn ðum ; r mn Þ ¼ 4eðum ; r mn Þ
rGBB 0
n
12
R ol
rGBB 0 R
6 ) ;
ð9Þ
^ m ; ^r mn Þ2 ; eðu^ m ; ^rmn Þ ¼ eGBB 1 v0 ðu 0
ð10Þ
^ m ; ^rmn Þ þ rGBB R ¼ r mn rGBB ðu ; 0
ð11Þ
u^ m
particle
u^ n r^mn
(n)
(m)
u^ m
dumbbell
particle
r^mn
(m)
(n)
Fig. 2. Schematic diagram of vectors ^rmn used in the Gay–Berne and the modified Gay–Berne potentials.
Xe in their papers. The model parameter eGBB used in the present r paper represents the ratio of the potential energy minima of the end and side configurations (Fig. 2b in Ref. [28]) and corresponds to ee /eS in Ref. [39]. Typical potential curves for these functions are available in Ref. [28]. We set the values of the model parameters to make the potential well depth smaller than that of the inter-particle case. The RLJ potential /RLJ ðrÞ is used to represent the excluded-volume effect among FENE dumbbells:
8 h i pffiffiffi < 4eBB rB 12 rB 6 for 0 < r<6 2r B 0 r r RLJ / ðrÞ ¼ pffiffiffi : BB e0 for rP6 2r B
ð13Þ
where eBB 0 denotes the strength of the potential and r is the distance between the two beads under consideration. 3. Numerical scheme 3.1. Brownian dynamics simulation The numerical procedure of BD simulation is the same as that used in our previous study [28] and hence we explain the procedure only briefly and omit its detail in the present paper. The motions of model particles and beads in dumbbells are numerically simulated by a BD simulation. Forces acting on a disk-like particle are the particle–particle interaction force calculated from the GB potential F GB , the particle–bead interaction force calculated from the MGB potential F MGB , the hydrodynamic drag force F V , and the Brownian force F B . When the acceleration term is neglected, the equation of motion for particle m is expressed by the force balance equation: MGB V B F GB ðmÞ þ F ðmÞ F ðmÞ þ F ðmÞ ¼ 0:
ð14Þ
Furthermore, the torque balance equation for particle m is given by:
1 2
rGBB ¼ ðr0 þ rB Þ;
ð12Þ
where ^ r mn is a unit vector indicating the direction from the centroid of a bead n to the center of particle m. The equations of the MGB model are not given in their entirety here because the details regarding the shape of this function have been described by Lintuvuori et al. [39] and the notation employed in the present paper is the same as that used in their paper: l and v0 . Terms with the superscript GB-B correspond to terms with the superscript GB-
MGB V B T GB ðmÞ þ T ðmÞ T ðmÞ þ T ðmÞ ¼ 0;
ð15Þ
where T denotes the torque and the superscripts are used in a similar way to the description of the force. The potential force F GB and the torque T GB are computed by ij ij GB ^ ^ @/GB =@ r and u ð@/ =@ r^i Þ, respectively. i i ij ij The hydrodynamic drag force and torque acting on a particle are computed using resistance functions for oblate spheroids [28]:
n o ^mu ^mu ^ m þ Y A ðI u ^ m Þ ðv m U Þ; F VðmÞ ¼ 3pgs r0 X A u
ð16Þ
4
T VðmÞ
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
n o ^mu ^mu ^ m þ Y C ðI u ^mÞ ¼ pgs r30 X C u ^mu ^ m Þ : D; ðxm XÞ pgs r30 Y H ðe u
ð17Þ
where I is the unit tensor, e is Eddington’s epsilon. v m and U are the velocity vector of a particle m and that of background flow, respectively. xm and X are the angular velocity vector of a particle m and that of background flow, respectively. X A ; Y A ; X C ; Y C , and Y H are resiatance functions that depend on the aspect ratio of oblate particles, and their concrete expressions are available in Ref. [28]. The Brownian force F B and the Brownian torque T B are assumed to be Gaussian [28], while they have different variances between the normal and the tangential directions to the rotation axis. The averages of both F B and T B in both normal and the tangential directions are zero, and their variances for the normal (?) and the tangential (k) directions are as follows:
6pgs r0 Y A kB T ; Dt A 6pgs r0 X kB T hF Bk ðDtÞ F Bk ðDtÞi ¼ : Dt
Fig. 3. Definitions of shear flow, coordinate system, and director angles /in and /out .
hF B? ðDtÞ F B? ðDtÞi ¼
U0* = 14
hT B? ðDtÞ T B? ðDtÞi ¼
ð19Þ
The equation of motion of a bead of a FENE dumbbell is expressed by the force balance equation. The forces acting on a bead are the bead–particle interacting force F MGB , the bead–bead interaction force computed from the RLJ potential F RLJ , the spring force F S computed from Eq. (1), the hydrodynamic drag force F VB computed as the Stokes drag force using the resistance coefficients fn and ff , and the Brownian force F BB . Consequently, the force balance equation for bead m is given by: RLJ VB BB S F MGB ðmÞ þ F ðmÞ F ðmÞ þ F ðmÞ þ F ðmÞ ¼ 0:
ð20Þ
Note that the friction coefficient of beads has a constant value in our previous study [28], while in the present simulation, it varies according to the type of beads, i.e., ’node’ or ’free’, and the coefficient for node beads is Z times larger than that for free beads. The Brownian force F BB is assumed to be Gaussian with zero mean and the variance of 2fa kB T=Dt (a=n or f) [28]. We rewrite the governing equations in a non-dimensionalized form using the representative length r0 , time sr = 3pgs r30 =ð4e0 Þ, and temperature T r = 4e0 =kB . The numerical simulation was performed based on the governing equations in the non-dimensional form. According to this formulation, the number density q, shear rate c_ , and angular velocity x are scaled as q =r30 q; c_ =sr c_ , and x =sr x, respectively. Here the variables with an asterisk denote dimensionless variables. The discrete equation of motion for beads in dimensionless form is expressed as follows:
rðmÞ ðt þ Dt Þ ¼ rðmÞ ðt Þ þ
o Dt r0 n MGB BB S F ðmÞ þ F RLJ ðmÞ þ F ðmÞ þ F ðmÞ Rf r B
þ Dt U ; Rf ¼
Z
for node beads
1 for free beads
Pd
ð18Þ
ð21Þ ð22Þ
where r ðmÞ is the position vector of bead m; Dt is the time interval. Solving the equations of motion, one obtains the velocity vector of both particles and beads and the angular velocity vector of particles. Then the position and orientation of both particles and dumbbells are updated according to the vectors obtained.
Probability
2pgs r30 Y C kB T ; Dt C 2pgs r30 X kB T hT Bk ðDtÞ T Bk ðDtÞi ¼ : Dt
1
U0* = 17 0.5
U0* = 20 U0* = 23
Pd Pa
0 0
0.5 *
1 *
Q /Q 0 Fig. 4. Probability functions for dissociation pffiffiffi pffiffiffiffiffiffiffiffiffiffi (P d ) and association (P a ) processes: sf ¼ 8 109 , d ¼ 3 8 102 , Q 0 ¼ 0:08, a0 ¼ 0:83, a1 ¼ 1, a2 ¼ 1, T ¼ 1, and Dt ¼ 0:1.
We investigate the shear flow of the system at several shear rates. The flow direction is x and the velocity gradient direction is y as shown in Fig. 3. The flow is assumed to be homogeneous and hence, in a computational region, the velocity vector of background flow U is denoted as c_ y ex , where ex is the base vector in the x direction. The Lees–Edwards moving boundary conditions [40] are used for the simulation of simple shear flows. 3.2. Model parameters The parameters of RN model are set as follows: sf =8 pffiffiffi 109 ; d =3 8 102 ; a0 =0.83, a1 =1, and a2 =1. Probability functions Pd and Pa used in the present simulation are plotted in Fig. 4. In the present simulation, the number of disk-like particles in a computational cell N and the number density of the particles q are 180 and 2, respectively. We chose the number of particles to obtain the results of simulation in a reasonable computational time. We preliminary performed computations at several shear rates for 320 particles and confirmed that the difference in computational results between N = 180 and N = 320 was within the acceptable range; for example the relative error in viscosity between the two cases was around 8% for c_ ¼ 102 . The number of dumbbells in a computational cell N d and the number density of the particles qd are 1350 and 15, respectively. The dimensionless time step Dt was chosen as it ranges between 5 109 and 5 108 . We preliminary confirmed for a few simulation conditions that the dependence of numerical results on Dt is small when the time step is shorter than 107 .
5
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
4. Results and discussion 4.1. Evaluation of flow-induced structures We analyzed the orientation behavior of particles and dumbbells in a system using the director d and the orientation order parameter S. The director is a unit vector indicating the average direction of particles or dumbbells. The director of particles is ^ m , and the director of dumbbells is calculated by the average of u calculated by the average of unit vectors Q ij =jQ ij j. In the present simulation, the behavior in shear flows is considered and d is represented by the in-plane angle /in and the out-of-plane angle /out that are defined as depicted in Fig. 3. The orientation order parameter S is defined by K 1X 3ðcos2 bi Þ 1 ; K i¼1 2
ð23Þ
where K is the number of particles (N) or dumbbells (N d ) and bi is the angle between the vector indicating the orientation of particle i or dumbbell i and the director d. The orientation order parameter S is zero when the system is in a completely random state and is unity when the system is in perfect alignment. Furthermore, we analyzed the system structures using the radial distribution function gðr Þ and the fraction of active and dangling dumbbells. The radial distribution function gðr Þ is expressed by
gðr Þ ¼
nðr ; t Þ ; 4pqr 2 Dr
ð24Þ
where nðr ; t Þ is the ensemble average of the number of particles located in a distance between r and r þ Dr from a targeted particle at time t . 4.2. Reversible network model Before analyzing the results for the polymer/clay system, we consider the results for a single phase system of the RN model. The dimensionless temperature T is 1. The entanglement strength U 0 U 0 =ðkB TÞ is varied as 14, 17, 20, and 23; The shear rate c_ is varied as 1, 4, 10, 4 10; 102 ; 4 102 ; 103 ; 4 103 ; 104 . The Peclet number Pe defined by
Pe
3pgS r 3B c_ c_ ¼ 4kB T 4T
ð25Þ
r0
and the Weissenberg number We defined by
We
2 Q 0
3pgS r B c_ ¼ 4H 4bT
rB
r0
0.5
10
0
10
2
10
4
Shear rate γ* Fig. 5. Active dumbbell fraction /a as a function of shear rate c_ at U 0 ¼ 14, 17, 20, and 23.
Fig. 5 shows the relation between the active dumbbell fraction /a and shear rate c_ . The fraction /a is higher for larger U 0 . /a is very small for U 0 ¼ 14 and 17; For U 0 ¼ 20 and 23, /a slightly depends on c_ at low shear rates, increases at c_ of around 102 , and remarkably decreases at c_ of around 103 . This phenomenon is relevant to the distribution of dumbbell length. In Fig. 6, the fraction of dumbbell length wd for the system of U 0 ¼ 23 at c_ ¼ 10, 102 ; 103 , and 104 are indicated. The dumbbell length distribution is relatively narrow at low shear rates; it is broad and the fraction of long dumbbells increases at high shear rates. At c_ =103 , strongly stretched dumbbells appear and P d rapidly increases as the dumbbell length approaches the maximum length, while P a increases just moderately as shown in Fig. 4. Consequently, the rapid decrease in /a occurs around c_ =103 as shown in Fig. 5. Fig. 7 shows the dependence of the in-plane angle /in on shear rate for U 0 ¼ 14, 17, 20, and 23. Error bars indicate the variance. At low shear rates, the variance is quite large and the degree of orientation is low as indicated in Fig. 8. Consequently, the system is in a nearly random state at low shear rates. At high shear rates, the dumbbells tend to be aligned in the flow direction and the angle /in approaches zero as shear rate increases. This tendency is more remarkable for larger U 0 because the fraction of active dumbbells, which are less affected by Brownian motion as compared to dangling dumbbells, increases with increasing U 0 as shown in Fig. 5.
0.03
3
rB
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0
c_
ð26Þ
for each c_ are listed in Table 1. Each calculation was performed until shear strain c_ t of 200 was imposed on a system. Table 1 Numerical conditions. The Peclet number Pe and the Weissenberg number We at each shear rate c_ .
c_
1
4
10
4 10
102
4 102
103
4 103
104
Pe We
0.002 0.0004
0.008 0.0016
0.02 0.004
0.08 0.016
0.2 0.04
0.8 0.16
2 0.4
8 1.6
20 4
Dumbbell length fraction ψd
S¼
1
Active dumbbell fraction φa
The FENE parameter b=HQ 20 =ðkB TÞ is set to 10 and Q 0 =Q 0 /r0 is pffiffiffiffiffiffiffiffiffiffi fixed to 0:08. The aspect ratio of oblate spheroid particles r a , which is defined by (minor axis length)/(major axis length), is fixed to 0.2. The values of the potential function parameters employed are: l = 1, m ¼ 0:625, and er ¼ 5 for the GB potential; l = 1, eGBB / 0 BB e0 =0.54, eGBB ¼ 5, and r / r ¼ 0:2 for the MGB potential; e B 0 r 0 / e0 ¼ 0:04 for the RLJ potential.
4
0
γ * = 10
0.03 3
γ * = 10 0 0.03
2
γ * = 10 0 0.03
γ * = 10 0 0.4
0.6
0.8
1
Dumbbell length Q*/ Q0* Fig. 6. Dumbbell length fraction wd for the system of U 0 =23 at c_ =10, 102 ; 103 , and 104 .
6
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
In−plane angle φ in (deg)
40
20
0
U0* = 14 U0* = 17 U0* = 20 U0* = 23
−20
−40
1
10
2
3
10
10
4
10
Shear rate γ* Fig. 7. Dependence of the in-plane angle /in on shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
Orientational order parameter Sd
1
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.5
0 0
10
2
10
4
10
Shear rate γ * Fig. 8. Orientation order parameter for dumbbells Sd as a function of shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
Fig. 8 shows the shear rate dependence of the orientation order parameter for dumbbells Sd . The degree of orientation increases consistently with shear rate. The orientation order parameter Sd for each U 0 is similar at high shear rates, whereas at low shear rates it is higher for larger U 0 . At very low shear rates, Sd is similar at low value for each U 0 because the system is in a nearly random state. At relatively low to moderate shear rates, Sd increases with increasing shear rate and is lower for smaller U 0 . In this region, the effect of Brownian motion is relatively large. Sd tends to decrease when the fraction of node beads, which is more affected by the Brownian motion, increases. Consequently, the degree of dumbbell orientation is lower for smaller U 0 , at which the active dumbbell fraction is smaller. On the other hand, at high shear rates, Sd takes a similar value again for each U 0 because dumbbells are well aligned in flows and the orientation degree is high. Furthermore, in cases of U 0 ¼ 20 and 23, Sd decreases between c_ ¼ 103 and 4 103 . This phenomenon is caused by the rapid decrease in the active dumbbell fraction in this region.
4.3. Polymer/clay systems We next consider the shear flow of polymer/clay systems. The dimensionless temperature T is 1. The value of U 0 is varied as 14, 17, 20, and 23. The shear rate c_ is varied as 1, 4, 10,
Fig. 9. Snapshots for the polymer/clay systems of U 0 ¼ 14 (a) and 23 (b) after shear strain of 200 is imposed at c_ ¼ 103 .
4 10; 102 ; 4 102 ; 103 ; 4 103 , and 104 . The definition of Pe and We is equivalent to the case of the RN model, and numerical conditions corresponding to each c_ are indicated in Table 1. Each calculation was performed until shear strain of 200 was imposed on a system. Fig. 9 shows snapshots for the polymer/clay systems of U 0 ¼ 14 (a) and 23 (b) after shear strain of 200 has been imposed at c_ ¼ 103 . Disk-like particles and free beads are depicted by disks and spheres, respectively. Node beads and free beads are colored by light and dark gray, respectively. Many active dumbbells are observed in case of U 0 ¼ 14, while less active dumbbells are seen for U 0 ¼ 23. We next analyze the shear-induced structure of the polymer/ clay system. Fig. 10 shows the dependence of the active dumbbell fraction /a on shear rate. We can see trends similar to the system of polymers described by the RN model. Because the dissociation probability P d is high at small U 0 , the number of dangling dumbbells is larger at smaller U 0 . Consequently, the active dumbbell fraction is small even at low shear rates for U 0 ¼ 14 and 17. On the other hand, /a keeps a large value over wide range of c_ for U 0 ¼ 20 and 23, while it rapidly decreases around c_ ¼ 104 because of the rapid increase in P d due to the increases in c_ .
7
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
40
In−plane angle φ in (deg)
Active dumbbell fraction φa
1
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.5
0
(dumbbells)
20
0 U0* = 14 U0* = 17 U0* = 20 U0* = 23
−20
−40 10
1
10
2
10
3
10
4
10
1
Fig. 10. Active dumbbell fraction /a as a function of shear rate c_ for the polymer/ clay system of U 0 =14, 17, 20, and 23.
U0* = 14 U0* = 23
4
10
1
0 0.03 3
γ * = 10
0 0.03 γ
*=
2
10
0 0.03 γ * = 10
0 0.2
0.4
10
3
10
4
0.6
0.8 *
Dumbbell length Q /
Fig. 12. Orientation angle /in for dumbbells as a function of shear rate c_ for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
Orientational order parameter Sd
Dumbbell length fraction ψd
0.03 γ
2
Shear rate γ*
Shear rate γ *
*=
10
1
Q0*
Fig. 11. Dumbbell length fraction wd for the polymer/clay systems of U 0 ¼ 14 and 23 at c_ ¼ 10, 102 ; 103 , and 104 .
Fig. 11 shows the fraction of dumbbell length wd at several shear rates. At c_ ¼ 10, the dumbbell length is distributed widely for small U 0 systems, in which dangling dumbbells are dominant. Dangling dumbbells are affected by Brownian motion strongly as compared to active dumbbells. Then the dumbbell length fluctuates more largely for smaller U 0 system. For U 0 ¼ 23, the distribution is narrow and has a peek around Q ¼ 0:6. In large U 0 systems, active dumbbells are dominant and hence the effect of Brownian motion is relatively small. Consequently, a narrow distribution is formed. In addition, at low shear rates, spring force can keep the dumbbell length be relatively short. At higher shear rate c_ ¼ 102 , the distribution for U 0 ¼ 23 widens because dumbbells are stretched by shear deformation then the fraction of long dumbbells increases. As the shear rate increases, the distribution tends to be broad. In the system of U 0 ¼ 23 at c_ ¼ 104 , the number of dumbbells stretched to a length close to their maximum length increases. As a result, dangling probability increases to cause the rapid decrease in /a . In Fig. 12, the orientation angle /in for dumbbells is indicated. Error bars indicate the variance. The average orientation angle ranges from 5° to 25° and slightly depends on shear rate. At low shear rates or for small U 0 , the variance is large because of the effect of Brownian motion. However, the variance is small as compared to the system of dumbbells only because the orientation of disk-like particles affect the dumbbell motion. The variance
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.5
0 1
10
2
10
3
10
4
10
Shear rate γ * Fig. 13. Orientation order parameter for dumbbells Sd as a function of shear rate c_ for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
decreases with increasing shear rate, and dumbbells tend to be aligned in the flow direction. These results mean that the existence of clay particles affect the polymer orientation behavior. In addition, no tumbling behavior is observed in this system although it appears in systems of disk-like particles [27,28]. It is probably because that the orientation of polymers restricts the synchronization of the rotation of each particle. As for the out-of-plane angle /out , it fluctuates around zero and the out-of-plane orientation is not observed in the present simulation. Fig. 13 shows the orientation order parameter Sd for dumbbells. For the system of U 0 ¼ 14 and 17, Sd increases consistently with c_ . On the other hand, for the system of U 0 ¼ 20 and 23, Sd also increases with increasing c_ below c_ ¼ 103 and rapidly decreases at c_ above 103 . It is due to the decrease in /a shown in Fig. 10. This phenomenon does not appear in a polymer/clay system with dumbbells having no association/dissociation process [28]. Next, we analyze the behavior of disk-like particles. The orientation angle indicated in Fig. 14 slightly depends on U 0 and takes angles that the rotation axis of a disk-like particle is nearly perpendicular to the flow direction. The variance of orientation angle for each condition was within 4°. The orientation degree of the particles, S, gradually increases with increasing shear rate as shown in Fig. 15. Although S takes a similar value at high shear rates, the orientation degree depends
8
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
U0* = 14
95
90
85 10
1
10
2
10
3
10
4
1 0
0
4
γ * = 10
3
γ * = 10
2
1 0 1
γ * = 10 0
Fig. 14. Orientation angle /in for disk-like particles as a function of shear rate c for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
γ * = 10
1
0
Shear rate γ *
1
2
Particle−particle distance r *
_
Fig. 16. Radial distribution function gðr Þ of disk-like particles in the polymer/clay systems of U 0 ¼ 14 and 23 at c_ ¼ 10, 102 ; 103 , and 104 .
apparent at higher shear rate c_ ¼ 103 for both U 0 ¼ 14 and 23. These results indicate flow-induced structures of disk-like particles are formed at those shear rates. At very high shear rate c_ ¼ 104 , a peak of the face-to-face configuration is remarkable. The results in S; /in , and the radial distribution function indicate that the shear flow induces a structure that particles are well aligned as their face normals are oriented in the velocity gradient direction with forming a face-to-face configuration.
1
Orientational order parameter S
U0* = 23
(face−to−face) (side−by−side)
U0* = 14 U0* = 17 U0* = 20 U0* = 23
Radial distribution function g
In−plane angle φ in (deg)
100
0.8
0.6 U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.4
10
1
10
2
10
3
10
4
Shear rate γ * Fig. 15. Orientation order parameter for disk-like particles S as a function of shear rate c_ for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
4.3.1. Rheology We investigate also rheological properties of the present model. The stress tensor can be estimated using data of the simulation of shear flows and then shear viscosity and normal stress differences are evaluated: The stress tensor r is described by
r ¼ rs þ rp þ rv þ re þ rm ; rs ¼ 2gs D;
rpab ¼ on U 0 over wide range of c_ . This result indicates that the orientation degree is affected by the polymer orientation, and is larger for larger U 0 , at which the degree of surrounding polymer melts are well aligned in the flow direction. Similarly to the results in our previous study [28], the orientation of polymers accelerates the alignment of particles to increase the degree of particle orientation. Dykes et al. [14] evaluated the shear-induced orientation for dispersion systems of montmorillonite or fluorohectrite clays in PDMS and PDPS copolymer matrix using a shear rate dependent anisotropy factor and orientation angle. The anisotropy factor takes a value between zero and unity, zero for a random state and unity for a perfect orientation. Clay particles tend to be oriented as the particle surface normal is perpendicular to the flow direction as shear rate increases. Furthermore, the anisotropy factor gradually increases with increasing shear rate. These results are qualitatively consistent with the prediction of the present simulation. Fig. 16 shows the radial distribution function gðr Þ for U 0 ¼ 14 and 23 at several shear rates. At relatively low shear rate c_ ¼ 10, distinct peak is not seen, while a broad peak around r ¼ 1 that corresponds to a side-by-side configuration (see inserted figures in Fig. 16) is seen. At c_ ¼ 102 , a peak emerges at r around 0.4 for U 0 ¼ 23. This peak corresponds to a face-to-face configuration. In addition, a broad peak around r ¼ 1 is also seen. The peaks of both the face-to-face and side-by-side configurations become
1 V 1 V
N1 X N X i¼1 j>i
ð28Þ
Nb N X 1X rPP F GB rPB F MGB ij ij ij ij a b a b V i¼1 j¼1
N Nb b 1X X
r BB ij
i¼1 j>i
ð27Þ
Nd 1X F LJ þ ðQ k Þa F Sk ij a b b V k¼1
ð29Þ
rv ¼ 2gs ufAhu^ u^ u^ u^ i : D þ Bðhu^ u^ i D þ D hu^ u^ iÞ þ CD þ FDr0 hu^ u^ ig; ð30Þ 2 r 1 ^u ^ i; re ¼ 3qkB T a2 hu ra þ 1
ð31Þ
2Nd kB T I; V
ð32Þ
rm ¼
PB BB where r PP ij ; r ij , and r ij are the unit vector from a particle i to a particle j, that from a bead j to a particle i, and that from a bead j to a bead i, respectively; Q k (=Q ij ) is the connection vector of a dumbbell k that consists of beads i and j, and F sk (=F sij ) is the spring force vector between the beads i and j; N b is the number of beads; The coefficients A; B; C, and F are functions of the aspect ratio and shape of an oblate particle [41,42]; u ¼ qpra r30 =6 is the volume fraction of oblate particles; Dr0 ¼ 3kB Tflnð2ra Þ 0:5g=ðpgs r30 Þ is the rotary diffusivity of oblate particles. Fig. 17 shows the relative viscosity gr defined by gr ¼ g=gs , where g is the shear viscosity of system. The relative viscosity shows shear-thinning properties, which are often observed for polymer/clay nanocomposites [6,8,11,14]. The
Relative viscosity ηr
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
the effect of polymers on the particle orientation in shear flows using the present model. The numerical results indicate that the disk-like particles were well aligned in flows and that the active/dangling transition of dumbbells affects the orientation behavior of the particles. Furthermore, it is found that the degree of particle orientation is higher when the intensity of entanglement of polymers U 0 is stronger because polymers with larger U 0 tend to be more aligned in the flow direction and affect the particle motion behavior. Furthermore, the present model predicts qualitatively reasonable rheological properties of polymer/clay nanocomposites.
U0* = 14 U0* = 17 U0* = 20 U0* = 23
5
9
4
3
Acknowledgments 1
2
10
3
10
4
10
10
This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI, Grant-in-Aid for Science Research (B), No. 20360084.
Shear rate γ * Fig. 17. Relative viscosity gr as a function of shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
Normal stress difference N1* , N2*
Appendix A. Supplementary material U0* = 14 U0* = 17 U0* = 20 U0* = 23
200
Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jnnfm.2012. 06.005.
N1*
References 100
0
N2*
10
1
10
2
10
3
10
4
Shear rate γ * Fig. 18. The first and second normal stress differences N 1 and N 2 as a function of shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
viscosity is higher for larger U 0 because the entanglement strength is larger and it makes the polymer matrix viscosity higher. Fig. 18 shows the first normal stress difference N 1 and the second normal stress difference N 2 . Both N 1 and N 2 change its sign at a certain shear rate. The first normal stress difference is negative at low shear rates and increases with increasing shear rate. On the other hand, N 2 decreases from positive to negative with increasing shear rate. In both cases, properties of disk-like particles is dominant at low shear rates and the contribution of polymers exceeds as shear rate increases. These changes in sign of the first normal stress difference have been reported for some polymer/clay nanocomposites [11]. Lin-Gibson et al. [11] suggested that negative N 1 was due to tumbling or continuous reorientation of clay particles. 5. Conclusion We proposed a Brownian dynamics simulation model for the analysis of the behavior of polymer/clay nanocomposites under flows. In the present model, polymer/clay nanocomposites are described by dispersion systems of disk-like particles in polymer melt. The particles are represented by oblate spheroid particles that interact with each other via the Gay–Berne potential, and the polymer melts are modeled using a reversible network model. We investigated the orientation behavior of disk-like particles and
[1] Y. Kojima, A. Usuki, M. Kawasumi, A. Okada, Y. Fukushima, T. Kurauchi, O. Kamigaito, Mechanical properties of nylon 6-clay hybrid, J. Mater. Res. 8 (1993) 1185–1189. [2] Y. Kojima, A. Usuki, M. Kawasumi, A. Okada, T. Kurauchi, O. Kamicaito, K. Kaji, Fine structure of nylon-6-clay hybrid, J. Polym. Sci. B 32 (1994) 625–630. [3] P.C. LeBaron, Z. Wang, T.J. Pinnavaia, Polymer-layered silicate nanocomposites: an overview, Appl. Clay Sci. 15 (1999) 11–29. [4] M. Okamoto, P.H. Nam, P. Maiti, T. Kotaka, N. Hasegawa, A. Usuki, A house of cards structure in polypropylene/clay nanocomposites under elongational flow, Nano Lett. 1 (2001) 295–298. [5] R.A. Vaia, E.P. Giannelis, Liquid crystal polymer nanocomposites: direct intercalation of thermotropic liquid crystalline polymers into layered silicates, Polymer 42 (2001) 1281–1285. [6] T.D. Fornes, P.J. Yoon, H. Keskkula, D.R. Paul, Nylon 6 nanocomposites: the effect of matrix molecular weight, Polymer 42 (2001) 9929–9940. [7] P.J. Yoon, T.D. Fornes, D.R. Paul, Thermal expansion behavior of nylon 6 nanocomposites, Polymer 43 (2002) 6727–6741. [8] A. Lele, M. Mackley, G. Galgali, C. Ramesh, In situ rheo-X-ray investigation of flow-induced orientation in layered silicate-syndiotactic polypropylene nanocomposite melt, J. Rheol. 46 (2002) 1091–1110. [9] S.S. Ray, M. Okamoto, Polymer/layered silicate nanocomposites: a review from preparation to processing, Prog. Polym. Sci. 28 (2003) 1539–1641. [10] A. Bafna, G. Beaucage, F. Mirabella, S. Mehta, 3D hierarchical orientation in polymer-clay nanocomposite films, Polymer 44 (2003) 1103–1115. [11] S. Lin-Gibson, H. Kim, G. Schmidt, C.C. Han, E.K. Hobbie, Shear-induced structure in polymer-clay nanocomposite solutions, J. Colloid Interface Sci. 274 (2004) 515–525. [12] T. Aubry, T. Razafinimaro, P. Médéric, Rheological investigation of the melt state elastic and yield properties of a polyamide-12 layered silicate nanocomposite, J. Rheol. 49 (2005) 425–440. [13] S.C. Tjong, Structural and mechanical properties of polymer nanocomposites, Mater. Sci. Eng. R 53 (2006) 73–197. [14] L.M.C. Dykes, J.M.T. Torkelson, W.R. Burghardt, R. Krishnamoorti, Shearinduced orientation in polymer/clay dispersions via in situ X-ray scattering, Polymer 51 (2010) 4916–4927. [15] A. Balazs, V. Ginzburg, F. Qui, G. Peng, D. Jasnow, Multi-scale model for binary mixtures containing nanoscopic particles, J. Phys. Chem. B 104 (2000) 3411– 3422. [16] V. Ginzburg, A. Balazs, Calculating phase diagrams for nanocomposites: the effect of adding end-functionalized chains to polymer/clay mixtures, Adv. Mater. 23 (2000) 1805–1809. [17] V. Ginzburg, C. Gibbons, F. Qiu, G. Peng, A. Balazs, Modeling the dynamic behavior of diblock copolymer/particle composites, Macromolecules 33 (2000) 6140–6147. [18] V. Ginzburg, F. Qiu, A. Balazs, Three-dimensional simulations of diblock copolymer/particle composites, Polymer 43 (2002) 461–466. [19] R. Thompson, V. Ginzberg, M. Matsen, A. Balazs, Block copolymer-directed assembly of nanoparticles: Forming mesoscopically ordered hybrid materials, Macromolecules 35 (2002) 1060–1071. [20] M. Vacatello, Monte Carlo simulations of polymer melts filled with solid nanoparticles, Macromolecules 34 (2001) 1946–1952.
10
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
[21] F. Starr, T. Schroder, S. Glotzer, Effects of nanoscopic filler on the structure and dynamics of a simulated polymer melt and the relationship to ultrathin films, Phys. Rev. E 64 (2001) 021802. [22] J.S. Smith, D. Bedrov, G.D. Smith, A molecular dynamics simulation study of nanoparticle interactions in a model polymer–nanoparticle composite, Compos. Sci. Technol. 63 (2003) 1599–1605. [23] A.P. Singh, A.D. Rey, Microstructure constitutive equation for discotic nematic liquid crystalline materials. Part I. Selection procedure and shear flow predictions, Rheol. Acta 37 (1998) 30–45. [24] A.P. Singh, A.D. Rey, Microstructure constitutive equation for discotic nematic liquid crystalline materials. Part II. Microstructure-rheology relations, Rheol. Acta 37 (1998) 374–386. [25] A.P. Singh, A.D. Rey, Effect of long-range elasticity and boundary conditions on microstructural response of sheared discotic mesophases, J. Non-Newtonian Fluid Mech. 94 (2000) 87–111. [26] E. Bertevas, X. Fan, R.I. Tanner, Simulation of the rheological properties of suspensions of oblate spheroidal particles in a Newtonian fluid, Rheol. Acta 49 (2010) 53. [27] T. Yamamoto, T. Suga, N. Mori, Brownian dynamics simulation of orientational behavior, flow-induced structure, and rheological properties of suspension of oblate spheroid particles under simple shear, Phys. Rev. E 72 (2005) 021509. [28] T. Yamamoto, H. Kasama, Brownian dynamics simulation of multiphase suspension of disc-like particles and polymers, Rheol. Acta 49 (2010) 573– 584. [29] J.G. Hernández Cifre, Th.M.A.O.M. Barenbrug, J.D. Schieber, B.H.A.A. van den Brule, Brownian dynamics simulation of reversible polymer networks under shear using a non-interacting dumbbell model, J. Non-Newtonian Fluid Mech. 113 (2003) 73–96. [30] H.R. Warner Jr., Kinetic theory and rheology of dilute suspensions of finitely extensible dumbbells, Ind. Eng. Chem. Fundamen. 11 (1972) 379–387. [31] A. Sunarso, T. Yamamoto, Prediction of rheological properties of associative polymeric fluids using a reversible network with non-interacting FENE
[32]
[33]
[34] [35]
[36]
[37]
[38]
[39]
[40] [41]
[42]
dumbbell model, Nihon Reoroji Gakkaishi, J. Soc. Rheol. Jpn 35 (2007) 129– 136. F. Tanaka, S.F. Edwards, Viscoelastic properties of physically crosslinked networks Part 1. Non-linear stationary viscoelasticity, J. Non-Newtonian Fluid Mech. 43 (1992) 247–271. B.H.A.A. van den Brule, P.J. Hoogerbrugge, Brownian dynamics simulation of reversible polymeric networks, J. Non-Newtonian Fluid Mech. 60 (1995) 303– 334. J.G. Gay, B.J. Berne, Modification of the overlap potential to mimic a linear site– site potential, J. Chem. Phys. 74 (1981) 3316–3319. G.R. Luckhurst, R.A. Stephens, R.W. Phippen, Computer simulation studies of anisotropic systems. XIX. Mesophases formed by the Gay–Berne model mesogen, Liq. Cryst. 8 (1990) 451–464. M.A. Bates, G.R. Luckhurst, Computer simulation studies of anisotropic systems. XXX. The phase behavior and structure of a Gay–Berne mesogen, J. Chem. Phys. 110 (1999) 7087–7108. N. Mori, R. Semura, K. Nakamura, Simple shear flows of suspensions of Brownian ellipsoids interacting via the Gay–Berne potential, Mol. Cryst. Liq. Cryst. Sci. 367 (2001) 445–454. N. Mori, H. Fujioka, R. Semura, K. Nakamura, Brownian dynamics simulations for suspension of ellipsoids in liquid crystalline phase under simple shear flows, Rheol. Acta 42 (2003) 102–109. J. Lintuvuori, M. Straka, J. Vaara, Nuclear magnetic resonance parameters of atomic xenon dissolved in Gay–Berne model liquid crystal, Phys. Rev. E 75 (2007) 031707. A.W. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C 5 (1972) 1921–1929. E.J. Hinch, L.G. Leal, The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles, J. Fluid Mech. 52 (1972) 683–712. R.G. Larson, The Structure and Rheology of Complex Fluids, Oxford Univ. Press, New York, 1999.
Contents lists available at SciVerse ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: http://www.elsevier.com/locate/jnnfm
Computational model for Brownian dynamics simulation of polymer/clay nanocomposites under flow Takehiro Yamamoto ⇑, Nobuhiro Kanda Department of Mechanical Engineering, Graduate School of Engineering, Osaka University, 2-1, Yamadaoka, Suita, Osaka 565-0871, Japan
a r t i c l e
i n f o
Article history: Received 17 January 2012 Received in revised form 6 June 2012 Accepted 7 June 2012 Available online 28 June 2012 Keywords: Polymer/clay nanocomposites Brownian dynamics simulation Gay–Berne potential Reversible network model Shear flow
a b s t r a c t A Brownian dynamics (BD) simulation model is proposed for simulating the behavior of polymer/clay nanocomposites under flows. In this model, polymer/clay nanocomposites are described by dispersion systems of disk-like particles in polymer melt. The particles are represented by oblate spheroids that interact with each other via the Gay–Berne potential, and the polymer melts are modeled using a reversible network model, in which polymer networks are represented by a set of FENE dumbbells having two states, active and dangling. The particle–bead and the bead–bead interactions are also considered. In the present study, BD simulations for the shear flow of this system was conducted to investigate the orientation behavior of disk-like particles and the effect of polymers on the particle orientation. The numerical results indicate that the disk-like particles are well aligned in flows and that the active/dangling transition of dumbbells affects the orientation behavior of the particles. Furthermore, it is found that the degree of particle orientation is higher when the intensity of entanglement of polymers U 0 is stronger because polymers with larger U 0 tend to be more aligned in the flow direction and affect the particle motion behavior. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Polymer/clay nanocomposites are remarkable materials because of their superior mechanical, thermal, and optical properties. In polymer/clay nanocomposites, nanoscale particles, typically 10–100 nm in size, are dispersed in a polymer matrix as fillers. Because of their high surface-to-volume ratio, nanoscale particles effectively modify the material properties of polymers as compared to usual polymer composites. The effects of nanoscale fillers on material properties, their orientation behavior, the flowinduced structure of the system, and rheological properties have been actively studied, e.g. [1–14]. Okamoto et al. [4] have found a house of cards structure in polypropylene/clay nanocomposites under uniaxial elongational flow by TEM observation, and reported that both strong strain-thickening and rheopexy are originated from the perpendicular alignment of the silicate layers to the stretching direction. Fornes et al. [6] have investigated the dependence of both mechanical properties and rheology on the molecular weight of matrix polymer using organoclay nanocomposites based on three different molecular weight grades of nylon 6. Lele et al. [8] have carried out measurements of both in situ X-ray diffraction and rheology for capillary flow of a system of organically modified hydrophobic montmorillonite clay and syndiotactic ⇑ Corresponding author. Fax: +81 6 6879 7308. E-mail address: [email protected] (T. Yamamoto). 0377-0257/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnnfm.2012.06.005
polypropylene. They investigated the relation between the flow-induced structure and rheology, and reported that yield behavior that might be linked to the clay orientation occurs at high shear rate and that the clay particles are easily oriented by shear to achieve a high degree of orientation after the yielding transition. Lin-Gibson et al. [11] have investigated the shear-induced structure of systems of polyethylene oxide (PEO) and synthetic hectorite clay Laponite using X-ray scattering, light scattering, optical microscopy, and rheometry. This system forms physical gel via the bridging of neighboring clay particles by polymers absorbed on the clay surface. They reported that this system showed the clay particle orientation with their surface normal parallel the vorticity direction in shear flows. Recently, Dykes et al. [14] have studied the shear-induced orientation for systems of montmorillonite or fluorohectrite clays and PDMS and PDPS copolymer matrix using in situ X-ray scattering and rheometry. They introduced an anisotropy factor, which represents the orientation degree of particles and is evaluated from data of small angle X-ray scattering measurements. In a moderately concentrated system, shear rate dependent anisotropy factor and orientation angle, and partial relaxation upon flow cessation were observed. Conversely, a more highly concentrated system showed less rate dependence in orientation state, and no orientation relaxation. Because the properties of nanocomposites depend on both the orientation and distribution of the particles, the analysis of the behavior of particles in flows is an important issue, and numerous
2
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
experimental investigations have been done as mentioned above. On the other hand, numerical simulation is an alternative effective method for the analysis of particle behaviors. Monte-Carlo and self-consistent field simulations have been used to numerically investigate polymer nanocomposites to determine the effects of nanoparticle–polymer interactions in diblock copolymer/nanoparticle mixtures [15–19]. Molecular dynamics simulations have also been used to analyze the effects of nanoparticle–polymer interactions [20–22]. Singh and Rey [23–25] have proposed constitutive equations for discotic nematic liquid crystals and analyzed the flow-induced structure of discotic mesophases using the constitutive equations. Bertevas et al. [26] have proposed an algorithm of simulations for suspension of oblate spheroidal particles in a Newtonian fluid. In their algorithm, the motion of each particle is described by Jeffery’s solution modified by the interactions between the particles. They investigated the relation between rheological properties and the orientation behavior of particles. Yamamoto and coworkers [27,28] have investigated both particle orientation behavior and flow-induced structure in suspensions of disk-like particles using Brownian dynamics (BD) simulation. Yamamoto and Kasama [28] applied a dispersion system of oblate spheroids and elastic dumbbells in Newtonian solvent to a model of polymer/clay nanocomposites. In the present study, we modified this model for the application to a dispersion system of clays in polymer melts. Polymer melts are modeled using a reversible network (RN) model, which was originally developed for associative polymers [29]. In a RN model, polymer network is represented by a set of numerous non-interacting elastic dumbbells, and hence the deformation of network topology is not directly computed. Consequently, this coarse graining representation for polymer network can save much computational cost as compared to a numerical simulation in which the deformation network structure is simulated. We propose this system as a model for dispersion systems of clay particles in polymer melts, and numerically simulated the shear flow of this system to investigate the flow-induced structure and the orientation behavior of both disk-like particles and polymers.
node bead association Pa
We first describe the computational modeling of polymer/clay nanocomposites employed in the present study. As a model of polymer/clay nanocomposites, we consider a dispersion system of disk-like particles in a polymer melt. In the present model, the disk-like particles are approximately represented by oblate spheroids, and the polymer melt is described by a RN model. We modeled polymer melts by modifying the RN model proposed by Hernández Cifre et al. [29], which has been developed for describing the behavior of associative polymers. In this model a polymer network is expressed by a set of active and dangling segments. The former is a chain that has both ends attached to the network, and the latter is a chain that is attached to the network only one end. Active and dangling segments are modeled by two separate sets of non-interacting elastic dumbbells. A dumbbell consists of two beads linked by a finitely extensible non-linear elongational (FENE) spring [30]. In a FENE chain, the spring force between beads i and j; ðF S Þij , is expressed by
ðF Þij ¼
1 ðQ ij =Q 0 Þ2
;
free bead
Fig. 1. Schematic diagram of the association/dissociation process in a reversible network model.
Both ends of an active dumbbell belong to a polymer network and hence the beads of active dumbbell are denoted by ‘node’ beads. A dangling dumbbell has a node bead and a ‘free’ bead, which represents the end of dumbbell that is not linked to a network (Fig. 1). The resistance coefficient of a bead is different between node and free beads. In the present simulation, the resistance coefficient ratio Z defined by Z=fn =ff is larger than unity. Here fn and ff denote the resistance coefficient for node beads and that for free beads, respectively. The resistance coefficient ff is calculated by ff =3pgs rB , where gs is the Newtonian contribution to polymer viscosity and r B is the bead diameter. The transition from an active dumbbell to a dangling dumbbell and the inverse transition are expressed using probability functions Pd and Pa , respectively (Fig. 1). The RN model has been originally developed for associative polymers. However, Sunarso and Yamamoto [31] adopted this model to non-associative polymer melts by modifying the association probability functions. We used this model to adopt the RN model to the numerical simulation for polymer melts. As for the dissociation process, we used a probability function Pd proposed by Hernández Cifre et al. [29] expressed by Eqs. (2)– (4). This function describes a probability that the transition from an active dumbbell to a dangling dumbbell occurs within time interval Dt:
2 Dt ; Pd ¼ 1 exp
ð2Þ
s
2
2.1. Models for polymer melts and clay particles
HQ ij
dangling dumbbell
3
2
H2 d Q 2 7 s ¼ s0 exp 6 4 n o2 5; 4U 0 kB T 1 ðQ =Q 0 Þ2
2. Computational model
S
Pd dissociation
active dumbbell
ð1Þ
where H is the spring constant, Q ij is the connector vector between the beads, Q ij is the length of vector Q ij , and Q 0 is the maximum length of the spring.
s0 ¼ sf exp
ð3Þ
U0 : kB T
ð4Þ
Here s is the lifetime of an associated bead, s0 is the lifetime at U=U 0 , and sf is a fundamental molecular time scale related to the thermal vibration. Tanaka and Edwards [32] estimated that sf is on the order of 1010 s. U 0 =ðkB TÞ indicates the strength of entanglement. d is the depth of entanglement potential well. Considering the numerical results of BD simulations by van den Brule and Hoogerbrugge [33], Hernández Cifre et al. [29] proposed an association probability function Pa , which indicates a probability that the transition from a dangling dumbbell to an active dumbbell occurs in Dt. The probability is an increasing function of the dumbbell length and is expressed by
" ( Pa ¼ 1 exp a0 þ a1
Q
)
1 ðQ =Q 0 Þ2
#
Dt :
ð5Þ
On the other hand, Sunarso and Yamamoto [31] proposed an association probability function for polymer melts:
a2 Q Pa ¼ 1 exp a0 þ a1 Dt ; Q0
ð6Þ
where a0 ; a1 , and a2 are model parameters. This model was proposed to avoid the shear-thickening behavior of a system by
3
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
restricting an abrupt increase in association probability. This model can predict both the strain-thickening at low elongation rates and the strain-thinning at high elongation rates, which are usually observed for polymer melts. We adopted this model to the present numerical simulation. In the present numerical simulations, the transition probability for each dumbbell is calculated and whether the transition occurs in a considering dumbbell is stochastically determined using a computationally created random number. For the dangling transition, it is also determined at random which node bead in a dumbbell changes to a free node. The effect of excluded volume is usually ignored in BD simulations of the RN model, whereas it is considered in our simulation because the excluded volume affects the interaction between polymers and clay particles. The excluded-volume effect is introduced using the repulsive Lennard–Jones (RLJ) potential as explained later. In the present simulation, disk-like particles are approximately represented by oblate spheroids. In the present system, we need to consider the interactions among the particles and polymer chains such as the particle–particle, particle–bead, and bead–bead interactions. These interactions are represented by potential functions. The Gay–Berne (GB) potential [34] is used to represent the particle–particle interaction. The potential function between particles m and n; /GB mn , is expressed by
12 6 r0 r0 ^ ^ ^ ^ ^ ^ ; /GB mn ðum ; un ; r mn Þ ¼ 4eðum ; un ; r mn Þ R R
ð7Þ
^m; u ^ n ; ^rmn Þ þ r0 ; R ¼ r mn rðu
ð8Þ
^ m is a unit vector indicating the orientation of the rotation where u axis of a spheroid particle m; ^ r mn is a unit vector representing the direction from the centroid of particle n to that of particle m, and r mn is the distance between particles m and n. r0 is the major axis of an oblate spheroid particle. In the present paper, we will leave out the explanation of the entire equations for the GB model because they are rather numerous and are available in the literatures, e.g. [34–38]. Those references also give details regarding the strength of the potential e and the inter-particle separation r. We employed the same notation used in previous studies [37,38] for the parameters of the GB potential: e0 ; er ; l, and m. As for oblate particles, the face-to-face configuration, in which a disk surface of a particle faces to that of another particle, is the most stable. A modified Gay–Berne (MGB) potential [39] is adopted for the potential between a disk-like particle m and a bead n. The equations of the MGB potential function /MGB mn are similar in form to those for the GB potential:
(
^ ^ ^ ^ /MGB mn ðum ; r mn Þ ¼ 4eðum ; r mn Þ
rGBB 0
n
12
R ol
rGBB 0 R
6 ) ;
ð9Þ
^ m ; ^r mn Þ2 ; eðu^ m ; ^rmn Þ ¼ eGBB 1 v0 ðu 0
ð10Þ
^ m ; ^rmn Þ þ rGBB R ¼ r mn rGBB ðu ; 0
ð11Þ
u^ m
particle
u^ n r^mn
(n)
(m)
u^ m
dumbbell
particle
r^mn
(m)
(n)
Fig. 2. Schematic diagram of vectors ^rmn used in the Gay–Berne and the modified Gay–Berne potentials.
Xe in their papers. The model parameter eGBB used in the present r paper represents the ratio of the potential energy minima of the end and side configurations (Fig. 2b in Ref. [28]) and corresponds to ee /eS in Ref. [39]. Typical potential curves for these functions are available in Ref. [28]. We set the values of the model parameters to make the potential well depth smaller than that of the inter-particle case. The RLJ potential /RLJ ðrÞ is used to represent the excluded-volume effect among FENE dumbbells:
8 h i pffiffiffi < 4eBB rB 12 rB 6 for 0 < r<6 2r B 0 r r RLJ / ðrÞ ¼ pffiffiffi : BB e0 for rP6 2r B
ð13Þ
where eBB 0 denotes the strength of the potential and r is the distance between the two beads under consideration. 3. Numerical scheme 3.1. Brownian dynamics simulation The numerical procedure of BD simulation is the same as that used in our previous study [28] and hence we explain the procedure only briefly and omit its detail in the present paper. The motions of model particles and beads in dumbbells are numerically simulated by a BD simulation. Forces acting on a disk-like particle are the particle–particle interaction force calculated from the GB potential F GB , the particle–bead interaction force calculated from the MGB potential F MGB , the hydrodynamic drag force F V , and the Brownian force F B . When the acceleration term is neglected, the equation of motion for particle m is expressed by the force balance equation: MGB V B F GB ðmÞ þ F ðmÞ F ðmÞ þ F ðmÞ ¼ 0:
ð14Þ
Furthermore, the torque balance equation for particle m is given by:
1 2
rGBB ¼ ðr0 þ rB Þ;
ð12Þ
where ^ r mn is a unit vector indicating the direction from the centroid of a bead n to the center of particle m. The equations of the MGB model are not given in their entirety here because the details regarding the shape of this function have been described by Lintuvuori et al. [39] and the notation employed in the present paper is the same as that used in their paper: l and v0 . Terms with the superscript GB-B correspond to terms with the superscript GB-
MGB V B T GB ðmÞ þ T ðmÞ T ðmÞ þ T ðmÞ ¼ 0;
ð15Þ
where T denotes the torque and the superscripts are used in a similar way to the description of the force. The potential force F GB and the torque T GB are computed by ij ij GB ^ ^ @/GB =@ r and u ð@/ =@ r^i Þ, respectively. i i ij ij The hydrodynamic drag force and torque acting on a particle are computed using resistance functions for oblate spheroids [28]:
n o ^mu ^mu ^ m þ Y A ðI u ^ m Þ ðv m U Þ; F VðmÞ ¼ 3pgs r0 X A u
ð16Þ
4
T VðmÞ
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
n o ^mu ^mu ^ m þ Y C ðI u ^mÞ ¼ pgs r30 X C u ^mu ^ m Þ : D; ðxm XÞ pgs r30 Y H ðe u
ð17Þ
where I is the unit tensor, e is Eddington’s epsilon. v m and U are the velocity vector of a particle m and that of background flow, respectively. xm and X are the angular velocity vector of a particle m and that of background flow, respectively. X A ; Y A ; X C ; Y C , and Y H are resiatance functions that depend on the aspect ratio of oblate particles, and their concrete expressions are available in Ref. [28]. The Brownian force F B and the Brownian torque T B are assumed to be Gaussian [28], while they have different variances between the normal and the tangential directions to the rotation axis. The averages of both F B and T B in both normal and the tangential directions are zero, and their variances for the normal (?) and the tangential (k) directions are as follows:
6pgs r0 Y A kB T ; Dt A 6pgs r0 X kB T hF Bk ðDtÞ F Bk ðDtÞi ¼ : Dt
Fig. 3. Definitions of shear flow, coordinate system, and director angles /in and /out .
hF B? ðDtÞ F B? ðDtÞi ¼
U0* = 14
hT B? ðDtÞ T B? ðDtÞi ¼
ð19Þ
The equation of motion of a bead of a FENE dumbbell is expressed by the force balance equation. The forces acting on a bead are the bead–particle interacting force F MGB , the bead–bead interaction force computed from the RLJ potential F RLJ , the spring force F S computed from Eq. (1), the hydrodynamic drag force F VB computed as the Stokes drag force using the resistance coefficients fn and ff , and the Brownian force F BB . Consequently, the force balance equation for bead m is given by: RLJ VB BB S F MGB ðmÞ þ F ðmÞ F ðmÞ þ F ðmÞ þ F ðmÞ ¼ 0:
ð20Þ
Note that the friction coefficient of beads has a constant value in our previous study [28], while in the present simulation, it varies according to the type of beads, i.e., ’node’ or ’free’, and the coefficient for node beads is Z times larger than that for free beads. The Brownian force F BB is assumed to be Gaussian with zero mean and the variance of 2fa kB T=Dt (a=n or f) [28]. We rewrite the governing equations in a non-dimensionalized form using the representative length r0 , time sr = 3pgs r30 =ð4e0 Þ, and temperature T r = 4e0 =kB . The numerical simulation was performed based on the governing equations in the non-dimensional form. According to this formulation, the number density q, shear rate c_ , and angular velocity x are scaled as q =r30 q; c_ =sr c_ , and x =sr x, respectively. Here the variables with an asterisk denote dimensionless variables. The discrete equation of motion for beads in dimensionless form is expressed as follows:
rðmÞ ðt þ Dt Þ ¼ rðmÞ ðt Þ þ
o Dt r0 n MGB BB S F ðmÞ þ F RLJ ðmÞ þ F ðmÞ þ F ðmÞ Rf r B
þ Dt U ; Rf ¼
Z
for node beads
1 for free beads
Pd
ð18Þ
ð21Þ ð22Þ
where r ðmÞ is the position vector of bead m; Dt is the time interval. Solving the equations of motion, one obtains the velocity vector of both particles and beads and the angular velocity vector of particles. Then the position and orientation of both particles and dumbbells are updated according to the vectors obtained.
Probability
2pgs r30 Y C kB T ; Dt C 2pgs r30 X kB T hT Bk ðDtÞ T Bk ðDtÞi ¼ : Dt
1
U0* = 17 0.5
U0* = 20 U0* = 23
Pd Pa
0 0
0.5 *
1 *
Q /Q 0 Fig. 4. Probability functions for dissociation pffiffiffi pffiffiffiffiffiffiffiffiffiffi (P d ) and association (P a ) processes: sf ¼ 8 109 , d ¼ 3 8 102 , Q 0 ¼ 0:08, a0 ¼ 0:83, a1 ¼ 1, a2 ¼ 1, T ¼ 1, and Dt ¼ 0:1.
We investigate the shear flow of the system at several shear rates. The flow direction is x and the velocity gradient direction is y as shown in Fig. 3. The flow is assumed to be homogeneous and hence, in a computational region, the velocity vector of background flow U is denoted as c_ y ex , where ex is the base vector in the x direction. The Lees–Edwards moving boundary conditions [40] are used for the simulation of simple shear flows. 3.2. Model parameters The parameters of RN model are set as follows: sf =8 pffiffiffi 109 ; d =3 8 102 ; a0 =0.83, a1 =1, and a2 =1. Probability functions Pd and Pa used in the present simulation are plotted in Fig. 4. In the present simulation, the number of disk-like particles in a computational cell N and the number density of the particles q are 180 and 2, respectively. We chose the number of particles to obtain the results of simulation in a reasonable computational time. We preliminary performed computations at several shear rates for 320 particles and confirmed that the difference in computational results between N = 180 and N = 320 was within the acceptable range; for example the relative error in viscosity between the two cases was around 8% for c_ ¼ 102 . The number of dumbbells in a computational cell N d and the number density of the particles qd are 1350 and 15, respectively. The dimensionless time step Dt was chosen as it ranges between 5 109 and 5 108 . We preliminary confirmed for a few simulation conditions that the dependence of numerical results on Dt is small when the time step is shorter than 107 .
5
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
4. Results and discussion 4.1. Evaluation of flow-induced structures We analyzed the orientation behavior of particles and dumbbells in a system using the director d and the orientation order parameter S. The director is a unit vector indicating the average direction of particles or dumbbells. The director of particles is ^ m , and the director of dumbbells is calculated by the average of u calculated by the average of unit vectors Q ij =jQ ij j. In the present simulation, the behavior in shear flows is considered and d is represented by the in-plane angle /in and the out-of-plane angle /out that are defined as depicted in Fig. 3. The orientation order parameter S is defined by K 1X 3ðcos2 bi Þ 1 ; K i¼1 2
ð23Þ
where K is the number of particles (N) or dumbbells (N d ) and bi is the angle between the vector indicating the orientation of particle i or dumbbell i and the director d. The orientation order parameter S is zero when the system is in a completely random state and is unity when the system is in perfect alignment. Furthermore, we analyzed the system structures using the radial distribution function gðr Þ and the fraction of active and dangling dumbbells. The radial distribution function gðr Þ is expressed by
gðr Þ ¼
nðr ; t Þ ; 4pqr 2 Dr
ð24Þ
where nðr ; t Þ is the ensemble average of the number of particles located in a distance between r and r þ Dr from a targeted particle at time t . 4.2. Reversible network model Before analyzing the results for the polymer/clay system, we consider the results for a single phase system of the RN model. The dimensionless temperature T is 1. The entanglement strength U 0 U 0 =ðkB TÞ is varied as 14, 17, 20, and 23; The shear rate c_ is varied as 1, 4, 10, 4 10; 102 ; 4 102 ; 103 ; 4 103 ; 104 . The Peclet number Pe defined by
Pe
3pgS r 3B c_ c_ ¼ 4kB T 4T
ð25Þ
r0
and the Weissenberg number We defined by
We
2 Q 0
3pgS r B c_ ¼ 4H 4bT
rB
r0
0.5
10
0
10
2
10
4
Shear rate γ* Fig. 5. Active dumbbell fraction /a as a function of shear rate c_ at U 0 ¼ 14, 17, 20, and 23.
Fig. 5 shows the relation between the active dumbbell fraction /a and shear rate c_ . The fraction /a is higher for larger U 0 . /a is very small for U 0 ¼ 14 and 17; For U 0 ¼ 20 and 23, /a slightly depends on c_ at low shear rates, increases at c_ of around 102 , and remarkably decreases at c_ of around 103 . This phenomenon is relevant to the distribution of dumbbell length. In Fig. 6, the fraction of dumbbell length wd for the system of U 0 ¼ 23 at c_ ¼ 10, 102 ; 103 , and 104 are indicated. The dumbbell length distribution is relatively narrow at low shear rates; it is broad and the fraction of long dumbbells increases at high shear rates. At c_ =103 , strongly stretched dumbbells appear and P d rapidly increases as the dumbbell length approaches the maximum length, while P a increases just moderately as shown in Fig. 4. Consequently, the rapid decrease in /a occurs around c_ =103 as shown in Fig. 5. Fig. 7 shows the dependence of the in-plane angle /in on shear rate for U 0 ¼ 14, 17, 20, and 23. Error bars indicate the variance. At low shear rates, the variance is quite large and the degree of orientation is low as indicated in Fig. 8. Consequently, the system is in a nearly random state at low shear rates. At high shear rates, the dumbbells tend to be aligned in the flow direction and the angle /in approaches zero as shear rate increases. This tendency is more remarkable for larger U 0 because the fraction of active dumbbells, which are less affected by Brownian motion as compared to dangling dumbbells, increases with increasing U 0 as shown in Fig. 5.
0.03
3
rB
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0
c_
ð26Þ
for each c_ are listed in Table 1. Each calculation was performed until shear strain c_ t of 200 was imposed on a system. Table 1 Numerical conditions. The Peclet number Pe and the Weissenberg number We at each shear rate c_ .
c_
1
4
10
4 10
102
4 102
103
4 103
104
Pe We
0.002 0.0004
0.008 0.0016
0.02 0.004
0.08 0.016
0.2 0.04
0.8 0.16
2 0.4
8 1.6
20 4
Dumbbell length fraction ψd
S¼
1
Active dumbbell fraction φa
The FENE parameter b=HQ 20 =ðkB TÞ is set to 10 and Q 0 =Q 0 /r0 is pffiffiffiffiffiffiffiffiffiffi fixed to 0:08. The aspect ratio of oblate spheroid particles r a , which is defined by (minor axis length)/(major axis length), is fixed to 0.2. The values of the potential function parameters employed are: l = 1, m ¼ 0:625, and er ¼ 5 for the GB potential; l = 1, eGBB / 0 BB e0 =0.54, eGBB ¼ 5, and r / r ¼ 0:2 for the MGB potential; e B 0 r 0 / e0 ¼ 0:04 for the RLJ potential.
4
0
γ * = 10
0.03 3
γ * = 10 0 0.03
2
γ * = 10 0 0.03
γ * = 10 0 0.4
0.6
0.8
1
Dumbbell length Q*/ Q0* Fig. 6. Dumbbell length fraction wd for the system of U 0 =23 at c_ =10, 102 ; 103 , and 104 .
6
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
In−plane angle φ in (deg)
40
20
0
U0* = 14 U0* = 17 U0* = 20 U0* = 23
−20
−40
1
10
2
3
10
10
4
10
Shear rate γ* Fig. 7. Dependence of the in-plane angle /in on shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
Orientational order parameter Sd
1
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.5
0 0
10
2
10
4
10
Shear rate γ * Fig. 8. Orientation order parameter for dumbbells Sd as a function of shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
Fig. 8 shows the shear rate dependence of the orientation order parameter for dumbbells Sd . The degree of orientation increases consistently with shear rate. The orientation order parameter Sd for each U 0 is similar at high shear rates, whereas at low shear rates it is higher for larger U 0 . At very low shear rates, Sd is similar at low value for each U 0 because the system is in a nearly random state. At relatively low to moderate shear rates, Sd increases with increasing shear rate and is lower for smaller U 0 . In this region, the effect of Brownian motion is relatively large. Sd tends to decrease when the fraction of node beads, which is more affected by the Brownian motion, increases. Consequently, the degree of dumbbell orientation is lower for smaller U 0 , at which the active dumbbell fraction is smaller. On the other hand, at high shear rates, Sd takes a similar value again for each U 0 because dumbbells are well aligned in flows and the orientation degree is high. Furthermore, in cases of U 0 ¼ 20 and 23, Sd decreases between c_ ¼ 103 and 4 103 . This phenomenon is caused by the rapid decrease in the active dumbbell fraction in this region.
4.3. Polymer/clay systems We next consider the shear flow of polymer/clay systems. The dimensionless temperature T is 1. The value of U 0 is varied as 14, 17, 20, and 23. The shear rate c_ is varied as 1, 4, 10,
Fig. 9. Snapshots for the polymer/clay systems of U 0 ¼ 14 (a) and 23 (b) after shear strain of 200 is imposed at c_ ¼ 103 .
4 10; 102 ; 4 102 ; 103 ; 4 103 , and 104 . The definition of Pe and We is equivalent to the case of the RN model, and numerical conditions corresponding to each c_ are indicated in Table 1. Each calculation was performed until shear strain of 200 was imposed on a system. Fig. 9 shows snapshots for the polymer/clay systems of U 0 ¼ 14 (a) and 23 (b) after shear strain of 200 has been imposed at c_ ¼ 103 . Disk-like particles and free beads are depicted by disks and spheres, respectively. Node beads and free beads are colored by light and dark gray, respectively. Many active dumbbells are observed in case of U 0 ¼ 14, while less active dumbbells are seen for U 0 ¼ 23. We next analyze the shear-induced structure of the polymer/ clay system. Fig. 10 shows the dependence of the active dumbbell fraction /a on shear rate. We can see trends similar to the system of polymers described by the RN model. Because the dissociation probability P d is high at small U 0 , the number of dangling dumbbells is larger at smaller U 0 . Consequently, the active dumbbell fraction is small even at low shear rates for U 0 ¼ 14 and 17. On the other hand, /a keeps a large value over wide range of c_ for U 0 ¼ 20 and 23, while it rapidly decreases around c_ ¼ 104 because of the rapid increase in P d due to the increases in c_ .
7
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
40
In−plane angle φ in (deg)
Active dumbbell fraction φa
1
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.5
0
(dumbbells)
20
0 U0* = 14 U0* = 17 U0* = 20 U0* = 23
−20
−40 10
1
10
2
10
3
10
4
10
1
Fig. 10. Active dumbbell fraction /a as a function of shear rate c_ for the polymer/ clay system of U 0 =14, 17, 20, and 23.
U0* = 14 U0* = 23
4
10
1
0 0.03 3
γ * = 10
0 0.03 γ
*=
2
10
0 0.03 γ * = 10
0 0.2
0.4
10
3
10
4
0.6
0.8 *
Dumbbell length Q /
Fig. 12. Orientation angle /in for dumbbells as a function of shear rate c_ for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
Orientational order parameter Sd
Dumbbell length fraction ψd
0.03 γ
2
Shear rate γ*
Shear rate γ *
*=
10
1
Q0*
Fig. 11. Dumbbell length fraction wd for the polymer/clay systems of U 0 ¼ 14 and 23 at c_ ¼ 10, 102 ; 103 , and 104 .
Fig. 11 shows the fraction of dumbbell length wd at several shear rates. At c_ ¼ 10, the dumbbell length is distributed widely for small U 0 systems, in which dangling dumbbells are dominant. Dangling dumbbells are affected by Brownian motion strongly as compared to active dumbbells. Then the dumbbell length fluctuates more largely for smaller U 0 system. For U 0 ¼ 23, the distribution is narrow and has a peek around Q ¼ 0:6. In large U 0 systems, active dumbbells are dominant and hence the effect of Brownian motion is relatively small. Consequently, a narrow distribution is formed. In addition, at low shear rates, spring force can keep the dumbbell length be relatively short. At higher shear rate c_ ¼ 102 , the distribution for U 0 ¼ 23 widens because dumbbells are stretched by shear deformation then the fraction of long dumbbells increases. As the shear rate increases, the distribution tends to be broad. In the system of U 0 ¼ 23 at c_ ¼ 104 , the number of dumbbells stretched to a length close to their maximum length increases. As a result, dangling probability increases to cause the rapid decrease in /a . In Fig. 12, the orientation angle /in for dumbbells is indicated. Error bars indicate the variance. The average orientation angle ranges from 5° to 25° and slightly depends on shear rate. At low shear rates or for small U 0 , the variance is large because of the effect of Brownian motion. However, the variance is small as compared to the system of dumbbells only because the orientation of disk-like particles affect the dumbbell motion. The variance
U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.5
0 1
10
2
10
3
10
4
10
Shear rate γ * Fig. 13. Orientation order parameter for dumbbells Sd as a function of shear rate c_ for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
decreases with increasing shear rate, and dumbbells tend to be aligned in the flow direction. These results mean that the existence of clay particles affect the polymer orientation behavior. In addition, no tumbling behavior is observed in this system although it appears in systems of disk-like particles [27,28]. It is probably because that the orientation of polymers restricts the synchronization of the rotation of each particle. As for the out-of-plane angle /out , it fluctuates around zero and the out-of-plane orientation is not observed in the present simulation. Fig. 13 shows the orientation order parameter Sd for dumbbells. For the system of U 0 ¼ 14 and 17, Sd increases consistently with c_ . On the other hand, for the system of U 0 ¼ 20 and 23, Sd also increases with increasing c_ below c_ ¼ 103 and rapidly decreases at c_ above 103 . It is due to the decrease in /a shown in Fig. 10. This phenomenon does not appear in a polymer/clay system with dumbbells having no association/dissociation process [28]. Next, we analyze the behavior of disk-like particles. The orientation angle indicated in Fig. 14 slightly depends on U 0 and takes angles that the rotation axis of a disk-like particle is nearly perpendicular to the flow direction. The variance of orientation angle for each condition was within 4°. The orientation degree of the particles, S, gradually increases with increasing shear rate as shown in Fig. 15. Although S takes a similar value at high shear rates, the orientation degree depends
8
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
U0* = 14
95
90
85 10
1
10
2
10
3
10
4
1 0
0
4
γ * = 10
3
γ * = 10
2
1 0 1
γ * = 10 0
Fig. 14. Orientation angle /in for disk-like particles as a function of shear rate c for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
γ * = 10
1
0
Shear rate γ *
1
2
Particle−particle distance r *
_
Fig. 16. Radial distribution function gðr Þ of disk-like particles in the polymer/clay systems of U 0 ¼ 14 and 23 at c_ ¼ 10, 102 ; 103 , and 104 .
apparent at higher shear rate c_ ¼ 103 for both U 0 ¼ 14 and 23. These results indicate flow-induced structures of disk-like particles are formed at those shear rates. At very high shear rate c_ ¼ 104 , a peak of the face-to-face configuration is remarkable. The results in S; /in , and the radial distribution function indicate that the shear flow induces a structure that particles are well aligned as their face normals are oriented in the velocity gradient direction with forming a face-to-face configuration.
1
Orientational order parameter S
U0* = 23
(face−to−face) (side−by−side)
U0* = 14 U0* = 17 U0* = 20 U0* = 23
Radial distribution function g
In−plane angle φ in (deg)
100
0.8
0.6 U0* = 14 U0* = 17 U0* = 20 U0* = 23
0.4
10
1
10
2
10
3
10
4
Shear rate γ * Fig. 15. Orientation order parameter for disk-like particles S as a function of shear rate c_ for the polymer/clay system of U 0 ¼ 14, 17, 20, and 23.
4.3.1. Rheology We investigate also rheological properties of the present model. The stress tensor can be estimated using data of the simulation of shear flows and then shear viscosity and normal stress differences are evaluated: The stress tensor r is described by
r ¼ rs þ rp þ rv þ re þ rm ; rs ¼ 2gs D;
rpab ¼ on U 0 over wide range of c_ . This result indicates that the orientation degree is affected by the polymer orientation, and is larger for larger U 0 , at which the degree of surrounding polymer melts are well aligned in the flow direction. Similarly to the results in our previous study [28], the orientation of polymers accelerates the alignment of particles to increase the degree of particle orientation. Dykes et al. [14] evaluated the shear-induced orientation for dispersion systems of montmorillonite or fluorohectrite clays in PDMS and PDPS copolymer matrix using a shear rate dependent anisotropy factor and orientation angle. The anisotropy factor takes a value between zero and unity, zero for a random state and unity for a perfect orientation. Clay particles tend to be oriented as the particle surface normal is perpendicular to the flow direction as shear rate increases. Furthermore, the anisotropy factor gradually increases with increasing shear rate. These results are qualitatively consistent with the prediction of the present simulation. Fig. 16 shows the radial distribution function gðr Þ for U 0 ¼ 14 and 23 at several shear rates. At relatively low shear rate c_ ¼ 10, distinct peak is not seen, while a broad peak around r ¼ 1 that corresponds to a side-by-side configuration (see inserted figures in Fig. 16) is seen. At c_ ¼ 102 , a peak emerges at r around 0.4 for U 0 ¼ 23. This peak corresponds to a face-to-face configuration. In addition, a broad peak around r ¼ 1 is also seen. The peaks of both the face-to-face and side-by-side configurations become
1 V 1 V
N1 X N X i¼1 j>i
ð28Þ
Nb N X 1X rPP F GB rPB F MGB ij ij ij ij a b a b V i¼1 j¼1
N Nb b 1X X
r BB ij
i¼1 j>i
ð27Þ
Nd 1X F LJ þ ðQ k Þa F Sk ij a b b V k¼1
ð29Þ
rv ¼ 2gs ufAhu^ u^ u^ u^ i : D þ Bðhu^ u^ i D þ D hu^ u^ iÞ þ CD þ FDr0 hu^ u^ ig; ð30Þ 2 r 1 ^u ^ i; re ¼ 3qkB T a2 hu ra þ 1
ð31Þ
2Nd kB T I; V
ð32Þ
rm ¼
PB BB where r PP ij ; r ij , and r ij are the unit vector from a particle i to a particle j, that from a bead j to a particle i, and that from a bead j to a bead i, respectively; Q k (=Q ij ) is the connection vector of a dumbbell k that consists of beads i and j, and F sk (=F sij ) is the spring force vector between the beads i and j; N b is the number of beads; The coefficients A; B; C, and F are functions of the aspect ratio and shape of an oblate particle [41,42]; u ¼ qpra r30 =6 is the volume fraction of oblate particles; Dr0 ¼ 3kB Tflnð2ra Þ 0:5g=ðpgs r30 Þ is the rotary diffusivity of oblate particles. Fig. 17 shows the relative viscosity gr defined by gr ¼ g=gs , where g is the shear viscosity of system. The relative viscosity shows shear-thinning properties, which are often observed for polymer/clay nanocomposites [6,8,11,14]. The
Relative viscosity ηr
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
the effect of polymers on the particle orientation in shear flows using the present model. The numerical results indicate that the disk-like particles were well aligned in flows and that the active/dangling transition of dumbbells affects the orientation behavior of the particles. Furthermore, it is found that the degree of particle orientation is higher when the intensity of entanglement of polymers U 0 is stronger because polymers with larger U 0 tend to be more aligned in the flow direction and affect the particle motion behavior. Furthermore, the present model predicts qualitatively reasonable rheological properties of polymer/clay nanocomposites.
U0* = 14 U0* = 17 U0* = 20 U0* = 23
5
9
4
3
Acknowledgments 1
2
10
3
10
4
10
10
This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI, Grant-in-Aid for Science Research (B), No. 20360084.
Shear rate γ * Fig. 17. Relative viscosity gr as a function of shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
Normal stress difference N1* , N2*
Appendix A. Supplementary material U0* = 14 U0* = 17 U0* = 20 U0* = 23
200
Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jnnfm.2012. 06.005.
N1*
References 100
0
N2*
10
1
10
2
10
3
10
4
Shear rate γ * Fig. 18. The first and second normal stress differences N 1 and N 2 as a function of shear rate c_ for U 0 ¼ 14, 17, 20, and 23.
viscosity is higher for larger U 0 because the entanglement strength is larger and it makes the polymer matrix viscosity higher. Fig. 18 shows the first normal stress difference N 1 and the second normal stress difference N 2 . Both N 1 and N 2 change its sign at a certain shear rate. The first normal stress difference is negative at low shear rates and increases with increasing shear rate. On the other hand, N 2 decreases from positive to negative with increasing shear rate. In both cases, properties of disk-like particles is dominant at low shear rates and the contribution of polymers exceeds as shear rate increases. These changes in sign of the first normal stress difference have been reported for some polymer/clay nanocomposites [11]. Lin-Gibson et al. [11] suggested that negative N 1 was due to tumbling or continuous reorientation of clay particles. 5. Conclusion We proposed a Brownian dynamics simulation model for the analysis of the behavior of polymer/clay nanocomposites under flows. In the present model, polymer/clay nanocomposites are described by dispersion systems of disk-like particles in polymer melt. The particles are represented by oblate spheroid particles that interact with each other via the Gay–Berne potential, and the polymer melts are modeled using a reversible network model. We investigated the orientation behavior of disk-like particles and
[1] Y. Kojima, A. Usuki, M. Kawasumi, A. Okada, Y. Fukushima, T. Kurauchi, O. Kamigaito, Mechanical properties of nylon 6-clay hybrid, J. Mater. Res. 8 (1993) 1185–1189. [2] Y. Kojima, A. Usuki, M. Kawasumi, A. Okada, T. Kurauchi, O. Kamicaito, K. Kaji, Fine structure of nylon-6-clay hybrid, J. Polym. Sci. B 32 (1994) 625–630. [3] P.C. LeBaron, Z. Wang, T.J. Pinnavaia, Polymer-layered silicate nanocomposites: an overview, Appl. Clay Sci. 15 (1999) 11–29. [4] M. Okamoto, P.H. Nam, P. Maiti, T. Kotaka, N. Hasegawa, A. Usuki, A house of cards structure in polypropylene/clay nanocomposites under elongational flow, Nano Lett. 1 (2001) 295–298. [5] R.A. Vaia, E.P. Giannelis, Liquid crystal polymer nanocomposites: direct intercalation of thermotropic liquid crystalline polymers into layered silicates, Polymer 42 (2001) 1281–1285. [6] T.D. Fornes, P.J. Yoon, H. Keskkula, D.R. Paul, Nylon 6 nanocomposites: the effect of matrix molecular weight, Polymer 42 (2001) 9929–9940. [7] P.J. Yoon, T.D. Fornes, D.R. Paul, Thermal expansion behavior of nylon 6 nanocomposites, Polymer 43 (2002) 6727–6741. [8] A. Lele, M. Mackley, G. Galgali, C. Ramesh, In situ rheo-X-ray investigation of flow-induced orientation in layered silicate-syndiotactic polypropylene nanocomposite melt, J. Rheol. 46 (2002) 1091–1110. [9] S.S. Ray, M. Okamoto, Polymer/layered silicate nanocomposites: a review from preparation to processing, Prog. Polym. Sci. 28 (2003) 1539–1641. [10] A. Bafna, G. Beaucage, F. Mirabella, S. Mehta, 3D hierarchical orientation in polymer-clay nanocomposite films, Polymer 44 (2003) 1103–1115. [11] S. Lin-Gibson, H. Kim, G. Schmidt, C.C. Han, E.K. Hobbie, Shear-induced structure in polymer-clay nanocomposite solutions, J. Colloid Interface Sci. 274 (2004) 515–525. [12] T. Aubry, T. Razafinimaro, P. Médéric, Rheological investigation of the melt state elastic and yield properties of a polyamide-12 layered silicate nanocomposite, J. Rheol. 49 (2005) 425–440. [13] S.C. Tjong, Structural and mechanical properties of polymer nanocomposites, Mater. Sci. Eng. R 53 (2006) 73–197. [14] L.M.C. Dykes, J.M.T. Torkelson, W.R. Burghardt, R. Krishnamoorti, Shearinduced orientation in polymer/clay dispersions via in situ X-ray scattering, Polymer 51 (2010) 4916–4927. [15] A. Balazs, V. Ginzburg, F. Qui, G. Peng, D. Jasnow, Multi-scale model for binary mixtures containing nanoscopic particles, J. Phys. Chem. B 104 (2000) 3411– 3422. [16] V. Ginzburg, A. Balazs, Calculating phase diagrams for nanocomposites: the effect of adding end-functionalized chains to polymer/clay mixtures, Adv. Mater. 23 (2000) 1805–1809. [17] V. Ginzburg, C. Gibbons, F. Qiu, G. Peng, A. Balazs, Modeling the dynamic behavior of diblock copolymer/particle composites, Macromolecules 33 (2000) 6140–6147. [18] V. Ginzburg, F. Qiu, A. Balazs, Three-dimensional simulations of diblock copolymer/particle composites, Polymer 43 (2002) 461–466. [19] R. Thompson, V. Ginzberg, M. Matsen, A. Balazs, Block copolymer-directed assembly of nanoparticles: Forming mesoscopically ordered hybrid materials, Macromolecules 35 (2002) 1060–1071. [20] M. Vacatello, Monte Carlo simulations of polymer melts filled with solid nanoparticles, Macromolecules 34 (2001) 1946–1952.
10
T. Yamamoto, N. Kanda / Journal of Non-Newtonian Fluid Mechanics 181–182 (2012) 1–10
[21] F. Starr, T. Schroder, S. Glotzer, Effects of nanoscopic filler on the structure and dynamics of a simulated polymer melt and the relationship to ultrathin films, Phys. Rev. E 64 (2001) 021802. [22] J.S. Smith, D. Bedrov, G.D. Smith, A molecular dynamics simulation study of nanoparticle interactions in a model polymer–nanoparticle composite, Compos. Sci. Technol. 63 (2003) 1599–1605. [23] A.P. Singh, A.D. Rey, Microstructure constitutive equation for discotic nematic liquid crystalline materials. Part I. Selection procedure and shear flow predictions, Rheol. Acta 37 (1998) 30–45. [24] A.P. Singh, A.D. Rey, Microstructure constitutive equation for discotic nematic liquid crystalline materials. Part II. Microstructure-rheology relations, Rheol. Acta 37 (1998) 374–386. [25] A.P. Singh, A.D. Rey, Effect of long-range elasticity and boundary conditions on microstructural response of sheared discotic mesophases, J. Non-Newtonian Fluid Mech. 94 (2000) 87–111. [26] E. Bertevas, X. Fan, R.I. Tanner, Simulation of the rheological properties of suspensions of oblate spheroidal particles in a Newtonian fluid, Rheol. Acta 49 (2010) 53. [27] T. Yamamoto, T. Suga, N. Mori, Brownian dynamics simulation of orientational behavior, flow-induced structure, and rheological properties of suspension of oblate spheroid particles under simple shear, Phys. Rev. E 72 (2005) 021509. [28] T. Yamamoto, H. Kasama, Brownian dynamics simulation of multiphase suspension of disc-like particles and polymers, Rheol. Acta 49 (2010) 573– 584. [29] J.G. Hernández Cifre, Th.M.A.O.M. Barenbrug, J.D. Schieber, B.H.A.A. van den Brule, Brownian dynamics simulation of reversible polymer networks under shear using a non-interacting dumbbell model, J. Non-Newtonian Fluid Mech. 113 (2003) 73–96. [30] H.R. Warner Jr., Kinetic theory and rheology of dilute suspensions of finitely extensible dumbbells, Ind. Eng. Chem. Fundamen. 11 (1972) 379–387. [31] A. Sunarso, T. Yamamoto, Prediction of rheological properties of associative polymeric fluids using a reversible network with non-interacting FENE
[32]
[33]
[34] [35]
[36]
[37]
[38]
[39]
[40] [41]
[42]
dumbbell model, Nihon Reoroji Gakkaishi, J. Soc. Rheol. Jpn 35 (2007) 129– 136. F. Tanaka, S.F. Edwards, Viscoelastic properties of physically crosslinked networks Part 1. Non-linear stationary viscoelasticity, J. Non-Newtonian Fluid Mech. 43 (1992) 247–271. B.H.A.A. van den Brule, P.J. Hoogerbrugge, Brownian dynamics simulation of reversible polymeric networks, J. Non-Newtonian Fluid Mech. 60 (1995) 303– 334. J.G. Gay, B.J. Berne, Modification of the overlap potential to mimic a linear site– site potential, J. Chem. Phys. 74 (1981) 3316–3319. G.R. Luckhurst, R.A. Stephens, R.W. Phippen, Computer simulation studies of anisotropic systems. XIX. Mesophases formed by the Gay–Berne model mesogen, Liq. Cryst. 8 (1990) 451–464. M.A. Bates, G.R. Luckhurst, Computer simulation studies of anisotropic systems. XXX. The phase behavior and structure of a Gay–Berne mesogen, J. Chem. Phys. 110 (1999) 7087–7108. N. Mori, R. Semura, K. Nakamura, Simple shear flows of suspensions of Brownian ellipsoids interacting via the Gay–Berne potential, Mol. Cryst. Liq. Cryst. Sci. 367 (2001) 445–454. N. Mori, H. Fujioka, R. Semura, K. Nakamura, Brownian dynamics simulations for suspension of ellipsoids in liquid crystalline phase under simple shear flows, Rheol. Acta 42 (2003) 102–109. J. Lintuvuori, M. Straka, J. Vaara, Nuclear magnetic resonance parameters of atomic xenon dissolved in Gay–Berne model liquid crystal, Phys. Rev. E 75 (2007) 031707. A.W. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C 5 (1972) 1921–1929. E.J. Hinch, L.G. Leal, The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles, J. Fluid Mech. 52 (1972) 683–712. R.G. Larson, The Structure and Rheology of Complex Fluids, Oxford Univ. Press, New York, 1999.