Interfacial Interaction Enhanced Rheological Behavior in PAM/CTAC/Salt Aqueous Solution—A Coarse-Grained Molecular Dynamics Study

Interfacial interactions within a multi-phase polymer solution play critical roles in processing control and mass transportation in chemical engineering. However, the understandings of these roles remain unexplored due to the complexity of the system. In this study, we used an efficient analytical method—a nonequilibrium molecular dynamics (NEMD) simulation—to unveil the molecular interactions and rheology of a multiphase solution containing cetyltrimethyl ammonium chloride (CTAC), polyacrylamide (PAM), and sodium salicylate (NaSal). The associated macroscopic rheological characteristics and shear viscosity of the polymer/surfactant solution were investigated, where the computational results agreed well with the experimental data. The relation between the characteristic time and shear rate was consistent with the power law. By simulating the shear viscosity of the polymer/surfactant solution, we found that the phase transition of micelles within the mixture led to a non-monotonic increase in the viscosity of the mixed solution with the increase in concentration of CTAC or PAM. We expect this optimized molecular dynamic approach to advance the current understanding on chemical–physical interactions within polymer/surfactant mixtures at the molecular level and enable emerging engineering solutions.


Bonded Interaction and Nonbonded Interaction
Bond interaction between two neighboring beads is determined by a harmonic potential [45] 2 0 where bond K and 0 r represent the spring constant and the equilibrium distance, respectively. For CTA + , the spring constant is 1250 KJ·mol −1 ·nm −2 and the equilibrium distance is 0.47 nm. The equilibrium distance for the ring of Saland PAM are 0.27 nm and 0.28 nm, respectively. The stiffness of a chain is described by a weak harmonic potential [45]: where ij σ represents the effective distance of approach between two beads and ij ε is the depth of the potential well. The potential is truncated and shifted at a cutoff distance =1.2 nm r , which can be found in many studies [50][51][52]. ij σ and ij ε are defined by the MARTINI force field [45].
The electrostatic interaction between the charged beads i and j is defined as [45]: 0 ( ) 4 i j el ij r ij q q U r r πε ε = (4) where i q and j q are the charges; 0 ε and r ε are the permittivity of the vacuum and the relative dielectric constant, respectively; and r ε is 15 [45].

Bonded Interaction and Nonbonded Interaction
Bond interaction between two neighboring beads is determined by a harmonic potential [45] U bond (r ij ) = where K bond and r 0 represent the spring constant and the equilibrium distance, respectively. For CTA + , the spring constant is 1250 KJ·mol −1 ·nm −2 and the equilibrium distance is 0.47 nm. The equilibrium distance for the ring of Sal − and PAM are 0.27 nm and 0.28 nm, respectively. The stiffness of a chain is described by a weak harmonic potential [45]: U angle (θ) = 1 2 K angle (cos θ − cos θ 0 ) 2 (2) where K angle represents the force constant and θ 0 represents the equilibrium bond angle. The force constant and equilibrium bond angle are K angle = 25 KJ·mol −1 and θ 0 = 180 • for aliphatic chains. The nonbonded sites i and j interact through a Lennard-Jones potential [45]: where σ ij represents the effective distance of approach between two beads and ε ij is the depth of the potential well. The potential is truncated and shifted at a cutoff distance r = 1.2 nm, which can be found in many studies [50][51][52]. σ ij and ε ij are defined by the MARTINI force field [45]. The electrostatic interaction between the charged beads i and j is defined as [45]: U el (r ij ) = q i q j 4πε 0 ε r r ij (4) where q i and q j are the charges; ε 0 and ε r are the permittivity of the vacuum and the relative dielectric constant, respectively; and ε r is 15 [45].

Simulations
The Large-scale Atomic/Molecular Massively Parallel Simulator(LAMMPS) package was used to calculate the shear viscosity of the salt-added polymer and surfactant solutions. PACKMOL software was used to obtain a random distribution of surfactants, salt, polymers, and water particles in a simulation cube with dimensions of 23.1 nm × 23.1 nm × 23.1 nm. The van der Waals Polymers 2020, 12, 265 4 of 17 interaction shifted from 0.9 to 1.2 nm, and the electrostatic interaction shifted from 0.0 nm to 1.2 nm. The particle-particle particle-mesh (PPPM) solver was used to calculate the long-range electrostatic interactions. Minimization, isothermal-isobaric ensemble(NPT) simulation, and Canonical ensemble(NVT) simulation were carried out at a temperature of 300 K via a Nose-Hoover thermostat and under a pressure of 1 bar prior to calculating the shear viscosity. The timestep was 20 fs and the NVT simulation lasted for 600 ns to reach equilibrium. The simulation time was sufficient for the system to get to an equilibrium state because the potential energy will not change after 400 ns. The potential energy curve as a function of time can be found in Figures S1 and S2 in the Supplementary Materials (SM). The shear viscosity was calculated from NEMD simulation using the equilibrium data by deforming the box. This means that we needed to obtain the equilibrium structure of the simulation system before calculating the viscosity. Then, the viscosity was calculated by a non-equilibrium MD (NEMD) simulation by changing the shape of the simulation box. More simulation details can be found in the Supplementary Materials.

