Difference Finite Scheme of Diffusive Equation derived from Lattice Boltzmann

Due to the similarities in the grid structure between the Lattice Boltzmann method and finite difference schemes, several studies have been conducted over the years to investigate how these methods can be directly related and interpreted [3, 4, 5, 6, 7, 8]. In this section, we develop a finite-difference interpretation for different Lattice Boltzmann formulations.

Interpretation for Diffusive Equation

In the remainder of this section, we develop a FDM interpretation of the diffusive LBM, following the approach proposed by Suga [5].

Given the LBM formulation for the diffusive equation:

(55)\[ f_i( x_{\alpha} + e_{i,\alpha} \delta t, t+\delta t) = f_i(x_{\alpha}, t) -\left( \frac{ f_{i} - f_{i}^{eq} }{ \tau } \right) \qquad \rightarrow \qquad f_i( x_{\alpha} + e_{i,\alpha} \delta t, t+\delta t) = f_i(x_{\alpha}, t) - \omega\left( f_{i} - f_{i}^{eq} \right), \]

where the equilibrium distribution function is defined by

(56)\[ f^{eq}_{i} = w_{i}\phi(x_{\alpha}, t), \]

and the equilibrium moments are given by

(57)\[ \displaystyle\sum_{i=0} f_{i}^{eq}=\phi, \quad \quad \displaystyle\sum_{i}f_{i}^{eq} e_{i,\alpha} =0 \quad \quad \textrm{and} \quad \quad \displaystyle\sum_{i}f_{i}^{eq} e_{i,\alpha} e_{i,\beta} =c_{s}^{2} \phi \delta_{\alpha\beta}. \]

Through Chapman–Enskog analysis, it is demonstrated that the first-order non-equilibrium moment describes the pressure gradient and, consequently, the average velocity:

\[ m^{neq}_{\alpha} = \displaystyle\sum_{i} e_{i,\alpha}\left( f_i-f_i^{eq}\right) = -(c_{s}^{2} \tau \delta_{t})\partial_{\alpha}\phi \quad \quad \textrm{and} \quad\quad \tau = \displaystyle\frac{\nu}{c_{s}^{2} \delta_{t}} +\frac{1}{2}. \]

FDM Interpretation

For simplification we redefine the notation for $f_i( x_{\alpha} + e_{i,\alpha} \delta t, t+\delta t)=f_{i,x+c_{i}}^{t+1}$. Reinterpreting the Eq. (55) in pre-streaming moment:

(58)\[ f_{i,x}^{t+1} = f_{i,x-c_{i}}^{t} - \omega \left( f_{i,x-c_{i}}^{t} - f_{i,x-c_{i}}^{eq,t} \right). \]

Summing the Eq. (58) over all lattice direction and replacing $f_{i}^{eq}$ by Eq. (56), we have:

\[ \displaystyle\sum_{i}\Big( f_{i,x}^{t+1} \Big) = \displaystyle\sum_{i}\Big( f_{i,x-c_{i}}^{t} - \omega \left( f_{i,x-c_{i}}^{t} - f_{i,x-c_{i}}^{eq,t} \right) \Big) \qquad \rightarrow \qquad \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega) f_{i,x-c_{i}}^{t} + \omega w_{i}\phi_{x-c_{i}}^{t} \Big), \]

using the Eq. (57) to replace $f_{i,x-c_{i}}^{t}$ in the above equation:

\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega) f_{i,x-c_{i}}^{t} + \omega w_{i}\phi_{x-c_{i}}^{t} \Big) \qquad \rightarrow \qquad \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t} - (1-\omega) \displaystyle\sum_{k\neq i} f_{k,x-c_{i}}^{t} \Big). \]

Spliting the terms in the above equation and using Eq. (58) for $f_{i,x-c_{i}}^{t}$:

\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) - (1-\omega)\displaystyle\sum_{i}\left( \displaystyle\sum_{k\neq i} f_{k,x-c_{i}}^{t} \right), \]

\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) - (1-\omega)\displaystyle\sum_{i}\Bigg( \displaystyle\sum_{k\neq i} \left( f_{k,x-c_{i}-c_{k}}^{t-1} - \omega \left( f_{k,x-c_{i}-c_{k}}^{t-1} - f_{k,x-c_{i}-c_{k}}^{eq,t-1} \right) \right) \Bigg), \]

doing the same initial procedure to the terms $f_{k,x-c_{i}}^{t-1}$:

\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) -(1-\omega)\displaystyle\sum_{i}\Bigg( \displaystyle\sum_{k\neq i} \left( f_{k,x-c_{i}-c_{k}}^{t-1} - \omega \left( f_{k,x-c_{i}-c_{k}}^{t-1} - f_{k,x-c_{i}-c_{k}}^{eq,t-1} \right) \right) \Bigg) \]

\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) - (1-\omega) \displaystyle\sum_{i}\Bigg( \displaystyle\sum_{k\neq i} \left( \omega w_{k} \phi_{x-c_{i}-c_{k}}^{t-1} + (1-\omega) f_{k,x-c_{i}-c_{j}}^{t-1} \right) \Bigg), \]

notice that $\displaystyle\sum_{i}\Bigg( \displaystyle\sum_{k\neq i} \left( \omega w_{k} \phi_{x-c_{i}-c_{k}}^{t-1} \right) \Bigg)=\displaystyle\sum_{i}\Bigg( \displaystyle\sum_{k\neq i} \left(w_{k}\right) \omega \phi_{x+c_{i}}^{t-1} \Bigg)$ and using Eq. (57) we have $\displaystyle\sum_{i}\Bigg( \displaystyle\sum_{k\neq i} \left( f_{k,x-c_{i}-c_{j}}^{t-1} \right) \Bigg)= \displaystyle\sum_{i}\left(\phi_{x+c_{i}}^{t-1} - f_{i,x+c_{i}}^{t-1}\right) $. So rewriting the above equation base in the observations:

(59)\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) - (1-\omega)\displaystyle\sum_{i}\Bigg( \left(1-\omega + \omega \displaystyle\sum_{k\neq i} \left(w_{k}\right) \right) \phi_{x+c_{i}}^{t-1} \Bigg) + (1-\omega)^{2}\displaystyle\sum_{i}\Bigg( f_{i,x+c_{i}}^{t-1} \Bigg). \]

Interpreting $f_{i,x+c_{i}}^{t-1}$ with Eq. (58):

\[ f_{i,x+c_{i}}^{t-1} = f_{i,x}^{t-2} - \omega \left( f_{i,x}^{t-2} - f_{i,x}^{eq,t-2} \right) , \]

and replaing it in Eq. (59), we obtained a four-level explicity scheme:

\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) - (1-\omega)\displaystyle\sum_{i}\Bigg( \left(1-\omega + \omega \displaystyle\sum_{k\neq i} \left(w_{k}\right) \right) \phi_{x+c_{i}}^{t-1} \Bigg) + (1-\omega)^{2}\displaystyle\sum_{i}\Bigg(f_{i,x}^{t-2} - \omega \left( f_{i,x}^{t-2} - f_{i,x}^{eq,t-2} \right) \Bigg), \]

(60)\[ \phi_{x}^{t+1} = \displaystyle\sum_{i}\Big( (1-\omega + \omega w_{i}) \phi_{x-c_{i}}^{t}\Big) - (1-\omega)\displaystyle\sum_{i}\Bigg( \left(1-\omega + \omega \displaystyle\sum_{k\neq i} \left(w_{k}\right) \right) \phi_{x+c_{i}}^{t-1} \Bigg) + (1-\omega)^{2} \phi_{x}^{t-2} \]

Expanding the sums:

\[\begin{split} \begin{array}{c} \phi_{x}^{t+1} = (1-\omega + \omega w_{0}) \phi_{x}^{t} + (1-\omega + \omega w_{1}) \phi_{x-1}^{t} + (1-\omega + \omega w_{2}) \phi_{x+1}^{t} - (1-\omega)^{2} \Big( \phi_{x}^{t-1} + \phi_{x-1}^{t-1} + \phi_{x+1}^{t-1} \Big) \\ -(1-\omega)\omega \Big( (w_{1}+w_{2})\phi_{x}^{t-1} + (w_{0}+w_{1})\phi_{x-1}^{t-1} + (w_{0}+w_{2})\phi_{x+1}^{t-1} \Big) + (1-\omega)^{2} \phi_{x}^{t-2} \end{array} \end{split}\]

1st Benchmark: Decay of Sinoidal Wave

Considering a 1D convection–diffusion equation with constant coefficients $$ \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2} \qquad x\in\mathbb{R},\ t\ge 0, $$

with sinusoidal initial condition

