Thermodynamic Scaling of the Shear Viscosity of Lennard-Jones Chains of Variable Rigidity

: In this work, the thermodynamic scaling framework has been used to emphasize the limitations of fully ﬂexible coarse grained molecular models to yield shear viscosity of real liquids. In particular, extensive molecular dynamics simulations have conﬁrmed that, while being reasonable to describe the viscosity of short normal alkanes, the fully ﬂexible Lennard-Jones and Mie chains force ﬁelds are inadequate to capture the density dependence of shear viscosity of medium to long alkanes. In addition, it is shown that such a weakness in terms of coarse grained molecular models can be readily quantiﬁed by using the thermodynamic scaling framework. As a simple alternative to these force ﬁelds, it is demonstrated that the insertion of a variable intramolecular rigidity in the Lennard-Jones chains model exhibits promising results to model medium to long chain-like real ﬂuids from both thermodynamic and viscosity points of view.


Introduction
The knowledge of thermophysical properties of fluids in both dilute and dense phase is essential to understand and predict their behavior in many natural and industrial systems. For instance, the shear viscosity of liquids is a key element in many industrial applications [1,2]. As a result, a lot of research teams are interested in developing a rigorous theory which could describe thermodynamic and transport properties accurately and simultaneously from microscopic information. Despite significant progress [3][4][5][6], such a theory is not yet available apart from the dilute regime. As a matter of fact, current theories usually describe only one type of properties accurately and are often property dependent. Moreover, engineering tools allow for the prediction of well-known fluids in classical conditions satisfyingly, but they are often revealed to be inadapted when used to deal with "new" fluid types and/or complex thermodynamic conditions.
Recently, different theoretical frameworks have appeared to connect transport properties and thermodynamic properties in a single formalism, such as friction theory [7,8] and all the variants of excess entropy scaling [9][10][11][12][13] and isomorph theory [14][15][16]. In particular, the so-called Thermodynamic Scaling (TS) model reveals itself to be one of the most promising approaches to describe the shear viscosity of moderate to dense fluids, i.e., including the liquid state, under a large variety of thermodynamic conditions [17][18][19][20][21]. This approach relies on the assumption that the shear viscosity is related to the ratio of the density, raised at an exponent γ, divided by the temperature. It reduces the problem to the definition of γ, which can be connected to microscopic interaction characteristics of the studied systems [15,22].The TS framework has been extensively applied to monoatomic fluids, but is, in principle applicable to any molecular fluids. However, complex behaviors are known to occur in chain-like fluids [23][24][25]. In particular, γ is known to increase with the chain length when dealing with model fluids such as Lennard-Jones (LJ) fully flexible chains [10], whereas it seems to take relatively low values for real long chains molecules, such as in squalane [26].
Thus, the purpose of this work is to shed some light on the apparent contradiction in terms of TS of the viscosity of chain-like molecules. This would enable the improvement of the description of the shear viscosity of molecular fluids composed of long chains within the thermodynamic scaling framework and would allow the development of more efficient coarse grained force fields using the methodology described in Figure 1. To do so, extensive molecular dynamics simulations have been performed on Lennard-Jones chains of variable intramolecular "rigidity" and the results were analyzed relatively to the behavior noticed in normal alkanes. The article is structured as follows: Section 2 describes the different models and numerical methods employed. In Section 3, the results of our investigations are exposed and improvements are suggested to face the limitations of the fully flexible LJ chains model.

