Multi-scale modelling of dilute viscoelastic liquids: Atomistic to mesoscale mapping of polymer solutions

A framework predicting the rheological (storage and loss moduli, first normal stress coefficient, and relaxation time) and transport (viscosity, diffusivity) properties of non-Newtonian dilute polymer solutions at mesoscales (e.g. ∼ ns to µ s) from the atomistic-scale molecular behaviour is presented. More specifically, the rheological behaviour differences of OCP and PMA polymer solutions in PAO-2 oil are simulated using both atomistic molecular dynamics (MD) and many-body dissipative particle dynamics (mDPD) within a temperature range of 313-373 K. The simulation methodology described is able to distinguish itself from the standard DPD model by accurately reproducing the shear-thinning with high sensitivity, as seen in the atomistic MD simulations at high shear rates (e.g. 10 8 − 10 13 s − 1 ). It is shown that the model is well-suited to compute properties such as first normal stress differences and relaxation times that are difficult to estimate at atomistic scales due to the low signal-to-noise ratio. Moreover, the Schmidt numbers ( > 10 3 ) are predicted with high accuracy when compared with the values from atomistic-scale simulations. The proposed model is able to predict relaxation times of dilute polymer mixtures that are difficult to be obtained using state-of-the-art rheometers. Finally, it is found that the terminal relaxation time of PMA polymer chain does not vary monotonically as a function of temperature, unlike in the case of OCP; this is significant for describing viscoelastic behaviour at macroscales where satisfactory constitutive equations are not available.


Introduction
Dilute viscoelastic liquids are expected to be used in applications such as direct cooling in electronic and automotive industries, and soft electronics manufacturing, among many other industrial and biomedical processes.These liquids are valued owing to their viscoelastic nature, dielectric nature and thermal stability [1,2,3,4].New combinations of polymers and solvents are required to design novel viscoelastic dilute polymer solutions [5,6].In this work, we focus our investigation to the dilute polymer solutions relevant to heat transfer applications.Given the technical challenges associated with experimental methods in the dilute polymer regime [7,8], and the requirement of satisfactory constitutive equations to model such liquids in continuum scale simulations, an alternative approach for modelling such liquids is required [9].
The fluid chemistry influences its thermodynamic (e.g.density), morphological and transport (e.g., diffusion coefficient) properties, and those can be modelled via atomistic molecular dynamics simulations (MD).United-atom models are utilised to understand the Rouse to reptation cross-over behaviour of short-chain concentrated polymer solutions, their steady-state shear and elongational flow behaviour [10,11,12,13].Still, macroscopic properties such as deformation response, storage and loss moduli, and relaxation times of many such liquids are challenging to be obtained using a full-atomistic MD or united-atom MD, given the expensive computational schemes for long-chain polymer solutions and low signal-to-noise ratio [14,15].Therefore, a mesoscale modelling routine is required to bridge the gap between the microscale behaviour and macroscopic observations [16].
An extensive body of work over the past three decades in the relevant literature is dedicated to mesoscale modelling via the Lagrangian mechanics [17,18,19,20,21,22,23].Polymer melts and their entanglement characteristics are studied along with their relaxation modes using Kremer-Grest models [24,25].Concurrently, the widely applied methodology is dissipative molecular dynamics (DPD) developed by Groot and Warren [26].DPD uses a Langevin form of equation where the repulsive, dissipative and random forces lead to the description of the system of interest.Several works in the field have concentrated on predicting rheological and interfacial properties of aqueous solutions with dispersed macromolecules [27,28,19].As a response to anomalies in behaviour under anisotropic forces, vapour-liquid coexistence, and to assimilate heat transfer characteristics, DPD models have been modified to develop many-body DPD (mDPD) [29,30], energyconserving DPD (eDPD) [31] and so on.However, a modelling scheme to reproduce viscoelastic behaviour among a wide range of liquids is still being developed.
The standard DPD model is shown to provide misleading results for viscosity (e.g.shear thickening) in the case of non-equilibrium DPD simulations of simple fluids with spherical particles [32].This is in addition to the fact that conventional DPD underestimates Schmidt numbers by three orders of magnitude [30].Several suggestions extending from Lees-Edwards boundary condition modifications to thermostat modifications have been prescribed in the literature to overcome this [33,34,35,36,37].The motivation for schemes to improve the prediction of non-Newtonian behaviour is seen even recently.For example, Zhao et al. [38] tried to use Gaussian process regression (GPR) based active learning to connect bulk rheology to microstructure dynamics as a way to reproduce experimental shear-thinning of viscoelastic liquids.Recent efforts are engaged in parametrising the mDPD models to reproduce non-Newtonian behaviour in fluids.Zhao et al. [39] showcased the computation of surface tension and first normal stress difference of polymer solutions.Jamali et al. [40] proposed a methodology to relate the compressibility of the fluid with the Flory-Huggins theory to take into consideration the mixing effects of different components to generalise the modelling approach.At the same time, a modelling strategy to consistently map all the relevant properties into real units with accurate scaling is vital, and there is a scope for a consistent approach to be introduced [36,41].Moreover, the chemical structure of the organic molecules which influences the order of magnitude of Schmidt numbers needs to be included in the scaling-up process as well.
The investigation here demonstrates the simulation of polyalphaolefin (PAO-2) oil and polymer solutions of two different polymer chains, namely, olefin co-polymer (OCP) and polydodecylemethacrylate (PMA) in PAO-2.The specified oil and polymer compounds are extensively used in industrial applications such as lubrication, and are proposed as candidates for cooling technologies [42,43,44,45,46].The radius of gyration of OCP and PMA polymer chains varies differently in PAO-2 as a function of temperature [47,48].These structural differences lead to differences in their rheological behaviour, as shown in the literature [49].
The mapping of the chemistries of the different polymer chains in the oil appropriately to mesoscale shows the ability of the simulation methodology to be a tool for simulating any such liquids.The study models realistic Schmidt numbers (> 10 3 ) in mesoscale as seen for such liquids in experiments [50].The Schmidt number matching is essential to capture the correct hydrodynamic interactions of the liquid [51].The polymer chains are distinguished by inducting an angular cosine potential into the mDPD model, whose parameters vary with the chemistry and temperature to qualitatively model the radius of gyration.The simulated mesoscale polymer solution is able to behave as described by Zimm's model [52] and is able to show diluteness.
The simulation carried out here is able to showcase the variation in mildly viscoelastic non-Newtonian liquids with shear via non-equilibrium mDPD.Moreover, the zero shear viscosity obtained via non-equilibrium mDPD agrees with equilibrium mDPD computations, validating the process.Additionally, properties such as first normal stress coefficient, storage and loss moduli, and relaxation times which influence fluid behaviour in continuum scales are computed.Therefore, the detailed description of the methodology provides a novel tool that can be extended to the simulation of any non-Newtonian polymer liquids with minimal tuning of the scaling units and mDPD model parametrisation.