\[ u(x,0)=A\sin(kx+\phi), \]

where $A$ is the initial amplitude, $k$ the wavenumber, $\phi$ a phase shift.

Then the analytical solution is $$ \boxed{u(x,t)=Ae^{-\nu k^{2}t}\sin\big(k(x-ct)+\phi\big) } $$

So:

Convection shifts the wave: $x \mapsto x-ct$ (wave travels at speed $c$).
Diffusion damps the amplitude exponentially: $A(t)=Ae^{-\nu k^{2}t}$.

If you use a cosine initial condition $u(x,0)=A\cos(kx+\phi)$, the solution is the same form: $$ u(x,t)=Ae^{-\nu k^{2}t}\cos\big(k(x-ct)+\phi\big). $$

LBM D1Q3 Code - FDM Code (Only Diffusive Case $c=0$)

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
A   = 1.0                         # Amplitude
phi = 0.0                         # Wave initial shift
k_w   = 1.0                       # Wave number
c   = 0.0                         # wave speed
nu  = 0.0066666666666666666*np.pi # Diffusive coefficient
m_waves = 1                       # Wave Period
L = m_waves * (2.0*np.pi / k_w)   # Domain Lentgh
t_end = 48.0                       # Final time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int16")  
cs=1.0/np.sqrt(3.0);
#------------------------------------- Parameters (numerical units) -------------------------------------------
Nx = 25 # Numerical Length
x  = np.linspace(0.0, L, Nx, endpoint=False,dtype="float64") # Numerical Domain
dx = L / Nx # Grid size
r=2**(1) # dx/dt relation
dt = dx/r
ce = c/r
nue = nu * dt / (dx*dx)
tau = 0.5 + nue / cs**2
omg=1.0/tau
nt = int(np.ceil(t_end / dt))
print(f"dx={dx:.4e}\t dt={dt:.4e}\t nt={nt:d}") # Print values for check
print(f"c_lbm={ce:.3f}, tau={tau:.4f}\t omega={omg:.4f}")
print(f"nt={nt:.3f}, nue={nue:.4f}")
#--------------------------------- Initialization - LBM - Numerical Arrays -----------------------------------------
u=np.zeros((Nx),dtype="float64")
u = A * np.sin(k_w*(x+dx/2) + phi)
f=np.zeros((3,Nx),dtype="float64")
fp=np.zeros((3,Nx),dtype="float64")
for k in range(0,3):
    f[k,:]=w[k]*(u[:]+cx[k]*ce*u[:]/cs**2)
    fp[k,:]=w[k]*(u[:]+cx[k]*ce*u[:]/cs**2)
#--------------------------------- Initialization - FDM - Numerical Arrays -----------------------------------------
ufm2=np.zeros((Nx),dtype="float64")    # u field at the step t-2
ufm=np.zeros((Nx),dtype="float64")    # u field at the step t-1
uf=np.zeros((Nx),dtype="float64")     # u field at the step t
ufp=np.zeros((Nx),dtype="float64")    # u field at the step t+1
ufm2 = A * np.sin(k_w*(x+dx/2) + phi)
ufm = A * np.sin(k_w*(x+dx/2) + phi)
uf = A * np.sin(k_w*(x+dx/2) + phi)
ufp = A * np.sin(k_w*(x+dx/2) + phi)
#--------------------------------- Initialization - Save Data Arrays -----------------------------------------------
snaps=5 # number of states saved over time including 0 and t_end
u_snaps = np.empty(snaps, dtype=object) # array field used to same data over time
uf_snaps = np.empty(snaps, dtype=object) # array field used to same data over time
snaps_id = np.linspace(0, nt, snaps, dtype="int16") # timesteps to take the field
snap_index = {sid: i for i, sid in enumerate(snaps_id)}
for i in range(snaps):
    u_snaps[i] = np.zeros((Nx), dtype="float64")
    uf_snaps[i] = np.zeros((Nx), dtype="float64")
u_snaps[0][:]=u[:]
uf_snaps[0][:]=uf[:]
#----------------------------------------- Maind Loop --------------------------------------------------------
for t in range(nt):
    #========================================= LBM - Solution =============================================
    #--------------------Collision----------------
    for k in range(0,3):
        fp[k,:]= f[k,:] - (f[k,:] - w[k]*(u[:]+cx[k]*(ce*u[:])/cs**2))/tau
    #-----------------streaming-------------------
    for k in range(0,3):
        f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
    #----------------------Macro------------------
    u[:]=f[0,:]+f[1,:]+f[2,:]
    #---------------------save-snaps--------------
    if t+1 in snaps_id:
        i = snap_index[t+1]
        u_snaps[i][:]=u[:]
    #========================================= FDM - Solution =============================================
    if (t<2) :
        ufm2[:]=ufm[:]     # Swap for update t-1 to t-2
        ufm[:]=uf[:]       # Swap for update t to t-1
        uf[:]=u[:]         # Swap for update t to t-1
    else :
        ufp[:]=( ( (1-omg+omg*w[1])*np.roll(uf[:],1,axis=0) + (1-omg+omg*w[0])*uf[:] + (1-omg+omg*w[2])*np.roll(uf[:],-1,axis=0) )
                 - (1-omg)**(2)*( np.roll(ufm[:], 1, axis=0) + ufm[:] + np.roll(ufm[:], -1, axis=0) )
                 - (1-omg)*omg*( (w[0]+w[1])*np.roll(ufm[:], 1, axis=0) + (w[2]+w[1])*ufm[:] + (w[2]+w[0])*np.roll(ufm[:], -1, axis=0) )
                 + (1.0-omg)**(2)*ufm2[:]  )
        ufm2[:]=ufm[:]     # Swap for update t-1 to t-2
        ufm[:]=uf[:]       # Swap for update t to t-1
        uf[:]=ufp[:]       # Swap for update t+1 to t
    if t+1 in snaps_id:
        i = snap_index[t+1]
        uf_snaps[i][:]=ufp[:]

../_images/cfebaa2c874f367e958242232ecc0de1e981727e6d2fa6ea6a30d0a1dfe32f64.png

2nd Benchmark: Temperature ramp

Analyzing the accuracy of diffusive equation solver, we apply it to a one-dimensional analytical case. The geometry is described by a domain of length $L$ initialized with a constant temperatura $T(x,t=0)=0$. The boundary conditions are given $T(x=0,t)=1$ and $T(x=L,t)=0$.

Steady-State Analytical Solution

The analytical solution of this problem is given by:

\[ T(x,t\rightarrow \infty)= -\frac{x}{L} + 1. \]

Transient Analytical Solution

For the 1D heat equation $\partial_t T=\nu\,\partial_{xx}T$ on $0<x<L$ with

\[ T(0,t)=1,\qquad T(L,t)=0, \]

the steady state is $T_{\text{ss}}(x)=1-\tfrac{x}{L}$. With $u=T-T_{\text{ss}}$ (so $u(0,t)=u(L,t)=0$), the solution is

\[ T(x,t)=1-\frac{x}{L}+\sum_{n=1}^{\infty} B_n\,e^{-\nu (n\pi/L)^2 t}\,\sin\!\Big(\frac{n\pi x}{L}\Big), \quad B_n=\frac{2}{L}\int_0^L\!\Big[T_0(\xi)-\Big(1-\frac{\xi}{L}\Big)\Big]\sin\!\Big(\frac{n\pi \xi}{L}\Big)\,d\xi. \]

Special case $T_0(x)=0$: $B_n=-\dfrac{2}{n\pi}$, so

\[ T(x,t)=1-\frac{x}{L}-\frac{2}{\pi}\sum_{n=1}^{\infty}\frac{1}{n}\,e^{-\nu (n\pi/L)^2 t}\,\sin\!\Big(\frac{n\pi x}{L}\Big). \]

Here’s a ready-to-run Python snippet (with a single plot of $T(x,t)$ vs $x$ for multiple times):

LBM D1Q3 Code - FDM Code

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
L = 1.0             # Length of the domain
nu = 0.1/2.0        # Diffusion coefficient
t_end = 8.0             # Total simulation time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
cs=1.0/np.sqrt(3.0);
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int8")  
#------------------------------------- Parameters (numerical units) -------------------------------------------
Nx=21               # Domain Length
dx = L / (Nx-1)     # Spatial step size
c=1.0*2**(2)        # c=dx/dt
dt=dx/c             # Time step value
nt = int(t_end / dt)    # Time step number
nue=nu/(c*dx)       # Diffusion coefficient
tau=nue/cs**2+0.5   # Relaxation Time
omg=1.0/tau
print(f"dx={dx:.4e}\t dt={dt:.4e}\t nt={nt:d}") # Print values for check
print(f"nue={nue:.4f}\t tau={tau:.4f}\t omega={omg:.4f}")
#--------------------------------- Initialization - LBM - Numerical Arrays -----------------------------------------
T = np.zeros((Nx),dtype="float64")
T[0]=1.0 # Initital Temperature Field
f = np.zeros((3,Nx),dtype="float64")
fp = np.zeros((3,Nx),dtype="float64")
for k in range(0,3):
    f[k,:]=w[k]*T[:]
    fp[k,:]=w[k]*T[:]