The Lennard-Jones Chains Fluid Model
The Lennard-Jones chains fluid model, which is used in this work, consists in describing a fluid molecule as a chain of N c tangent spheres with bond lengths constrained [27]. The non-bonded intramolecular and intermolecular interactions between spheres are described by the Lennard-Jones fluid potential [28], i.e., where ε is the potential strength, σ is the sphere "diameter", r is the center to center distance between the two considered spheres and r c is the cutoff radius (taken equal to 2.5σ in this work).
In addition, to include a variable intramolecular "rigidity" of the chains, a simple harmonic bending potential [29] has been introduced as: where θ 0 is the equilibrium angle (θ 0 = π) between three adjacent spheres and k is a stiffness coefficient. It is important to note that a fully flexible chain is described by k = 0 whereas a fully rigid chain is represented by k → ∞. Furthermore, the classical LJ dimensionless units are used in what follows: where ρ = N T V is the number density, N T is the total number of spheres in the simulation box (i.e., N T = N c * N mol , with N mol the number of molecules and N c the chain length), V is the volume of the simulation box, M is the mass of a sphere, P is the pressure, T is the temperature, k B is the Boltzmann constant, and η r is reduced viscosity introduced in [17].

Numerical Methods and Parameters
An extensive set of non-equilibrium molecular dynamics simulations [28][29][30] have been performed by means of the homemade code Transpore, already validated and employed in previous studies and fluid types [31][32][33]. The velocity Verlet algorithm and the RATTLE method [34] are used to integrate the equations of motion. Moreover, classical periodic boundaries with Verlet neighbor's lists are applied [29], and a Berendsen thermostat is used to maintain the aimed temperature [35]. The error bars have been calculated via the sub-block method [28].
The shear viscosity was computed for various thermodynamic states and estimated by using the NEMD method of Müller-Plathe [36], which consists in performing a momentum exchange between molecules in the central and the edge parts of the simulations domain in order to be compatible with periodic boundary conditions. More precisely, this method [31,36] is applied with a subdivision of the simulation box in N s = 24 slabs, and the momentum exchanges are carried out between the edge slabs, 1 and N s , and the central slabs, N s 2 and N s 2 + 1. The fluid velocity in momentum exchange slabs and their neighbors are discarded to calculate the shear rate. Then, the viscosity is estimated as the shear stress divided by the shear rate.
The shear viscosity of the LJ chains fluid has been computed thanks to numerous NEMD simulations on a specific amount of particles, with various rigidities (0 ≤ k * ≤ 600), chain lengths (1 ≤ N c ≤ 16), densities (0.6 ≤ ρ * ≤ 1.1), temperatures (1 ≤ T * ≤ 5) and, through a minimum of 1.5 × 10 7 non-equilibrium time steps at the steady state and a reduced time step δt * = 0.002. A minimum of 350 molecules is required to have sufficient molecules in each slab of the simulation domain. Hence, 1500 molecules are used for N c = 1 and 350 molecules are employed for N c = 16.
It is worth pointing out that an adapted exchange frequency must be chosen during the NEMD simulations to avoid shear thinning [31,37,38]. As a matter of fact, these quantities are closely linked until a threshold value from which the shear viscosity remains constant whatever the exchange frequency (the Newtonian linear-response regime). Hence, its estimation requires the performance of several simulations with the same parameters except for the exchange frequency. However, this approach results in being fastidious all the more so as the value threshold of the shear rate is liable to depend on the rigidity or length of the chains (see Figure 2). As an alternative, given the fact that shear thinning is related to local fluid reorganization, it can be interesting to calculate the first and second normal stresses differences : N 1 = τ xx − τ yy and N 2 = τ yy − τ zz which are equal to zero in an isotropic fluid. The normal stress τ ll is calculated from the following formula: where l = x (resp y, z), v * i,l is the l component of the velocity of particle i, r * ij,l is the l component of the center to center distance between particle i and j and F * ij,l is the l component of the interaction forces between particle i and j.  As exhibited in Figure 2, the shear thinning persists until a threshold value from which η r remains constant for all rigidities. Nevertheless, this threshold (Newtonian plateau) depends on the rigidity of the chain, namely 0.001, which corresponds to an exchange frequency of about 300. Thus, to prevent shear thinning, an exchange frequency of 500 has been chosen to produce data. Similarly, the results provided in Figure 3 show that the differences of the normal stresses decrease with the shear rate and become close to zero at values which are consistent with those estimated from Figure 2. These results clearly confirm the existence of a correlation between the shear thinning threshold and the non-null differences of the normal stresses. As a consequence, N 1 and N 2 have been computed during all simulations to assess the validity of the evaluation of η r .  Figure 3. Shear rate effect on the differences of the normal stresses for ρ * = 0.95, T * = 3, N c = 4. (a) Normalized differences of the normal stress N 1 versus shear rate, (b) Normalized differences of the normal stress N 2 versus shear rate.