Solvent modelling
The simulation is designed in such a manner so as to map the important system parameters such as the characteristic length, mass and time for the temperature range of interest from atomistic MD to mDPD.Here, the solvent PAO-2 being investigated is an oil that has high Schmidt numbers (Sc) of the order of 10 3 − 10 4 .The goal is to find the appropriate mDPD parameter-set such that they reproduce Sc (as computed by Equation 13).This ensures that the relative motion of the particles due to momentum transport and mass transport manifest in the same manner at the larger length-and time-scales.At the same time, it must be ensured that the Sc is derived at the same thermodynamic conditions such as temperature and pressure (consistent with unit mapping) as obtained at the atomistic scale simulations.The fundamentals of the mDPD modelling are described in the Supporting Information.To provide maximum control on the modelling of the polymer solutions of interest, the work involves the tuning of five parameters, namely ρ, A ij , B ij , γ ij and r D (see Equations S1 to S7).A simulation framework is shown in Figure 1 for the mDPD simulation at temperature, T * = 313 K and pressure, p * = 1 atm with ' * ' from here on indicating quantities whose values are in real units.PAO-2 molecule has a molecular weight (M * w ) of 282 g mol −1 .
By considering the density of PAO-2 at the temperature of 313 K obtained from atomistic MD (ρ * = 780 kg m −3 ), we compute the molecular volume as 600.34 Å3 .The simulation procedure starts with a number density, ρ > 6 and local density-dependent conservative force cut-off, r d = 0.75 and a timestep, dt = 0.01.The reduced density range is in agreement with the description in the literature [53,39].The equation of state for an mDPD fluid has been approximated as follows: where α = 0.101, c = 4.16 and d = 19 for an organic fluid of the nature of PAO-2 [30,53].This simplifies the problem of mDPD model parametrisation to two parameters, namely A ii and B ii .To find their exact solutions, another equation is required.By differentiating the Equation 1, we can relate the reduced units isothermal compressibility, κ to these parameters The isothermal compressibility of the solvent from the atomistic MD (κ * ) is used in this Here, N m = 1 is the number of molecules of PAO-2 solvent being represented by a single mDPD bead and R is the universal gas constant.The κ * values agree with the experimental results available in the literature [46].The initial values of ρ and r d are used to solve the above equations to obtain A ii and B ii .The solvent system with these A ii and B ii is simulated.Since Equation 1 is an approximated equation of state, the value of B ii is optimised such that the virial pressure simulated oscillates around the pressure as computed in Table 1.The mDPD scaling units for the rest of the temperatures are provided in the Supporting Information from Table S9 to S11.Following this, simulations are carried out for maximum stable values of dissipative force parameters at this dt, γ max = 18.0 and r D,max = 2.0.If the Sc mDPD < Sc atomistic , ρ is increased and the above steps are repeated.If Sc mDPD > Sc atomistic , γ ij and r D are optimised till Sc mDPD = Sc atomistic .Figure 2 shows the typical size difference between the simulation boxes described in this study, with the mesoscale box ensuring a larger ensemble-averaging of properties.