#--------------------------------- Initialization - FDM - Numerical Arrays -----------------------------------------
Tfm2 = np.zeros((Nx),dtype="float64")    # Temperature field at the step t-2
Tfm = np.zeros((Nx),dtype="float64")    # Temperature field at the step t-1
Tf = np.zeros((Nx),dtype="float64")     # Temperature field at the step t
Tfp = np.zeros((Nx),dtype="float64")    # Temperature field at the step t+1
Tfm2[0]=1.0                              # Initital Temperature Field
Tfm[0]=1.0                              # Initital Temperature Field
Tf[0]=1.0                               # Initital Temperature Field
Tfp[0]=1.0                              # Initital Temperature Field
#--------------------------------- Initialization - Save Data Arrays -----------------------------------------
snaps=5 # number of states saved over time including 0 and t_end
T_snaps = np.empty(snaps, dtype=object) # array field used to same data over time
Tf_snaps = np.empty(snaps, dtype=object) # array field used to same data over time
snaps_id = np.linspace(0, nt, snaps, dtype="int16") # timesteps to take the field
snap_index = {sid: i for i, sid in enumerate(snaps_id)}
for i in range(snaps):
    T_snaps[i] = np.zeros((Nx), dtype="float64")
    Tf_snaps[i] = np.zeros((Nx), dtype="float64")
T_snaps[0][:]=T[:]
Tf_snaps[0][:]=Tf[:]    
#----------------------------------------- Maind Loop --------------------------------------------------------
for t in range(nt):
    #========================================= LBM - Solution =============================================
    #--------------------Collision----------------
    for k in range(0,3):
        fp[k,:]=f[k,:] - (f[k,:] - w[k]*T[:])/tau
    #-----------------streaming-------------------
    for k in range(0,3):
        f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
    #-----------------Boundaries-----------------------
    f[1,0]= 1.0 - f[0,0]-f[2,0]
    f[2,Nx-1]= 0.0 - f[0,Nx-1]-f[1,Nx-1]
    #----------------------Macro------------------
    T[:]=f[0,:]+f[1,:]+f[2,:]
    #---------------------save-snaps--------------
    if t+1 in snaps_id:
        i = snap_index[t+1]
        T_snaps[i][:]=T[:]
    #========================================= FDM - Solution =============================================
    if (t<2) :
        Tfm2[:]=Tfm[:]     # Swap for update t-1 to t-2
        Tfm[:]=Tf[:]     # Swap for update t to t-1
        Tf[:]=T[:]     # Swap for update t to t-1
    else :
        Tfp[:]=( ( (1-omg+omg*w[1])*np.roll(Tf[:],1,axis=0) + (1-omg+omg*w[0])*Tf[:] + (1-omg+omg*w[2])*np.roll(Tf[:],-1,axis=0) )
                 - (1-omg)**(2)*( np.roll(Tfm[:], 1, axis=0) + Tfm[:] + np.roll(Tfm[:], -1, axis=0) )
                 - (1-omg)*omg*( (w[0]+w[1])*np.roll(Tfm[:], 1, axis=0) + (w[2]+w[1])*Tfm[:] + (w[2]+w[0])*np.roll(Tfm[:], -1, axis=0) )
                 + (1.0-omg)**(2)*Tfm2[:]  )
        Tfp[0]=1.0
        Tfp[Nx-1]=0.0
        Tfm2[:]=Tfm[:]     # Swap for update t-1 to t-2
        Tfm[:]=Tf[:]     # Swap for update t to t-1
        Tf[:]=Tfp[:]     # Swap for update t+1 to t
    #---------------------save-snaps------------
    if t+1 in snaps_id:
        i = snap_index[t+1]
        Tf_snaps[i][:]=Tfp[:]

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

# ---- Modes and series -------------------------------------------------------
def modes_dirichlet_dirichlet(L, N):
    n = np.arange(1, N+1)        # n = 1,2,...,N
    k = n * np.pi / L            # eigenvalues mu_n
    return n, k

def project_coeffs_dirichlet(T0_func, L, N, grid_pts=4000):
    """
    B_n = (2/L) ∫( T0(x) - T_ss(x) ) sin(nπx/L) dx, with T_ss = 1 - x/L.
    """
    x = np.linspace(0.0, L, grid_pts, endpoint=True)
    Tss = 1.0 - x / L
    u0 = T0_func(x) - Tss
    n, k = modes_dirichlet_dirichlet(L, N)          # (N,), (N,)
    S = np.sin(k[:, None] * x[None, :])             # (N, M)
    B = (2.0 / L) * np.trapz(u0[None, :] * S, x, axis=1)
    return B, n, k

def T_series_dirichlet(x, t, L, nu, B, k):
    """
    T(x,t) = T_ss(x) + sum_n B_n e^{-nu k_n^2 t} sin(k_n x), T_ss = 1 - x/L.
    """
    x = np.asarray(x)
    Tss = 1.0 - x / L
    coeff = B * np.exp(-nu * (k**2) * t)            # (N,)
    add = np.sum(coeff[:, None] * np.sin(k[:, None] * x[None, :]), axis=0)
    return Tss + add

# ---- Parameters -------------------------------------------------------------
L   = 1.0
nu  = 0.1/2.0 
N   = 200
x   = np.linspace(0, L, 600)

# ---- Case A: initially cold rod T0(x)=0 (closed-form coefficients) ----------
n, k = modes_dirichlet_dirichlet(L, N)
B_cold = -2.0 / (n * np.pi)   # exact coefficients for T0 ≡ 0

times = [0.0, 2.0, 4.0, 6.0, 8.0]

#--------------------------------------- Plot Results -------------------------------------------------------
xl = np.linspace(0, L, Nx)
times = np.linspace(0, t_end, snaps, dtype="float64") # timesteps to take the field
plt.figure(figsize=(7, 5))
colors = plt.cm.viridis(np.linspace(0, 1, snaps))
for i in range(snaps):
    T_xt = T_series_dirichlet(x, times[i], L, nu, B_cold, k)
    plt.plot(x, T_xt, "--", color=colors[i], label=f"Exact $t = {t:.2f}$")
    plt.plot(xl, T_snaps[i], "s:", color=colors[i], lw=1, label=f"LBM  $t={(snaps_id[i]/nt)*t_end :.2f}$",fillstyle='none')
    plt.plot(xl, Tf_snaps[i], "o:", color=colors[i], lw=1, label=f"FDM  $t={(snaps_id[i]/nt)*t_end :.2f}$",fillstyle='none')
plt.xlabel(r"$x$",fontsize=18)
plt.ylabel(r"$\phi(x,t)$",fontsize=18)
plt.legend(ncol=3,fontsize=8)
plt.xlim(-0.05,1)
plt.ylim(-0.2,1.01)
plt.grid(True)
plt.tight_layout()
plt.show()

../_images/3c0acc32afd9c73d31376812bbe500b2c90785d7d5b4ac070c34314ac8d3f3a0.png

Fourth-order Diffusive Equation

Expanding the sums:

\[\begin{split} \begin{array}{c} \phi_{x}^{t+1} = (1-\omega + \omega w_{0}) \phi_{x}^{t} + (1-\omega + \omega w_{1}) \phi_{x-1}^{t} + (1-\omega + \omega w_{2}) \phi_{x+1}^{t} - (1-\omega)^{2} \Big( \phi_{x}^{t-1} + \phi_{x-1}^{t-1} + \phi_{x+1}^{t-1} \Big) \\ -(1-\omega)\omega \Big( (w_{1}+w_{2})\phi_{x}^{t-1} + (w_{0}+w_{1})\phi_{x-1}^{t-1} + (w_{0}+w_{2})\phi_{x+1}^{t-1} \Big) + (1-\omega)^{2} \phi_{x}^{t-2} \end{array} \end{split}\]

Considering $$ \alpha_{1}=\Omega + a\omega, \quad \alpha_{2}=\Omega + (1-2a)\omega, \quad \beta_{1}=-\Omega(\Omega + (1-a)\omega), \quad \textrm{and} \quad \beta_{2}=-\Omega(\Omega + 2a\omega), $$

