Dynamic behavior of dual cross-linked nanoparticle networks under oscillatory shear

We first consider the dynamic response of the system to oscillatory shear deformations for the case of a small amplitude of the applied shear, setting ${{\gamma }_{0}}=0.09$. In figure 2, we plot the shear force and average number of bonds in the network as functions of the shear strain, and show that the response of the network depends on both the energy of the labile bonds and frequency of oscillation. In particular, depending on the values of $U_{0}^{(l)}$ and $\omega$, the force–strain Lissajous curves form a straight line (see figures 2(a) and (b)) or an ellipse (see figures 2(b)–(d)). The Lissajous curves indicate that the PGN network exhibits linear viscoelastic behavior at ${{\gamma }_{0}}=0.09$. The tilt of an ellipse (slope of a straight line) is proportional to the storage modulus $G^{\prime}$, which characterizes the elasticity of the system. The width of an ellipse is proportional to the loss modulus $G^{\prime\prime}$, which is a measure of the dissipation in the system due to a lag between the applied strain and the resulting force.

Zoom In Zoom Out Reset image size

At a given frequency, $\omega$, the samples with the stronger labile bonds exhibit a larger tilt of the force–strain Lissajous curves, indicating that they have higher storage moduli (see figures 2(a)–(d)). The latter behavior can be attributed to the increase in the number of bonds with the increase in bond strength (see figures 2(e)–(h)). Furthermore, figure 2 shows that for both the weaker and stronger labile bonds, the number of labile bonds and the storage moduli (tilt of the Lissajous curves) increase with an increase in the frequency of oscillation, $\omega$. Finally, the dissipation in the system is seen to depend on $\omega$, as indicated by the variations in the width of the Lissajous curves in figures 2(a)–(d).

The storage and loss moduli as functions of the oscillation frequency (at the shear amplitude of ${{\gamma }_{0}}=0.09$) are shown in figures 3(a) and (b), respectively. To calculate the moduli $G^{\prime}$ and $G^{\prime\prime}$, we first determine the amplitude of the force, ${{F}_{0}}$, and the phase difference, $\delta$, by fitting the shear force $F$, which is obtained from the strain-controlled simulations at $\gamma ={{\gamma }_{0}}{\rm sin} \;(\omega {\mkern 1mu} t)$, to the equation $F={{F}_{0}}{\rm sin} \;(\omega t+\delta )$. Then, the moduli are found as:

$$\begin{eqnarray}&&G^{\prime} ={{F}_{0}}\;\gamma _{0}^{-1}{\rm cos} \;(\delta ),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&G^{\prime\prime} ={{F}_{0}}\;\gamma _{0}^{-1}{\rm sin} \;(\delta ).\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

The storage and loss moduli could also be determined directly from the Lissajous curves [15].

The storage modulus of the network shows two distinct features. First, as noted above, $G^{\prime}$ increases with an increase in the energy of the labile bonds $U_{0}^{(l)}$, so that it is consistently greater at $U_{0}^{(l)}=37{{k}_{B}}T$ than at $33{{k}_{B}}T$ (see figure 3(a)). Second, $G^{\prime}$ exhibits a transition from a lower plateau to a higher plateau value with an increase in frequency for both the $U_{0}^{(l)}$ considered here. The increase in $U_{0}^{(l)}$ shifts the transition to the lower frequencies (figure 3(a)). On the other hand, the loss modulus $G^{\prime\prime}$ as a function of $\omega$ exhibits distinctly different behavior at the weaker and stronger labile bonds (see figure 3(b)). Namely, $G^{\prime\prime}$ increases monotonically at $U_{0}^{(l)}=33{{k}_{B}}T$ (curve 1 in figure 3(b)), whereas at $U_{0}^{(l)}=37{{k}_{B}}T$, the loss modulus exhibits a peak (see curve 2 figure 3(b)). The latter peak is located within the frequency range where the storage modulus of the system, $G^{\prime}$, exhibits a transition between the two plateaus (compare curves 2 in figures 3(a) and (b)).