Viscosity Modeling and Thermodynamic Scaling
In dense fluids (see [17]), like liquids, the physical quantity ρ γ T behaves as an invariant of the reduced shear viscosity η r i.e., η r = f ( ρ γ T ) (see [39] for more details) where f is an unknown function and γ is a parameter depending on the characteristics of the considered fluid. Indeed, the γ parameter is a measure of the relative impact of density and temperature on shear viscosity that is closely related to the repulsion shape between particles in monoatomic fluids [21]. However, as shown in [10,17], this thermodynamic scaling approach needs to be adapted to deal with low to moderate density conditions by considering the reduced residual shear viscosity, η r res = η r − η r 0 , instead of the reduced shear viscosity η r , where η r 0 is the reduced zero density viscosity which can be estimated from the relation provided in [32]. This can be understood because the shear viscosity results from the sum of two contributions. The first contribution, which describes the translations of the molecules, is equivalent to the zero density viscosity η 0 and can be evaluated from the kinetic theory [40]. The second contribution, which corresponds to the collisions and interactions between molecules, is the residual viscosity. Thus, the translation term prevails at low densities whereas the collision/interaction term is dominant at high densities, as in the system studied here.
To determine the γ parameters, the reduced residual viscosity has been correlated using the following empirical relation: where b 1 , b 2 , b 3 , b 4 , b 5 are adjusted numerical parameters and X = ρ * γ T * . It is important to note that Equation (9) ensures that the reduced residual viscosity tends to zero when X tends to zero as it should. Moreover, the parameters γ and b i are determined by minimizing the average absolute deviation between the estimated reduced residual viscosity in Equation (9) and the value deduced from the molecular dynamics simulations.