Results and Discussion
In this study, the NEMD simulation was employed to understand the interactions in the PAM and CTAC mixture, and the associate impact on the rheology behavior of the PAM and CTAC containing aqueous solutions with added salt (PAM/CTAC/Salt multiphase aqueous solutions). In our simulation, there were two series of PAM/CTAC/Salt multiphase aqueous solutions. One was the PAM based solution, which had a given concentration of PAM, but variable concentrations of CTAC from 0.1 to 0.3 mol·L −1 . Another was the CTAC based solution with a fixed CTAC concentration, but variable concentrations of PAM from 2.96 × 10 −3 to 2.02 × 10 −2 mol·L −1 . NaSal is used as a counter-ion salt to form stable micelle structures with CTAC. In this work, the mole ratio of NaSal to CTAC (R) was 0.8 in both PAM based solutions and CTAC based solutions. Surfactants can self-assembly into micelles with various shapes such as spherical, wormlike, and branched micelles. As mentioned in the introduction, surfactants are efficient in drag reduction. The reason why surfactants can be used as drag reducing additives is that surfactant micelles can form shear induced structures (SIS). When the shear induced structures are destroyed, the solution will lose its drag reduction efficiency. The change in micelle structures is also called shear induced phase transition (SPS) [53,54]. The SPS phenomenon also leads to the shear thickening phenomenon in rheology behavior, which means that the evolution of structures can influence the rheology behavior during the shear process. In our study, changes in the concentrations of CTAC and/or PAM in the mixed solutions could influence the initial structures of the mixtures and further influence the rheology behavior of the solutions. Thus, the structure evolution (phase transition) and the rheology behavior can be studied.

Study of Shear Viscosity for PAM/CTAC/Salt Multiphase Aqueous System
The shear viscosities of the PAM/CTAC/Salt multiphase aqueous were analyzed at shearing rates ranging from 7.79 × 10 8 to 1.73 × 10 11 s −1 . The viscosity results of the PAM based solutions in Figure 2a show an initial plateau under low shear rates (from 7.79 × 10 8 to 1.3 × 10 9 s −1 ). When the shear rates increase from 1.3 × 10 9 to 1.73 × 10 11 s −1 , all samples showed a typical shear thinning phenomenon. Thus, we think that the point when the shear rate is 1.3 × 10 9 s −1 is considered to be the turning point; before this point is the plateau zone, and after this point is the shear thinning zone. The steady state shear viscosity data agreed well with the reported experimental results of the pure surfactant solutions and pure polymer solutions because both results had the plateau and the shear thinning zone [55][56][57]. The shear viscosity data of the CTAC based solutions in Figure 2b also show an initial plateau at low shear rates and shear thinning in the high shear rates region. The above results indicate that even with various concentrations of CTAC or PAM, the trend of the change in the shear viscosity was similar in the PAM/CTAC/Salt multiphase aqueous system. We have to point out that the shear rates used in our simulation were much larger than those applied in the experiments. This is because the calculation area is small in molecular simulation, and the shear rates applied in the experiments were not large enough to access reliable data in the molecular simulations. If the shear rate is not enough, the calculation noise will be very large and the statistical data cannot be accepted because the viscosity value will become a negative number, which has no physical significance. It can be found in the published studies [47,50,52,58] that the shear rates in the molecular simulation were far above the experimental data, but the simulation results shown qualitatively agreed with experimental results because the simulation curves and the experimental curves had similar trends (both the plateau and the shear thinning zone). However, if the shear rate is too large, the random motion of the sites will also lead to no physical significance viscosity value because the data cannot be well analyzed. Thus, the reliable shear rates have a range for each simulation case, for our simulation, the reliable range was from 7.79 × 10 8 to 1.73 × 10 11 s −1 and the large shear rates make sense. The corresponding inner structures can give us valuable information in understanding the rheology behavior of polymer/surfactant solutions.

Medium and High Shear Rates
The viscosity values of the PAM/CTAC/Salt solutions, as a function of simulation time, are shown in Figure 3a. The data of the PAM based solutions revealed a similar trend at a medium shear rate (1.3× 10 10 s −1 ) and a high (8.65× 10 10 s −1 ) shear rate. Herein, the data can be regarded as an equilibrium value when the value of the viscosity fluctuates ≤2% and last for a simulation time of ≥2 ns. It can be observed in Figure 3a that there was a set of data points at the beginning of the simulation, followed by a steady region at a later stage of simulation. The viscosity tended to reach its stable region earlier at the high shear rate than under the medium shear rate. More information can be found in Supplementary Materials Figure S1. We have to point out that the shear rates used in our simulation were much larger than those applied in the experiments. This is because the calculation area is small in molecular simulation, and the shear rates applied in the experiments were not large enough to access reliable data in the molecular simulations. If the shear rate is not enough, the calculation noise will be very large and the statistical data cannot be accepted because the viscosity value will become a negative number, which has no physical significance. It can be found in the published studies [47,50,52,58] that the shear rates in the molecular simulation were far above the experimental data, but the simulation results shown qualitatively agreed with experimental results because the simulation curves and the experimental curves had similar trends (both the plateau and the shear thinning zone). However, if the shear rate is too large, the random motion of the sites will also lead to no physical significance viscosity value because the data cannot be well analyzed. Thus, the reliable shear rates have a range for each simulation case, for our simulation, the reliable range was from 7.79 × 10 8 to 1.73 × 10 11 s −1 and the large shear rates make sense. The corresponding inner structures can give us valuable information in understanding the rheology behavior of polymer/surfactant solutions.