It is instructive to plot the moduli $G^{\prime}$ and $G^{\prime\prime}$ as functions of the normalized frequency $\omega /{{\omega }_{0}}$ (see figures 3(c) and (d), respectively), where ${{\omega }_{0}}=2\pi k_{0r}^{\left( l \right)}$ is the characteristic frequency, which depends on the energy of labile bonds $U_{0}^{(l)}$ through the rupture rate constant, $k_{0r}^{\left( l \right)}$. The storage modulus $G^{\prime}$ as a function of $\omega /{{\omega }_{0}}$ shows the same behavior for the networks containing the weaker and stronger labile bonds. Namely, $G^{\prime}$ exhibits a lower plateau regime at $\omega /{{\omega }_{0}}<1$, a higher plateau regime at $\omega /{{\omega }_{0}}>100$, and a smooth transition between the latter two regimes at $1<\omega /{{\omega }_{0}}<100$. Figure 3(c) indicates, therefore, that at low strain amplitudes, the network elasticity (storage modulus) as a function of frequency is controlled by the rupture rate $k_{0r}^{\left( l \right)}$ of the labile bonds. The latter conclusion cannot be made, however, for the frequency-dependent dissipation in the network. Figure 3(d) shows that the loss modulus $G^{\prime\prime}$ as a function of $\omega /{{\omega }_{0}}$ behaves quite differently at $U_{0}^{(l)}=33{{k}_{B}}T$ and $37{{k}_{B}}T$. The reason for the observed behavior of $G^{\prime\prime}$ is associated with dissipation due to the particle–solvent friction, which does not depend on the energy of the labile bonds and is linearly proportional to $\omega$.

The effect of the particle–solvent friction is demonstrated in figure 4, where we plot the storage and loss moduli as functions of $\omega /{{\omega }_{0}}$ at three values of the mobility, which is inversely proportional to friction. Figure 4(a) shows that the storage modulus $G^{\prime}$ is unaffected by the variations in the value of the mobility. In contrast, the loss modulus $G^{\prime\prime}$ exhibits notable changes as the mobility is varied. Within the high frequency region (see figure 4(b)), $G^{\prime\prime}$ decreases with an increase in the mobility, i.e., with a decrease in the particle–solvent friction. Namely, the effect of friction on $G^{\prime\prime}$ becomes noticeable at $\omega >0.1$ for both the weaker and stronger labile bonds. It is worth noting that at $U_{0}^{(l)}=33{{k}_{B}}T$, the behavior of $G^{\prime\prime} (\omega /{{\omega }_{0}})$ at the mobility of $2.5\mu$ indicates the onset of a maximum (at $\omega /{{\omega }_{0}}\approx 6$) that is not observed in the samples with lower mobility.

Zoom In Zoom Out Reset image size

To gain greater insight into the simulation results, we consider a simple model, which captures all the salient features of the viscoelastic behavior arising from the dynamics of the labile bonds. In our simulations, the nanoparticles in the PGN network form a hexagonal lattice at equilibrium (as observed in recent experimental studies on a comparable system [16]). We note that the hexagonal lattice is a stable configuration and alternate initial configurations of square and face centered square lattices spontaneously transform to the hexagonal lattice during equilibration. Hence, we consider the system of three particles shown in figure 5 as a representative cell of the lattice shown in figure 1, and determine the response of this cell to a strain-controlled oscillatory shear deformation. Specifically, particles 1 and 2 are fixed in space, whereas particle 3 exhibits oscillations in the horizontal direction, as indicated in figure 5. The position of particle 3 in the vertical direction is held constant at ${{Y}_{0}}={{3}^{1/2}}{{\Delta }_{0}}/2\ \;$ and the horizontal coordinate changes in time as $X(t)={{Y}_{0}}\;{{\gamma }_{0}}\;{\rm sin} (\omega \;t)$, where ${{\Delta }_{0}}$ and ${{\gamma }_{0}}$ are the equilibrium center-to-center distance between the particles and shear amplitude, respectively. The time-dependent distance between the particles in pairs 1–3 and 2–3 is calculated as $R(t)={{Y}_{0}}{{[1+{{({{3}^{-1/2}}\pm {{\gamma }_{0}}\;{\rm sin} \left( \omega \;t \right))}^{2}}]}^{1/2}}$, where the signs '+' and '−'correspond to the pairs 1–3 and 2–3, respectively. Within each pair, the dynamics of the rupture and formation of the labile bonds is described by the following master equation [17] for the probability of finding $n\geqslant 1$ labile bonds ${{p}_{n}}$:

$$\begin{eqnarray}&&\frac{{\rm d}{{p}_{n}}}{{\rm d}t}=-\left[ n\;{{k}_{r}}(R)+{{k}_{f}}(R,n) \right]{{p}_{n}}+(n+1){{k}_{r}}(R){{p}_{n+1}}+\;{{k}_{f}}(R,n-1){{p}_{n-1}}.\end{eqnarray} \tag{ 8 }$$

Here, ${{k}_{r}}(R)$ is the rate coefficient for the rupture of labile bonds defined in section 2 and ${{k}_{f}}(R,n)={{[{{N}_{{\rm max} }}(R)-n]}^{2}}{{P}_{c}}(R)\;k_{0f}^{(l)}$ is the rate coefficient for the bond formation (cf equation (5)). The probability of having no labile bonds between a pair of particles, ${{p}_{0}}$, is found from the normalization condition to be ${{p}_{0}}=1-\mathop \sum \nolimits_{n=1}^{\infty } {{p}_{n}}$. It is worth recalling that the inter-particle distance $R$ in equation (8) is a function of time $t$ (as described above).