where $\Omega=(1-\omega)$ and $a=w_{1}=w_{2}$. Replacing in the Eq. (60):

\[ \begin{array}{c} \phi_{x}^{t+1} = \alpha_{1} \phi_{x-1}^{t} + \alpha_{2} \phi_{x}^{t} + \alpha_{1} \phi_{x+1}^{t} + \beta_{1}\phi_{x-1}^{t-1} + \beta_{2}\phi_{x}^{t-1} + \beta_{1}\phi_{x+1}^{t-1} + (1-\omega)^{2} \phi_{x}^{t-2} \end{array} \]

To recover the partial differetial equation described by the above equation, we expand each $\phi$ term in Taylor series:

\[ \phi_{x}^{t\pm 1} = \phi \pm \Delta_{t}\frac{\partial\phi}{\partial t} + \frac{\Delta_{t}^{2}}{2}\frac{\partial^{2}\phi}{\partial t^{2}} \pm \frac{\Delta_{t}^{3}}{6}\frac{\partial^{3}\phi}{\partial t^{3}} + \frac{\Delta_{t}^{4}}{24}\frac{\partial^{4}\phi}{\partial t^{4}} \pm \mathcal{O}(\Delta t)^{4}, \]

\[ \phi_{x\pm 1}^{t} = \phi \pm \Delta_{x}\frac{\partial\phi}{\partial x} + \frac{\Delta_{x}^{2}}{2}\frac{\partial^{2}\phi}{\partial x^{2}} \pm \frac{\Delta_{x}^{3}}{6}\frac{\partial^{3}\phi}{\partial x^{3}} + \frac{\Delta_{x}^{4}}{24}\frac{\partial^{4}\phi}{\partial x^{4}} \pm \mathcal{O}(\Delta x)^{4}, \]

\[ \phi_{x}^{t - 2} = \phi - 2\Delta_{t}\frac{\partial\phi}{\partial t} + \frac{4\Delta_{t}^{2}}{2}\frac{\partial^{2}\phi}{\partial t^{2}} - \frac{8\Delta_{t}^{3}}{6}\frac{\partial^{3}\phi}{\partial t^{3}} + \frac{16\Delta_{t}^{4}}{24}\frac{\partial^{4}\phi}{\partial t^{4}} \pm \mathcal{O}(\Delta t)^{4}, \]

\[ \phi_{x\pm 1}^{t-1} = \phi \pm \Delta_{x}\frac{\partial\phi}{\partial x} - \Delta_{t}\frac{\partial\phi}{\partial t} + \frac{\Delta_{x}^{2}}{2}\frac{\partial^{2}\phi}{\partial x^{2}} \mp \Delta_{t}\Delta_{x} \frac{\partial^{2}\phi}{\partial x\partial t} + \frac{\Delta_{t}^{2}}{2}\frac{\partial^{2}\phi}{\partial t^{2}} \pm \frac{\Delta_{x}^{3}}{6}\frac{\partial^{3}\phi}{\partial x^{3}} \mp \frac{\Delta_{x}^{2}\Delta_{t}}{2}\frac{\partial^{3}\phi}{\partial x^{2}\partial t} \pm \frac{\Delta_{x}\Delta_{t}^{2}}{2}\frac{\partial^{3}\phi}{\partial x\partial t^{2}} - \frac{\Delta_{t}^{3}}{6}\frac{\partial^{3}\phi}{\partial t^{3}} + \mathcal{O}(\Delta x \Delta t)^{3}. \]

Replacing the expanded terms and simplifying the term (this link provide a sympy code that reach the above equation Sympy code for FDM expansion of LBM diffusive equation), we have:

\[ 0=\Delta_{t} \omega^{2} \frac{\partial\phi}{\partial t} + \Delta_{x}^{2} a \omega \left(\omega - 2\right) \frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\Delta_{t}^{2} \omega \left(4 - 3 \omega\right) }{2} \frac{\partial^{2}\phi}{\partial t^{2}} - \Delta_{t} \Delta_{x}^{2} \left(a \omega^{2} - a \omega - \omega + 1\right) \frac{\partial^{3}\phi}{\partial x^{2}\partial t} + \frac{\Delta_{x}^{4} a \omega \left(\omega - 2\right) }{12}\frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \Delta_{t} \omega^{2} \frac{\partial\phi}{\partial t} = - \Delta_{x}^{2} a \omega \left(\omega - 2\right) \frac{\partial^{2}\phi}{\partial x^{2}} - \frac{\Delta_{t}^{2} \omega \left(4 - 3 \omega\right) }{2} \frac{\partial^{2}\phi}{\partial t^{2}} + \Delta_{t} \Delta_{x}^{2} \left(a \omega^{2} - a \omega - \omega + 1\right) \frac{\partial^{3}\phi}{\partial x^{2}\partial t} - \frac{\Delta_{x}^{4} a \omega \left(\omega - 2\right) }{12}\frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial\phi}{\partial t} = \frac{a \left(2 - \omega \right)}{\omega} \frac{\Delta_{x}^{2}}{\Delta_{t}} \frac{\partial^{2}\phi}{\partial x^{2}} + \frac{ \left(3 \omega - 4\right) \Delta_{t}}{2\omega} \frac{\partial^{2}\phi}{\partial t^{2}} + \frac{\left(a \omega^{2} - a \omega - \omega + 1\right)\Delta_{x}^{2}}{\omega^{2}} \frac{\partial^{3}\phi}{\partial x^{2}\partial t} + \frac{ a \left(2 - \omega\right) }{12\omega} \frac{\Delta_{x}^{4}}{\Delta_{t}} \frac{\partial^{4}\phi}{\partial x^{4}} \]

where $\nu=\frac{a \left(2 - \omega \right)}{\omega} \frac{\Delta_{x}^{2}}{\Delta_{t}}$ is the same term recovered in the Chapman-Enskog Analysis. To enforce fourth-order $\Delta_{x}$ convergence rate, employ the relations

\[ \frac{\partial\phi}{\partial t} = \nu \frac{\partial^{2}\phi}{\partial x^{2}}, \qquad \frac{\partial^{2}\phi}{\partial t^{2}} = \nu^{2} \frac{\partial^{4}\phi}{\partial x^{4}}, \qquad \frac{\partial^{3}\phi}{\partial x^{2}\partial t} = \nu \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial\phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{a \left(2 - \omega \right)}{\omega} \frac{ \left(3 \omega - 4\right) \nu\Delta_{x}^{2}}{2\omega} \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\left(a \omega^{2} - a \omega - \omega + 1\right)\nu\Delta_{x}^{2}}{\omega^{2}} \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\nu\Delta_{x}^{2}}{12} \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial\phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{ a\left(6 \omega - 8 - 3\omega^{2} + 4 \omega \right) \nu\Delta_{x}^{2}}{2\omega^{2}} \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\left(a \omega^{2} - a \omega - \omega + 1\right)\nu\Delta_{x}^{2}}{\omega^{2}} \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\omega^{2}\nu\Delta_{x}^{2}}{12\omega^{2}} \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial\phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\nu\Delta_{x}^{2}}{2\omega^{2}}\left( a\left(10 \omega - 8 - 3\omega^{2} \right) + 2\left(a \omega^{2} - a \omega - \omega + 1\right) + \frac{\omega^{2}}{6} \right) \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial\phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\nu\Delta_{x}^{2}}{2\omega^{2}}\left( a\left(8 \omega - 8 - \omega^{2} \right) + 2\left(1 - \omega \right) + \frac{\omega^{2}}{6} \right) \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial\phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} - \frac{\nu\Delta_{x}^{2}}{2\omega^{2}}\left( a\left(\omega^{2} - 8 \omega + 8 \right) - \frac{\omega^{2} - 12\omega + 12 }{6} \right) \frac{\partial^{4}\phi}{\partial x^{4}} + \mathcal{O}(\Delta x)^{4} \]

../_images/46430f4b0a93c85e1b1516f1720dedad8668b0d547142476fc289ba710f3f807.png

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
A   = 1.0                         # Amplitude
phi = 0.0                         # Wave initial shift
k_w   = 1.0                       # Wave number
c   = 0.0                         # wave speed
nu  = 0.066666666666666666*np.pi # Diffusive coefficient
m_waves = 1                       # Wave Period
L = m_waves * (2.0*np.pi / k_w)   # Domain Lentgh
t_end = 12.0                       # Final time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int16")  
cs=1.0/np.sqrt(3.0);
#--------------------------------- Initialization - Save Data Arrays -----------------------------------------------
cases=4
Nx0 = np.array([10, 20, 40, 80],dtype="int64")
r0 = 2**(0)*np.array([2.0**(1), 2.0**(2), 2.0**(3), 2.0**(4)],dtype="float64")
u_cases = np.empty(cases, dtype=object) # array field used to same data over time
for i in range(cases):
    u_cases[i] = np.zeros((Nx0[i]), dtype="float64")
