Dear LAMMPS users,
This is a follow up question from a discussion about the temp/deform compute from a few weeks ago. I have seen the same phenomenon as described by the original questioner. I wanted to pick up on the point below that the absolute particle velocity v_i and the streaming velocity v_si are uncorrelated (note v_i = v_si + v_fluctuating). Is it right that this is only the case for homogeneous flow, i.e. simple shear in periodic domain?
Hi Ling,
I think this can be rationalised (but I’m not sure!).
Let v_i be the velocity of atom i, and v_si be its streaming velocity.
T_K = (alpha / N) * sum(v_i ** 2) = alpha * <v_i**2> (this is the
temperature c_myTemp)
T_D = (alpha / N) * sum([v_i - v_si] ** 2) = alpha * <(v_i - v_si)**2>
(this is the temperature c_stream)
Expanding the average,
T_D = alpha * <v_i2> + alpha * <v_si2> - 2 * alpha * <v_i * v_si>
If we then think about the difference T_D - T_K:
T_D - T_K = alpha * <v_si**2> - 2 * alpha * <v_i * v_si>
If v_i and v_si are uncorrelated (I think they are), then the second
term disappears, and
T_D - T_K = alpha * <v_si ** 2>, which is positive.
Chris,
The point is moot. T_D is the variance of the velocity, T_K is the second central moment E(x^2), which if you remember from stat class Var(x)=E(x^2)-E(x)^2, hence T_D-T_K=-alpha*<v_si>^2. If there is a streaming component, the difference between these two statistics is non-zero. Hope this helped.
Just for clarity, note that v_si is a constant = <v_i>.
Thanks for the reply, much appreciated - I follow your argument here, but now we have (T_D - T_K) < 0, which is what we might expect, but not what we (Ling+colleague) observe, and not what Niall demonstrated before.
The arguments there can be cleaned up as well, starting from equation 4 in Niall’s e-mail.
T_D - T_K = alpha * <v_si**2> - 2 * alpha * <v_i * v_si>
noting that v_si=<v_i>
=alpha(<v_i>^2-2<v_i><v_i>)
=-alpha(<v_i>^2)
Same result
I haven’t gone through their e-mails, but you can believe the math above.
To answer your question in your first e-mail, on physical grounds, the temperature is not typically strongly correlated to the streaming pieces given high enough temperatures and low enough Knudsen numbers (the collision integral, which models scattering from collisions, is larger than the advection term in the Boltzmann equation). They can, however, in certain systems(e.g. granular) become correlated through physical processes like viscous heating.
Ok thanks for clearing that up. Do you have an explanation, then, for the phenomenon described by Ling in the original post (fact that in simple shear, LAMMPS will give a higher temperature with compute temp/deform than with compute/temp)? Have you experienced this?
Sorry Chris, I post process all of my kinetic data with code I’ve written, as such I’m not intimately familiar with the computes. In addition, I’m a granular guy, so I don’t have a ton of knowledge on the thermostatting protocols as my ensembles are not of constant temperature.
One phenomena, of which I am familiar with, is structure effects on temperature calcs. For example, in simple shear of cohesive granular systems below a certain shear rate we see the formation of linear clusters. In such calcs, once the cluster length scale diverges to the box length, the velocity statistics are no longer Maxwellian and instead have long tails (commonly referred to as excess kurtosis). These systems behave like glasses, where movement is constrained and relaxation happens in jerky rearrangements. As such, the second moment, i.e. temperature, becomes unrepresentatively large because these are extremal statistics rather than the normal statistics we are used to seeing in more dilute fluids and crystalline solids. (Note that this secondary rise in temperature also happens in dense granular systems)
This is a digression, and I am certainly not indicating that this is what is occurring. There may be many reasons such calculations are different.
It’s simple. If the fluid is in fact not moving
with the box as it is deformed then compute temp
can give a small temperature (e.g. 0 if the atoms aren’t moving).
But compute temp/deform will subtract a large deformation
velocity and then compute the temp which will
thus be large. This is a misuse of compute temp/deform,
so it’s not surprising that it gives a bogus answer.
Steve
Thanks for the reply Steve - I agree with your point, but I have also checked the velocity profile of the particles to confirm that they are flowing in the way that I expect based on the fix deform. They are, and yet I still see this anomaly with the temperature calculations.
Chris,
Could you send me your input script. I’d like to look over for things like how long you’re running these, volume fraction, spring constants, etc. Since you’re doing granular, it could be something as simple as computing temperature before stationary statistics are achieved.
Hi Chris,
Yes, I found the same as you have stated: the velocity profile was linear in my case. But later I changed my definition for temperature in the shearing case, i.e., I used rather compute temp/partial than compute temp/deform, to thermostat only the non-streaming velocity direction.
As I am not using a global thermostat, so this definition works for my system.
But as for explaining the higher temperature than it is supposed to be, I still have no clue yet. I am not sure if resorting to another way of defining the temperature would be suitable for your case.
I would like to know the reason as well.
Regards,
Ling
Hi Ling,
My concern with using temp/partial to remove a streaming velocity is that you will also be neglecting the fluctuating velocity that exists in the direction that you are neglecting.
For example, if you remove a streaming velocity in the z direction using:
compute newT flow temp/partial 1 1 0
I think your temperature calculation, or thermostating (if applicable), will also be neglecting the z-component of the fluctuating velocity, which may or may not be appropriate for your case.
Also, try outputting your particle velocities in the dump file and calculating the “temp/deform” yourself by removing the streaming velocity based on the linear velocity profile you obtain - I think you will find the temp/deform to be lower than the full temperature, as you would expect.
The problem you are seeing in LAMMPS may be to do with the way the deformation rate is used to remove a streaming velocity from the particles. Suppose you have a 10x10x10 box, and you shear using fix deform at a shear rate of 1. You probably find a velocity profile that is linear from -5 at the “bottom” of the box, 0 at the box centre, and 5 at the “top” of the box.
When deforming the box, however, you can see from the box bound data in the dump file that the deformation happens at the “top” of the box, meaning the “top” surface will have a “velocity” of 10, and the “bottom” surface will have a “velocity” of 0. So I think the linear velocity used by temp/deform is this one (going from 0 to 10), while the actual velocity you observe goes from -5 to 5. So if you subtract one from the other, you will not end up with just the fluctuating part of the velocity, you will get something completely different.
This would explain the problem you are seeing, but I am not sure I am correct - I hope someone will be able to rationalise this.
Chris
Chris,
I understand that you are going coarse. But its probably important to address actual systems you will simulate if there truly is a problem. You may require more time but you can rule out any number of things that could go wrong for this simple system due to poorly selected parameters. You have the advantage that this system has been simulated by many people in the granular community using many different simulation engines.
Chris and Lin appear to be seeing the same thing.
Namely a system flowing with the fix deform profile,
but compute temp/deform is giving a higher temp
than compute temp. I don’t think that is possible
(unless there is a bug).
I assume you are not using fix deform remap x,
which as the compute temp/deform doc page
explains is a mistake (and will give a WARNING):
IMPORTANT NOTE: Fix deform has an option for remapping either atom coordinates or velocities to the changing simulation box. When using this compute in conjunction with a deforming box, fix deform should NOT remap atom positions, but rather should let atoms respond to the changing box by adjusting their own velocities (or let fix deform remap the atom velocities, see it’s remap option). If fix deform does remap atom positions, then they appear to move with the box but their velocity is not changed, and thus they do NOT have the streaming velocity assumed by this compute. LAMMPS will warn you if fix deform is defined and its remap setting is not consistent with this compute.
It is also possible you “see” the atoms flowing
but they do not have the same profile as
the deforming box, in which case compute temp/deform
will give anomalous answers (as will compute temp).
You can try compute temp/profile which should give
the same answer as compute temp/deform, if (and only
if) the flow coincides with the box deformation.
Failing all that, I suggest you define a group
of 5-10 atoms, dump their velocities, and use
comptue temp/deform and compute temp on
that group. Then verify by hand that the 2
computes are giving correct answers.
Steve
Steve,
Thanks for the reply.
Yes I always use remap v
Using temp/profile, I obtain the same result for temperature/pressure/KE as when I dump the velocity components myself, remove streaming velocity and then calculate in matlab.
Therefore the problem is that the true flow velocity does not correspond to the box deformation. As I tried to explain previously, for a strain rate of say 10, the true velocity profile will be linear from -5 to +5, whereas the box deformation rate will be linear from 0 to 10. So when temp/deform subtracts the deformation rate from the true velocity profile, we are not left with the fluctuating part about a zero mean.
I think this is the source of the problem.
Chris
As I tried to explain previously, for a strain rate of say 10, the true velocity profile will be linear from -5 to +5, whereas the box
deformation rate will be linear from 0 to 10. So when temp/deform subtracts the deformation rate from the true velocity profile, we >are not left with the fluctuating part about a zero mean.
I missed that. How are you getting an actual flow from -5 to +5 when you use
fix deform which, as you say, should induce 0 to 10 ??
Steve
Chris,
Have you subtracted the mean streaming velocity from these particles to obtain -5 to +5. I’ve run many homogeneous shear problems with gran/hooke/history. I never obtain what you are describing.
Eric, Steve,
I have not done any post-processing to obtain the -5 to 5 velocity profile.
Here is an input script that demonstrates the difference between temp, temp/deform and temp/profile. Please see the data output from the fix ave/spatial, or even look at VX in the dump file, to see the actual velocity profile going from -5 to 5.
atom_style sphere
boundary p p p
newton off
communicate single vel yes
region reg prism 0 10 0 10 0 10 0 0 0 units box
create_box 1 reg
lattice sc 1
create_atoms 1 box
neighbor 0.2 bin
neigh_modify delay 0
pair_style gran/hooke/history 20000 NULL 20 NULL 0.5 1
pair_coeff * *
timestep 0.00011374
velocity all create 10 4928459 dist gaussian
fix 1 all nve/sphere
thermo_style custom step
thermo 50
run 100000
velocity all zero linear
velocity all zero angular
compute 1 all temp
compute 2 all temp/deform
compute 3 all temp/profile 1 0 0 y 10
thermo_style custom step c_1 c_2 c_3
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
fix 2 all deform 1 xy erate 1 remap v
dump id all custom 440 A1.dump id vx vy vz omegax omegay omegaz y
thermo 440
fix 3 all ave/spatial 1 100 440 y lower 0.05 vx units reduced file A1.velocityprofile
run 44000
There is something odd going on in your input script, this system should not be doing this. There is probably something glaringly obvious. I am getting the same thing on my 2 year old version.
Chris,
I’m not sure what you’ve stumbled upon. I toyed around with your script making a few changes for consistency with my older scripts, it is attached. For some reason if you give your box an initial velocity gradient, it grows to match the shear rate symmetrically about the initial gradient. E.g. If I give the initial velocity difference at the bottom and top to be from 0 to 5, the ending difference is -2.5 and 7.5. If you give it to be 0 and 10 the system behaves as I would expect it to.
in.ChrisTest (1.18 KB)