Zoom In Zoom Out Reset image size

The system of equations generated from equation (8) was truncated at $n=17$ and solved numerically for the two pairs of particles using the Mathematica™ software. The model parameters used in these calculations were the same as used in the simulations of the array of particles. Upon solving equation (8), the value of attractive force between two particles due to bonded polymer arms was calculated according to equation (3) at the number of bonds ${{N}_{b}}=1+\mathop \sum \nolimits_{n=1}^{17} n\;{{p}_{n}}$. The latter equation means that the link consists of one permanent bond and an average number of labile bonds given by $\left\langle n \right\rangle =\sum n n\;{{p}_{n}}$. Finally, the values of $G^{\prime}$ and $G^{\prime\prime}$ were determined using the numerical Fourier transform applied to the X component of the total force acting on particle 3.

In figures 6(a) and (b), we compare the behavior of the storage and loss moduli as functions of $\omega /{{\omega }_{0}}$ for the three particle model (figure 5(a)) and for the array of particles (figure 1) at $U_{0}^{(l)}=33$ and $37{{k}_{B}}T$, and the strain amplitude ${{\gamma }_{0}}=0.09$. Only a qualitative comparison is possible since, in our simulations, the magnitudes of the calculated modulus is based on the force required to shear the array of particles and, hence, depends on the size of the array. To make these comparisons, the moduli $G^{\prime}$ and $G^{\prime\prime}$ were normalized by ${{G}_{0}}$, which corresponds to the plateau value of the storage modulus at the lower frequencies, i.e., ${{G}_{0}}=G^{\prime} (\omega /{{\omega }_{0}}\ll 1)$. The normalized storage moduli obtained from the three-particle model (solid lines) and from simulations of the array of particles (symbols) are in quantitative agreement at both values of the labile bond energies considered here (see figure 6(a)). Note that the three-particle model takes into account only the processes of rupture and formation of the labile bonds between the individual pairs of PGNs. The latter result (figure 6(a)) indicates, therefore, that at low shear strains, the individual PGN pairs contribute independently to the elastic properties of the material, i.e., there are no collective effects.

Zoom In Zoom Out Reset image size

Figure 6(b) shows the normalized loss moduli for the three-particle model (solid lines) and PGN array (symbols). For the three-particle model, $G^{\prime\prime}$ exhibits a peak at the value of $\omega /{{\omega }_{0}}$ where $G^{\prime}$ exhibits an inflection point within the transition zone between the lower and higher plateaus (figure 6(a)). Such behavior is typical for cross-linked viscoelastic polymeric systems [18]. In networks of cross-linked PGNs, the viscoelastic behavior is caused by the rupture and formation of the labile bonds, and the dissipation is maximal when the characteristic times of shear deformation and of bond rupture are comparable $(\omega /{{\omega }_{0}}\tilde{\ }1)$. Notably, for $U_{0}^{(l)}=37{{k}_{B}}T$ (the stronger bond), the peak in the loss modulus for the three-particle model coincides with the peak for the PGN array (figure 6(b)), indicating the predictive ability of the simplistic three-particle model. The peak in $G^{\prime\prime}$ at $U_{0}^{(l)}=33{{k}_{B}}T$ is not observed for the cross-linked array of PGN because it is masked by dissipation due to the particle–solvent friction. (We will return to this figure to discuss the behavior in the nonlinear regime, which is captured by figures 6(c)–(f).)

The three-particle model served to highlight the importance of the breaking and making of labile bonds in the viscoelastic response of the system to the oscillatory shear. Figure 7 shows how the number of labile bonds depends on changes in the scaled shear frequency $\omega /{{\omega }_{0}}$ for $U_{0}^{(l)}=37{{k}_{B}}T$. (Recall that ${{\omega }_{0}}$ depends on the rupture rate constant for the bond.) In this figure, we show the difference in number of labile bonds, $\Delta n$, between each bonded pair of PGNs in the array at the maximal $(\gamma ={{\gamma }_{0}})$ and zero $\left( \gamma =0 \right)$ shear strains. The variations $\Delta n$ are shown for three values of $\omega /{{\omega }_{0}}$ (figures 7(a)–(c)). The respective frequencies $\omega$ in figures 7(a)–(c) are the same as those in figures 2(e), (f) and (h), which show the Lissajous curves for the number of the labile bonds per particle, $(N_{b}^{(l)}/{{N}_{{\rm PGN}}})$. At the lowest frequency $\omega /{{\omega }_{0}}\approx 0.66$ (figure 7(a)), the oscillatory shear causes the greatest variation in the number of labile bonds between individual pairs of PGNs. The variation in $\Delta n$ decreases with an increase in the frequency of oscillation, so that at the highest frequency $\omega /{{\omega }_{0}}\approx 660$, only a few pairs of PGNs exhibit a change in the number of labile bonds in the course of deformation from $\gamma =0$ to $\gamma ={{\gamma }_{0}}$ (figure 7(c)).

