regarding the computation of global pressure during NEMD simulation

LAMMPS LAMMPS Mailing List Mirror
1

Dear steve,
As some part of this question is not related to lammps technical question, so please allow me for such post.

Dear all,

If one system have C-O-M velocity ( as the consequence of shock generation method,basically momentum mirror method ) how can I calculate global pressure, where lammps thermo variable and compute stress/atom both would be affected by com velocity in the kinetic energy contribution term. Is there any idea how to do that ?

From the previous discussions regarding that by Ray and Oscar G, I got some idea basically how to compute stress locally that by using the velocity (after subtracting of the com velocity ), it was like-
Stress_xx = [ Sum(Force_xx of i) / (LyLz) ] + [ Sum(mass_i * Vx* Vx) / (LxLyLz)]. Here my question is, as a shock phenomenon is beyond hydro-static case so there xx, yy and zz comp of the stress tensor should differ significantly, so if I do average it out, would it be safe to say that this is the pressure, locally for a region ? From there, would it make sense if I do the above computation assuming the whole system as one group ?

I’ve done the above calculation for a local pressure (assuming a region as a group ) but I doubt about the result, may be I have done some wrong in implementation within lammps, here it is what I’ve done —
compute “property/atom fx” , “compute reduce/region fx” and used “vcm()” to calculate com velocity of a region. The Lx, Ly,Ly and mass was computed for a region as readily available in lammps. Is there anything which I have meshed up ?

One thing regarding shock physics:

If I calculate P from momentum conservation i.e. P = density_unshocked * Us*Up so it shocked signify that it would be the pressure at the shock-precursor/front, is it wrong ? If so it should be well merged with the computed local pressure ( by the above calculation) just behind the shock front if I safely track the shock front by Mott-Smith density profile or by any means, is it so ?

If I wish to study the (P Vs T), or (P Vs density) profile so if I compute this variables by averaging over a fixed region, spatially averaged fixed volume /fixed no of particles bins, in a defined region whatever may be my sample length along shock, would it be correct ?

Thanks a lot

Rajdeep Behera
IOP, Bhubneswar

2

Hi,

Your post is really long and i'm not really a good reader. It will
take me time to digest the information and give an answer. However,
this discussion thread gives me a great opportunity to Post My
Tutorial about Shock/quasi-isentropic/uniaxial compression using
LAMMPS (please look at the attachment). The tutorial still needs a
lot of refinement , but It may be useful and a fine addition to the
lammps documentation and the SHOCK package?

The Tutorial is basically the recompilation of the ideas and feedback
obtained from the LAMMPS's community . I just tried to assemble and
organize the ideas into and Script. Hopefully people will continue to
post ideas , and I will keep adding up information .

I wanted to wait a little longer to Post this tutorial , because there
are some thing I still need to review and some basic concepts are not
clear to me e.g the center of mass substraction,etc

A salute
Oscar Guerrero

PS: @Rajdeep , please tell me if you find it useful and provide feedback

main.pdf (853 KB)

3

Dear Oscar,

Firstly, it’s a great enthusiastic effort from you, thanks a lot for that. Specially for beginners, like me this tutorial would be a great help to understand the capability and usage of different fixes in Lammps.

Now, I’m going through the whole and initially it seems to me that you have done one great job. After having a total look, I would come to you, if any thing you miss.

The Shock portion I read till now, it seems to me that you have used flyer plate method so you need not to over come ice-cube problem, is it so ?
Then the figures you added, those figures are centro-symmetry analysis results, and that you have done in AtomEye ?

Again a special thanks

4

Rajdeep,

I will try my best to answer your questions below, but anyone please feel free to chime in and/or correct me if I am wrong.

Thanks,
Ray

Dear steve,
As some part of this question is not related to lammps technical question, so please allow me for such post.

Dear all,

If one system have C-O-M velocity ( as the consequence of shock generation method,basically momentum mirror method ) how can I calculate global pressure, where lammps thermo variable and compute stress/atom both would be affected by com velocity in the kinetic energy contribution term. Is there any idea how to do that ?

Pressure, in thermodynamics, is a state function, which describes the equilibrium state of a system. Therefore a global pressure in non-equilibrium MD simulation is not meaningful.

From the previous discussions regarding that by Ray and Oscar G, I got some idea basically how to compute stress locally that by using the velocity (after subtracting of the com velocity ), it was like-
Stress_xx = [ Sum(Force_xx of i) / (LyLz) ] + [ Sum(mass_i * Vx* Vx) / (LxLyLz)]. Here my question is, as a shock phenomenon is beyond hydro-static case so there xx, yy and zz comp of the stress tensor should differ significantly, so if I do average it out, would it be safe to say that this is the pressure, locally for a region ? From there, would it make sense if I do the above computation assuming the whole system as one group ?