Medium and High Shear Rates
The viscosity values of the PAM/CTAC/Salt solutions, as a function of simulation time, are shown in Figure 3a. The data of the PAM based solutions revealed a similar trend at a medium shear rate (1.3 × 10 10 s −1 ) and a high (8.65 × 10 10 s −1 ) shear rate. Herein, the data can be regarded as an equilibrium value when the value of the viscosity fluctuates ≤2% and last for a simulation time of ≥2 ns. It can be observed in Figure 3a that there was a set of data points at the beginning of the simulation, followed by a steady region at a later stage of simulation. The viscosity tended to reach its stable region earlier at the high shear rate than under the medium shear rate. More information can be found in Supplementary Materials Figure S1. Figure 3b,c show the structure of the mixture under medium shear rates (1.3× 10 10 s −1 in Figure 3a) and high shear rates (8.65 × 10 10 s −1 in Figure 3a), respectively. It can be seen that the CTAC micelles and PAM polymer chains orientate along the flow direction under the shear force and form a stretched band at medium and high shear rates. The alignments of the mixture along the flow direction could   Figure  3a) and high shear rates (8.65 × 10 10 s −1 in Figure 3a), respectively. It can be seen that the CTAC micelles and PAM polymer chains orientate along the flow direction under the shear force and form a stretched band at medium and high shear rates. The alignments of the mixture along the flow direction could lead to the reduction in viscosity of the mixture solution, and the saturated viscosity value whilst the phase interactions became stable.

Low Shear Rates for PAM Based Solutions
The volume fraction of CTAC was larger than that of PAM in the PAM based solutions, which means that the rheological behavior of the mixture will be significantly influenced by surfactants, especially under low shear rates. As shown in Figure 4a, typical "two peaks" curves were obtained with a CTAC concentration of 0.3 mol·L −1 (more curves are available in the Supplementary Materials Figure S2), which has been previously reported in many surfactant solution experiments [59,60]. After the two peaks, the viscosity curve finally reached a steady state after 9 ns, and the plateau lasted until 11 ns, which can be easily noted in Figure 4a. Next, the time dependent viscosities of the PAM based solutions under low shear rates are shown in the small figure in Figure 4a, and the time dependent viscosities of the PAM based solutions under low shear rates are shown in Figure 4a. Figure 4b shows a multi-stage evolution for the viscosity and is a theoretical construct summarization of Figure 4a. Color scheme: purple, hydrophilic part of CTAC; cyan and blue, hydrophobic tail of CTAC; magenta, aromatic ring of Sal -; khaki, charge group of Sal − ; yellow, Na + ; dark green, Cl − ; red, the acidammide group of PAM; green, carbon backbone of PAM. Water sites (reddish) and antifreeze water sites (grey) were hidden for clear display. The applied flow direction is indicated by a red arrow in (b) and (c).