Polymer modelling
Two different polymer chemistries are simulated in this work, namely an olefin copolymer (OCP) composed of ethylene and propylene monomers, and poly dodecylmethacrylate (PMA) that has a branched structure at the atomistic scale with a carboxylate function group.In the case of OCP, the monomer consists of short alkene molecules that make it chemically similar  to the PAO-2 solvent molecules.Therefore, in mDPD simulation, the solvent-polymer bead interactions can be treated as similar to that of solvent-solvent bead interactions.In other words, A ij = A ii and B ij = B ii .This leads to the treatment of PAO-2 oil as a theta solvent in case of the OCP polymer in the mDPD simulations.In a theta solvent, the repulsive and attractive interactions between the monomer beads of a polymer chain are balanced [54].The bonded interactions of the polymer chain are derived from FENE potential, and angle stiffness is provided by a cosine force-field (see Equations S8 to S11).The parameters for the simulation of OCP solution are shown in Table 2. Since OCP is composed of two types of monomers, i.e., ethylene and propylene, the differentiation of those has been reflected in the use of k angle,1 and k angle,2 , respectively.After obtaining the initial values from Everaers et al. [55], the final values of the two parameters are optimised so as to agree with the radius of gyration (R g ) values computed from the atomistic scale MD.
On the other hand, PMA has a different monomer structure and therefore, the polymersolvent bead interactions must consider the mixing effects.In mesoscale simulations, such effects are modelled by considering the Flory-Huggins theory [56].Under this theory, the magnitude of the attractive interaction parameters in the conservative force in Equation S2 are computed as where χ ij is called the Flory-Huggins interaction parameter that is defined as where V is the volume of the simulation box.Here, δ i is computed using the cohesive energy density (CED) of a species i from their non-bonded energies obtained from the relationship Here, E k nb is the ensemble average of the non-bonded energies of k molecules of species i and k i=1 E iso nb is the sum of non-bonded energies of k individual molecules in an isolated simulation box of species i.From the atomistic scale simulations, the total pairwise interaction energies using a simulation box of pure PAO-2 or PMA are used to compute E k nb .Similarly, the E iso nb is computed using the total pairwise interaction energies in the same simulation box retaining a single molecule of PAO-2 or PMA.The parameters obtained are shown in  for the mDPD simulation of the PMA polymer solution at 313 K.A similar exercise in case of OCP polymer chain shows that δ i = δ j , confirming the hypothesis to treat the OCP polymer solution as a theta solution.The parameters for the rest of the temperatures are shown in the Supporting Information (see Tables S6 to S8).At the same time, the modelling of the polymer chain used in the study is complete only with the correct parameters of angle stiffness, k angle .It is tuned to reproduce the variation of R g as a function of temperature.

Simulation details
The time-integration of the mDPD equations is done using a modified velocity-verlet algorithm using the λ parameter of 0.5 [26].A typical periodic simulation box consists of 10989 solvent beads.The ρ equal to 7.9 obtained at 313 K is used for the modelling and simulation of the mDPD system at other temperatures.Polymer chains of different lengths between 11 to 77 beads are simulated.The upper limit of the chain length simulated here ensures that the radius of gyration is less than or equal to half of the simulation box length ≈ 11.94 mDPD units, thereby preventing any artefacts.The computation of equilibrium properties and rheological analysis of polymer solutions have been reported for a single OCP or PMA chain of a length of 11 beads representing a molar mass of 3102 g mol −1 and a concentration of ∼ 1000 ppm (0.1%) by weight.Additionally, an OCP polymer solution corresponding to a similar concentration in atomistic simulation of 5.9% by weight (= 203 polymer chains) is also simulated for the sake of validation.
The reported results are from five independent simulations of 10 7 timesteps or 1.7 µs in real units (see Table 1).The non-equilibrium simulation is carried out by box deformation with equations of motion expressed in SLLOD form [57], making it conceptually equivalent to Lees-Edwards periodic boundary conditions [58].The simulated solvent systems are compared to the results obtained from the simulations of original DPD model using parameters described by Groot and Warren [26] at ρ = 3 and ρ = 7.9.The former ρ value is used in the original DPD model, whereas the latter is used to compare with the mesoscale simulations carried out in this study.
The atomistic MD simulation details including the force-field parameters (Table S1 to S5, and Figure S1) are provided in the Supporting Information.In brief, the simulations are done using a box consisting of 170 PAO-2 solvent molecules in which a single polymer chain is dispersed.The simulation uses OPLS-AA force field to account for the polymer and solvent interactions.The mesoscale and atomistic simulations are done using the LAMMPS software with source-code modified in-house to implement the mDPD modelling scheme.