Zoom In Zoom Out Reset image size

As noted with respect to the discussion of figure 6(b), the dissipation in the system (value of the loss modulus) reaches a maximum value when $\omega /{{\omega }_{0}}\tilde{\ }1$. The rupture and formation of the labile bonds do not contribute to the loss modulus at both the lower and higher frequencies (see the solid lines in figure 6(b)), although the variations $\Delta n$ are drastically different at these two extremes (compare figures 7(a) and (c)). At the lower frequency, the labile bonds break and reform faster than the progression of the shear deformation. Therefore, no hysteresis is observed in the number of bonds in the system, as evident from figure 2(e) (see curve 2), and hence, there is no contribution to dissipation. In contrast, at the higher shear frequency, the rate of rupture of the labile bonds is much lower than the shear rate. The number of labile bonds does not change noticeably in the course of an oscillation cycle (see curve 2 in figure 2(h)), so again there is no contribution to dissipation. The rupture and reformation of the labile bonds contribute to the loss modulus when the bond rupture rate and shear rate are comparable, so that the number of labile bonds exhibits hysteretic behavior (see loops in curve 2 in figure 2(f)).

As demonstrated above, the viscoelastic properties of the PGN networks that are caused by rupture and formation of the labile bonds are governed by the bond rupture rate constant, $k_{0r}^{\left( l \right)}$, which depends on the bond energy. The viscoelastic behavior of the PGN networks cannot, however, be considered as a result of a single relaxation process characterized by the relaxation time $\tau$ ${\mkern 1mu} \tilde{\ }{\mkern 1mu} 1/k_{0r}^{\left( l \right)}$. This is evident from figure 8, which shows Cole−Cole plots for the PGN networks with labile bonds of energy $U_{0}^{(l)}=33{{k}_{B}}T$ (plot 1) and $37{{k}_{B}}T$ (plot 2). In a Cole–Cole plot, the properly normalized frequency-dependent loss and storage moduli are plotted against each other. The normalization is performed as indicated in figure 8, where ${{G}_{0}}$ is the low frequency plateau value of the storage modulus, and $\Delta G$ is the difference between the high and low frequency plateau values of $G^{\prime}$. It is known that if viscoelastic relaxation is characterized by a single relaxation time, then the Cole–Cole plot forms a semi-circle indicated by curve 3 (dashed line) in figure 8 [19, 20]. In contrast, figure 8 shows that the Cole–Cole plots for the viscoelastic moduli obtained from both the results of simulations (symbols) and the three-particle model (solid curves) lie well below the semicircle. The only exceptions are the symbols on the right that correspond to the simulations at the higher frequency of shear oscillations where the particle-solvent friction dominates (see above).

Zoom In Zoom Out Reset image size

The features of viscoelastic behavior revealed through the Cole–Cole plots (figure 8) are characteristic of rheologically complex materials, for which the stress relaxation exhibits a stretched exponential dependence on time and hence, the spectrum of relaxation times is continuous [19, 20]. The viscoelasticity of rheologically complex materials could be described in terms of the customary linear springs and dashpots only if an infinite, hierarchical system of springs and dashpots is considered [19, 20]. It is rather remarkable that the simple three-particle model, which takes into account only the breakage and formation of bonds (see figure 5 and equation (8)), is capable of reproducing the characteristic features of the rheologically complex behavior of the PGN network as seen in figure 8. Hence, we can conclude that the relaxation processes in this system are controlled mostly by the kinetics of bond rupture and formation. In order to demonstrate how the bond kinetics results in a spectrum of relaxation times, in figures 9(a)–(c) we plot the evolution of the average number of bonds, $\left\langle n \right\rangle$, obtained for the three-particle model through solving equation (8) at the different values of $\omega /{{\omega }_{0}}$ indicated in the figures. The lines labelled by ${{n}_{13}}$ and ${{n}_{23}}$ correspond to the average number of bonds between particles 1 and 3 and particles 2 and 3, respectively (see figure 5 and figures 9(a)–(c)). The values of ${{n}_{13}}$ and ${{n}_{23}}$ exhibit step-wise variations at the frequencies of $\omega /{{\omega }_{0}}\approx 0.66$ and 6.6 (see figures 9(a), (b), respectively), and saw-tooth kinetics is observed at the higher frequency of $\omega /{{\omega }_{0}}\approx 66$ (figure 9(c)). Such a highly nonharmonic response of the average numbers of bonds ${{n}_{13}}$ and ${{n}_{23}}$ to the sinusoidal shear deformations cannot be represented by a combination of a few Fourier harmonics. As the attractive force between two particles is proportional to the average number of bonds, the stress relaxation exhibits a spectrum of relaxation times, which is a characteristic feature of the rheologically complex behavior.