for case in range(0,cases):
    #------------------------------------- Parameters (numerical units) -------------------------------------------
    Nx=Nx0[case]  # Numerical Length
    x  = np.linspace(0.0, L, Nx, endpoint=False,dtype="float64") # Numerical Domain
    dx = L / Nx # Grid size
    r = 1*r0[case] # dx/dt relation
    dt = dx/r
    ce = c/r
    nue = nu * dt / (dx*dx)
    tau = 0.5 + nue / cs**2
    omg=1.0/tau
    nt = int(np.ceil(t_end / dt))
    print(f"dx={dx:.4e}\t dt={dt:.4e}\t nt={nt:d}\t tn={nt*dt:.2f}") # Print values for check
    print(f"c_lbm={ce:.3f}, tau={tau:.4f}\t omega={omg:.2f}")
    print(f"nt={nt:.3f}, nue={nue:.4f}")
    #--------------------------------- Initialization - LBM - Numerical Arrays -----------------------------------------
    u=np.zeros((Nx),dtype="float64")
    u = A * np.sin(k_w*(x+dx/2) + phi)
    f=np.zeros((3,Nx),dtype="float64")
    fp=np.zeros((3,Nx),dtype="float64")
    for k in range(0,3):
        f[k,:]=w[k]*(u[:]+cx[k]*ce*u[:]/cs**2)
        fp[k,:]=w[k]*(u[:]+cx[k]*ce*u[:]/cs**2)
    #----------------------------------------- Maind Loop --------------------------------------------------------
    for t in range(nt):
        #========================================= LBM - Solution =============================================
        #--------------------Collision----------------
        for k in range(0,3):
            fp[k,:]= f[k,:] - (f[k,:] - w[k]*(u[:]+cx[k]*(ce*u[:])/cs**2))/tau
        #-----------------streaming-------------------
        for k in range(0,3):
            f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
        #----------------------Macro------------------
        u[:]=f[0,:]+f[1,:]+f[2,:]
    #---------------------save-field-cases--------------
    u_cases[case][:]=u[:]

../_images/9a9d0fedd0c5d98771916c2dcbbb63c063c33e0b831df1afa5aedf24fbfc5b34.png

../_images/a1b62e754403dd35e7c6e1a60c65923a6d25558da918decbfbb0dacdf8f3b64d.png

Extension for Higher-Dimensions and Anisotropy

Suga [9] also extend the approach for 2 and 3D cases considering anisotropic scenarios. Notice that different second order errors will appear with other dimensions:

\[\begin{split} \begin{array}{ll} 1D-Isotropic: \qquad & \displaystyle\frac{\partial\phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} - \textcolor{red} {\frac{\nu\Delta_{x}^{2}}{2\omega^{2}}\left( a\left(\omega^{2} - 8 \omega + 8 \right) - \frac{\omega^{2} - 12\omega + 12 }{6} \right) \frac{\partial^{4}\phi}{\partial x^{4}} } + \mathcal{O}(\Delta x)^{4}\\ 2D-Isotropic: \qquad & \displaystyle\frac{\partial\phi}{\partial t} = \nu\left(\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\partial^{2}\phi}{\partial y^{2}}\right) - \textcolor{red} { \Delta_{x}^{2} R1 \left( \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\partial^{4}\phi}{\partial y^{4}} \right) } - \textcolor{red} { \Delta_{x}^{2} R12 \left( \frac{\partial^{4}\phi}{\partial x^{2}\partial y^{2}} \right) } + \mathcal{O}(\Delta x)^{4}\\ 2D-Anisotropic: \qquad & \displaystyle\frac{\partial\phi}{\partial t} = \nu_{x}\frac{\partial^{2}\phi}{\partial x^{2}} + \nu_{y}\frac{\partial^{2}\phi}{\partial y^{2}} - \textcolor{red} { \Delta_{x}^{2} R1 \frac{\partial^{4}\phi}{\partial x^{4}} + \Delta_{x}^{2} R2\frac{\partial^{4}\phi}{\partial y^{4}} } - \textcolor{red} { \Delta_{x}^{2} R12 \left( \frac{\partial^{4}\phi}{\partial x^{2}\partial y^{2}} \right) } + \mathcal{O}(\Delta x)^{4}\\ 3D-Isotropic: \qquad & \displaystyle\frac{\partial\phi}{\partial t} = \nu\left(\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\partial^{2}\phi}{\partial y^{2}} + \frac{\partial^{2}\phi}{\partial z^{2}}\right) - \textcolor{red} { \Delta_{x}^{2} R1 \left( \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\partial^{4}\phi}{\partial y^{4}} \right) } - \textcolor{red} { \Delta_{x}^{2} R12 \left( \frac{\partial^{4}\phi}{\partial x^{2}\partial y^{2}} + \frac{\partial^{4}\phi}{\partial x^{2}\partial z^{2}} + \frac{\partial^{4}\phi}{\partial y^{2}\partial z^{2}} \right) } + \mathcal{O}(\Delta x)^{4}\\ 3D-Anisotropic: \qquad & \displaystyle\frac{\partial\phi}{\partial t} = \nu_{x}\frac{\partial^{2}\phi}{\partial x^{2}} + \nu_{y}\frac{\partial^{2}\phi}{\partial y^{2}} + \nu_{z}\frac{\partial^{2}\phi}{\partial z^{2}} - \textcolor{red} { \Delta_{x}^{2} \left( R1\frac{\partial^{4}\phi}{\partial x^{4}} + R2\frac{\partial^{4}\phi}{\partial y^{4}} + R3\frac{\partial^{4}\phi}{\partial y^{4}} \right) } - \textcolor{red} { \Delta_{x}^{2} \left( R12\frac{\partial^{4}\phi}{\partial x^{2}\partial y^{2}} + R13\frac{\partial^{4}\phi}{\partial x^{2}\partial z^{2}} + R23\frac{\partial^{4}\phi}{\partial y^{2}\partial z^{2}} \right) } + \mathcal{O}(\Delta x)^{4} \end{array} \end{split}\]

Interpretation for Diffusive Equation in a MRT-Framework

