I noticed that in the ExternalForceDynamics.hh file there is something that may be wrong. Consider the following line :
00283 Array<T,Descriptor::d> jCorrected(j+invOmega*Descriptor::fullRho(rhoBar) * force);
Actually, in the Shanchen multiphase model, this class is not used but due to some reasons, I wanted to utilize it. As you see,a function of density is multiplied to force, but it seems incorrect. Let’s consider the other part of this file. In the line 290 velocity is calculated as follows:
uSqr += util::sqr(j[iD]invRho + 0.5force[iD]);
Which there is not any density beside force term. In fact, there is a contradiction between these two lines. Actually, this line (290) is correct and line 283 should be rewritten to
Array<T,Descriptor::d> jCorrected(j+invOmega * force);
In addition, to apply the EDM model, I used this class which with the default code it diverges, but when I make change as above. The results are pretty well.
Please tell me if there is any special consideration that I do not know.