Results and discussion
The atomistic scale modelling of the viscoelastic liquids shows unique behaviour for each polymer solution when computing the radius of gyration as a function of temperature.The mDPD models aim to reproduce these fundamental features such that they can be used to derive different rheological properties at higher scales.In this respect, we initially validate the mDPD model to represent the diluteness of the polymer solution, following which the hydrodynamic radii variation with respect to the temperatures is compared between atomistic and mDPD simulations.It is followed by the computation of equilibrium properties showing the Schmidt number comparison at various temperatures.The analysis is concluded with the computation of the rheological properties of the two different polymer solutions to showcase the similarities and differences in their viscoelastic nature.

Validation of diluteness of polymer solution
If a model is expected to replicate a dilute polymer solution, it has to generally agree with Zimm's polymer theory [52].In essence, it establishes that the polymer chains in the solution do not interact with each other leading to entanglements.In turn, the morphology and rheology exhibited by the polymer chains are functions of their chain lengths.For example, the theory states that in a dilute solution with a theta solvent, where R g is the radius of gyration of the polymer chain (as computed by Equation 11), D poly is the self-diffusion coefficient of the polymer chain in the solution, N is the number of monomer beads and ν is the so called Flory exponent ≈ 0.53.The methodologies to compute the two quantities are described in later sections.In this work, the OCP polymer in PAO-2 acts as a theta solution.As shown in Figure 3, R g increases as a function of N − 1 and D poly decreases as a function of N .The green curves show equations as per Zimm's model in Equation 7. The data points fit well with the curve indicating the dilute state of the polymer solution.Additionally, the computation of the mean-squared displacement of the central monomer of the polymer chain (g 1 (t)) [59] has the following power-law sequence as a function of time Figure 3c shows the three distinct temporal regimes of g 1 (t) for the different polymer chain lengths, agreeing with the Zimm's theory [60].The obtained values of τ Z for different chain lengths agree with the relationship τ Z ∝ N 3ν .However, increasing N beyond 33 leads to deviation from the Zimm relationships of D poly and τ Z , indicating the role of the local stiffness of the chain in sustaining the intermediate regime longer than expected, as described by Hinczewski et al. [61].
In addition to Zimm's relationship, the primary characteristic to verify if a mesoscale fluid acts as a dilute mixture is to establish the ability of solvent molecules to provide effective hydrodynamic interactions to screen the different monomer beads of the polymer [62].To verify if the model reported here is able to provide the required hydrodynamic screening, the computation of the variable is done, where b is the average bond length of the polymer chain and r H is the hydrodynamic radius of the polymer bead computed from the Stoke-Einstein relationship Here, D monomer is the self-diffusion coefficient of the polymer bead [63].As mentioned in the literature, the value of h * should be >= 0.25 to consider the simulated system as a dilute solution.From Figure 3d, it is seen that h * value is above 0.25 at all simulated conditions for the model with the values increasing from 0.53 to 1.017 with increase in the OCP polymer chain length at 313 K.

Structural characteristics
The basic structural characteristic of a polymer chain in a solution is its hydrodynamic radius defined as where M is the mass of a molecule, m i is the mass and r i is the position vector of an atom/bead i of a molecule consisting of n atoms/beads, and r com is the centre of mass of a molecule.Therefore, we require the polymer behaviour at mesoscale to match the one shown at atomistic scale MD as a function of temperature.
As observed in Figure 4, the radius of gyration of OCP in PAO-2 does not vary as a function of temperature at atomistic scale.The comparison of the end-to-end chain length distribution of OCP (see Figure S4) shows a similar mean value in both the scales, even though the distribution is broader in the mesoscale.This can be attributed to the longer runtime in the mesoscale ensuring the chain relaxation to be complete.On the other hand, the R g of PMA polymer chain in PAO-2 increases as a function of temperature.The viscosity index and shear response of the solutions have been demonstrated to be influenced by the variation of R g as a function of temperature [64,65].The PMA polymer chain is reported to be a better viscosity index improver (VII) than OCP in olefin oils because of how its R g increases as a function of temperature.Hence, the mDPD model here is expected to reproduce the same trend as in the case of the atomistic MD.
The behaviour of the atomistic MD is replicated in the mDPD model by tuning the angle stiffness energy, k angle at every temperature.k angle limits the movements of the bonded string of beads.As OCP chain is shown to have a constant R g with temperature, the k angle is given the same values at all temperatures of the mDPD model.The values of R g oscillate around a mean of 1.90 in mDPD units with an increase in temperature.On the other hand, as the PMA chain increases its R g values with temperature, the k angle values are tuned to reproduce the behaviour (see Figure S5 in the Supporting Information).For that, the R g value of mDPD simulation at 313 K is treated as the benchmark, and the k angle is increased at every temperature such that the ratio of R g at the simulated temperature and 313 K is the same for atomistic MD and mDPD.As shown in Figure 4a, this procedure is shown to increase the R g as a function of different temperatures of the mDPD simulation of PMA polymer in PAO-2 with the mean values increasing from 1.83 at 313 K to 2.17 at 373 K.While ensuring the R g mapping, the intramolecular bond distribution of the polymer  chain is verified by computing the pair distribution function g(r) polymer bond = n(r) 4π ρ bead r 2 dr .
For comparison, the bond distribution in the atomistic scale is computed from the sum of all the distances between the bonded atoms of the same monomeric mass as the mDPD bead (see Figure 4b).The distances between the OCP monomer beads are close to 16 Å with the difference of the peaks in the two scales being approximately equal to 1 Å.Similarly, PMA monomer beads are found at a similar distance, with the mDPD scale having beads closer than the equivalent bonded distance between atoms at the atomistic scale.Further comparison of the dihedral distribution function [66] of the polymer chains in both the scales is shown in Figure S6 in the Supporting Information.It shows that the specific peaks of the distribution plateaus with increase in coarse-graining.The decrease in specificity of the peaks of intramolecular features in the implemented modelling scheme has been attributed to the insufficient accounting of entropy in recent literature [67,68].An attempt to match these structural distributions for smaller segments of alkanes is done by Reith et al. [69].The radial distribution function (g(r) solvent−polymer ) between solvent and polymer beads is shown in Figure S7.The impact of the structural differences of g(r) solvent−polymer in atomistic and mesoscale modelling on the dynamics is looked into in a later section.

