Non-equilibrium initial conditions

I am trying to initialize my PDF fields as shown in Skordos 1993[/url] and [url=http://www.springerlink.com/content/p832u0u402761261/]Caiazzo 2005, with the correct shear rates.
In particular I am trying to recreate the Taylor-Green vortex presented as a test case for initial conditions in Skordos '93.
However, my errors are not decreasing as they should. In fact I get larger errors than with the pure equilibrium initialization.
I wrote a D2Q9 Matlab code, and the perturbation from the shear stresses is calculated by a function,


function [f1] = CalcFirstOrderPerturbation(dxux, dxuy, dyux, dyuy)
    for dir=1:9
        f1(:,:,dir) = -tau*3*weight(dir) * ...
            (cx(dir)^2 * dxux + cy(dir)^2 * dyuy + ...
            cx(dir)*cy(dir)*(dxuy + dyux) );
    end
end

where dxux is the partial derivative d[sub]x[/sub]u[sub]x[/sub] and so on…
The function was working great in a 2D-Couette Flow where you have only dyux and all other derivatives vanish. In this more complicated case it does not seem to work at all.
I would be grateful for any hints!

Hi,

Your formula looks correct to me (except that your off-equilibrium population should additionally be multiplied by the density).

Do you also compute the initial value of the pressure, and initialize the density accordingly, to something else than just a constant? This is crucial, and has a much larger impact than the off-equilibrium populations.

Cheers,
Jonas

Hi, Jonas!
Thank you for the quick reply!
Yes, for the Taylor-Green scenario I am initializing the density according to the analytic pressure, like


Rho = 1 + 1/csqr * dt^2 / dx^2 *P;

where P is the pressure field for the vortex.


I am now getting the expected results! Multiplying the non-equilibrium initialization by the density as you suggested, however, seems to spoil my results: Actually, if I do that, I get the same results as in the case of the pure equilibrium initial conditions.