I’m trying to use Guo dynamics to integrate the force term calculated in ShanChen data processor. An easy modification was made in origin ShanChen data processor in Palabos. And I made a simulation to test this code.(Code is written below) A 3d droplet was placed in the center of simulation domain, with the density get by Maxwell construction. However, the simulation result seemed wrong. The radius of droplet become larger in the beginning stage (approximately 1000 time steps). Has anyone met this problem while using this method?
Here is the code for ShanChen data processor:
plint nx = domain.getNx() + 2; // Include a one-cell boundary
plint ny = domain.getNy() + 2; // Include a one-cell boundary
plint nz = domain.getNz() + 2; // Include a one-cell boundary
plint offsetX = domain.x0 - 1;
plint offsetY = domain.y0 - 1;
plint offsetZ = domain.z0 - 1;
ScalarField3D<T> psiField(nx, ny, nz);
// Compute density and momentum on every site and store result in external scalars;
// furthermore, evaluate the interaction potential Psi and store it into a ScalarField.
// Envelope cells are included, because they are needed to compute the interaction potential
// in the following. Note that the value of the momentum is stored temporarily only, as
// it is corrected later on to include corrections due to the interaction potential.
for (plint iX = domain.x0 - 1; iX <= domain.x1 + 1; ++iX) {
for (plint iY = domain.y0 - 1; iY <= domain.y1 + 1; ++iY) {
for (plint iZ = domain.z0 - 1; iZ <= domain.z1 + 1; ++iZ) {
// Get "intelligent" value of density through cell object, to account
// for the fact that the density value can be user-defined, for example
// on boundaries.
Cell<T, Descriptor> &cell = lattice.get(iX, iY, iZ);
T rho = cell.computeDensity();
// Evaluate potential function psi.
psiField.get(iX - offsetX, iY - offsetY, iZ - offsetZ) = Psi->compute(rho);
// Store density into the corresponding external scalar.
*cell.getExternal(densityOffset) = rho;
}
}
}
// Compute the interparticle forces, and store they by means of a
// velocity correction in the external velocity field.
for (plint iX = domain.x0; iX <= domain.x1; ++iX) {
for (plint iY = domain.y0; iY <= domain.y1; ++iY) {
for (plint iZ = domain.z0; iZ <= domain.z1; ++iZ) {
Array<T, D::d> rhoContribution;
rhoContribution.resetToZero();
// Compute the term \sum_i ( t_i psi(x+c_i,t) c_i )
for (plint iPop = 0; iPop < D::q; ++iPop) {
plint nextX = iX + D::c[iPop][0];
plint nextY = iY + D::c[iPop][1];
plint nextZ = iZ + D::c[iPop][2];
T psi = psiField.get(nextX - offsetX, nextY - offsetY, nextZ - offsetZ);
for (int iD = 0; iD < D::d; ++iD) {
rhoContribution[iD] += D::t[iPop] * psi * D::c[iPop][iD];
}
}
// Computation and storage of the final momentum, including tho momentum
// difference due to interaction potential and the external force.
Cell<T, Descriptor> &cell = lattice.get(iX, iY, iZ);
T *momentum = cell.getExternal(momentumOffset);
for (int iD = 0; iD < D::d; ++iD) {
// Initialize force contribution with force from external fields if there
// is any, or with zero otherwise.
T forceContribution = T();
// Add interaction term.
T psi = psiField.get(iX - offsetX, iY - offsetY, iZ - offsetZ);
forceContribution -= G * psi * rhoContribution[iD];
*(cell.getExternal(forceOffset) + iD) = forceContribution;
}
}
}
}
Here is the main code for simulation:
T lx = 1., ly = 1., lz = 1., gasRho = 0.032, liquidRho = 0.279, wallRho = 0.12, uMax = 0.1,
Cv = 1.;
plint resolution = 150;
Array<T, 3> bodyForce(0., 0., 0);
ShanChenParamsBase<T> Params(
lx, ly, lz, resolution, liquidRho, gasRho, wallRho, uMax, Cv, bodyForce);
plint nx = Params.getNx();
plint ny = Params.getNy();
plint nz = Params.getNz();
MultiBlockLattice3D<T, NSDESCRIPTOR> nsLattice(
nx, ny, nz, new GuoExternalForceBGKdynamics<T, NSDESCRIPTOR>(1.));
nsLattice.periodicity().toggleAll(true);
nsLattice.toggleInternalStatistics(false);
programSetup(nsLattice, Params);
for (plint iT = 0; iT < iMax+1; iT++) {
if (iT % frequency == 0) {
writeVTKs(nsLattice, iT, Params);
pcout << iT << endl;
}
nsLattice.collideAndStream();
}
}
ps: I have changed the functional: computeVelocity()