Computing viscosity vs shear rate for shear-thinning fluid with deform option

User discussions
1

GROMACS version: 2023.4
GROMACS modification: No

Dear GROMACS community (especially @MichelePellegrino and @hess),

I am currently working on obtaining the viscosity vs. shear rate for a non-Newtonian fluid exhibiting shear-thinning behavior. I want to ensure that my approach follows best practices and can yield reproducible and reliable results.

Current Workflow

My workflow involves performing NVT MD simulations using the deform option to achieve specific shear rates, ranging from 1 s⁻¹ to 100 s⁻¹. My system is enclosed in a cubic simulation box with a side length of 5 nm. For example, to achieve a shear rate of 1 s⁻¹, I set deform = 0.0 0.0 0.0 5.0e-12 0.0 0.0, where the deformation velocity in the xy-plane is calculated using shear_rate * box_length. Velocities are initiated using the deform-init-flow option.

For obtaining the viscosity at each shear rate, I am following the method outlined in Shear viscosity coefficient of aqueous glycerol from equilibrium Molecular Dynamics simulations (Zenodo link). Specifically, I run 20 independent MD simulations starting from different configurations, taken close to the average density obtained from an NPT simulation (P = 1 atm). My system is an organic solvent with a salt, simulated at 300 K. Experimentally, the fluid’s viscosity is known to decrease from approximately 100 cP at low shear rates to around 5 cP at shear rates above 30 s⁻¹.

Questions & Concerns

  1. Is this workflow appropriate for obtaining the apparent viscosity vs. shear rate trend?
  2. Is running 20 independent simulations sufficient to get a reliable estimate, or should I aim for more, perhaps 60 simulations as suggested in the reference method?
  3. Regarding the shear stress integral vs. time:
    • In my simulations, the log-log plot of the shear stress integral reaches a linear regime only after long simulation times (e.g., >15 ns for low shear rates). Is this behavior expected for non-newtonian fluids?
    • Should I simply fit the linear portion to obtain the viscosity, or am I missing some critical consideration? My assumption is that the non-linear portion represents the disruption of the fluid’s local structure before reaching steady state, and as such should be not included in the fit.
  4. Could the size of the simulation box (5 nm) affect my ability to capture the correct shear-thinning behavior, or is it adequate for this study?

I have thoroughly reviewed the literature and forum posts, but I haven’t come across any step-by-step guidance regarding this topic, so I apologize if any of these questions seem trivial.

Additionally, I have experimented with the cosine acceleration method for calculating viscosity vs shear rate trend. However, I noticed that at increasing shear rates, the viscosity decreases linearly. From what I understand, and also according to literature in this regard, this method should primarily be used to estimate the zero-shear viscosity by extrapolating the linear trend back to zero shear rate for newtonian fluids. Can anyone confirm if this interpretation is correct?

Future Contribution

If I can establish a reliable workflow, I would be happy to contribute an in-depth tutorial on this topic with a practical example, complete with scripts for automating the necessary steps.

Any insights or advice would be greatly appreciated!

Thank you for your time and help,

Mattia

2

Hi Mattia, I’ll try to answer you questions the best I can.

For obtaining the viscosity at each shear rate, I am following the method outlined in Shear viscosity coefficient of aqueous glycerol from equilibrium Molecular Dynamics simulations (Zenodo link).

Why are you following the procedure for equilibrium simulations, when you are using the deform option? The equilibrium simulation workflow uses Einstein relations to compute viscosity. It doesn’t relate the shear rate to the off-diagonal stresses, but rather uses the integral of the fluctuations of the off-diagonal stress components to estimate the ‘zero-shear-rate’ viscosity. Hence, it cannot really be used to investigate non-Newtonian fluids. The non-equilibrium approaches (being cosine perturbation or deform) should in principle converge faster, so probably 60 replicates are not needed. But you should not use -eviscoi: with deform, just compute the off-diagonal stress in the deformation direction with gmx energy and relate it to the imposed deformation rate. This I hope answers both point 1. and 2.