Zoom In Zoom Out Reset image size

It is important to note that the average numbers of bonds ${{n}_{13}}$ and ${{n}_{23}}$ shown in figures 9(a)–(c) are obtained by solving the master equation, equation (8), and in this way, we account for stochastic effects. The contribution of the stochastic effects can be estimated by calculating the value of relative deviation $\varepsilon ={{(\langle {{n}^{2}}\rangle /{{\langle n\rangle }^{2}}-1)}^{1/2}}$; the stochastic effects can be neglected at $\varepsilon \ll 1$. Figures 9(d)–(f) show the time-dependent values of $\varepsilon$ calculated for the pairs of particles 1–3 and 2–3 and denoted as ${{\varepsilon }_{13}}$ and ${{\varepsilon }_{23}}$, respectively, at the same frequencies of shear oscillation as in figures 9(a)–(c). The stochastic effects cannot be neglected as $\varepsilon \leqslant 0.1$ at the frequencies of $\omega /{{\omega }_{0}}\approx 0.66$ and 66 (figures 9(d) and (f), respectively), and $\varepsilon \leqslant 0.2$ at the intermediate frequency $\omega /{{\omega }_{0}}\approx 6.6$ (figure 9(e)). It is worth noting that the frequency $\omega /{{\omega }_{0}}\approx 6.6$, at which the stochastic effects are especially important, corresponds to the peak in the loss modulus (see figure 6(b)).

The response of the PGN networks to oscillatory shear deformations is nonlinear under sufficiently high shear amplitudes. In the nonlinear regime, the force–strain Lissajous curves deviate from the elliptical shape. Figure 10 shows Lissajous plots obtained from the simulations of the PGN network at shear frequencies of $\omega /{{\omega }_{0}}\approx 0.3$ (figures 10(a)–(c)) and 3 (figures 10(d)–(f)) and shear amplitudes from ${{\gamma }_{0}}=0.15$ to 0.45. Here, we fixed the labile bond energy at $U_{0}^{(l)}=33{{k}_{B}}T$. The shape of the Lissajous curves changes dramatically as the shear amplitude ${{\gamma }_{0}}$ is increased, and that the shape distortion is more pronounced at the lower relative frequency of oscillation $\omega /{{\omega }_{0}}$. Furthermore, the shear force (marked in gray in figure 10) exhibits noticeable fluctuations, which increase with an increase in the amplitude of shear ${{\gamma }_{0}}$.

Zoom In Zoom Out Reset image size

The nonlinear response of the material to the large amplitude oscillatory shear (LAOS), $\gamma (t)={{\gamma }_{0}}{\rm sin} \;(\omega \;t)$, can be analyzed by representing the time dependent shear force $F(t)$ as a Fourier series, i.e.,

$$\begin{eqnarray}&&F(t)={{\gamma }_{0}}\mathop \sum \limits_{n\,:\,{\rm odd}} \left\{ G_{n}^{^{\prime} }\left( \omega ,{{\gamma }_{0}} \right){\rm sin} \left( n\;\omega \;t \right)+G_{n}^{^{\prime\prime} }\left( \omega ,{{\gamma }_{0}} \right){\rm cos} \left( n\;\omega \;t \right) \right\}.\end{eqnarray} \tag{ 9 }$$

Here, $G_{n}^{^{\prime} }$ and $G_{n}^{^{\prime\prime} }$ are the Fourier coefficients, which depend on the shear frequency $\omega$ and amplitude ${{\gamma }_{0}}$. It is more useful, however, to consider the force $F$ as a function of the instantaneous values of shear strain $\gamma$ and shear strain rate $\dot{\gamma }$, and then decompose $F$ into the elastic part ${{F}_{e}}(\gamma )$, which depends only on the strain, and the strain rate-dependent viscous contribution ${{F}_{v}}(\dot{\gamma })$:

$$\begin{eqnarray}&&F(\gamma ,\dot{\gamma })={{F}_{e}}(\gamma )+{{F}_{v}}(\dot{\gamma }).\end{eqnarray} \tag{ 10 }$$

The above decomposition can be performed in a straightforward manner if the elastic and viscous contributions are represented by Chebyshev polynomials of the first kind [21, 22]:

$$\begin{eqnarray}&&{{F}_{e}}\left( \gamma \right)={{\gamma }_{0}}\mathop \sum \limits_{n\,:\,{\rm odd}} {{e}_{n}}(\omega ,{{\gamma }_{0}})\;{{T}_{n}}(\gamma /{{\gamma }_{0}}),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{F}_{v}}\left( {\dot{\gamma }} \right)={{\dot{\gamma }}_{0}}\mathop \sum \limits_{n\,:\,{\rm odd}} {{v}_{n}}(\omega ,{{\gamma }_{0}})\;{{T}_{n}}(\dot{\gamma }/{{\dot{\gamma }}_{0}}),\end{eqnarray} \tag{ 12 }$$

