Hello Tim

First of all thank you for your interest in my problem.

The way I have defined my numerical error was as the sum (over all computational nodes) of the absolute value of the difference between the numerical solution and the analytical one over the absolute value of the sum of the analytical solution, i.e. I used a relative error with the conventional L1 norm definition, as it was used for instance for He and Luo in the paper: “Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation” (1997).

The results that I am reporting were obtained with the model that considers the body force to be Fi=wi*ei*F/cs^2. My implementation with the Guo model for the forcing term yields equivalent results…

So as you suggested me, here it is the relation between tau and the error for Re=10, LY=30 (channel width) and LX=15 (channel length) at steady state is as follows

tau Error Relat

0.51 0.0019225

0.55 0.0019304

0.60000 0.0018953

0.70000 0.0016857

0.80000 1.29450E-03

0.90000 7.29980E-04

0.99000 8.02250E-05

0.99900 8.09370E-06

0.99990 8.10080E-07

0.99999 8.10E-08

1.00000 5.32550E-14

1.00001 8.10170E-08

1.00010 8.10240E-07

1.00100 8.10950E-06

1.01000 8.18040E-05

1.10000 8.87810E-04

1.50000 5.86270E-03

1.70000 9.08980E-03

2.00000 1.46570E-02

2.50000 2.52710E-02

3.00000 3.66480E-02

The error increases as tau value departs from 1 (with the exception of tau=0.51).

Keeping tau=1 and still with LY=30 and LX=15 the steady state solutions yield the following error values as function of Reynolds:

Re Error Relat

1 9.0978E-13

5 2.3235E-13

10 5.32550E-14

15 2.4712E-14

20 3.8433E-14

25 1.1738E-14

30 4.5003E-14

35 2.8743E-14

40 1.286E-14

45 1.2921E-14

50 2.352E-14

60 2.931E-14

70 0.061668

80 0.36923

90 NaN

100 NaN

Note that for Re=70 the Mach is above 0.58. However for Re=60 Mach is 0.53 and the numerical solution is equal to the analytical one to machine accuracy.

Using LY=30, LX=15 but now tau=0.9 the Mach vs Error evolution is:

0.5 0.020547

0.4 0.017771

0.3 0.0099619

0.2 0.0050707

0.1 0.0013894

0.09 0.0011588

0.08 0.00091038

0.07 0.00068389

0.06 0.0005326

0.05 0.00036605

This results in a 1.79 accuracy. If I stick to Mach numbers below 0.1 the slope is 1.94 which is much closer to 2. In my opinion this result is acceptable.

First I thought that the problem was related to the forcing term, so I removed it and checked the results.

So, setting the body force to zero and using just the pressure boundary conditions to drive the flow (for Re=10, LY=30 and LX=15) I will have a similar behavior, i.e. the minimum error is not for the minimum tau (where I have the minimum Mach number) but now for tau=0.9!:

tau Error Relat

0.51 0.0019097

0.55 0.0018667

0.60000 0.0017678

0.70000 0.0014308

0.80000 0.00091209

0.90000 0.0002201

1.00000 0.00063735

1.10000 0.0016526

1.50000 0.0071377

1.70000 0.01062

2.00000 0.016571

2.50000 0.027825

3.00000 0.039844

Doing exactly the same thing, i.e. using only pressure gradients to move the flow, but setting constant tau=1, Re=10 I have:

Ma Error Relat

0.5 0.019608

0.4 0.013699

0.3 0.0077519

0.2 0.0029499

0.1 0.00079936

0.09 0.00063735

0.08 0.00052002

0.07 0.00038565

0.06 0.00028337

0.05 0.00019996

This gives a slop of 2.0155 which agrees very well with the expected 2nd order accuracy.

So my code produces logical results for tau=1 if no body force is used or if a body force is used but tau is different from 1.

If you want I can send you my codes (they are simple and short MATLAB codes)…

Again thank you for your help.

Regards

Goncalo Silva

P.S. I don’t know if this is relevant for the understanding of my problem but, to have a poiseuille flow without a developing region, i.e. to have a linear pressure drop from the inlet to the outlet I defined my body force as:

g=0.5*12.*Re*kvisc.^2./((LY-2).^3) (0.5 factor is just to account that this body force will only contribute for the flow movement in 50%)

and my outlet pressure drop value as:

rho_outlet=(rho_inlet-3.*12.**Re*kvisc.^2.(rho0).*(LX-1)./((LY-2).^3)) + 3.*g(1).*rho0.*(LX-1) (in my case rho_inlet was simply set to 1).

So my body force/pressure difference between inlet/outlet changes whenever I change Re, tau, Ma, etc…