Fully Flexible LJ Chains Model for Alkanes
This subsection is devoted to the evaluation the ability of the fully flexible LJ chains model to describe the viscosity of normal alkanes within the TS framework. As previously shown by Galliero [10], this model seems to be a good compromise to describe the thermophysical properties of alkanes, in so far as it allows one to obtain a good estimation of the thermodynamic properties, whatever the chain length, and a good description of the shear viscosity for small chains. However, it seems to fail to represent the viscosity of long normal alkanes correctly, as already pointed out in [41]. To emphasize such a limitation, a simple metric is to look at the value of the γ parameter of the coarse grained molecular model relatively to that of the corresponding normal alkane in the TS framework. So, in a first step, the correlation between the shear viscosity η r res of the LJ chain fluid and the term ρ * γ T * has been verified. As exhibited in Figure 4, the adjusted parameter γ allows for the representation of η r res as a monotonously increasing function of ρ γ T , whatever the chain length. This confirms that the fully flexible LJ chains model allows for the identification of a thermodynamic scaling for the viscosity as η r res = f ( ρ γ T ), i.e., a master curve for each chain length, over a wide range of thermodynamic conditions covering liquid and dense supercritical phases.  In addition, the analysis of the behavior of the parameter γ with the chain length N c reveals a monotonously increasing evolution until reaching a threshold from N c = 4 (see Figure 5), which is consistent with previous results [10], even if this result is counterintuitive. Indeed, an increase in γ should correspond to a stronger repulsion between the molecule, as this parameter is a measure of the repulsive steepness potential between the fluid molecules [21]. So, for fully flexible LJ chains, one could have expected a softening of the repulsion with the increase of the chain length when N c is larger than 2. This unexpected behavior is probably related to a fundamentally different behavior, in the thermodynamic scaling framework, of atomic and polyatomic fluids. To better identify the qualitative difference between the TS of Lennard-Jones chains and real molecular fluids, our NEMD results are compared with the numerous data correlated in the N IST reference database [42] for several normal alkanes, namely : Methane, Heptane, Decane and n-dodecane, and with experimental data on: n-octodecane [43] and squalane [26]. For each compound, an efficient thermodynamic scaling has been found, i.e., the reduced residual viscosity can be represented by a function of ρ γ alkane T , which confirms once again the robustness of the thermodynamic scaling framework to deal with viscosities of dense fluids. The estimated γ alkane for the studied hydrocarbons and the γ Nemd evaluated from simulations are compared in Figure 6. In order to compare these results qualitatively, the equivalent LJ chain length of real fluids, N c , has to be defined as a LJ chain of a coarse grained model, i.e., one sphere corresponds to some real atoms. To do so, the equivalent LJ chain length has been estimated from the number of carbons, n c of a given alkane, by using the following often-used relation which can be found in the literature using the LJ chain force field to describe thermophysical properties [44]: As shown in Figure 6, γ Nemd and γ Alkanes for short chains not only have the same evolution with the chain length (monotonously increasing) but also exhibit quite similar values. Unfortunately, their behavior strongly differs from N c = 4, i.e., normal decane, since γ Nemd remains constant whereas γ Alkanes decreases with N c . As a consequence, γ Nemd differs more with γ Alkanes when N c increases.
These observations confirm that the fully flexible LJ model is adequate to obtain a reliable estimation of the viscosity for small chains but highlight strong limitations to describe the viscosity density dependence for longer chains, typically for normal decane and longer molecules. As an alternative to the fully flexible LJ chains model, one option is to use the fully flexible Mie chains model [33,45,46]. Such force field implies an additional molecular parameter (repulsion exponent) compared to the LJ chain force field but exhibits improved results on predicting thermophysical properties of many fluids [33]. Interestingly, it appears that a TS for the viscosity can be established when using the fully flexible Mie chains force field. However, the γ behavior along with N c remains similar to that of a fully flexible LJ chain (see Figure 7), i.e., monotonously increasing with the chain length whatever the repulsion exponent. This important result indicates that using a Mie potential instead of a LJ one, i.e., a more complex interaction potential, could not solve the viscosity density scaling problem for long chains. Another alternative could be to model long normal alkanes by partially rigid chains, as proposed by Galliero [41]. As a matter of fact, a fully flexible LJ chain model was used in order to try to scale the viscosity of the n-heptane and n-decane along vapor-liquid equilibrium line. However, strong viscosity deviations were observed for low temperatures in the liquid phase, when compared with experiments data. Hence, introducing a weak intramolecular rigidity in the fully LJ chains model enabled one to obtain an accurate estimation of the viscosity at low temperatures, validating the idea of introducing "rigidity" in the molecular model but only on the vapor-liquid equilibrium line. Using the thermodynamic scaling framework in a more systematic way, as done in the following, could allow for the confirmation of these preliminary results over a wider range of thermodynamic conditions and systems.