Equilibrium properties
The motion of particles in a fluid occurs simultaneously due to the momentum transport as well as molecular diffusive motion due to thermal fluctuations.Schmidt number that compares the two kinds of motion is computed as where ν is the kinematic viscosity and D is the self-diffusion coefficient.Here, the momentum diffusivity (ν) values at zero shear rate are calculated from the dynamic viscosity (η 0 ) computed via the Green-Kubo formulation [70].At the same time, the mass diffusivity (D) values are computed via Einstein's relation [71].In accordance with this modelling approach, the Schmidt number has to match at both the atomistic and mesoscale levels.Sc values in the atomistic scale MD match well with with experiments reported in the literature [46].Figure 5 shows the comparison of the PAO-2 solvent Schmidt numbers (Sc solvent ) in atomistic and mDPD simulations at different temperatures under equilibrium conditions.The Sc solvent value decreases as a function of temperature.As seen, the Sc solvent ≈ 40000 at 313 K, before it undergoes a drastic reduction in magnitude with increasing temperature leading to a value ≈ 3000 at 373 K.After carrying out the mDPD modelling following the steps described in Section IV A, the Sc values in the mesoscale simulation agree with that of the atomistic MD at all temperatures, thereby validating the solvent physics.
The relative motion of the OCP polymer chain with respect to the solvent motion is measured by the computation of Sc poly = η 0,solution − η 0,solvent ρD poly .
Here, Sc poly is the Schmidt number of the polymer chain with η solution and η 0,solvent representing dynamic viscosity of polymer solution and pure solvent at zero shear rate, respectively.
Figure 6 shows the Sc poly as a function of temperature.As seen in Figure 6a, Sc poly values for atomistic MD are higher than the Sc poly values obtained in the case of mDPD.Sc poly decreases from 30000 at 313 K to 5100 at 373 K in the atomistic MD, whereas it decreases from 26800 to 867 for the same temperature range in mDPD simulation.In order to verify the veracity of the model, mDPD simulation at the exact polymer concentration of 5.9% by weight as in the case of atomistic MD simulation is carried out.Consequently, Figure 6b shows that Sc poly values are in good agreement expressing unit consistency of polymer concentrations in atomistic and mesoscale models.