Low Shear Rates for PAM Based Solutions
The volume fraction of CTAC was larger than that of PAM in the PAM based solutions, which means that the rheological behavior of the mixture will be significantly influenced by surfactants, especially under low shear rates. As shown in Figure 4a, typical "two peaks" curves were obtained with a CTAC concentration of 0.3 mol·L −1 (more curves are available in the Supplementary Materials Figure S2), which has been previously reported in many surfactant solution experiments [59,60]. After the two peaks, the viscosity curve finally reached a steady state after 9 ns, and the plateau lasted until 11 ns, which can be easily noted in Figure 4a The multi-stage evolution can be approximately divided into five stages: oscillating shear-thickening (Stage I), shear-thinning (Stage II), oscillating sharp shear-thickening (Stage III), oscillating adjustment (Stage IV), and oscillating plateau (Stage V).  To further understand the relationship between the structural evolution and the changes in shear viscosity, the continuous phase changes for the mixed solution in Figure 4a were analyzed at a shear rate of 1.08 × 10 9 s −1 . A series of snapshots (Figures 5a-g) were taken from the CGMD simulation to visualize the time dependent multi-stage phase evolution. In Figure 4b, from a to b, there was an increase in viscosity, but the increase only lasted for a short time. From Stage I in Figures 5a,b, we can see that the structure of the mixture had no significant changes, apart from the slight stretching under shearing (top view). The oscillating increase in viscosity at Stage I can be attributed to the increase of viscosity due to the extension of the PAM molecules [61,62] and the elastic response of the resistance of disentanglement action of surfactant micelles [59]. At Stage II, the viscosity decreased with time and the stage lasted for about 0.8 ns (Figure 4b, from b to c). From the corresponding snapshots (Figures 5b,c from the left view), we observed that the structures were slightly stretched, and the structures became smaller, indicating the alignment of structures with flow direction, which is likely to induce a slight decrease in the viscosity [52].  To further understand the relationship between the structural evolution and the changes in shear viscosity, the continuous phase changes for the mixed solution in Figure 4a were analyzed at a shear rate of 1.08 × 10 9 s −1 . A series of snapshots (Figure 5a-g) were taken from the CGMD simulation to visualize the time dependent multi-stage phase evolution. In Figure 4b, from a to b, there was an increase in viscosity, but the increase only lasted for a short time. From Stage I in Figure 5a,b, we can see that the structure of the mixture had no significant changes, apart from the slight stretching under shearing (top view). The oscillating increase in viscosity at Stage I can be attributed to the increase of viscosity due to the extension of the PAM molecules [61,62] and the elastic response of the resistance of disentanglement action of surfactant micelles [59]. At Stage II, the viscosity decreased with time and the stage lasted for about 0.8 ns (Figure 4b, from b to c). From the corresponding snapshots (Figure 5b,c from the left view), we observed that the structures were slightly stretched, and the structures became smaller, indicating the alignment of structures with flow direction, which is likely to induce a slight decrease in the viscosity [52].  To further understand the relationship between the structural evolution and the changes in shear viscosity, the continuous phase changes for the mixed solution in Figure 4a were analyzed at a shear rate of 1.08 × 10 9 s −1 . A series of snapshots (Figures 5a-g) were taken from the CGMD simulation to visualize the time dependent multi-stage phase evolution. In Figure 4b, from a to b, there was an increase in viscosity, but the increase only lasted for a short time. From Stage I in Figures 5a,b, we can see that the structure of the mixture had no significant changes, apart from the slight stretching under shearing (top view). The oscillating increase in viscosity at Stage I can be attributed to the increase of viscosity due to the extension of the PAM molecules [61,62] and the elastic response of the resistance of disentanglement action of surfactant micelles [59]. At Stage II, the viscosity decreased with time and the stage lasted for about 0.8 ns (Figure 4b, from b to c). From the corresponding snapshots (Figures 5b,c from the left view), we observed that the structures were slightly stretched, and the structures became smaller, indicating the alignment of structures with flow direction, which is likely to induce a slight decrease in the viscosity [52].  In the following stage (Figure 4b, from c to d), the viscosity increased with time and this part lasted for about 1 ns. In this stage, we can see from the snapshots in Figure 5 that the interconnected networks consisting of PAM and CTAC were very stretched, both the CTAC micelles and PAM polymer chains were elongated, and the original tight structures were unfolded (Figure 5c,d). The expenditure of the contact area with fluids leads to an increase in viscosity [63]. Coupled with the alignment of structures to the shear direction (from the left view in Figure 5c,d), the viscosity showed an oscillation increase. Next, a reduction in viscosity appeared at Stage IV (Figure 4b, from d to e), where the reduction lasted for a long time, about 7 ns. The viscosity decreased very slowly with time. We believe that this reduction was due to the breakage of aggregates and the further alignment of structures. In the top view images in Figure 5d,e, we observed some micelles of small size as well as a reduction in the projection area of the aggregate structures. In Stage V, the shear viscosity reached a constant value independent of the time duration (Figure 4, from f to g), the viscosity reached a plateau regime at 9 ns, and the stage lasted for 2 ns. From Figure 5f to 5g, this plateau can be attributed to the dynamical equilibrium between inner structures (Figure 5f,g).

Low Shear Rates for CTAC Based Solutions
Compared with the PAM based solutions, the viscosity curves (Figure 6a) of the CTAC based solutions showed less complexity with only one peak in Figure 6b. The curve in Figure 6b is a theoretical construct tendency for the viscosity curves of CTAC based solutions under low shear rates (Figure 6a). The volume fraction of the PAM in the CTAC based solutions was much larger than that in the PAM based solutions due to the increase in the concentration of PAM in the CTAC based solutions while the volume fraction of the PAM was close to or larger than that of the CTAC micelles. Thus, PAM becomes the main factor that influences the rheology of the CTAC based solutions. CTAC; cyan and blue, hydrophobic tail of CTAC; magenta, aromatic ring of Sal − ; khaki, charge group of Sal -; red, the acidammide group of PAM; green, carbon backbone of PAM. Na + (yellow), Cl − (dark green), water sites (reddish), and antifreeze water sites (grey) were hidden.
In the following stage (Figure 4b, from c to d), the viscosity increased with time and this part lasted for about 1 ns. In this stage, we can see from the snapshots in Figure 5 that the interconnected networks consisting of PAM and CTAC were very stretched, both the CTAC micelles and PAM polymer chains were elongated, and the original tight structures were unfolded (Figures 5c,d). The expenditure of the contact area with fluids leads to an increase in viscosity [63]. Coupled with the alignment of structures to the shear direction (from the left view in Figures 5c,d), the viscosity showed an oscillation increase. Next, a reduction in viscosity appeared at Stage IV (Figure 4b, from d to e), where the reduction lasted for a long time, about 7 ns. The viscosity decreased very slowly with time. We believe that this reduction was due to the breakage of aggregates and the further alignment of structures. In the top view images in Figures 5d,e, we observed some micelles of small size as well as a reduction in the projection area of the aggregate structures. In Stage V, the shear viscosity reached a constant value independent of the time duration (Figure 4, from f to g), the viscosity reached a plateau regime at 9 ns, and the stage lasted for 2 ns. From Figure 5f to 5g, this plateau can be attributed to the dynamical equilibrium between inner structures (Figures 5f,g).