Extending the approach Developed above for now considering a sistems of multiple relaxatio time. The present development follow sequence presente in Lin et al. [10].

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 1 & 1\\0 & 1 & -1\\-2 & 1 & 1\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}\phi_{x}^{t} \left(w_{0} + w_{1} + w_{2}\right)\\\phi_{x}^{t} \left(w_{1} - w_{2}\right)\\\phi_{x}^{t} \left(- 2 w_{0} + w_{1} + w_{2}\right)\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}s_{0} & 0 & 0\\0 & s_{1} & 0\\0 & 0 & s_{2}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{1}{3} & 0 & - \frac{1}{3}\\\frac{1}{3} & \frac{1}{2} & \frac{1}{6}\\\frac{1}{3} & - \frac{1}{2} & \frac{1}{6}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{s_{0} \left(f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3} - \frac{s_{2} \left(- 2 f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3}\\\frac{s_{0} \left(f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3} + \frac{s_{1} \left(f_{1,x+c_{i}}^t - f_{2,x+c_{i}}^t\right)}{2} + \frac{s_{2} \left(- 2 f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{6}\\\frac{s_{0} \left(f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3} - \frac{s_{1} \left(f_{1,x+c_{i}}^t - f_{2,x+c_{i}}^t\right)}{2} + \frac{s_{2} \left(- 2 f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{6}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}f_{0,x+c_{i}}^t\\f_{1,x+c_{i}}^t\\f_{2,x+c_{i}}^t\end{matrix}\right] = \left[\begin{matrix}f_{0,x+c_{i}}^t - \frac{s_{0} \left(- \phi_{x-c_{i}}^{t} w_{0} - \phi_{x-c_{i}}^{t} w_{1} - \phi_{x-c_{i}}^{t} w_{2} + f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3} + \frac{s_{2} \left(2 \phi_{x-c_{i}}^{t} w_{0} - \phi_{x-c_{i}}^{t} w_{1} - \phi_{x-c_{i}}^{t} w_{2} - 2 f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3}\\f_{1,x+c_{i}}^t - \frac{s_{0} \left(- \phi_{x-c_{i}}^{t} w_{0} - \phi_{x-c_{i}}^{t} w_{1} - \phi_{x-c_{i}}^{t} w_{2} + f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3} - \frac{s_{1} \left(- \phi_{x-c_{i}}^{t} w_{1} + \phi_{x-c_{i}}^{t} w_{2} + f_{1,x+c_{i}}^t - f_{2,x+c_{i}}^t\right)}{2} - \frac{s_{2} \left(2 \phi_{x-c_{i}}^{t} w_{0} - \phi_{x-c_{i}}^{t} w_{1} - \phi_{x-c_{i}}^{t} w_{2} - 2 f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{6}\\f_{2,x+c_{i}}^t - \frac{s_{0} \left(- \phi_{x-c_{i}}^{t} w_{0} - \phi_{x-c_{i}}^{t} w_{1} - \phi_{x-c_{i}}^{t} w_{2} + f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{3} + \frac{s_{1} \left(- \phi_{x-c_{i}}^{t} w_{1} + \phi_{x-c_{i}}^{t} w_{2} + f_{1,x+c_{i}}^t - f_{2,x+c_{i}}^t\right)}{2} - \frac{s_{2} \left(2 \phi_{x-c_{i}}^{t} w_{0} - \phi_{x-c_{i}}^{t} w_{1} - \phi_{x-c_{i}}^{t} w_{2} - 2 f_{0,x+c_{i}}^t + f_{1,x+c_{i}}^t + f_{2,x+c_{i}}^t\right)}{6}\end{matrix}\right]\end{split}\]

\[\displaystyle f_{0,x}^{t+1} = \phi_{x}^{t} s_{2} w_{0} - f_{0,x}^t s_{2} + f_{0,x}^t\]

\[\displaystyle f_{1,x}^{t+1} = - \frac{\phi_{x-1}^{t} s_{2} w_{0}}{2} + \frac{f_{0,x-1}^t s_{2}}{2} - \frac{f_{1,x-1}^t s_{1}}{2} + f_{1,x-1}^t + \frac{f_{2,x-1}^t s_{1}}{2}\]

\[\displaystyle f_{2,x}^{t+1} = - \frac{\phi_{x+1}^{t} s_{2} w_{0}}{2} + \frac{f_{0,x+1}^t s_{2}}{2} + \frac{f_{1,x+1}^t s_{1}}{2} - \frac{f_{2,x+1}^t s_{1}}{2} + f_{2,x+1}^t\]

\[\displaystyle \phi_{x}^{t+1} = f_{0,x}^t + f_{1,x-1}^t + f_{2,x+1}^t + \frac{s_{1} \left(f_{1,x+1}^t - f_{1,x-1}^t - f_{2,x+1}^t + f_{2,x-1}^t\right)}{2} + s_{2} w_{0} \left(- \frac{\phi_{x+1}^{t}}{2} - \frac{\phi_{x-1}^{t}}{2} + \phi_{x}^{t}\right) + s_{2} \left(\frac{f_{0,x+1}^t}{2} + \frac{f_{0,x-1}^t}{2} - f_{0,x}^t\right)\]

Reinterpreting the above equation, we can derive

(61)\[ f_{0,x}^{t}+f_{1,x+1}^{t}+f_{2,x-1}^{t}=\phi^{t-1}_{x} \]

\[\displaystyle \phi_{x}^{t+1} = f_{0,x}^t + f_{1,x-1}^t + f_{2,x+1}^t + \frac{s_{1} \left(\phi_{x}^{t-1} - f_{0,x}^t - f_{1,x-1}^t - f_{2,x+1}^t\right)}{2} + s_{2} w_{0} \left(- \frac{\phi_{x+1}^{t}}{2} - \frac{\phi_{x-1}^{t}}{2} + \phi_{x}^{t}\right) + s_{2} \left(\frac{f_{0,x+1}^t}{2} + \frac{f_{0,x-1}^t}{2} - f_{0,x}^t\right)\]

To repalce the term $f_{0,x}^{t}+f_{1,x-1}^{t}+f_{2,x+1}^{t}$, in the above equation we use the $\sum_{i}f_{i,x}^{t}=\phi_{x}^{t}$:

\[\displaystyle f_{0,x}^t + f_{1,x-1}^t + f_{2,x+1}^t = \phi_{x+1}^{t} + \phi_{x-1}^{t} - \phi_{x}^{t-1} - f_{0,x+1}^t - f_{0,x-1}^t + 2 f_{0,x}^t\]

\[\displaystyle \phi_{x}^{t+1} = \phi_{x+1}^{t} \left(- \frac{s_{1}}{2} - \frac{s_{2} w_{0}}{2} + 1\right) + \phi_{x-1}^{t} \left(- \frac{s_{1}}{2} - \frac{s_{2} w_{0}}{2} + 1\right) + \phi_{x}^{t-1} \left(s_{1} - 1\right) + \phi_{x}^{t} s_{2} w_{0} + \left(f_{0,x+1}^t + f_{0,x-1}^t - 2 f_{0,x}^t\right) \left(\frac{s_{1}}{2} + \frac{s_{2}}{2} - 1\right)\]

\[\displaystyle f_{0,x+1}^t + f_{0,x-1}^t - 2 f_{0,x}^t = - 2 \phi_{x}^{t} + f_{0,x+1}^t + f_{0,x-1}^t + 2 f_{1,x}^t + 2 f_{2,x}^t\]

Replace the distribution functions in above relation by your collisional form

\[\displaystyle f_{0,x+1}^t + f_{0,x-1}^t - 2 f_{0,x}^t = - 2 \phi_{x}^{t} + f_{0,x+1}^{t-1} + f_{0,x-1}^{t-1} + 2 f_{1,x-1}^{t-1} + 2 f_{2,x+1}^{t-1} + s_{1} \left(f_{1,x+1}^{t-1} - f_{1,x-1}^{t-1} - f_{2,x+1}^{t-1} + f_{2,x-1}^{t-1}\right)\]

Replacing in the above relation the terms $f_{1,x-1}^{t-1}$ and $f_{2,x+1}^{t-1}$ usgin Eq. (61):

\[\displaystyle f_{0,x+1}^t + f_{0,x-1}^t - 2 f_{0,x}^t = 2 \phi_{x+1}^{t-1} + 2 \phi_{x-1}^{t-1} - 2 \phi_{x}^{t} - f_{0,x+1}^{t-1} - f_{0,x-1}^{t-1} - 2 f_{1,x+1}^{t-1} - 2 f_{2,x-1}^{t-1} + s_{1} \left(- \phi_{x+1}^{t-1} - \phi_{x-1}^{t-1} + f_{0,x+1}^{t-1} + f_{0,x-1}^{t-1} + 2 f_{1,x+1}^{t-1} + 2 f_{2,x-1}^{t-1}\right)\]

\[\displaystyle f_{0,x+1}^t + f_{0,x-1}^t - 2 f_{0,x}^t = - \phi_{x+1}^{t-1} \left(s_{1} - 2\right) - \phi_{x-1}^{t-1} \left(s_{1} - 2\right) - 2 \phi_{x}^{t} + \left(s_{1} - 1\right) \left(2 \phi_{x}^{t-2} + f_{0,x+1}^{t-1} + f_{0,x-1}^{t-1} - 2 f_{0,x}^{t-1}\right)\]

Replacing the above relation in $\phi_{x}^{t+1}$ term:

\[\displaystyle \phi_{x}^{t+1} = \phi_{x}^{t-1} \left(s_{1} - 1\right) - \phi_{x}^{t} \left(s_{1} - s_{2} w_{0} + s_{2} - 2\right) + \left(- \frac{\phi_{x+1}^{t-1}}{2} - \frac{\phi_{x-1}^{t-1}}{2}\right) \left(s_{1}^{2} + s_{1} s_{2} - 4 s_{1} - 2 s_{2} + 4\right) + \left(- \frac{\phi_{x+1}^{t}}{2} - \frac{\phi_{x-1}^{t}}{2}\right) \left(s_{1} + s_{2} w_{0} - 2\right) + \left(s_{1} - 1\right) \left(s_{1} + s_{2} - 2\right) \left(\phi_{x}^{t-2} + \frac{f_{0,x+1}^{t-1}}{2} + \frac{f_{0,x-1}^{t-1}}{2} - f_{0,x}^{t-1}\right)\]

Applying a shift in time for the Eq. XX.

\[\displaystyle \phi_{x}^{t+1} = \frac{\phi_{x+1}^{t-1} s_{1} s_{2} w_{0}}{2} - \frac{\phi_{x+1}^{t-1} s_{1} s_{2}}{2} + \frac{\phi_{x+1}^{t-1} s_{1}}{2} - \frac{\phi_{x+1}^{t-1} s_{2} w_{0}}{2} + \phi_{x+1}^{t-1} s_{2} - \phi_{x+1}^{t-1} - \frac{\phi_{x+1}^{t} s_{1}}{2} - \frac{\phi_{x+1}^{t} s_{2} w_{0}}{2} + \phi_{x+1}^{t} + \frac{\phi_{x-1}^{t-1} s_{1} s_{2} w_{0}}{2} - \frac{\phi_{x-1}^{t-1} s_{1} s_{2}}{2} + \frac{\phi_{x-1}^{t-1} s_{1}}{2} - \frac{\phi_{x-1}^{t-1} s_{2} w_{0}}{2} + \phi_{x-1}^{t-1} s_{2} - \phi_{x-1}^{t-1} - \frac{\phi_{x-1}^{t} s_{1}}{2} - \frac{\phi_{x-1}^{t} s_{2} w_{0}}{2} + \phi_{x-1}^{t} - \phi_{x}^{t-1} s_{1} s_{2} w_{0} + \phi_{x}^{t-1} s_{1} + \phi_{x}^{t-1} s_{2} w_{0} - \phi_{x}^{t-1} + \phi_{x}^{t-2} s_{1} s_{2} - \phi_{x}^{t-2} s_{1} - \phi_{x}^{t-2} s_{2} + \phi_{x}^{t-2} + \phi_{x}^{t} s_{2} w_{0} - \phi_{x}^{t} s_{2} + \phi_{x}^{t}\]

\[\displaystyle \phi_{x+1}^{t-1} \left(\frac{s_{1} s_{2} w_{0}}{2} - \frac{s_{1} s_{2}}{2} + \frac{s_{1}}{2} - \frac{s_{2} w_{0}}{2} + s_{2} - 1\right) + \phi_{x+1}^{t} \left(- \frac{s_{1}}{2} - \frac{s_{2} w_{0}}{2} + 1\right) + \phi_{x-1}^{t-1} \left(\frac{s_{1} s_{2} w_{0}}{2} - \frac{s_{1} s_{2}}{2} + \frac{s_{1}}{2} - \frac{s_{2} w_{0}}{2} + s_{2} - 1\right) + \phi_{x-1}^{t} \left(- \frac{s_{1}}{2} - \frac{s_{2} w_{0}}{2} + 1\right) + \phi_{x}^{t-1} \left(- s_{1} s_{2} w_{0} + s_{1} + s_{2} w_{0} - 1\right) + \phi_{x}^{t-2} \left(s_{1} s_{2} - s_{1} - s_{2} + 1\right) + \phi_{x}^{t} \left(s_{2} w_{0} - s_{2} + 1\right)\]

Fourth-order MRT Diffusive Equation

To recover the partial differetial equation with Fourth-order of convergence error, we expand each $\phi$ term in Taylor series (this link provide a sympy code that reach the above equation Sympy code for FDM expansion of LBM diffusive equation up to Fourth-order):

\[ 0=\frac{\Delta_{t}^{2} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2} \frac{\partial^{2}\phi}{\partial t^{2}} + \frac{\Delta_{t} \Delta_{x}^{2} \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right)}{2} \frac{\partial^{3}\phi}{\partial x^{2}\partial t} + \Delta_{t} s_{1} s_{2} \frac{\partial \phi}{\partial t} + \frac{\Delta_{x}^{4} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right)}{24} \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\Delta_{x}^{2} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right)}{2} \frac{\partial^{2}\phi}{\partial x^{2}} \]

\[ \Delta_{t} s_{1} s_{2} \frac{\partial \phi}{\partial t} = - \frac{\Delta_{x}^{2} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right)}{2} \frac{\partial^{2}\phi}{\partial x^{2}} -\frac{\Delta_{t}^{2} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2} \frac{\partial^{2}\phi}{\partial t^{2}} - \frac{\Delta_{t} \Delta_{x}^{2} \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right)}{2} \frac{\partial^{3}\phi}{\partial x^{2}\partial t} - \frac{\Delta_{x}^{4} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right)}{24} \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial \phi}{\partial t} = \frac{\Delta_{x}^{2} (1-w_{0})\left(2 - s_{1}\right)}{2\Delta_{t} s_{1} } \frac{\partial^{2}\phi}{\partial x^{2}} - \frac{\Delta_{t} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2 s_{1} s_{2}} \frac{\partial^{2}\phi}{\partial t^{2}} - \frac{\Delta_{x}^{2} \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right)}{2 s_{1} s_{2}} \frac{\partial^{3}\phi}{\partial x^{2}\partial t} + \frac{\Delta_{x}^{4} (1-w_{0})\left(2 - s_{1}\right)}{24\Delta_{t} s_{1}} \frac{\partial^{4}\phi}{\partial x^{4}} \]

Considering diffusive coeficient $\nu=(1-w_{0})\left(\frac{1}{s_{1}} - \frac{1}{2}\right)\frac{\Delta_{x}^{2}}{\Delta_{t}}$, we can rewrite the equation

\[ \frac{\partial \phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} -\frac{\Delta_{t} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2 s_{1} s_{2}} \frac{\partial^{2}\phi}{\partial t^{2}} - \frac{ \Delta_{x}^{2} \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right)}{2 s_{1} s_{2}} \frac{\partial^{3}\phi}{\partial x^{2}\partial t} + \frac{\nu\Delta_{x}^{2}}{12} \frac{\partial^{4}\phi}{\partial x^{4}} \]

Using the relations:

\[ \frac{\partial\phi}{\partial t} = \nu \frac{\partial^{2}\phi}{\partial x^{2}}, \qquad \frac{\partial^{2}\phi}{\partial t^{2}} = \nu^{2} \frac{\partial^{4}\phi}{\partial x^{4}}, \qquad \frac{\partial^{3}\phi}{\partial x^{2}\partial t} = \nu \frac{\partial^{4}\phi}{\partial x^{4}}, \]

we can reduce the equation in terms of $\frac{\partial^{4}\phi}{\partial x^{4}}$:

\[ \frac{\partial \phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} - \frac{\Delta_{x}^{2} (1-w_{0})\left(2 - s_{1}\right)}{2 s_{1} }\frac{\nu \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2 s_{1} s_{2}} \frac{\partial^{4}\phi}{\partial x^{4}} - \frac{ \Delta_{x}^{2} \nu \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right)}{2 s_{1} s_{2}} \frac{\partial^{4}\phi}{\partial x^{4}} + \frac{\nu\Delta_{x}^{2}}{12} \frac{\partial^{4}\phi}{\partial x^{4}} \]

\[ \frac{\partial \phi}{\partial t} = \nu\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\Delta_{x}^{2} \nu \left(3 s_{1}^{2} s_{2} w_{0} - 2 s_{1}^{2} s_{2} - 6 s_{1}^{2} w_{0} - 18 s_{1} s_{2} w_{0} + 12 s_{1} s_{2} + 12 s_{1} w_{0} + 12 s_{2} w_{0} - 12 s_{2}\right)}{12 s_{1}^{2} s_{2}} \frac{\partial^{4}\phi}{\partial x^{4}} \]