Rheology
For a model in the Lagrangian framework to be used for compounds with flow properties different from Newtonian fluids, the rheological behaviour must be reproduced in scales relevant to continuum mechanics.For different polymer chemistries, the dilute solutions have different rheological behaviour at the atomistic scale with respect to the changes in viscosity as a function of shear rate (i.e., shear-thinning) at different temperatures.However, it is essential to see if the shear-thinning behaviour itself is replicated while using the mesoscale models as well or not.Moreover, having an accurate mesoscale model helps in deriving properties that are not possible to be computed at the atomistic scale due to low signal-tonoise ratio and slow dynamics.
The solvent and the two different polymer solutions undergo a Couette flow simulation.The flow is simulated by deforming the periodic box in the y-direction and the velocity gradient is computed in the z-direction of the simulation box.The η( γ) as a function of shear rate ( γ) is computed from the relationship where τ yz is the off-diagonal component of the stress-tensor.In the atomistic scale, it is clearly seen that the solvent and polymer solution shear-thins at very high shear rates (refer to Figure 7).A distinct behaviour with respect to shear-thinning is observed in case of OCP and PMA in their respective solutions at different temperatures by fitting a Carreau-model of the form where η 0 is the Newtonian viscosity that is obtained from the extrapolation of the data to γ = 0, τ s indicates the time constant related to shear-thinning and m is the strain-rate sensitivity coefficient.
The PMA solution starts shear-thinning at a lower shear rate (τ s = 2.07 ns) than the OCP solution (τ s = 0.73 ns) at 313 K. Similarly, at a higher temperature of 373 K, the shearthinning behaviour is observed at lower shear rates for PMA (τ s = 0.80 ns) than what is observed in case of the pure solvent (τ s = 0.10 ns) and OCP polymer solution (τ s = 0.14 ns).At the same time, the mechanism of shear-thinning enhancement in case of PMA and OCP solution is different.While the OCP polymer chain leads to an increase in the hydrodynamic radius of PAO-2 on average, the presence of PMA does not impact the overall radius of PAO-2 molecules (see Figure S2 in the Supporting Information).The enhanced shear-thinning by PMA can be attributed to ease of shear-direction alignment of branched PMA polymer compared to the linear OCP polymer.It is essential to reproduce the shear-thinning behaviour at such high shear rates for these mildly viscoelastic liquids via mesoscale model while proposing a universal methodology that can be extended to strongly viscoelastic medium as well.An oligomer solvent itself, such as the one used in this work, can show shear-thinning better than solvents such as water.This is attributed to the chains aligning in the direction of deformation [46].Figure 8 shows the dynamic viscosity (η) of the PAO-2 solvent and DPD solvent as a function of shear rate (0.001 ≤ γ ≤ 1.0 in mDPD units).The DPD solvents show no viscosity variation with the imposed shear rate and follow more of a Newtonian character within the range, as observed in literature [72].This is attributed to the inherent nature of the DPD equations that cannot capture a significant dependence of η with γ [73].However, the PAO-2 solvent in this work is able to reproduce the shear-thinning behaviour as seen in the atomistic scale.Oils or hydrocarbons can shear-thin at shear rates well below those solvents represented by conventional DPD.Additionally, the model clearly shows that with NEMD, the zero-shear viscosity obtained by the Carreau-model agrees with the zero-shear viscosity obtained via the Green-Kubo methodology in equilibrium simulations of PAO-2.This underlines the consistency of the modelling technique described in this study.Figure 9 shows the viscosity as a function of shear rate for the PAO-2 solvent and the two polymer solutions of OCP and PMA polymer chains.With the addition of the polymer chain, the resultant polymer solutions of OCP and PMA are shown to have a higher viscosity than that of the solvent at lower shear rates.For the concentration of OCP and PMA polymers simulated in mesoscale, they are not adding any more favourable shear-thinning than that as seen in case of the PAO-2 solvent.
Even though shear-thinning at simulated shear rates depicts the mild non-Newtonian behaviour, the characteristic that is essential to prove that a liquid is viscoelastic is the first normal stress coefficient (ψ 1 ) computed as The horizontal lines show the zero-shear viscosity computed using Green-Kubo method under equilibrium conditions.The dashed and dash-dotted lines show the Carreau-model fit.
Here, τ yy is the diagonal component of stress tensor in the direction of flow and τ zz is the component of stress tensor in the direction of the velocity gradient in the Couette flow being examined.For a dilute viscoelastic liquid, the value of ψ 1 has been demonstrated to be greater than zero [74].It is reported that the greater the magnitude of ψ 1 , the greater the fluid tendency to augment instabilities leading to vortex roll-up.In the continuum-scale simulations, the prominent viscoelastic models such as Oldroyd-3-constant and second-orderfluid models show the relationship between the magnitude of ψ 1 and the magnitude of vortical motion [75,76].Thus, the influence of fluid rheology on flow and heat-transfer applications becomes evident.
Figure 10 shows ψ 1 as a function of shear rate in case of a DPD solvent, the simulated PAO-2 solvent and the polymer solutions.As observed, the DPD solvent has physically insignificant values of ψ 1 <≈ 0 at all values of shear rate [77].Similarly, the simulated PAO-2 solvent shows a Newtonian behaviour at low shear rate.With a higher shear rate, the ψ 1 obtains positive values for the simulated solvent.However, with the addition of a polymer chain, ψ 1 is positive at all shear rates simulated.It confirms the viscoelastic nature achieved with the addition of polymer chains in the solvent [78].
The individual contribution of the viscous and elastic components towards the viscoelasticity can be obtained with the computation of storage (G ′ ) and loss moduli (G ′′ ) from which is fitted with with a sum of exponential functions with M relaxation modes.Figure 11a shows the G(t) values smoothened by doing a running average between 0.9t and 1.1t, as described by Sen et al. [79], and the corresponding fits for OCP and PMA polymer solutions at 313 K and 373 K.The long scattered tail of G(t) (t > 1000 mDPD units) emerging due to numerical precision issues at small magnitudes in the long time regime [80] is further smoothened by taking a weighted average [81].G ′ and G ′′ are subsequently estimated from the Fourier transformation of the Equation 19such that and Figure 11b shows G ′ and G ′′ of the different polymer solutions as functions of the frequency ω at the two different temperatures.As observed, G ′ and G ′′ values decrease with increase in temperature.In other words, the elasticity and viscosity decrease with temperature.The ω values at which the slope of G ′ and G ′′ change increase with temperature.This is clearly demonstrated by computing as depicted in Figure S8 in the Supporting Information.The data shows consistency with the results from Green-Kubo formulation for η 0 .At the same time, the value of η * (ω) decreases at a higher ω with increase in temperature.
The final aspect to verify for the present modelling approach is about how the relaxation time provided by the mDPD model behaves as a function of temperature.In the atomistic MD simulation, the longest relaxation time, i.e., terminal relaxation time is computed using the formula where with R(t) depicting the end-to-end vector (r 1 − r l ) of the polymer chain having l monomers [59].As shown in Figure S3 in the Supporting Information, the τ term values of OCP polymer chain decrease with increasing temperature.However, τ term values of PMA could not be computed due to the extremely slow decay of C(t) at the atomistic scale and hence, they are not reported.The computation of the relaxation time using the same procedure is performed for the case of the mDPD simulations (see Figure S9 in the Supporting Information) and the results are shown in Figure 12a.The comparison of the longest relaxation time, τ Z obtained from Equation 19 is performed against τ term in Table S12.The values of τ Z are showcased to be in the range of 0.3τ term to 0.4τ term depending on the temperature and the nature of the polymer.In literature, these variations are attributed to the fact that C(t) decays very slowly compared to G(t), given the requirement of true decorrelation in space for C(t) to reach zero [79].Furthermore, in Figure 12a, the relaxation time decreases in case of OCP polymer chain from 92.7 ns to 29.9 ns as a function of temperature.However, the relaxation time of PMA shows a different trend wherein the values decrease from 74.9 ns to 28.6 ns in the temperature range of 313 K -353 K after which they increase to 40.4 ns at 373 K.This is a possible trend since the expansion of PMA chain with increasing temperature compensates for the additional thermal energy at 373 K, leading to slower dynamics in the solution.The veracity of the claim is supported by the argument that relaxation time is proportional to the physical quantity [η] × η s T according to the Zimm model [82].Here, [η] is the intrinsic viscosity of the polymer chain, η s is the solvent viscosity and T is the temperature.The same proportionality can be simplified in this case as follows: Figure 12b shows that the product in the right-hand side of Equation 25 decreases with increasing temperature in case of OCP.However, for PMA, the value increases beyond 353 K for the modelling strategy implemented here.This in turn could lead to the explanation of how OCP in PAO solvent acts merely as a thickener whereas PMA in the same solvent acts as a viscosity index improver (VII) [83,47].
Figure 13a shows the terminal relaxation time values as a function of shear rate.The values in the reduced mDPD units, and the corresponding values in the real units using the scaling parameters in Table 1 are depicted in the figure.In case of the OCP polymer chain, the τ term decreases with an increase in the shear rate at both 313 K and 373 K. Additionally, it is observed that beyond γ = 0.1, the τ term values drop significantly at both the temperature values.At the same time, PMA polymer chain also shows a decrease in τ term as a function of shear rate.The reduction in τ term values as a function of shear rate is evident from the Zimm model relationship as described above.However, the structural reason is explored in this work.By computing the ratio of the y-component of squared end-to-end distance (R ey ) of the polymer chains to the total squared end-to-end distance (R e ) such that i.e., the fraction of molecular stretching aligned in the direction of the shear (y-direction) is observed.As shown in Figure 13b, the S alignment ≈ 0.33 at zero shear rate.As the shear rate increases, the polymer chains become aligned in the direction of shear as depicted by S alignment values increasing from 0.45 to 0.9.This in turn leads to the drop in viscosity of the polymer and consequently, the relaxation time as well.