Low Shear Rates for CTAC based Solutions
Compared with the PAM based solutions, the viscosity curves (Figure 6a) of the CTAC based solutions showed less complexity with only one peak in Figure 6b. The curve in Figure 6b is a theoretical construct tendency for the viscosity curves of CTAC based solutions under low shear rates (Figure 6a). The volume fraction of the PAM in the CTAC based solutions was much larger than that in the PAM based solutions due to the increase in the concentration of PAM in the CTAC based solutions while the volume fraction of the PAM was close to or larger than that of the CTAC micelles. Thus, PAM becomes the main factor that influences the rheology of the CTAC based solutions. The structure evolution for the mixture at a shear rate of 8.65 × 10 8 s −1 from Figure 6a was captured and presented in Figure 7 to analyze the time dependent phase evolution within the mixture. The CTAC micelles exist as a minority phase, which is mostly wrapped by PAM polymer chains. When shear force is applied, the PAM polymer chains first endure most of the shear load and then extend, leading to the increase in viscosity (Figure 6b, from a to b). This stage only lasts for a The structure evolution for the mixture at a shear rate of 8.65 × 10 8 s −1 from Figure 6a was captured and presented in Figure 7 to analyze the time dependent phase evolution within the mixture. The CTAC micelles exist as a minority phase, which is mostly wrapped by PAM polymer chains. When shear force is applied, the PAM polymer chains first endure most of the shear load and then extend, leading to the increase in viscosity (Figure 6b, from a to b). This stage only lasts for a short time when compared to the other two stages. The resistance of the disentanglement of CTAC micelles also contributes to the increment in viscosity. This can be traced from the top view in Figure 7a,b. Then, the aggregates, especially CTAC micelles, gradually align with the flow (Figure 7b-d), resulting in a continuous decrease in both viscosity and stabilization (Figure 6b, from b to c). At the last stage, the viscosity value reached a plateau regime (Figure 6b, from c to d); from the snapshots in Figure 7c,d, it can be observed that the structures have reached a dynamic equilibrium and no structural changes appear.
Polymers 2020, 12, 265 9 of 17 short time when compared to the other two stages. The resistance of the disentanglement of CTAC micelles also contributes to the increment in viscosity. This can be traced from the top view in Figures  7a,b. Then, the aggregates, especially CTAC micelles, gradually align with the flow (Figures 7b-d), resulting in a continuous decrease in both viscosity and stabilization (Figure 6b, from b to c). At the last stage, the viscosity value reached a plateau regime (Figure 6b, from c to d); from the snapshots in Figures 7c and 7d, it can be observed that the structures have reached a dynamic equilibrium and no structural changes appear. Comparing the CTAC based solution and PAM based solution, it can be found that the viscosity is more influenced by the structural evolution of the micelles and the viscosity curves are complicated when CTAC plays a dominant role in rheology. When the concentration of PAM increases, the curves become much smoother due to the long chains of PAM. The polymer chains can absorb and transfer the shear stress to the elastic energy during an enduring and gradual process [47].

Relation between Shear Rates and Characteristic Time
After calculating the shear viscosity versus time, we now compare the relation between shear rates and characteristic time. As some of the simulation results data were fluctuant, where it is difficult to calculate the characteristic time; here, we selected one group of data as an example. We used the "plateau" time as the characteristic time to describe the dynamics of the solutions, which has been used in many other studies [64,65]. The plateau time t is defined at which viscosity reaches the dynamically stable value in the plateau region under the shear rate. It has been found that the relation between characteristic time and shear rate is consistent with the power low dependence of the form ∝ γ in experiments [64,65]. Here are our simulation results.
We chose a solution with CCTAC = 0.3 mol·L −1 and CPAM = 2.96 × 10 −3 mol·L −1 as an example. In the CTAC based solutions, surfactant micelles are the minority and wrapped by PAM polymer chains (PAM molecules) in the PAM-dominated solution, while the PAM polymer chains are the backbones of the aggregates. The excess PAM polymer chains are connected by insufficient surfactant micelles, and PAM chains cannot form strong network structures by themselves. When the aggregates are stretched and unfolded, the structures are divergent and almost have no capability in improving the fluid retention capacity of the solution, and the viscosity no longer increases. Comparing the CTAC based solution and PAM based solution, it can be found that the viscosity is more influenced by the structural evolution of the micelles and the viscosity curves are complicated when CTAC plays a dominant role in rheology. When the concentration of PAM increases, the curves become much smoother due to the long chains of PAM. The polymer chains can absorb and transfer the shear stress to the elastic energy during an enduring and gradual process [47].

