In Sarkar’s article “Thermal Conductivity Computation of Nanofluids by Equilibrium Molecular Dynamics Simulation: Nanoparticle Loading and Temperature Effect” as attached, he said "We calculated the average enthalpy as the sum of the average kinetic energy, potential energy and average virial terms per particle of each species."
So I add the enthalpy term into “compute_heat_flux.cpp” like this:
Am I right? thanks
*****************************************************************
/ ---------------------------------------------------------------------- */
void ComputeHeatFlux::compute_vector()
{
invoked_vector = update->ntimestep;
// invoke 3 computes if they haven’t been already
if (!(c_ke->invoked_flag & INVOKED_PERATOM)) {
c_ke->compute_peratom();
c_ke->invoked_flag |= INVOKED_PERATOM;
}
if (!(c_pe->invoked_flag & INVOKED_PERATOM)) {
c_pe->compute_peratom();
c_pe->invoked_flag |= INVOKED_PERATOM;
}
if (!(c_stress->invoked_flag & INVOKED_PERATOM)) {
c_stress->compute_peratom();
c_stress->invoked_flag |= INVOKED_PERATOM;
}
// heat flux vector = jc[3] + jv[3] + jh[3]
// jc[3] = convective portion of heat flux = sum_i (ke_i + pe_i) v_i[3]
// jv[3] = virial portion of heat flux = sum_i (stress_tensor_i . v_i[3])
// jh[3] = enthalpy portion of heat flux
// normalization by volume is not included
double ke = c_ke->vector_atom;
double pe = c_pe->vector_atom;
double **stress = c_stress->array_atom;
double **v = atom->v;
int mask = atom->mask;
int nlocal = atom->nlocal;
int type = atom->type;
int ntypes = atom->ntypes;
int itype;
double jc[3] = {0.0,0.0,0.0};
double jv[3] = {0.0,0.0,0.0};
double jh[3] = {0.0,0.0,0.0};
double eng;
int nn[ntypes+1]; // total atom number for each type
double hh[ntypes+1]; // total enthalpy for each type
int data1[ntypes+1]; // aim to be used to compute total atom number for each type
double data2[ntypes+1]; // aim to be used to compute total enthalpy for each type
for (itype = 1; itype <= ntypes; itype++) {
data1[itype] = 0;
data2[itype] = 0.0;
}
// convert jv from stressvolume to energy units via nktv2p factor
double nktv2p = force->nktv2p;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
itype = type[i];
data1[itype] += 1;
eng = pe[i] + ke[i];
jc[0] += engv[i][0];
jc[1] += engv[i][1];
jc[2] += engv[i][2];
jv[0] -= (stress[i][0]*v[i][0] + stress[i][3]*v[i][1] + stress[i][4]*v[i][2]) / nktv2p;
jv[1] -= (stress[i][3]*v[i][0] + stress[i][1]*v[i][1] + stress[i][5]*v[i][2]) / nktv2p;
jv[2] -= (stress[i][4]*v[i][0] + stress[i][5]*v[i][1] + stress[i][2]*v[i][2]) / nktv2p;
data2[itype] -= (stress[i][0] + stress[i][1] + stress[i][2]) / nktv2p;
data2[itype] += eng;
}
}
// sum across all procs
MPI_Allreduce(data1,nn,ntypes+1,MPI_INT,MPI_SUM,world);
MPI_Allreduce(data2,hh,ntypes+1,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
itype = type[i];
jh[0] -= v[i][0]*hh[itype]/nn[itype];
jh[1] -= v[i][1]*hh[itype]/nn[itype];
jh[2] -= v[i][2]*hh[itype]/nn[itype];
}
}
// sum across all procs
// 1st 3 terms are total heat flux excluding enthalpy portion
// 2nd 3 terms are just conductive portion
// 3nd 3 terms are total heat flux including enthalpy portion
double data[9] = {jc[0]+jv[0],jc[1]+jv[1],jc[2]+jv[2],jc[0],jc[1],jc[2],jc[0]+jv[0]+jh[0],jc[1]+jv[1]+jh[1],jc[2]+jv[2]+jh[2]};
MPI_Allreduce(data,vector,9,MPI_DOUBLE,MPI_SUM,world);
}
Xiaoliang
|