Salient features of the model
Simultaneous implementation of the SLLOD thermostat and the mDPD thermostat is functioning appropriately.The consistency in the zero-shear viscosity results computed via equilibrium and non-equilibrium simulations acts as proof of it.The mesoscale model of the liquids of interest is essentially that of a droplet having a diameter of 0.02 µm.
The simulation results reported here are intended to reproduce the major features such as molecular structure and rheology, as seen in atomistic simulation studies using a mesoscale model.Properties such as viscosity, relaxation time and first normal stress coefficient are reported for a shear rate range in real units in order to be compared with the atomistic-scale MD.The consistency of the results in the range, as seen in both scales, indicates their mild viscoelasticity.Increasing the polymeric concentration makes the non-Newtonian behaviour more pronounced, as shown in Figure 14, with the shear-thinning of 5.9% by weight polymer solution starting at a lower shear rate than the 0.1% polymer solution, an observation that is consistent in the real units as well.At the same time, it needs to be pointed out that the mesoscale model here shows a slower rate of shear thinning than the atomistic scale model.
The effect could emerge from the increase in structural ordering of the beads in mesoscale compared with the atomistic scale (see Figure S7 in the Supporting Information).So it will be vital to improve the model to account for these variations by optimising the excluded volume further.

Conclusions
The responses of terminal relaxation time, viscosity and first normal stress difference as a function of shear rates are vital to determine the distinct aspects of a fluid's non-Newtonian behaviour.To predict those, a framework for non-Newtonian dilute polymer solutions at mesoscale has been developed, capturing the structural differences of OCP and PMA polymers in PAO-2 solvent.The order of magnitude differences in storage and loss moduli between fluids simulated by both mDPD and DPD simulations establishes the high sensitivity of the model to identify viscoelastic behaviour.It is achieved by a combination of density-dependent conservative force parameters, and systematically-derived dissipative force  cut-off and angle stiffness energies.As a result, the reported relationships between the mDPD (reduced) units and real (SI) units make the model accessible for direct comparison with results from macroscopic experiments of industrially important viscoelastic liquids.By the application of the developed methodology, distinct aspects of the behaviour of the different solutions at mesoscales are demonstrated, such as the relationship of terminal relaxation time with respect to temperature.Furthermore, the model is able to predict the viscoelasticity introduced with the addition of polymer chains, even at a small concentration of 0.1% by weight, as shown by the computation of the first normal stress coefficient, and storage and loss moduli.