Relation between Shear Rates and Characteristic Time
After calculating the shear viscosity versus time, we now compare the relation between shear rates and characteristic time. As some of the simulation results data were fluctuant, where it is difficult to calculate the characteristic time; here, we selected one group of data as an example. We used the "plateau" time as the characteristic time to describe the dynamics of the solutions, which has been used in many other studies [64,65]. The plateau time t is defined at which viscosity reaches the dynamically stable value in the plateau region under the shear rate. It has been found that the relation between characteristic time and shear rate is consistent with the power low dependence of the form t ∝ γ n in experiments [64,65]. Here are our simulation results.
Under low shear rates (shear rates equal and smaller than 1.3 × 10 9 s −1 , Figure 8a), the four curves were similar to each other, thus we used one plateau time t 1 to represent the four curves, and the corresponding shear rate was 1.3 × 10 9 s −1 , which was the turning point (shear thinning starts). Similarly, four other characteristic time points can be obtained from Figure 8b. When the shear rate is 1.73 × 10 11 s −1 , it is hard to distinguish the plateau time point, thus we did not use this data. Figure 9 shows the relation between the characteristic time and the shear rates. Under low shear rates (shear rates equal and smaller than 1.3 × 10 9 s −1 , Figure 8a), the four curves were similar to each other, thus we used one plateau time t1 to represent the four curves, and the corresponding shear rate was 1.3 × 10 9 s −1 , which was the turning point (shear thinning starts). Similarly, four other characteristic time points can be obtained from Figure 8b. When the shear rate is 1.73 × 10 11 s −1 , it is hard to distinguish the plateau time point, thus we did not use this data. Figure  9 shows the relation between the characteristic time and the shear rates.  From the figure, we can note that the relation between characteristic time and shear rate is consistent with the power law dependence of the form ∝ γ , which coincides with the experimental studies [65,66]. The exponent n is −0.94, which shows an agreement with Ma's experimental results [65].

Shear Viscosity Versus the Concentration of CTAC for Polymer/Surfactant Solutions.
We next summarize the shear viscosity as a function of CTAC concentration for all PAM based solutions in Supplementary Materials Figure S5. The results showed that with the increase in the concentration of CTAC, the shear viscosity increased and then decreased, indicating good agreement with the reported experimental results [67][68][69]. Here, we analyzed the influence of CTAC concentration in a mixture solution by using the initial equilibrated morphology data (before shear action). According to Wang's published work [70], the shapes of the micelles can be classified into three certain types based on the length ratio of the long axis to the short axis of the micelle. In the present work, when the ratio was around 1, the micelle was defined as a spherical micelle; when the Under low shear rates (shear rates equal and smaller than 1.3 × 10 9 s −1 , Figure 8a), the four curves were similar to each other, thus we used one plateau time t1 to represent the four curves, and the corresponding shear rate was 1.3 × 10 9 s −1 , which was the turning point (shear thinning starts). Similarly, four other characteristic time points can be obtained from Figure 8b. When the shear rate is 1.73 × 10 11 s −1 , it is hard to distinguish the plateau time point, thus we did not use this data. Figure  9 shows the relation between the characteristic time and the shear rates.  From the figure, we can note that the relation between characteristic time and shear rate is consistent with the power law dependence of the form ∝ γ , which coincides with the experimental studies [65,66]. The exponent n is −0.94, which shows an agreement with Ma's experimental results [65].

Shear Viscosity Versus the Concentration of CTAC for Polymer/Surfactant Solutions.
We next summarize the shear viscosity as a function of CTAC concentration for all PAM based solutions in Supplementary Materials Figure S5. The results showed that with the increase in the concentration of CTAC, the shear viscosity increased and then decreased, indicating good agreement with the reported experimental results [67][68][69]. Here, we analyzed the influence of CTAC concentration in a mixture solution by using the initial equilibrated morphology data (before shear action). According to Wang's published work [70], the shapes of the micelles can be classified into three certain types based on the length ratio of the long axis to the short axis of the micelle. In the present work, when the ratio was around 1, the micelle was defined as a spherical micelle; when the From the figure, we can note that the relation between characteristic time and shear rate is consistent with the power law dependence of the form t ∝ γ n , which coincides with the experimental studies [65,66]. The exponent n is −0.94, which shows an agreement with Ma's experimental results [65].

