双松弛时间（Two-relaxation-time）格子 Boltzmann 方法是多松弛时间（Multi-Relaxation-Time）LBM 的一种特殊形式。它仅具有两个松弛时间，分为对称和非对称两部分，同时保持了与 LBGK 模型一致的简单算法和计算效率。
对于 t t t x \boldsymbol{x} x c i \boldsymbol{c}_i c i f i ( x , t ) f_i(\boldsymbol{x}, t) f i ( x , t ) Δ t \Delta t Δ t
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − 1 τ + [ f i + − f i ( e q ) + ] − 1 τ − [ f i − − f i ( e q ) − ] + F i Δ t (1) \begin{aligned}
f_i(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t + \Delta t) - f_i(\boldsymbol{x}, t) &=
-\frac{1}{\tau^{+}} \left[ f_i^{+} - f_i^{(\mathrm{eq})+} \right] \\
&\quad -\frac{1}{\tau^{-}} \left[ f_i^{-} - f_i^{(\mathrm{eq})-} \right] \\
&\quad + \mathbb{F}_i \Delta t
\end{aligned}
\tag{1}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − τ + 1 [ f i + − f i ( e q ) + ] − τ − 1 [ f i − − f i ( e q ) − ] + F i Δ t ( 1 )
其中：
τ + \tau^{+} τ + τ − \tau^{-} τ − ν \nu ν Λ = ( τ + − 1 2 ) ( τ − − 1 2 ) \Lambda = (\tau^{+} - \frac{1}{2}) (\tau^{-} - \frac{1}{2}) Λ = ( τ + − 2 1 ) ( τ − − 2 1 ) Λ \Lambda Λ Λ = 1 / 4 \Lambda=1/4 Λ = 1 / 4
τ + = ν c s 2 + 1 2 , τ − = 2 Λ 2 τ + − 1 + 1 2 \tau^{+} = \frac{\nu}{c_s^2} + \frac{1}{2} ,\quad
\tau^{-} = \frac{2 \Lambda}{2 \tau^{+} - 1} + \frac{1}{2}
τ + = c s 2 ν + 2 1 , τ − = 2 τ + − 1 2 Λ + 2 1
f − i ( x , t ) f_{-i}(\boldsymbol{x}, t) f − i ( x , t ) f − i ( e q ) ( x , t ) f_{-i}^{(\mathrm{eq})}(\boldsymbol{x}, t) f − i ( e q ) ( x , t ) c − i = − c i \boldsymbol{c}_{-i}=-\boldsymbol{c}_{i} c − i = − c i
f i + f_i^{+} f i + f i − f_i^{-} f i − 对称 部分和非对称 部分，即
f i + = 1 2 ( f i ( x , t ) + f − i ( x , t ) ) , f i − = 1 2 ( f i ( x , t ) − f − i ( x , t ) ) f_i^{+} = \frac{1}{2}\left( f_i(\boldsymbol{x}, t) + f_{-i}(\boldsymbol{x}, t) \right) ,\quad
f_i^{-} = \frac{1}{2}\left( f_i(\boldsymbol{x}, t) - f_{-i}(\boldsymbol{x}, t) \right)
f i + = 2 1 ( f i ( x , t ) + f − i ( x , t ) ) , f i − = 2 1 ( f i ( x , t ) − f − i ( x , t ) )
同理， f i ( e q ) + f_i^{(\mathrm{eq})+} f i ( e q ) + f i ( e q ) − f_i^{(\mathrm{eq})-} f i ( e q ) − f i ( e q ) f_i^{(\mathrm{eq})} f i ( e q )
f i e q = w i ρ [ 1 + c i ⋅ u c s 2 + ( c i ⋅ u ) 2 2 c s 4 − u 2 2 c s 2 ] , f i ( e q ) + = 1 2 ( f i ( e q ) ( x , t ) + f − i ( e q ) ( x , t ) ) = w i ρ [ 1 + ( c i ⋅ u ) 2 2 c s 4 − u 2 2 c s 2 ] , f i ( e q ) − = 1 2 ( f i ( e q ) ( x , t ) − f − i ( e q ) ( x , t ) ) = w i ρ c i ⋅ u c s 2 \begin{aligned}
f^{\mathrm{eq}}_i &= w_i \rho \left[ 1 + \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{u^2}{2 c_s^2} \right] ,\\
f_i^{(\mathrm{eq})+} &= \frac{1}{2}\left( f_i^{(\mathrm{eq})}(\boldsymbol{x}, t) + f_{-i}^{(\mathrm{eq})}(\boldsymbol{x}, t) \right) =w_i \rho \left[ 1 + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u})^2}{2 c_s^4} - \frac{u^2}{2 c_s^2} \right] ,\\
f_i^{(\mathrm{eq})-} &= \frac{1}{2}\left( f_i^{(\mathrm{eq})}(\boldsymbol{x}, t) - f_{-i}^{(\mathrm{eq})}(\boldsymbol{x}, t) \right) = w_i \rho \frac{\boldsymbol{c}_i \cdot \boldsymbol{u}}{c_s^2}
\end{aligned}
f i e q f i ( e q ) + f i ( e q ) − = w i ρ [ 1 + c s 2 c i ⋅ u + 2 c s 4 ( c i ⋅ u ) 2 − 2 c s 2 u 2 ] , = 2 1 ( f i ( e q ) ( x , t ) + f − i ( e q ) ( x , t ) ) = w i ρ [ 1 + 2 c s 4 ( c i ⋅ u ) 2 − 2 c s 2 u 2 ] , = 2 1 ( f i ( e q ) ( x , t ) − f − i ( e q ) ( x , t ) ) = w i ρ c s 2 c i ⋅ u
下面记 t i = w i / c s 2 t_i = {w_i} / {c_s^2} t i = w i / c s 2 F i + = 1 2 ( F i + F − i ) F_{i}^{+} = \frac{1}{2} (\mathbb{F}_i + \mathbb{F}_{-i}) F i + = 2 1 ( F i + F − i ) F i − = 1 2 ( F i − F − i ) F_{i}^{-} = \frac{1}{2} (\mathbb{F}_i - \mathbb{F}_{-i}) F i − = 2 1 ( F i − F − i )
Postma 等的研究对几种主流的源项模型进行了研究。不过他原文里列出的源项方程我还有些看不懂，我就照着Bawazeer等的文章 抄了两个简单的源项：
（1）Buick and Greated :
F i = ( 1 − 1 2 τ ) t q c i ⋅ F \mathbb{F}_i = \left(1 - \frac{1}{2 \tau} \right) t_q \boldsymbol{c}_{i} \cdot
\boldsymbol{F}
F i = ( 1 − 2 τ 1 ) t q c i ⋅ F
转写为正负两项后表示为：
F i + = 0 , F i − = ( 1 − 1 2 τ − ) t i c i ⋅ F F_{i}^{+} = 0, \quad
F_{i}^{-} = \left(1 - \frac{1}{2 \tau^{-}} \right) t_i \boldsymbol{c}_{i} \cdot
\boldsymbol{F}
F i + = 0 , F i − = ( 1 − 2 τ − 1 ) t i c i ⋅ F
（2）Guo et al. :
F i = ( 1 − 1 2 τ ) w i ⋅ [ c i − u c s 2 + ( c i ⋅ u ) c i c s 4 ] ⋅ F \mathbb{F}_{i} = \left(1 - \frac{1}{2 \tau} \right) w_i \cdot
\left[\frac{\boldsymbol{c}_i - \boldsymbol{u}}{c_s^2} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u}) \boldsymbol{c}_i}{c_s^4}\right] \cdot \boldsymbol{F}
F i = ( 1 − 2 τ 1 ) w i ⋅ [ c s 2 c i − u + c s 4 ( c i ⋅ u ) c i ] ⋅ F
转写为正负两项后表示为：
F i + = ( 1 − 1 2 τ + ) t i ⋅ [ − u + ( c i ⋅ u ) c i c s 2 ] ⋅ F F i − = ( 1 − 1 2 τ − ) t i c i ⋅ F \begin{aligned}
F_i^{+} =& \left(1 - \frac{1}{2 \tau^{+}} \right) t_i \cdot \left[-\boldsymbol{u} + \frac{(\boldsymbol{c}_i \cdot \boldsymbol{u}) \boldsymbol{c}_i}{c_s^2}\right] \cdot \boldsymbol{F} \\
F_i^{-} =& \left(1 - \frac{1}{2 \tau^{-}} \right) t_i \boldsymbol{c}_i \cdot \boldsymbol{F} \\
\end{aligned}
F i + = F i − = ( 1 − 2 τ + 1 ) t i ⋅ [ − u + c s 2 ( c i ⋅ u ) c i ] ⋅ F ( 1 − 2 τ − 1 ) t i c i ⋅ F
宏观密度 ρ \rho ρ j \boldsymbol{j} j
ρ = ∑ i = 0 b − 1 f i , j = ρ u = ∑ i = 0 b − 1 f i c i + F 2 \rho = \sum\limits_{i=0}^{b-1} f_i ,\quad
\boldsymbol{j} = \rho \boldsymbol{u} = \sum\limits_{i=0}^{b-1} f_i \boldsymbol{c}_{i} + \frac{\boldsymbol{F}}{2}
ρ = i = 0 ∑ b − 1 f i , j = ρ u = i = 0 ∑ b − 1 f i c i + 2 F
(1) 计算 f i + f_i^{+} f i + f i − f_i^{-} f i − f i ( e q ) + f_i^{(\mathrm{eq})+} f i ( e q ) + f i ( e q ) − f_i^{(\mathrm{eq})-} f i ( e q ) −
(2) 计算碰撞步
f i ( p c ) ( x , t ) = f i ( x , t ) − 1 τ + [ f i + − f i ( e q ) + ] − 1 τ − [ f i − − f i ( e q ) − ] + F i Δ t (2) f_i^{(\mathrm{pc})}(\boldsymbol{x}, t) = f_i(\boldsymbol{x}, t)
-\frac{1}{\tau^{+}} \left[ f_i^{+} - f_i^{(\mathrm{eq})+} \right]
-\frac{1}{\tau^{-}} \left[ f_i^{-} - f_i^{(\mathrm{eq})-} \right] + \mathbb{F}_i \Delta t
\tag{2}
f i ( p c ) ( x , t ) = f i ( x , t ) − τ + 1 [ f i + − f i ( e q ) + ] − τ − 1 [ f i − − f i ( e q ) − ] + F i Δ t ( 2 )
(3) 计算流动步（与常规LBM一致）
(4) 更新宏观量。
将式(1)整理，可得式(3)：
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − ( 1 2 τ + + 1 2 τ − ) [ f i ( x , t ) − f i ( e q ) ] − ( 1 2 τ + − 1 2 τ − ) [ f − i ( x , t ) − f − i ( e q ) ] + F i Δ t (3) \begin{aligned}
f_i(\boldsymbol{x}+\boldsymbol{c}_i \Delta t, t + \Delta t) - f_i(\boldsymbol{x}, t) &= - \left(\frac{1}{2 \tau^{+}} + \frac{1}{2 \tau^{-}}\right) [f_i(\boldsymbol{x}, t) - f^{(\mathrm{eq})}_i] \\
&\quad -\left(\frac{1}{2 \tau^{+}} - \frac{1}{2 \tau^{-}}\right) [f_{-i}(\boldsymbol{x}, t) - f^{(\mathrm{eq})}_{-i}] \\
&\quad + \mathbb{F}_i \Delta t
\end{aligned}
\tag{3}
f i ( x + c i Δ t , t + Δ t ) − f i ( x , t ) = − ( 2 τ + 1 + 2 τ − 1 ) [ f i ( x , t ) − f i ( e q ) ] − ( 2 τ + 1 − 2 τ − 1 ) [ f − i ( x , t ) − f − i ( e q ) ] + F i Δ t ( 3 )
所以，当 τ + = τ − \tau^{+} = \tau^{-} τ + = τ − Λ = ν 2 c s 4 \Lambda = \frac{\nu^2}{c_s^4} Λ = c s 4 ν 2 τ + \tau^{+} τ +
下文为方便起见，记 w 1 = ( 1 τ + + 1 τ − ) / 2 \mathtt{w}_1 = (\frac{1}{\tau^{+}} + \frac{1}{\tau^{-}})/2 w 1 = ( τ + 1 + τ − 1 ) / 2 w 2 = ( 1 τ + − 1 τ − ) / 2 \mathtt{w}_2 = (\frac{1}{\tau^{+}} - \frac{1}{\tau^{-}})/2 w 2 = ( τ + 1 − τ − 1 ) / 2
先定义：
f = [ f 0 ( x , t ) , . . . , f 8 ( x , t ) ] T \boldsymbol{f} = [f_0(\boldsymbol{x}, t), ..., f_8(\boldsymbol{x}, t)]^\mathrm{T} f = [ f 0 ( x , t ) , . . . , f 8 ( x , t ) ] T f ( p c ) = [ f 0 ( p c ) ( x , t ) , . . . , f 8 ( p c ) ( x , t ) ] T \boldsymbol{f}^{(\mathrm{pc})} = [f_0^{(\mathrm{pc})}(\boldsymbol{x}, t), ..., f_8^{(\mathrm{pc})}(\boldsymbol{x}, t)]^\mathrm{T} f ( p c ) = [ f 0 ( p c ) ( x , t ) , . . . , f 8 ( p c ) ( x , t ) ] T f ( e q ) = [ f 0 ( e q ) ( x , t ) , . . . , f 8 ( e q ) ( x , t ) ] T \boldsymbol{f}^{(\mathrm{eq})} = [f_0^{(\mathrm{eq})}(\boldsymbol{x}, t), ..., f_8^{(\mathrm{eq})}(\boldsymbol{x}, t)]^\mathrm{T} f ( e q ) = [ f 0 ( e q ) ( x , t ) , . . . , f 8 ( e q ) ( x , t ) ] T
则 TRT 方程的碰撞步在D2Q9模型下的矩阵形式为：
f ( p c ) = [ H ] ( f − f ( e q ) ) (4) \boldsymbol{f}^{(\mathrm{pc})} = [\bold{H}] (\boldsymbol{f} - \boldsymbol{f}^{(\mathrm{eq})})
\tag{4}
f ( p c ) = [ H ] ( f − f ( e q ) ) ( 4 )
其中矩阵 [ H ] [\bold{H}] [ H ]
[ H ] = [ w 1 + w 2 w 1 w 2 w 1 w 2 w 2 w 1 w 2 w 1 w 1 w 2 w 1 w 2 w 2 w 1 w 2 w 1 ] [\bold{H}] = \begin{bmatrix}
\mathtt{w}_1+\mathtt{w}_2 &&&&&&&\\
& \mathtt{w}_1 & & \mathtt{w}_2 &&&&&\\
&& \mathtt{w}_1 & & \mathtt{w}_2 &&&&\\
& \mathtt{w}_2 & & \mathtt{w}_1 &&&&&\\
&& \mathtt{w}_2 & & \mathtt{w}_1 &&&&\\
&&&&& \mathtt{w}_1 & & \mathtt{w}_2 &\\
&&&&&& \mathtt{w}_1 & & \mathtt{w}_2 \\
&&&&& \mathtt{w}_2 & & \mathtt{w}_1 &\\
&&&&&& \mathtt{w}_2 & & \mathtt{w}_1
\end{bmatrix}
[ H ] = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ w 1 + w 2 w 1 w 2 w 1 w 2 w 2 w 1 w 2 w 1 w 1 w 2 w 1 w 2 w 2 w 1 w 2 w 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤
对比 MRT 方程的碰撞步在D2Q9模型下的矩阵形式：
f ( p c ) = [ M − 1 ] [ S ] [ M ] ( f − f ( e q ) ) (5) \boldsymbol{f}^{(\mathrm{pc})} = [\bold{M}^{-1}][\bold{S}][\bold{M}] (\boldsymbol{f} - \boldsymbol{f}^{(\mathrm{eq})})
\tag{5}
f ( p c ) = [ M − 1 ] [ S ] [ M ] ( f − f ( e q ) ) ( 5 )
其中 [ S ] [\bold{S}] [ S ] [ M ] [\bold{M}] [ M ]
[ M ] = [ 1 1 1 1 1 1 1 1 1 − 4 − 1 − 1 − 1 − 1 2 2 2 2 4 − 2 − 2 − 2 − 2 1 1 1 1 0 1 0 − 1 0 1 − 1 − 1 1 0 − 2 0 2 0 1 − 1 − 1 1 0 0 1 0 − 1 1 1 − 1 − 1 0 0 − 2 0 2 1 1 − 1 − 1 0 1 − 1 1 − 1 0 0 0 0 0 0 0 0 0 1 − 1 1 − 1 ] [\bold{M}] = \begin{bmatrix}
1&1&1&1&1&1&1&1&1\\
-4&-1&-1&-1&-1&2&2&2&2\\
4&-2&-2&-2&-2&1&1&1&1\\
0&1&0&-1&0&1&-1&-1&1\\
0&-2&0&2&0&1&-1&-1&1\\
0&0&1&0&-1&1&1&-1&-1\\
0&0&-2&0&2&1&1&-1&-1\\
0&1&-1&1&-1&0&0&0&0\\
0&0&0&0&0&1&-1&1&-1
\end{bmatrix}
[ M ] = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1 − 4 4 0 0 0 0 0 0 1 − 1 − 2 1 − 2 0 0 1 0 1 − 1 − 2 0 0 1 − 2 − 1 0 1 − 1 − 2 − 1 2 0 0 1 0 1 − 1 − 2 0 0 − 1 2 − 1 0 1 2 1 1 1 1 1 0 1 1 2 1 − 1 − 1 1 1 0 − 1 1 2 1 − 1 − 1 − 1 − 1 0 1 1 2 1 1 1 − 1 − 1 0 − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤
求解方程 [ H ] = [ M − 1 ] [ S ] [ M ] [\bold{H}] = [\bold{M}^{-1}][\bold{S}][\bold{M}] [ H ] = [ M − 1 ] [ S ] [ M ]
[ S ] = [ M ] [ H ] [ M − 1 ] = d i a g ( w 1 + w 2 , w 1 + w 2 , w 1 + w 2 , w 1 − w 2 , w 1 − w 2 , w 1 − w 2 , w 1 − w 2 , w 1 + w 2 , w 1 + w 2 ) \begin{aligned}
[\bold{S}] =& [\bold{M}][\bold{H}][\bold{M}^{-1}] \\
=&\mathrm{diag}(
\mathtt{w}_1+\mathtt{w}_2,
\mathtt{w}_1+\mathtt{w}_2,
\mathtt{w}_1+\mathtt{w}_2,\\&\quad
\mathtt{w}_1-\mathtt{w}_2,
\mathtt{w}_1-\mathtt{w}_2,
\mathtt{w}_1-\mathtt{w}_2,\\&\quad
\mathtt{w}_1-\mathtt{w}_2,
\mathtt{w}_1+\mathtt{w}_2,
\mathtt{w}_1+\mathtt{w}_2
)
\end{aligned}
[ S ] = = [ M ] [ H ] [ M − 1 ] d i a g ( w 1 + w 2 , w 1 + w 2 , w 1 + w 2 , w 1 − w 2 , w 1 − w 2 , w 1 − w 2 , w 1 − w 2 , w 1 + w 2 , w 1 + w 2 )
可见 TRT 模型就是 MRT 模型的特例，将 MRT 中所使用的诸多松弛系数减少至到 2 个（1 τ + \frac{1}{\tau^{+}} τ + 1 1 τ − \frac{1}{\tau^{-}} τ − 1