Figure 1 :
Figure 1: Schematic of the modelling workflow of the solvent for mDPD simulations.

Figure 2 :
Figure 2: Representative simulation boxes of PAO-2 solvent in the atomistic scale consisting of 170 molecules and in the mesoscale consisting of 10980 molecules at 313 K.

Figure 3 :
Figure 3: (a) Radius of gyration of OCP polymer chain as a function of chain length, N − 1 at T = 313 K. (b) The self-diffusion coefficient of the OCP polymer chain as a function of chain length, N at T = 313 K. (c) g 1 (t) as a function of time showcasing different regimes for different chain lengths.Here, the horizontal dashed lines show the value of R 2 g and the vertical dashed lines show the corresponding values of τ Z for different values of N .(d) h * as a function of chain length, N at a temperature, T = 313 K.

Figure 4 :
Figure 4: (a) Radius of gyration as a function of temperature, T for OCP and PMA solutions in atomistic and mDPD simulations.(b) Intramolecular bond distribution of OCP and PMA polymer chains at 313 K in the atomistic and mDPD simulations.The schematic here compares the equivalent bonded distance of the beads in mesoscale with that in atomistic scale.

Figure 5 :
Figure 5: Schmidt number of solvent, Sc solvent as a function of temperature, T .

Figure 6 :
Figure 6: Comparison of OCP polymer Schmidt number (Sc poly ) as a function of temperature: (a) Different concentrations simulated in atomistic MD (5.9% by weight) and mDPD (0.001% by weight).(b) Same concentration simulated in atomistic MD and mDPD (5.9% by weight).

Figure 7 :
Figure 7: Viscosity (η( γ)) as a function of shear rate ( γ) at (a) 313 K and (b) 373 K for pure solvent and polymer solution systems simulated via atomistic MD.The symbols with the error bars show the simulated dynamic viscosity values, η( γ) at different shear rates, γ.The lines depict the Carreau-model fittings.

Figure 8 :
Figure 8: Shear viscosity, η as a function of shear rate, γ for solvent molecules computed in the mesoscale simulations.The horizontal lines show the zero-shear viscosity computed using Green-Kubo method under equilibrium conditions and the dashed line shows the Carreau-model fit.

Figure 9 :
Figure 9: Shear viscosity, η as a function of shear rate, γ at 313 K and 373 K in the mesoscale simulations.The horizontal lines show the zero-shear viscosity computed using Green-Kubo method under equilibrium conditions.The dashed and dash-dotted lines show the Carreau-model fit.

Figure 12 :
Figure 12: (a) Terminal relaxation time, τ term as a function of temperature, T at zero shear rate for OCP and PMA solutions.(b) Sc poly D poly • Sc solvent D solvent as a function of temperature (refer to Equation25).

Figure 13 :
Figure 13: (a) Terminal relaxation time, τ term as a function of shear rate, γ at 313 K. (b) S alignment of polymer chains as a function of shear rate at 313 K.The data on the y-axis shows the S alignment values at zero shear rate.

Figure 14 :
Figure 14:The ratio of shear viscosity, η( γ)/η 0 as a function of shear rate, γ at 313 K for OCP polymer solution with polymer concentrations of 0.1 % and 5.9% by weight.

Table 1 :
The properties of the PAO-2 system in reduced units (mDPD units) and corresponding values in real units at 313 K. (Note: N a is the Avogadro number)

Table 2 :
Various parameters of mDPD simulations at 313 K in reduced units (mDPD units).