Shear Viscosity Versus the Concentration of CTAC for Polymer/Surfactant Solutions
We next summarize the shear viscosity as a function of CTAC concentration for all PAM based solutions in Supplementary Materials Figure S5. The results showed that with the increase in the concentration of CTAC, the shear viscosity increased and then decreased, indicating good agreement with the reported experimental results [67][68][69]. Here, we analyzed the influence of CTAC concentration in a mixture solution by using the initial equilibrated morphology data (before shear action). According to Wang's published work [70], the shapes of the micelles can be classified into three certain types based on the length ratio of the long axis to the short axis of the micelle. In the present work, when the ratio was around 1, the micelle was defined as a spherical micelle; when the ratio was larger than 1 and there was no branch in the micelle, the micelle was defined as a wormlike micelle; when a branch exists, the micelle was a branched micelle. Figure 10 shows the representative morphologies of the micelle and branch node. The representative morphology was selected from the simulation results. From the figure, the spherical micelle had the smallest size; two or more spherical micelles coalesced and then formed a wormlike micelle; two wormlike micelles or one spherical micelle and one wormlike micelle coalesced and then a branch node exists and branched micelle forms, the branched micelle had the largest size [51]. Thus, if the concentration of CTAC is low, only small sized spherical micelles and wormlike micelles form; if the concentration increases, large wormlike micelles, and branched micelles will form. The concentration will influence the morphology of the micelle, which can further lead to different viscosity behaviors in rheology. Figure 11 shows the morphology proportion, number of branch nodes, and viscosity as a function of CTAC concentration in PAM based solutions.
Polymers 2020, 12, 265 11 of 17 ratio was larger than 1 and there was no branch in the micelle, the micelle was defined as a wormlike micelle; when a branch exists, the micelle was a branched micelle. Figure 10 shows the representative morphologies of the micelle and branch node. The representative morphology was selected from the simulation results. From the figure, the spherical micelle had the smallest size; two or more spherical micelles coalesced and then formed a wormlike micelle; two wormlike micelles or one spherical micelle and one wormlike micelle coalesced and then a branch node exists and branched micelle forms, the branched micelle had the largest size [51]. Thus, if the concentration of CTAC is low, only small sized spherical micelles and wormlike micelles form; if the concentration increases, large wormlike micelles, and branched micelles will form. The concentration will influence the morphology of the micelle, which can further lead to different viscosity behaviors in rheology. Figure  11 shows the morphology proportion, number of branch nodes, and viscosity as a function of CTAC concentration in PAM based solutions.  In our simulation, the volume fraction of PAM was lower than that of CTAC in the PAM based solutions, thus the viscosity of the mixture is likely to be affected by the rheology property of CTAC. In addition, the CTAC micelles have a diverse morphological adaptivity and can self-assemble into spherical micelles, wormlike micelles, branched micelles, network micelles, and so on. Those morphological changes also have a large impact on the viscosity of solution.
It can be seen in Figure 11 that CTAC micelles show a continuous morphological transition among the spherical shape, wormlike shape, and branched shape when the concentration of CTAC ratio was larger than 1 and there was no branch in the micelle, the micelle was defined as a wormlike micelle; when a branch exists, the micelle was a branched micelle. Figure 10 shows the representative morphologies of the micelle and branch node. The representative morphology was selected from the simulation results. From the figure, the spherical micelle had the smallest size; two or more spherical micelles coalesced and then formed a wormlike micelle; two wormlike micelles or one spherical micelle and one wormlike micelle coalesced and then a branch node exists and branched micelle forms, the branched micelle had the largest size [51]. Thus, if the concentration of CTAC is low, only small sized spherical micelles and wormlike micelles form; if the concentration increases, large wormlike micelles, and branched micelles will form. The concentration will influence the morphology of the micelle, which can further lead to different viscosity behaviors in rheology. Figure  11 shows the morphology proportion, number of branch nodes, and viscosity as a function of CTAC concentration in PAM based solutions.  In our simulation, the volume fraction of PAM was lower than that of CTAC in the PAM based solutions, thus the viscosity of the mixture is likely to be affected by the rheology property of CTAC. In addition, the CTAC micelles have a diverse morphological adaptivity and can self-assemble into spherical micelles, wormlike micelles, branched micelles, network micelles, and so on. Those morphological changes also have a large impact on the viscosity of solution.
It can be seen in Figure 11 that CTAC micelles show a continuous morphological transition among the spherical shape, wormlike shape, and branched shape when the concentration of CTAC Figure 11. Results for the morphology proportion, number of branch nodes, and viscosity as a function of CTAC concentration.
In our simulation, the volume fraction of PAM was lower than that of CTAC in the PAM based solutions, thus the viscosity of the mixture is likely to be affected by the rheology property of CTAC. In addition, the CTAC micelles have a diverse morphological adaptivity and can self-assemble into spherical micelles, wormlike micelles, branched micelles, network micelles, and so on. Those morphological changes also have a large impact on the viscosity of solution.
It can be seen in Figure 11 that CTAC micelles show a continuous morphological transition among the spherical shape, wormlike shape, and branched shape when the concentration of CTAC increased (more morphology images are in Supplementary Materials Figure S6). This transition can explain the viscosity changes in Figure 11. From Figure 11, we can observe that with the increase in the concentration of CTAC, CTA + ions can form larger sized micelles such as wormlike micelles and branched micelles, and small size micelles (spherical micelles) disappear. Furthermore, the larger micelle proportion also increased with the increase in CTAC concentration. In detail, when the CTAC concentration is from 0.1 to 0.2 mol·L −1 , only spherical micelles and wormlike micelles exist, and the proportion of wormlike micelles increases with the increase in CTAC concentration, which indicates that the micelle size increases with the increase in CTAC concentration. Since the wormlike micelles can be more effective in resisting the shear action and suppress the fluctuation of the flow field than spherical ones, so the viscosity increases with the increase in proportion of wormlike micelle. However, branched micelles appear when the CTAC concentration further increases. Although a high CTAC concentration leads to the formation of large sized micelles (branched micelles and wormlike micelles), the formation of branched micelles can reduce the viscosity of the solutions as branches provide a sliding mechanism for stress relaxation, as has been reported [50,70,71], leading to a reduction in the viscosity of the solutions. When the CTAC concentration is further increased from 0.25 to 0.3 mol·L −1 , both the branch nodes and the proportion of branched micelles increase. The stress relaxation is faster in a solution with a higher degree of branching [71], thus the viscosity has a further decrease. It should be noted that the concentration of micelles also increases during the period, but it is still insufficient to compensate the reduction effect on viscosity induced by the morphological change of CTAC micelles.