where ${{\dot{\gamma }}_{0}}=\omega \;{{\gamma }_{0}}$ is the strain rate amplitude. The Chebyshev polynomials ${{T}_{n}}(x)$ are defined by the following recursive relation [23]:

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{T}_{0}}(x) & = & 1, \\ {{T}_{1}}(x) & = & x, \\ {{T}_{n+1}}(x) & = & 2x{{T}_{n}}(x)-{{T}_{n-1}}(x). \\ \end{array}\end{eqnarray} \tag{ 13 }$$

The first harmonic coefficients in the Fourier series in equation (9) provide a measure of the average (over one cycle of oscillation) values of the storage and loss moduli $G^{\prime}$ and $G^{\prime}$, respectively, which are related to the first Chebyshev coefficients in equations (11) and (12) as follows:

$$\begin{eqnarray}&&G^{\prime} \equiv G_{1}^{^{\prime} }(\omega ,{{\gamma }_{0}})={{e}_{1}}(\omega ,{{\gamma }_{0}}),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&G^{\prime\prime} \equiv G_{1}^{^{\prime\prime} }(\omega ,{{\gamma }_{0}})=\omega \;{{v}_{1}}(\omega ,{{\gamma }_{0}}).\end{eqnarray} \tag{ 15 }$$

To determine the average storage and loss moduli from the results of the simulations, we employ equations (10)–(12), in which the first four terms of the Chebyshev expansion (n = 1, 3, 5 and 7) are used to fit the shear force over 20 consecutive cycles of oscillation. Fitting was performed using the Mathematica™ software. The Lissajous curves obtained by the fitting procedure are shown in figure 10 by the thick colored lines. Finally, the values of $G^{\prime}$ and $G^{\prime\prime}$ are calculated according to equations (14) and (15) followed by averaging over eight independent runs.

Figures 11(a)–(d) show the Lissajous curves obtained from the Chebyshev analysis of the shear force at various values of the strain amplitude ${{\gamma }_{0}}$ and scaled frequency $\omega /{{\omega }_{0}}$. The Lissajous curves showing the average number of labile bonds versus the shear strain in the respective samples are presented in figures 11(e)–(h). The plots for $\omega /{{\omega }_{0}}\approx 0.3$ and 3 are obtained from the results of simulations at the energy of labile bonds $U_{0}^{(l)}=33{{k}_{B}}T$, and the plots for $\omega /{{\omega }_{0}}\approx 3.3$ and 33 correspond to $U_{0}^{(l)}=37{{k}_{B}}T$. The Lissajous curves in figures 11(a)–(d) clearly show that the nonlinearity of the force response to the oscillatory shear deformations (indicated by deviation of the curves from elliptical shape) increase with the strain amplitude, and are more pronounced at the lower values of $\omega /{{\omega }_{0}}$. Notably, a strain softening is observed at low frequencies, i.e., the shear force decreases with an increase in shear strain at $\gamma >0.2$ (see figures 11(a) and (b)).

Zoom In Zoom Out Reset image size

The number of labile bonds in the system also depends on both the amplitude ${{\gamma }_{0}}$ and relative frequency $\omega /{{\omega }_{0}}$ of oscillations (figures 11(e)–(h)). At $\omega <{{\omega }_{0}}$, during a cycle of shear oscillation, the variation in the number of labile bonds increases with the shear strain amplitude ${{\gamma }_{0}}$ (figure 11(e)). At a sufficiently high frequency of oscillation, an increase in ${{\gamma }_{0}}$ leads to an increase in the number of labile bonds, but the number of bonds shows very little variation over an oscillation cycle (figure 11(h)). At $\omega \tilde{\ }{{\omega }_{0}}$, $N_{b}^{(l)}$ exhibits hysteretic behavior, as indicated by the looped Lissajous curves in figures 11(f) and (g), and both the average number of bonds and the range of variation increase with ${{\gamma }_{0}}$.

The average storage and loss moduli extracted from the results of simulations in the LAOS regime are shown in figures 12(a) and (b), respectively. Comparison of figures 3 and 12 reveals that both $G^{\prime}$ and $G^{\prime\prime}$ for the PGN networks exhibit very similar behavior as functions of $\omega /{{\omega }_{0}}$ in the linear and nonlinear regimes. (Recall that ${{\omega }_{0}}$ depends on the energy of the labile bonds.) An increase in the shear amplitude is seen to affect mostly the average storage modulus, as $G^{\prime}$ decreases consistently with an increase in the shear amplitude (figure 12(a)). The decrease in $G^{\prime}$ is attributed to the strain softening, which occurs at the shear strains $\gamma >0.2$ as discussed above. In contrast, variations in the amplitude of the strain produce very little effect on the dissipative properties of the PGN networks characterized by the average loss modulus $G^{\prime\prime}$, as a comparison of figures 3(d) and 12(b) shows.