Making the term $\left(3 s_{1}^{2} s_{2} w_{0} - 2 s_{1}^{2} s_{2} - 6 s_{1}^{2} w_{0} - 18 s_{1} s_{2} w_{0} + 12 s_{1} s_{2} + 12 s_{1} w_{0} + 12 s_{2} w_{0} - 12 s_{2}\right)$ equal to zero e can reach fourth order approximation. For $w_{0}=2/3$ (weight specific from the lattice D1Q3) the relation reduce to $4s_{1}^{2} -8s_{1}+4s_{2}$:

../_images/38ac13614e70ba0e929efd8746c1cb19c718d8edc1e506f411b002b79d3e0886.png

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
A   = 1.0                         # Amplitude
phi = 0.0                         # Wave initial shift
k_w   = 1.0                       # Wave number
c   = 0.0                         # wave speed
nu  = 0.066666666666666666*np.pi # Diffusive coefficient
m_waves = 1                       # Wave Period
L = m_waves * (2.0*np.pi / k_w)   # Domain Lentgh
t_end = 12.0                       # Final time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int16")  
cs=1.0/np.sqrt(3.0);
#--------------------------------- Initialization - Save Data Arrays -----------------------------------------------
cases=4
Nx0 = np.array([10, 20, 40, 80],dtype="int64")
r0 = 2**(0)*np.array([2.0**(1), 2.0**(2), 2.0**(3), 2.0**(4)],dtype="float64")
u_cases = np.empty(cases, dtype=object) # array field used to same data over time
for i in range(cases):
    u_cases[i] = np.zeros((Nx0[i]), dtype="float64")
