How to compute spatial resolved heat flux (per chunk) in LAMMPS?

LAMMPS LAMMPS General Discussion
1

Hai all,

I am working on shock waves, and I need the spatially resolved heat flux rate for calculating the entropy production.
The comment “compute heat/flux” generates a global vector output, whereas I need to calculate the spatial heat flux rate along the shock wave direction.
Is there any built-in way to calculate the heat flux rate over chunks (using chunk/atom or fix ave/chunk), or is the correct approach to reconstruct the energy current from per-atom quantities (ke/atom, pe/atom, stress/atom) and use “compute reduce/chunk” to sum it in each bin?

Thank you in advance.

2

Hi @ajputhoor,

Note that this only considers the case where you use the group all. You can define several sub regions in which you can apply a compute heat/flux each. However, the Green-Kubo formalism supposes equilibrium and takes time to converge. I’m not sure this is a correct method to apply to entropy generation via shock wave. I think only a search in the literature will tell you which method to use in such simulations.

1 Like
3

Thank you for the clarification. I have come across a paper where entropy production is evaluated using local heat flux and viscous shear stress, computed with an in-house Fortran code and validated against LAMMPS for a shock wave. I would like to know if there is any built-in way in LAMMPS to compute spatially resolved heat flux. I will also look into applying compute heat/flux to subregions. Thank you for your guidance.

4

Any feature available in LAMMPS is documented in the manual. If it is not mentioned, it is not available, and you have to either collect sufficient raw data and then post-process it, or add the feature you need as C++ code yourself. LAMMPS is designed to be easy to modify. In fact, a large part of LAMMPS is code contributed by LAMMPS users that needed a feature that was not yet available.

5

You can find the equation to calculate heat flux in the documentation (compute heat/flux). I ran some non-equilibirum simulations of Argon in a nanochannel and used the following method to get spatial heat flux (Heat flux per bin). The results were satisfactory and if applied for the whole simulation box, it gave me the exact same result given by the compute heat/flux command.

Note: I also used the compute heat/flux command to verify my code gives the same results. For my case, heat flux calculated using compute heat flux ( v_Jx) and the binning method (v_hfX) gave the exact same value.
Also, variable ‘Cx2’ is divided by 1.6/10^6 for correct unit conversion.

Heat Flux profile along X-direction

compute pe Argon pe/atom
compute SA Argon stress/atom NULL virial
compute flux Argon heat/flux ke pe SA

variable Jx equal c_flux[1]/${vol_ar}

variable Cx1 atom (c_ke+c_pe)*vx
compute Jcx Argon reduce sum v_Cx1

variable Cx2 atom -(c_SA[1]*vx+c_SA[4]*vy+c_SA[5]*vz)/1.6/10^6
compute Jvx Argon reduce sum v_Cx2

variable hfX equal (c_Jcx+c_Jvx)/{vol_ar} variable hf atom (v_Cx1+v_Cx2)/{binvol}

fix hf Argon ave/time {Nevery} {Nrepeat} {Nfreq} v_Jx v_hfX file heat.flux fix fluxprofile Argon ave/chunk {Nevery} {Nrepeat} {Nfreq} Bins v_hf norm none file heat.profile

Regards,
Wazih