Semi-Rigid LJ Chains Model for Alkanes
The first aim of this subsection is to study the TS for the shear viscosity on semi-rigid LJ chains. The second goal is to identify the capability of such coarse grained force field to describe transport properties of medium and long alkanes satisfyingly. To do so, abundant NEMD simulations have been achieved on semi-rigid LJ chains fluids for several rigidities (from the fully flexible case, k * = 0, to a nearly fully rigid case, k * = 400) and for different thermodynamic conditions to estimate its shear viscosity. However, the phase diagram of such semi-rigid chains is somewhat more complex than the one of fully flexible chains, especially for the most rigid molecule. Indeed, the coexistence of isotropic and nematic phases [47] in some thermodynamic conditions, detected by the computation of the normal stresses differences in this work, has led to constrain the computations to a certain range of chains lengths when studying the influence of the rigidity on the parameter γ (N c = 3, N c = 4 and N c = 6). Figure 8 describes the reduced residual viscosity η r res along with the term ρ * γ T * for N c = 4 and for different rigidities. Interestingly, it is possible to find a function which correlates η r res with ρ * γ T * for all tested rigidities taken separately. As a consequence, for a given chain, a TS for the viscosity, and the corresponding γ value, can be established regardless of the rigidity. γ=6.7, k*=0 γ=5.4, k*=5 γ=5.0, k*=10 γ=5.2, k*=50 γ=5.9, k*=100 γ=6.9, k*=400 Having ensured that the TS applies to semi-rigid chains, the influence of the rigidity along with the parameter γ has been studied . It reveals a non-monotonous evolution, with a minimum around k * = 10, as exhibited in Figure 9. This point is extremely interesting as it emphasizes that increasing the rigidity of a chain may induce a decrease of the parameter γ over a certain range of values. Such a behavior is somewhat complex to analyze as one could have expected an increase in γ as increasing the internal rigidity of a molecule should lead to a steeper effective repulsion between the molecules. This confirms that the extension of the TS framework results analysis from atomic to polyatomic fluids is not straightforward and may reveal some subtle behaviors. In a second step, the behavior of the parameter γ has been evaluated on a weakly rigid chain (k * = 5), when its chain length increases. Interestingly, it appears that γ of a weakly rigid chain exhibits a non-monotonous evolution with chain length similar to that obtained on normal alkanes as shown in Figure 10, which describes the behavior of different γ (γ k * =0 , γ k * =5 and γ Alkanes ) with N c . This figure shows more particularly that γ k * =5 increases with N c until a maximum value around N c = 3 and then decreases. Thus, even if preliminary and only qualitative, all these results highlight the fact that a semi-rigid LJ chains model is more suitable to better represent, from the shear viscosity point of view, alkanes' medium and long chains compared to a fully flexible LJ chains model or even a Mie chains force field. These results also confirm that the way internal degrees of freedom are dealt with in coarse grained force field is crucial, especially if one is interested in transport properties such as viscosity.

Conclusions
In this work, the thermodynamic scaling framework has been applied to emphasize the limitations of fully flexible molecular models, often used as coarse grained force fields in molecular simulations, to provide shear viscosity of long real molecules such as normal alkanes. More precisely, our investigations based on extensive non-equilibrium molecular dynamics simulations have confirmed that the fully flexible LJ chains model is well suited to provide thermodynamic and transport properties of small normal alkanes. This is emphasized by a similar behavior of the parameters γ Alkane and γ Nemd , which characterize the density dependence of the viscosity in the TS framework, with the chain length until N c ≈ 4 or equivalently until n-decane. However, an increasing deviation has been observed between γ Alkane and γ Nemd along with N c for long chains. More precisely, the former decreases with the chain length above N c ≈ 4 whereas the latter increases with the chain length. Such a finding clearly indicates that the fully flexible Lennard-Jones chain force field is inadequate to describe accurately the shear viscosity of long alkanes in dense phases, including liquid state.
In addition, a more complex fully flexible force field, based on the Mie chains model, has been tested to overcome the weaknesses of the fully flexible LJ chains model, but it ends up having the same limits. This last result points out that the main limitations of these coarse grained force fields probably arise from their fully flexible nature, i.e., the way internal degrees of freedom are dealt with.
Thus, as another alternative in terms of coarse grained force field, while keeping a simple representation, a variable intramolecular "rigidity" has been inserted into the LJ chains model, i.e., the properties of semi-rigid LJ chains force field have been simulated with variable rigidities (from fully flexible to fully rigid). As a matter of fact, the results have revealed that this rigidity influences the behavior of the parameter γ to such a point that it evolves non-monotonously with the chain length and exhibits a similar behavior to that of γ Alkane . As a consequence, including rigidity in a coarse grained force field seems to be an appropriate and simple option to capture qualitatively the shear viscosity density dependence of normal alkanes, whatever their size.
All the previous results also confirm that the way in which internal degrees of freedom are included in coarse grained force fields is crucial to correctly model simultaneously equilibrium and transport properties of dense fluids composed of long chain-like real fluids.

Data Availability Statement:
The data generated for this study are available on request from the corresponding author.