I have to say no and no. No to the first question since pressure is also only meaningful when the system is under hydrostatic. No to the second question because the system is not in equilibrium, so defining the whole system as one group is unreasonable. Under these circumstances, it is only reasonable to define local stress tensors.

I’ve done the above calculation for a local pressure (assuming a region as a group ) but I doubt about the result, may be I have done some wrong in implementation within lammps, here it is what I’ve done —
compute “property/atom fx” , “compute reduce/region fx” and used “vcm()” to calculate com velocity of a region. The Lx, Ly,Ly and mass was computed for a region as readily available in lammps. Is there anything which I have meshed up ?

Please note that Vx in the above equation is (V_i - V_COM). You may also try compute stress/atom suggested by Oscar in his nice Intro. I personally calculate stress tensors using post processing codes I wrote to analyze dump files, so I am not able to give you suggestions on how to calculate these using a series of LAMMPS commands.

One thing regarding shock physics:

If I calculate P from momentum conservation i.e. P = density_unshocked * Us*Up so it shocked signify that it would be the pressure at the shock-precursor/front, is it wrong ? If so it should be well merged with the computed local pressure ( by the above calculation) just behind the shock front if I safely track the shock front by Mott-Smith density profile or by any means, is it so ?

If you wish to estimate the pressure of a shocked system by using P = rho_0 * U_S * U_P, then you have one more condition to satisfy, in addition to 1) thermodynamic equilibrium, 2) hydrostatic. That is the system has to satisfy the Rankine-Hugoniot condition. Therefore this P can be estimated only when the shocked-system is in steady-state. For a region that is far, or far far far, behind the shock front, you may be able to compare the P estimated by this equation and the local stress tensors, provided the conditions are satisfied.

If I wish to study the (P Vs T), or (P Vs density) profile so if I compute this variables by averaging over a fixed region, spatially averaged fixed volume /fixed no of particles bins, in a defined region whatever may be my sample length along shock, would it be correct ?

Conceptually correct, although P is stress not pressure.

If you wish to obtain the P-T or P-V (same as P-rho) Hugoniot planes, then you may check out several Hugoniot state fixes, such as fix nphug and fix msst.

Cheers,
Ray

5

Dear Oscar & Ray,
In case of shock, you have used wall/piston so as per my experience it should give the whole system some translation due to center-of-mass velocity as well. And you have used stress/atom directly to get pressure but this pressure should contain com velocity, it should be removed, that was my quary today. However, will it be affect the whole if one use wall/reflect instead of wall/piston ?

6

I don’t use compute stress/atom so I can’t your first Q.

I personally use velocity set to the whole system and fix wall/reflect - so called momentum mirror.

But the above method is just the same as a flyer plate. Once the flyer plate impacts the system and the system starts to drift along the impact direction, you have COM that has to be removed when estimating local properties.

Ray

7

I have checked stress/atom to get local stress within Lammps, but it significantly differs from the experimental and as well as existing literature. Oscar, have you tally the results you get from this wall/piston and stress/atom with the expected value from literature, please make some comments about your experience.

@Ray,
The existing literature is very confusing regarding this pressure, or you may say local stress. There is many many of them have calculated averaged stress tensor and called them “pressure” and not only that they have shown hugoniot curves taking this “averaged stress (may be with com velocity)” as “pressure”. I also agree with you com velocity should be subtracted. Ray would you be so pleased to share your post-processing codes to calculate this local stress, so that we can understand what we have missed during this calculation. Because in this Stress_xx = [ Sum(Force_xx of i) / (LyLz) ] + [ Sum(mass_i * Vx* Vx)

I have subtracted com velocity from computed Vx, in this way

compute vx all property/atom vx
variable vxorgn atom c_vx-(gmask(reg1)*vcm(reg1,x))
so it should be correct, is there any wrong ?

If you share the code it would be a great help.

Thanks a lot

8

I agree with all of Ray’s helpful comments. I will add a slight qualification to the answer about P = rho_0 * U_S * U_P. This will never be exactly true, either in a NEMD shock simulation or in a shock experiment, but it is often a useful approximation. The biggest issue is non-hydrostaticity. There are obvious reasons why that will not be satisfied in some cases (elastic compression, strength). But, for strong shocks where Pzz greatly exceeds both the hugoniot elastic limit and the material strength, it is a good approximation. Independent of hydrostaticity, we can consider instead the relation:

Pzz = rho_0 * U_S * u_p

This is satisfied exactly when the wave is steady or the thickness of the transient region vanishes. In many cases, both of these things are good approximations. In summary, I think the above expression is an excellent first estimate of pressure. At the very least, it provides something to check other more complicated and error-prone estimates against.

Aidan

9

Hi,

Thank you for your valuable comments. I hope to digest this information and I am really going to take some time to read this thread. I will try to modify my tutorial and make some updates.

A salute
Oscar Guerrero