for case in range(0,cases):
    #------------------------------------- Parameters (numerical units) -------------------------------------------
    Nx=Nx0[case]  # Numerical Length
    x  = np.linspace(0.0, L, Nx, endpoint=False,dtype="float64") # Numerical Domain
    dx = L / Nx # Grid size
    r = 1*r0[case] # dx/dt relation
    dt = dx/r
    ce = c/r
    nue = nu * dt / (dx*dx)
    tau1 = 0.5 + nue / cs**2
    tau2=1.0/((1.0/tau1)*(2.0-(1.0/tau1)))
    nt = int(np.ceil(t_end / dt))
    print(f"dx={dx:.4e},\t dt={dt:.4e}\t nt={nt:d},\t tn={nt*dt:.2f}") # Print values for check
    print(f"c_lbm={ce:.3f},\t tau1={tau1:.4f},\t tau2={tau2:.4f}")
    print(f"s1={1.0/tau1:.4f},\t s2={1.0/tau2:.4f}")
    print(f"nt={nt:.3f}\t, nue={nue:.4f}")
    #--------------------------------- Initialization - LBM - Numerical Arrays -----------------------------------------
    u=np.zeros((Nx),dtype="float64")
    u = A * np.sin(k_w*(x+dx/2) + phi)
    mx=np.zeros((Nx),dtype="float64")
    mxx=np.zeros((Nx),dtype="float64")
    f=np.zeros((3,Nx),dtype="float64")
    feq=np.zeros((3,Nx),dtype="float64")
    fp=np.zeros((3,Nx),dtype="float64")
    for k in range(0,3):
        f[k,:]=w[k]*(u[:]+cx[k]*ce*u[:]/cs**2)
        fp[k,:]=w[k]*(u[:]+cx[k]*ce*u[:]/cs**2)
    #----------------------------------------- Init Loop --------------------------------------------------------
    for t in range(nt*10):
        #========================================= LBM - Solution =============================================
        #--------------------Collision----------------
        for k in range(0,3):
            feq[k,:]= w[k]*(u[:]+cx[k]*(ce*u[:])/cs**2)
        mx = np.einsum('i,ix->x', cx, f-feq)
        mxx = np.einsum('i,ix->x', 3.0*cx*cx-2.0, f-feq)
        for k in range(0,3):
            # fp[k,:]= f[k,:] - w[k]*( cx[k]*mx[:]/(tau1*cs**2) + (cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2*tau2) )
            fp[k,:]= feq[k,:] + w[k]*( (1.0-1.0/tau1)*cx[k]*mx[:]/cs**2 + (1.0-1.0/tau2)*(cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2) )
        #-----------------streaming-------------------
        for k in range(0,3):
            f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
    #----------------------------------------- Maind Loop --------------------------------------------------------
    for t in range(nt):
        #========================================= LBM - Solution =============================================
        #--------------------Collision----------------
        for k in range(0,3):
            feq[k,:]= w[k]*(u[:]+cx[k]*(ce*u[:])/cs**2)
        mx = np.einsum('i,ix->x', cx, f-feq)
        mxx = np.einsum('i,ix->x', 3.0*cx*cx-2.0, f-feq)
        for k in range(0,3):
            # fp[k,:]= f[k,:] - w[k]*( cx[k]*mx[:]/(tau1*cs**2) + (cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2*tau2) )
            fp[k,:]= feq[k,:] + w[k]*( (1.0-1.0/tau1)*cx[k]*mx[:]/cs**2 + (1.0-1.0/tau2)*(cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2) )
        #-----------------streaming-------------------
        for k in range(0,3):
            f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
        #----------------------Macro------------------
        # u[:]=f[0,:]+f[1,:]+f[2,:]
        u=np.einsum('ix->x', f)
    #---------------------save-field-cases--------------
    u_cases[case][:]=u[:]

../_images/53f862192818947da5d4e47e03b24a998a01ac6d5817b73b05e90e7dff81ae4a.png

Sixth-order MRT Diffusive Equation

Extending up to Sixth-order of convergence error (this link provide a sympy code that reach the above equation Sympy code for FDM expansion of LBM diffusive equation up to Sixth-order):

\[ \frac{\partial \phi}{\partial t} = \nu \frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\Delta_{x}^{2} \nu R^{4th}}{12 s_{1}^{2} s_{2}} \frac{\partial^{4}\phi}{\partial x^{4}} - \frac{\Delta_{x}^{4} \nu R^{6th} }{24 s_{1}^{3} s_{2}} \frac{\partial^{6}\phi}{\partial x^{6}}, \]

where $R^{4th}$ and $R^{6th}$ are the terms in fouth and sixth-order error, respectively, given by:

\[ R^{4th}=\left(3 s_{1}^{2} s_{2} w_{0} - 2 s_{1}^{2} s_{2} - 6 s_{1}^{2} w_{0} - 18 s_{1} s_{2} w_{0} + 12 s_{1} s_{2} + 12 s_{1} w_{0} + 12 s_{2} w_{0} - 12 s_{2}\right) \]

\[ R^{6th}=\left(7 s_{1}^{3} s_{2} w_{0}^{2} - 13 s_{1}^{3} s_{2} w_{0} + 6 s_{1}^{3} s_{2} - 6 s_{1}^{3} w_{0}^{2} + 12 s_{1}^{3} w_{0} - 5 s_{1}^{3} - 34 s_{1}^{2} s_{2} w_{0}^{2} + 67 s_{1}^{2} s_{2} w_{0} - 32 s_{1}^{2} s_{2} + 30 s_{1}^{2} w_{0}^{2} - 60 s_{1}^{2} w_{0} + 28 s_{1}^{2} + 52 s_{1} s_{2} w_{0}^{2} - 104 s_{1} s_{2} w_{0} + 52 s_{1} s_{2} - 48 s_{1} w_{0}^{2} + 96 s_{1} w_{0} - 48 s_{1} - 24 s_{2} w_{0}^{2} + 48 s_{2} w_{0} - 24 s_{2} + 24 w_{0}^{2} - 48 w_{0} + 24\right) \]

Making the term $\left(3 s_{1}^{2} s_{2} w_{0} - 2 s_{1}^{2} s_{2} - 6 s_{1}^{2} w_{0} - 18 s_{1} s_{2} w_{0} + 12 s_{1} s_{2} + 12 s_{1} w_{0} + 12 s_{2} w_{0} - 12 s_{2}\right)$ equal to zero e can reach fourth order approximation. For $w_{0}=2/3$ (weight specific from the lattice D1Q3) the relation reduce to $4s_{1}^{2} -8s_{1}+4s_{2}$:

Trying to make both $R^{4th}$ and $R^{6th}$ null for lattice D1Q3 ($w_{0}=2/3$), we have the relation:

\[ R^{4th}=4s_{1}^{2} -8s_{1}+4s_{2}=0 \qquad \rightarrow \qquad s_{2}=s_{1}(2-s_{1}) \]

\[ R^{6th} = \frac{s_{1}^{3}}{3} + \frac{4 s_{1}^{2}}{3} - \frac{16 s_{1}}{3} + s_{2} \left(\frac{4 s_{1}^{3}}{9} - \frac{22 s_{1}^{2}}{9} + \frac{52 s_{1}}{9} - \frac{8}{3}\right) + \frac{8}{3}=0 \qquad \rightarrow \qquad s_{2} = \displaystyle\frac{-\left( \frac{s_{1}^{3}}{3} + \frac{4 s_{1}^{2}}{3} - \frac{16 s_{1}}{3} + \frac{8}{3}\right) }{\left(\frac{4 s_{1}^{3}}{9} - \frac{22 s_{1}^{2}}{9} + \frac{52 s_{1}}{9} - \frac{8}{3}\right)} = \frac{3(2-s_{1})(s_{1}^{2} +6s_{1} -4)}{2(2s_{1}^{3} - 11s_{1}^{2} + 26s_{1} - 12)} \]

../_images/3ec7f6ccd9996079c878642705612fb0f32c89a2b98be0d56ec49f917c4b0be9.png

Difference Finite Scheme of Diffusive Equation derived from Lattice Boltzmann

Interpretation for Diffusive Equation

FDM Interpretation

1st Benchmark: Decay of Sinoidal Wave

LBM D1Q3 Code - FDM Code (Only Diffusive Case \(c=0\))

2nd Benchmark: Temperature ramp

Steady-State Analytical Solution

Transient Analytical Solution

LBM D1Q3 Code - FDM Code

Fourth-order Diffusive Equation

Extension for Higher-Dimensions and Anisotropy

Interpretation for Diffusive Equation in a MRT-Framework

Fourth-order MRT Diffusive Equation

Sixth-order MRT Diffusive Equation