Dear All,
I have a question about using SLLOD equations to generate a shear flow. According to “Statistical Mechanics of Nonequilibrium Liquids” Evans&Morriss Eq. 6.39 and Tuckerman et.al J. Chem. Phy 106 5615-5621 Eq. 21, the acceleration of the particle has a term from the time derivative of the strain rate. This term will contribute when the strain rate is time-dependent, for example, the switch-on of a constant strain rate at t=0 or a oscillatory shear implied by fix deform wiggle. My question is in fix nvt/sllod code, it looks like there isn’t such the term when computing the velocity of the particle. The buildup of the flow with expecting strain rate is only because of the constraint from the thermostat. I modified the code to remove the thermostat, the expecting strain rate will not be formed in a few timesteps as SLLOD algorithm suggested. The evolution of the velocity profile is almost the same when using nvt/sllod or just using nvt with modified thermostat, this result indicate even I use nvt/sllod, the flow is still boundary driven other than field driven. So the nvt/sllod doesn’t make any difference when generating a shear flow.
Another question is if someone wants to implement a profile unbiased thermostat (compute temp/profile) on the SLLOD equations, the velocity computing in the code might be pragmatic at high shear rate. Since in the code, the momenta is peculiar to the streaming velocity computed by compute temp/profile, however, the momenta should only peculiar to flow you imposing not the outcome flow.
Any comment on these questions will be very helpful. Thanks in advance.
Best
Frank
I'm not clear on what you are asking. The relevant LAMMPS code is
in fix_nvt_sllod::nh_vt_temp(). Each atom's velocity is adjusted as:
temperature->remove_bias(i,v[i]);
vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
vdelu[2] = h_two[2]*v[i][2];
v[i][0] = v[i][0]*factor_eta - dthalf*vdelu[0];
v[i][1] = v[i][1]*factor_eta - dthalf*vdelu[1];
v[i][2] = v[i][2]*factor_eta - dthalf*vdelu[2];
temperature->restore_bias(i,v[i]);
h_two = Hrate * Hinv, so there is a strain-rate dependent term.
Are you saying there should be an extra term if the strain
rate is not constant?
The t = 0 case you mention is a special case, since when strain
is switched on, the derivative of the strain rate is momentarily
infinite.
The remove_bias() and restore_bias() calls are what effectively
bring the system rapidly to the streaming deformation velocity.
I'm CCing Pieter who did the initial SLLOD implementation in LAMMPS.
He may wish to comment.
Steve
Dear Steve,
Thanks a lot for your reply. Yes, I mean if the strain rate is time dependent, there should be a extra term when computing acceleration as a=F/m+xd(Del u)/dt, where a is the particle acceleration in laboratory frame and Del u is the strain rate. However, in v[i][0] = v[i][0]factor_eta - dthalfvdelu[0] for example, the sllod correction term "dthalfvdelu[0]" and remove_bias, restore_bias essentially only change the moving frame momenta back to the laboratory frame momenta. The time-dependent strain rate contribution is not included.
This is my understanding, not sure whether it’s correct. Eagerly need comments.
Best wishes,
Frank
The practical effect of removing the bias, thermostatting the remainder, and
restoring the bias is to rapidly bring the velocity profile of the fluid to be
equal to the deforming box profile. Initially, if the fluid isn't moving, then
when the bias is removed, the fluid is far from the desired (small) temperature,
so the thermostat makes a big, rapid change in the fluid, bringing it closer
to the streaming profile after the bias is added back in. This continues until
the fluid is basically at the streaming profile, plus a small thermal
noise around it.
So in tests I have seen, it does quite well.
As for the acceleration term, I am guessing it is normally a quite
small correction,
since you don't want to be using fix deform (at least with SLLOD) in a way that
induces large derivatives of the strain rate. In other words, you shouldn't be
trying to drive your fluid that rapidly else you will have other
problems and difficulty
measuing any kind of equilibrated property, e.g. viscosity.
Let's see if Pieter weighs in.
Steve
Dear Steve,
Thanks a lot you for your test.I have done exact the same test. My point is, in the current code,
the buildup of the streaming profile is because of the thermostat as you test. However, in
original SLLOD algorithm, even there isn’t a thermostat, the algorithm should still be able
to build a streaming velocity fast.
The test I have done is that after removing the thermostat by modifying the code, the SLLOD
equations (fix nvt/sllod) evolves the system almost the same as the common Newtonian equations
(fix nvt), which is not I expected.
I think the issue comes in the acceleration. In the SLLOD equations, the second equation “dp/dt=F-pDel u"
gives the acceleration of the particle. However, as lots of literature mentioned, the momenta here is
the peculiar momenta with respect to the streaming velocity. So, this acceleration is the acceleration
in the moving reference with the streaming velocity. If the acceleration is converting back to the laboratory
reference, the term "pDel u” is canceled and “q*d(Del u)/dt” appears.
My question is that in the velocity Verlet integrator, when computing the increment of the velocity, whether
the laboratory reference acceleration or the moving reference acceleration should be used.
Best
Frank