How to change to density ratio in the example of rayleigh-taylor instability

I found that there is an example of multi-component example in palabos, which solve the rayleigh-taylor instability.
In that code the initial density of heavy fluid and light fluid are both set to be 1.0.
Does this means the density ratio is 1? I wonder how to change the density ratio?
can I just change the initial density of the light fluid to be 0.5 to obtain the density ratio 2?


you could apply different external body forces to different fluids (and I’m sure that the example you are looking at does the same thing). Doing that, you would “reproduce” the effect of different densities between the two fluids. For example you could apply a zero body force to the lighter fluid (initially placed at the bottom of your domain) and a non-zero force “looking” downward for the heavy fluid (initially place above). Now you can play with different values of this body force in order to get the dimensionless number in which you are interested about.

I hope it helps


In the original setup of the Rayleigh-Taylor (2D), density of Heavyfluid is 1 at the upper half of the model and almostzero at the lower half. Lightfluid density is has an exact opposite configuration.
what is the significance of the Almostzero density in the simulation results?
Are we assuming that the position of these lattices ( Heavyfluid and Lightfluid) are constant throughout the simulation? (I guess we don’t). But, what is the criteria for changing the configuration of these lattices?
I looked at the source code and this is the implementation of Shan-Chen Model in Palabos.

 for (plint iSpecies=0; iSpecies<numSpecies; ++iSpecies) {
                Cell<T,Descriptor>& cell = lattices[iSpecies]->get(iX,iY);
                T *momentum = cell.getExternal(momentumOffset);
                for (int iD = 0; iD < D::d; ++iD) {
                    momentum[iD] = uTot[iD];
                    // Initialize force contribution with force from external fields if there
                    //   is any, or with zero otherwise.
                    T forceContribution = getExternalForceComponent(cell, iD);
                    // Then, add a contribution from the potential of all other species.
                    for (plint iPartnerSpecies=0; iPartnerSpecies<numSpecies; ++iPartnerSpecies) {
                        if (iPartnerSpecies != iSpecies) {
                            forceContribution -= G * rhoContribution[iPartnerSpecies][iD];
                    momentum[iD] += invOmega[iSpecies]*forceContribution;
                    // Multiply by rho to convert from velocity to momentum.
                    momentum[iD] *= *cell.getExternal(densityOffset);

Theoretically, value of G is non-zero where we are located close to an interface between different phases; however, in this code,

 if (iPartnerSpecies != iSpecies) {
                            forceContribution -= G * rhoContribution[iPartnerSpecies][iD];

for many nodes, value of G could be non-zero!

I would appreciate it if you could help me with these questions.


I guess I found the answer to the question about the implementation of ShanChen on Palabos.
In order to handle non-zero values of G close to the interface, they defined regions with very low density (Almostzero density in Rayleigh-Taylor.cpp) which causes

forceContribution -= G * rhoContribution[iPartnerSpecies][iD];

to give values close to zero for forceContribution at points away from the interface.