Shear Viscosity Versus the Concentration of PAM for Polymer/Surfactant Solutions
The shear viscosities of CTAC based solutions are summarized as a function of PAM concentration in Supplementary Materials Figure S7, and the corresponding structure snapshots are in Figure S8. The group at a shear rate of 8.65 × 10 8 was selected for analyzing the influence of PAM concentration in a mixture solution ( Figure 12). The initial equilibrated morphology data (before shear action) were used to analyze the micelle shape and the morphology proportion. With the increase in PAM concentration, the shear viscosity of solutions first decreased and then increased. In the polymer/surfactant solutions, the PAM polymer chains normally work as "bridges" and bind surfactant micelles together [47]. When we increased the concentration of PAM from 2.96 × 10 −3 to 4.45 × 10 −3 mol·L −1 , the interaction between the polymers and surfactants became stronger and most of the micelles changed from wormlike micelles to branched micelles. The existence of branched micelles led to the reduction in viscosity for the mixed solutions, thus the viscosity of the CTAC based solutions first decreased with the increase in concentration of PAM. increased (more morphology images are in Supplementary Materials Figure S6). This transition can explain the viscosity changes in Figure 11. From Figure 11, we can observe that with the increase in the concentration of CTAC, CTA + ions can form larger sized micelles such as wormlike micelles and branched micelles, and small size micelles (spherical micelles) disappear. Furthermore, the larger micelle proportion also increased with the increase in CTAC concentration. In detail, when the CTAC concentration is from 0.1 to 0.2 mol·L −1 , only spherical micelles and wormlike micelles exist, and the proportion of wormlike micelles increases with the increase in CTAC concentration, which indicates that the micelle size increases with the increase in CTAC concentration. Since the wormlike micelles can be more effective in resisting the shear action and suppress the fluctuation of the flow field than spherical ones, so the viscosity increases with the increase in proportion of wormlike micelle. However, branched micelles appear when the CTAC concentration further increases. Although a high CTAC concentration leads to the formation of large sized micelles (branched micelles and wormlike micelles), the formation of branched micelles can reduce the viscosity of the solutions as branches provide a sliding mechanism for stress relaxation, as has been reported [50,70,71], leading to a reduction in the viscosity of the solutions. When the CTAC concentration is further increased from 0.25 to 0.3 mol·L −1 , both the branch nodes and the proportion of branched micelles increase. The stress relaxation is faster in a solution with a higher degree of branching [71], thus the viscosity has a further decrease. It should be noted that the concentration of micelles also increases during the period, but it is still insufficient to compensate the reduction effect on viscosity induced by the morphological change of CTAC micelles.

Shear Viscosity Versus the Concentration of PAM for Polymer/Surfactant Solutions
The shear viscosities of CTAC based solutions are summarized as a function of PAM concentration in Supplementary Materials Figure S7, and the corresponding structure snapshots are in Figure S8. The group at a shear rate of 8.65 × 10 8 was selected for analyzing the influence of PAM concentration in a mixture solution ( Figure 12). The initial equilibrated morphology data (before shear action) were used to analyze the micelle shape and the morphology proportion. With the increase in PAM concentration, the shear viscosity of solutions first decreased and then increased. In the polymer/surfactant solutions, the PAM polymer chains normally work as "bridges" and bind surfactant micelles together [47]. When we increased the concentration of PAM from 2.96 × 10 −3 to 4.45 × 10 −3 mol·L −1 , the interaction between the polymers and surfactants became stronger and most of the micelles changed from wormlike micelles to branched micelles. The existence of branched micelles led to the reduction in viscosity for the mixed solutions, thus the viscosity of the CTAC based solutions first decreased with the increase in concentration of PAM. The viscosity then increased slightly due to the further increase in the concentration of PAM (which increased to 7.41 × 10 −3 mol·L −1 ), but was still a relatively low value due to the continuous The viscosity then increased slightly due to the further increase in the concentration of PAM (which increased to 7.41 × 10 −3 mol·L −1 ), but was still a relatively low value due to the continuous formation of branched micelles in the CTAC based solution. When the concentration of PAM was higher than 7.41 × 10 −3 mol·L −1 , the volume fraction of PAM was larger than that of CTAC, which means that polymers dominate the rheology behavior and the polymer chains become the framework of the mixed structures. Thus, the morphology of the surfactant micelles is limited due to the polymer framework. Branched micelles disappear and the size of the micelles decreases, and small wormlike micelles form instead. The disappearance of branched micelles leads to the increase of the viscosity. Meanwhile, the increase in the concentration of PAM also results in the increase in viscosity. The viscosity change caused by PAM can compensate for the reduction effect on viscosity induced by the morphological change of CTAC micelles from wormlike micelles to spherical micelles, and leads to the further increase in viscosity.

Conclusions
We used a NEMD based numerical approach to understand the chemical-physical interactions within multi-phase PAM/CTAC/NaSal solutions and their impacts on macroscopic rheological behaviors. By analyzing the relationship between shear viscosity and shear rates, a plateau followed by a shear thinning zone was observed, which is in good agreement with the experimental results. We also discovered the "two peaks" phenomenon for the time by simulation after investigating the time dependent shear viscosity and the corresponding inner structural evolution between different phases under different shear rates in the PAM based solutions. The viscosity curve of the CTAC based multi-phase solution became smoother than that of the PAM based (CTAC dominated) solutions due to the fact that PAM chains absorb and transfer the shear energy into elastic energy by extending/stretching the polymer chains during an enduring and gradual process. The relation between characteristic time and shear rate is consistent with the power law. After looking into the relationship between shear viscosity and the concentration of the solutes, an explicit reduction in viscosity was found because of the phase transition of CTAC micelles from spherical micelles into branched micelles.