In my simulations, the log-log plot of the shear stress integral reaches a linear regime only after long simulation times (e.g., >15 ns for low shear rates). Is this behavior expected for non-newtonian fluids?

If the liquid is decently viscous (like 100cP at zero shear rate), then I am not really surprised. But again, since you are using Einstein relations, I would say this is not related to any non-Newtonian behavior as far as I can tell.

Should I simply fit the linear portion to obtain the viscosity, or am I missing some critical consideration? My assumption is that the non-linear portion represents the disruption of the fluid’s local structure before reaching steady state, and as such should be not included in the fit.

Again, this relates to the Einstein/equilibrium-MD workflow, which is not the one you would want to use for deform. But the answer is yes: you need to fit the linear regime, ideally for time going to infinity (which is never really possible, not just due to the simulation length being finite, but most importantly because the variance of the integral also increases with time).

Could the size of the simulation box (5 nm) affect my ability to capture the correct shear-thinning behavior, or is it adequate for this study?

My gut feeling is that this depends on the liquid you are simulating. How many solvation cells do you have? 5 nm should be fine for most simple liquids, I would imagine.

Additionally, I have experimented with the cosine acceleration method for calculating viscosity vs shear rate trend. However, I noticed that at increasing shear rates, the viscosity decreases linearly. From what I understand, and also according to literature in this regard, this method should primarily be used to estimate the zero-shear viscosity by extrapolating the linear trend back to zero shear rate for newtonian fluids. Can anyone confirm if this interpretation is correct?

Yes and no, in the sense that you should probably see a decrase in viscosity for very high shear rates (at least thi is what I’ve also experienced), but non-equilibrium methods, like cosine or deform, are the ones to use to study the viscosity of non-Newtonian fluids, since you can actually explore the shear rate vs. viscosity relation.

If I can establish a reliable workflow, I would be happy to contribute an in-depth tutorial on this topic with a practical example, complete with scripts for automating the necessary steps.

That would be great! I was planning to make a tutorial on how to compute the shear viscosity coefficient for a LJ liquid, we can maybe communicate about that.

3

Hello Michele,

Thank you for your kind feedback!

Regarding the use of the procedure for equilibrium simulations, this mainly stems from my limited understanding of how molecular interactions contribute to rheological phenomena and the physical chemistry involved. After looking into this more and considering your pointers, I now realize that this is not the right method for use with NEMD.

Just to clarify, with NEMD and the deform option, all I need to do to get the viscosity is to track the average of the three shear stress off-diagonal components to assess when it stabilizes, and once it does, divide it by the shear rate. Is this correct?

For the box size, I agree that it should be large enough to capture the structures responsible for the rheological phenomena. In my specific case, I have the formation of a stable complex between the cation and the solvent, which occupies a sphere of approximately 10 Å. Therefore, a 5 nm box should be more than sufficient!

Regarding my comment on the cosine-acceleration method, what I’ve observed is the same phenomenon reported in this article: Zhao, L. et al. (2008), Journal of Chemical Physics, 129(14). They show that, even for common Newtonian solvents, using the cosine acceleration method leads to a linear decrease in viscosity instead of a constant value as expected. They justify this by observing a formation of uneven density waves in the liquid, causing some areas to be denser and others less dense, which apparently makes the fluid appear more “slippery” than it really is.

That would be great! I was planning to make a tutorial on how to compute the shear viscosity coefficient for an LJ liquid. We can maybe communicate about that.

I would be more than happy to contribute! I will try to confirm that, following your suggestions for NEMD simulation, I can obtain at least a somewhat physically consistent result. Then, I’ll reach out to your KTH address, if that works for you.

Thank you again for your precious feedback!

Kind regards,

Mattia

4

You’re welcome.

Only the one corresponding to the deformation direction (two actually, since the stress tensor is symmetric).