Zoom In Zoom Out Reset image size

To further demonstrate the similarity in the viscoelastic response of the materials in the linear and nonlinear regimes, we compare the results of simulations for the LAOS deformations with the calculations obtained from the three-particle model (figure 5), which accurately described the linear case (figures 6(a), (b)). Figures 6(c)–(f) present the normalized average storage and loss moduli determined from the simulations (symbols) and from the model calculations (solid lines) at the strain amplitudes of ${{\gamma }_{0}}$ = 0.15 and 0.45. At ${{\gamma }_{0}}$ = 0.15, $G^{\prime}$ obtained from the three-particle model is in a good agreement with the data from the simulations of the network (figure 6(c)), and the agreement is fairly good at the larger strain amplitude ${{\gamma }_{0}}$ = 0.45 (figure 6(e)). Furthermore, similar to the linear regime, the three-particle calculations reproduce the peak observed in the loss modulus $G^{\prime\prime}$ for the PGN network containing the stronger labile bonds (compare curves 2 in figure 6(b) and in figures 6(d), (f)). Therefore, we can conclude that the viscoelastic relaxation in the LAOS regime is governed by the same processes as in the linear case.

We now utilize the three-particle model (figure 5) to discuss the effect of strain softening observed in the LAOS regime (see figure 11(a)). Figure 13 indicates that the strain softening effect can be explained through geometric arguments. Figures 13(a) and (b) show schematically the forces acting on particle 3 in the system of three PGNs at the shear strains $\gamma =0$ and 0.45, respectively, with ${\bf F}_{{\rm att}}^{p}$ and ${\bf F}_{{\rm rep}}^{p}$ being the respective attractive and repulsive force within the pair $p$, $p=(13)$ or (23). At $\gamma =0$, the repulsive forces acting on particle 3 within pairs 1–3 and 2–3 compensate each other, whereas the sum of the attractive forces is small but non-zero because ${{n}_{13}}\ne {{n}_{23}}$ due to the hysteretic behavior of the average number of bonds (figure 13(a)). With the shear deformation, the individual forces change their magnitude and direction. As seen in figure 13(b), the repulsive and attractive forces in pair 1–3 become more oriented towards the direction of shear (the horizontal direction) but decrease in magnitude. In contrast, the forces in pair 2–3 increase in magnitude, but become more oriented in the normal to the shear direction (figure 13(b)). As a result, the total force in the direction of the applied shear is small at $\gamma =0.45$ as seen in figure 13(b), i.e., strain softening takes place. The shear, ${{F}_{s}}$, and normal, ${{F}_{n}}$, components of the force as functions of the oscillatory shear strain (the Lissajous curves) obtained for the three-particle model at $\omega /{{\omega }_{0}}\approx 0.3$ and ${{\gamma }_{0}}$ = 0.45 are shown in figure 13(c). At $\gamma >0.2$, the magnitude of the shear force ${{F}_{s}}$ decreases with an increase of shear strain (curve 1), whereas the magnitude of the normal force ${{F}_{n}}$ continues to increase (curve 2). Notably, the three-particle model reproduces the value of shear strain at the onset of shear softening $\left( \gamma \approx 0.2 \right)$ observed in the simulations of the crosslinked array of PGNs (as can be seen by comparing curve 3 in figure 11(a) and curve 1 in figure 13(c)).

Zoom In Zoom Out Reset image size

Finally, it is important to note that the PGN networks could exhibit noticeable structural changes in the LAOS regime depending on the frequency of shear strain oscillation $\omega$ relative to the characteristic frequency ${{\omega }_{0}}=2\pi k_{0r}^{\left( l \right)}$, which, in turn, depends on the energy of the labile bonds. Figure 14 shows snapshots of the simulations at the higher strain amplitude ${{\gamma }_{0}}=0.45$ that correspond to the moments of time when the shear strain is $\gamma$ = 0 and 0.45. Only labile bonds between the particles are shown in figure 14, and the number of bonds in a link is indicated by the thickness of the line. Thus, the thicker lines indicate more labile bonds in the link. Figures 14(a) and (b) show the configuration of the labile bonds at the lower relative frequency $\omega /{{\omega }_{0}}\approx 0.3$ for the shear strain of $\gamma$ = 0 and 0.45, respectively; figures 14(a) and (b) were obtained at $U_{0}^{(l)}=33{{k}_{B}}T$. Figure 14(b) shows that at $\gamma$ = 0.45, the distribution of the labile bonds is non-uniform in the system, i.e., some PGN pairs have many labile bonds between them (thick lines) and some have none (voids). Nevertheless, when $\gamma =0$, the particles return to the hexagonal arrangement (compare figures 1 and 14(a)) because the labile bonds rupture and reform sufficiently rapidly at the given frequency of oscillation $\omega <{{\omega }_{0}}$.

Zoom In Zoom Out Reset image size

The PGN network exhibits quite different structural features in the opposite case of $\omega >{{\omega }_{0}}$, as shown in figures 14(c) and (d) for the configuration of the labile bonds at $\omega /{{\omega }_{0}}\approx 33$; $U_{0}^{(l)}=37{{k}_{B}}T$ in figures 14(c) and (d). Since the newly formed labile bonds do not have enough time to rupture during one cycle of oscillation, the particles start aggregating owing to the increased attractive forces, and the latter leads to the formation of cavities after several cycles of deformation (figures 14(c) and (d)).

Note that in figure 14, the opposite cases of lower ($\omega /{{\omega }_{0}}<1$ in figures 14(a) and (b)) and higher ($\omega /{{\omega }_{0}}>1$ in figures 14(c) and (d)) relative frequency were obtained by varying both the frequency $\omega$ and the energy of the labile bonds $U_{0}^{(l)}$. At a given value of $\omega /{{\omega }_{0}}$, the configurations of labile bonds look similar in the cases of the weaker and stronger labile bonds (not shown here but see the discussion below).

In the LAOS regime, the local distribution of shear strain within a deformed PGN network becomes spatially non-uniform if the relative frequency $\omega /{{\omega }_{0}}$ is sufficiently low. The value of shear strain in the vicinity of a particle is calculated by averaging the shear strains determined in the two triangular elements to which the particle belongs in the undeformed state. In figure 15(a), particle 1 belongs to the elements (123) and (145) formed by the neighboring particles 1, 2, 3 and 1, 4, 5, respectively. The shear strains in the two elements are determined as ${{\gamma }_{(123)}}=({{x}_{2}}+{{x}_{3}}-2{{x}_{1}}){{({{y}_{2}}+{{y}_{3}}-2{{y}_{1}})}^{-1}}$ and ${{\gamma }_{(145)}}=(2{{x}_{1}}-{{x}_{4}}-{{x}_{5}}){{(2{{y}_{1}}-{{y}_{4}}-{{y}_{5}})}^{-1}}$, where $\{{{x}_{n}},{{y}_{n}}\}$ are the coordinates of the particle $n$ $\left( n=1,2,\ldots ,5 \right)$. Then, the shear strain localized at the particle 1 is calculated as ${{\gamma }_{1}}=\;({{\gamma }_{(123)}}+{{\gamma }_{(145)}})/2$. Figure 15(b) shows the histograms of the local shear strain in the PGN network having the labile bonds of the energy $U_{0}^{(l)}=33{{k}_{B}}T$ obtained in the LAOS regime at the relative frequencies of $\omega /{{\omega }_{0}}\approx 3$ (top histogram) and 0.3 (bottom histogram) and the strain amplitude ${{\gamma }_{0}}=0.45$. The histograms correspond to the moments of time when the shear strain is maximal $\left( \gamma =0.45 \right)$. At the higher frequency of $\omega /{{\omega }_{0}}\approx 3$, the distribution of shear strain is bell-shaped and centered around the value of the imposed strain $\gamma =0.45$ (figure 15(b)). In contrast, at the lower frequency of $\omega /{{\omega }_{0}}\approx 0.3$, the distribution of strain is rather featureless but remarkably wide, revealing the presence of areas where the local shear strain is lower $\left( \gamma \approx 0.28{\mkern 1mu} <{\mkern 1mu} 0.45 \right)$ and higher $\left( \gamma \approx 0.52{\mkern 1mu} >{\mkern 1mu} 0.45 \right)$ than the imposed strain (figure 15(b)). The observed non-uniformity of the local strains is evidence of local rearrangements of particles that occur because the labile bonds between the neighboring PGNs break in the course of deformation if the deformation is sufficiently slow (low frequency). If the deformation is fast (high frequency), the breakage of labile bonds during deformation is less prominent so the shear strain is distributed more uniformly within the sample. Note that this is the relative frequency of oscillations $\omega /{{\omega }_{0}}$ that controls the effect. It is seen in figures 15(c) and (d) that the distributions of strain in the samples having weaker $(U_{0}^{(l)}=33{{k}_{B}}T)$ and stronger $(U_{0}^{(l)}=37{{k}_{B}}T)$ labile bonds are quite similar if the relative frequency of oscillations is the same ($\omega /{{\omega }_{0}}\approx 3$ and 33 in figures 15(c) and (d), respectively).

Zoom In Zoom Out Reset image size