格子Boltzmann方法(lattice Boltzmann method, LBM)已用于模拟流体流动^[1], 并扩展到模拟非线性演化方程(nonliner evolution equations, NLEEs), 如各向异性色散方程^[2], 变系数对流扩散方程^[3], Burgers方程^[4], Korteweg-de Vries (KdV)方程^[5], 广义Boussinesq方程^[6]和变系数Gardner方程^[7]等.

与根据时间和空间直接离散控制方程的传统宏观数值方法不同, 格子Boltzmann方法(LBM)是一种介观的数值计算方法^[8]. LB方程基于一个基本的离散速度动力学方程, 它包括粒子碰撞和迁移基本过程, 这两个过程可以清楚地还原出各种宏观物理现象. 为了将微观LB输运方程与宏观流体动力学方程联系起来, 文献[9]使用Chapmann-Enskog多尺度展开技术, 把Knudsen数$ \epsilon $(假设平均自由程$ \mathit{l} $与宏观特征长度L的比值很小)作为展开参数, 对LB方程进行逐阶展开. 文献[10]使用直接Taylor展开方法有效地恢复出Navier-Stokes方程和非线性对流扩散方程, 并指出直接Taylor展开是比Chapmann-Enskog多尺度展开更有效的方法.

基于LBM的NLEEs的数值研究一般考虑常系数NLEEs, 在物理环境中, 当考虑到介质的不均匀性和边界的不均一性时, 变系数非线性偏微分方程比常系数偏微分方程更能反映实际情况. 文献[11], 着重于一个具有时变系数的广义Gardner方程, 方程(1.1)被用来模拟弱非线性长波在KdV型介质中的传播, 这种介质的特征是色散和非线性系数会发生变化, 可以描述大气阻塞现象.

(1.1) $ \begin{equation} u_{t} -6f_{0}(t) uu_{x} +f_{2}(t) u_{xxx} -f_{3}(t) u_{x} +f_{4} (t)(Cu+xu_{x} ) = 0, \end{equation} $

其中, $ u $是关于空间$ x $和时间$ t $的波幅函数, $ C $是非零常数, $ f_{0} (t), f_{1} (t) $和$ f_{2} (t) $分别对应于二次非线性和三次非线性和色散系数, $ -f_{3} (t)+f_{4} (t)x $和$ f_{4} (t)C $分别表示耗散项和阻尼项.

文献[12]研究了一类变系数的复合Korteweg-de Vries-Burgers(vc-cKdVB)方程, 它是一种在固体材料、等离子体、流体等领域广泛使用的物理模型.

(1.2) $ \begin{equation} u_{t} +l(t)uu_{x} +m(t)u^{2p}u_{x} +n(t)u_{xx} -r(t)u_{xxx} +s(t)u_{x} = 0, \end{equation} $

其中, 参数$ p $是非负整数, 系数$ l (t), m (t), n (t), r (t) $和$ s (t) $是关于时间$ t $的解析函数, 下标表示偏导数. 在方程(1.2)中, 由于非线性项$ u^{2p}u_{x} $项为偶数次幂, 本文考虑更一般的情形, 将非线性项$ u^{2p}u_{x} $推广到$ u^{p}u_{x} $, 利用LBM模拟非线性项为任意次幂的偏微分方程, 即该文研究的变系数KdV-Burgers方程. 方程(1.3)第二项、第三项为平流项, 等号右边的项依次为扩散项、反扩散项或超扩散项.

(1.3) $ \begin{equation} u_{t} + l(t)uu_{x} +m(t)u^{p} u_{x} = n(t)u_{xx} -r(t)u_{xxx}. \end{equation} $

本文的主要内容如下: 在第2节中, 利用Taylor展开和Chapmann-Enskog多尺度展开技术, 证明了带有修正项的格子Boltzmann模型与变系数KdV-Burgers宏观方程是相容的, 并由宏观方程的系数构造了平衡态分布函数和修正函数. 在第3节中, 利用MATLAB进行数值模拟, 给出不同时刻变系数方程的时空演化图和误差分析图, 并对模型的空间精度和时间精度进行数值分析. 第4节对本文研究的结果进行了总结和展望.

带有修正项的离散分布函数的LB方程$ f_{i} $为

(2.1) $ \begin{equation} f_{i}\left(x+c_{i} \Delta t, t+\Delta t\right)-f_{i}(x, t) = -\frac{1}{\tau}\left[f_{i}(x, t)-f_{i}^{e q}(x, t)\right]+\Delta t h_{i}(x, t), \end{equation} $

其中, $ c_{i} $为离散晶格速度, $ \tau $为松弛时间, $ h_{i} $为修正项, $ f_{i}^{eq} $为平衡态分布函数, $ \left |c\right | $为特征晶格速度, 假设为1.

对(2.1)式进行多元Taylor展开, 并保留到$ O(\epsilon ^{4} ) $项, 可以得到

(2.2) $ \begin{eqnarray} &&\Delta t\left(\frac{\partial}{\partial t}+c_{i} \frac{\partial}{\partial x}\right) f_{i}+\frac{\Delta t^{2}}{2}\left(\frac{\partial}{\partial t}+c_{i} \frac{\partial}{\partial x}\right)^{2} f_{i}+\frac{\Delta t^{3}}{6}\left(\frac{\partial}{\partial t}+c_{i} \frac{\partial}{\partial x}\right)^{3} f_{i}+O(\epsilon ^{4} )\\& = &-\frac{1}{\tau}\left(f_{i}-f_{i}^{e q}\right)+\Delta t h_{i}. \end{eqnarray} $

利用Chapmann-Enskog多尺度展开技术对时间和空间导数, 分布函数$ f_{i} $, 修正函数$ h_{i} $展开

(2.3) $ \begin{equation} \frac{\partial}{\partial x} = \epsilon \frac{\partial}{\partial x_{1}}+O\left(\epsilon ^{2}\right), \end{equation} $

(2.4) $ \begin{equation} \frac{\partial}{\partial t} = \epsilon \frac{\partial}{\partial t_{1}}+\epsilon ^{2} \frac{\partial}{\partial t_{2}}+\epsilon ^{3} \frac{\partial}{\partial t_{3}}+O\left(\epsilon ^{4}\right), \end{equation} $

(2.5) $ \begin{equation} f_{i} = f_{i}^{(0)}+\epsilon f_{i}^{(1)}+\epsilon ^{2} f_{i}^{(2)}+\epsilon ^{3} f_{i}^{(3)}+O\left(\epsilon ^{4}\right), \end{equation} $

(2.6) $ \begin{equation} h_{i} = \epsilon h_{i}^{(1)}+O\left(\epsilon ^{2}\right), \end{equation} $

其中, $ \epsilon $是Knudsen数, 定义为$ \epsilon = l/L $, $ l $是平均自由程, $ L $是特征长度, $ \epsilon $可取为时间步长$ \Delta t $^[13]. 把(2.3)–(2.6)式代入方程(2.2), 我们可以得到关于$ \epsilon $的0到3阶微分方程

(2.7) $ \begin{equation} O\left(\epsilon ^{0}\right):-\frac{1}{\tau}\left(f_{i}^{(0)}-f_{i}^{e q}\right) = 0, \quad \mbox{ 即 }\ f_{i}^{(0)} = f_{i}^{e q}. \end{equation} $

(2.8) $ \begin{equation} O\left(\epsilon ^{1}\right):\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right) f_{i}^{(0)} = -\frac{1}{\tau \Delta t} f_{i}^{(1)}+h_{i}^{(1)}, \end{equation} $

(2.9) $ \begin{eqnarray} &&O\left(\epsilon ^{2}\right): \frac{\partial}{\partial t_{2}} f_{i}^{(0)}+\Delta t\left(\frac{1}{2}-\tau\right)\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{2} f_{i}^{(0)} \\&&+\Delta t \tau \left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right) h_{\alpha}^{(1)} = -\frac{1}{\tau \Delta t} f_{i}^{(2)}, \end{eqnarray} $

(2.10) $ \begin{eqnarray} &&O\left(\epsilon ^{3}\right): \frac{\partial}{\partial t_{3}} f_{i}^{(0)}+\Delta t(1-2 \tau) \frac{\partial}{\partial t_{2}}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right) f_{i}^{(0)} \\&& +\Delta t^{2} \left(\tau^{2}-\tau+\frac{1}{6}\right)\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{3} f_{i}^{(0)}\\&& +\Delta t^{2} \left(\frac{\tau}{2}-\tau^{2}\right)\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{2} h_{i}^{(1)} +\Delta t \tau \frac{\partial}{\partial t_{2}} h_{i}^{(1)} = -\frac{1}{\tau \Delta t} f_{i}^{(3)}. \end{eqnarray} $

与一般的LBM相似, 把宏观物理量$ u $定义为该点的局部粒子分布函数的和, 即

(2.11) $ \begin{equation} u(x, t) = \sum\limits_{i} f_{i}(x, t), \end{equation} $

由质量守恒定律, $ f_{i} $和$ f_{i}^{eq} $满足

(2.12) $ \begin{equation} u(x, t) = \sum\limits_{i} f_{i}(x, t) = \sum\limits_{i} f_{i}^{(e q)}(x, t), \end{equation} $

由(2.7)式, 有

(2.13) $ \begin{equation} \sum\limits_{i} f_{i}^{(0)} = u, \quad \sum\limits_{i} f_{i}^{(n)} = 0, \quad n = 1, 2, \cdots. \end{equation} $

由(2.6)式, 有

(2.14) $ \begin{equation} \sum\limits_{i} h_{i} = \epsilon \sum\limits_{i} h_{i}^{(1)}, \end{equation} $

(2.15) $ \begin{equation} \sum\limits_{i} c_{i} h_{i} = \epsilon \sum\limits_{i} c_{i} h_{i}^{(1)}, \end{equation} $

(2.16) $ \begin{equation} \sum\limits_{i} c_{i}^{2} h_{i} = \epsilon \sum\limits_{i} c_{i}^{2}h_{i}^{(1)}. \end{equation} $

为了恢复方程(1.3), 对平衡分布函数和修正函数做如下约束

(2.17) $ \begin{equation} \sum\limits_{i} c_{i} f_{i}^{(0)} = 0, \end{equation} $

(2.18) $ \begin{equation} \sum\limits_{i} c_{i}^{2} f_{i}^{(0)} = \lambda_{1}n\left (t \right )u, \end{equation} $

(2.19) $ \begin{equation} \sum\limits_{i} c_{i}^{3} f_{i}^{(0)} = \lambda_{2}r\left (t \right )u. \end{equation} $

(2.20) $ \begin{equation} \sum\limits_{i} h_{i} = 0, \end{equation} $

(2.21) $ \begin{equation} \sum\limits_{i} c_{i} h_{i} = \lambda _{3} l (t) u^{2}+\lambda _{4} m (t)u^{p+1}, \end{equation} $

(2.22) $ \begin{equation} \sum\limits_{i} c_{i}^{2} h_{i} = 0. \end{equation} $

对方程(2.7)两边关于$ i $求和, 有

(2.23) $ \begin{equation} \frac{\partial u}{\partial t_{1}}+\frac{\partial}{\partial x_{1}}\left(\sum\limits_{i} c_{i} f_{i}^{(0)}\right) = \sum\limits_{i} h_{i}^{(1)} = 0, \ \mbox{ 即 } \quad \frac{\partial u}{\partial t_{1}} = 0. \end{equation} $

对方程(2.8)两边关于$ i $求和, 有

(2.24) $ \begin{equation} \frac{\partial u}{\partial t_{2}}+\Delta t\left(\frac{1}{2}-\tau\right) \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{2} f_{i}^{(0)}+\Delta t \tau \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right) h_{i}^{(1)} = 0, \end{equation} $

其中

(2.25) $ \begin{eqnarray} \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{2} f_{i}^{(0)} & = &\frac{\partial^{2}}{\partial t_{1}^{2}}\left(\sum\limits_{i} f_{i}^{(0)}\right)+2 \frac{\partial^{2}}{\partial t_{1} \partial x_{1}}\left(\sum\limits_{i} c_{i} f_{i}^{(0)}\right)+\frac{\partial^{2}}{\partial x_{1}^{2}}\left(\sum\limits_{i} c_{i}^{2} f_{i}^{(0)}\right) \\& = &\frac{\partial^{2}}{\partial x_{1}^{2}}\left(\sum\limits_{i} c_{i}^{2} f_{i}^{(0)}\right) = \lambda_{1} n(t)\frac{\partial ^{2} u}{\partial x_{1}^{2}}, \end{eqnarray} $

(2.26) $ \begin{eqnarray} \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right) h_{i}^{(1)} & = &\frac{\partial}{\partial t_{1}}\left(\sum\limits_{i} h_{i}^{(1)}\right)+\frac{\partial}{\partial x_{1}}\left(\sum\limits_{i} c_{i} h_{i}^{(1)}\right) \\ & = &\frac{\partial}{\partial x_{1}}\left[\lambda _{3} l (t) u^{2}+\lambda _{4}m (t)u^{p+1} \right] \\ & = &2 \lambda _{3} l (t)u\frac{\partial u}{\partial x_{1} } +(p+1)\lambda _{4}m(t) u^{p} \frac{\partial u}{\partial x_{1} }, \end{eqnarray} $

把(2.25)式, (2.26)式代入(2.24)式中

(2.27) $ \begin{equation} \frac{\partial u}{\partial t_{2}}+\frac{\tau }{\epsilon } [2 \lambda _{3}\Delta t l(t) u \frac{\partial u}{\partial x_{1} }+(p+1)\lambda _{4}\Delta t n (t)u^{p}\frac{\partial u}{\partial x_{1} } ] +\Delta t(\frac{1}{2} -\tau )\lambda _{1} n(t)\frac{\partial^{2} u}{\partial x_{1}^{2} } = 0. \end{equation} $

对方程(2.9)两边关于$ i $求和, 有

(2.28) $ \begin{eqnarray} &&\frac{\partial u}{\partial t_{3}}+\Delta t^{2}\left(\tau^{2}-\tau+\frac{1}{6}\right) \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{3} f_{i}^{(0)}\\&& +\Delta t^{2}\left(\frac{\tau}{2}-\tau^{2}\right) \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{2} h_{i}^{(1)} = 0, \end{eqnarray} $

其中

(2.29) $ \begin{eqnarray} \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{3} f_{i}^{(0)} & = &\frac{\partial^{3}}{\partial t_{1}^{3}}\left(\sum\limits_{i} f_{i}^{(0)}\right)+3 \frac{\partial^{3}}{\partial x_{1} \partial t_{1}^{2}}\left(\sum\limits_{i} c_{i} f_{i}^{(0)}\right) \\&&+3 \frac{\partial^{3}}{\partial x_{1}^{2} \partial t_{1}}\left(\sum\limits_{i} c_{i}^{2} f_{i}^{(0)}\right)+\frac{\partial^{3}}{\partial x_{1}^{3}}\left(\sum\limits_{i} c_{i}^{3} f_{i}^{(0)}\right) \\& = &3\lambda _{1} \frac{\partial^{3} }{\partial x_{1}^{2} \partial t_{1} }[n(t) u] +\lambda _{2} r(t)\frac{\partial^{3} u}{\partial x_{1}^{3} }, \end{eqnarray} $

(2.30) $ \begin{eqnarray} \sum\limits_{i}\left(\frac{\partial}{\partial t_{1}}+c_{i} \frac{\partial}{\partial x_{1}}\right)^{2} h_{i}^{(1)}& = & \frac{\partial^{2}}{\partial t_{1}^{2}}\left(\sum\limits_{i} h_{i}^{(1)}\right)+2 \frac{\partial^{2}}{\partial x_{1} \partial t_{1}}\left(\sum\limits_{i} c_{i} h_{i}^{(2)}\right) +\frac{\partial^{2}}{\partial x_{1}^{2}}\left(\sum\limits_{i} c_{i}^{2} h_{i}^{(1)}\right) \\& = &\frac{2}{\epsilon }\frac{\partial^{2} }{\partial x_{1} \partial t_{1}}[\lambda _{3} l (t)u^{2} +\lambda _{4} m (t)u^{p+1} ], \end{eqnarray} $

把(2.29)式, (2.30)式代入(2.28)式, 得到

(2.31) $ \begin{eqnarray} && \frac{\partial u}{\partial t_{3}}+\Delta t^{2}(\tau^{2}-\tau+\frac{1}{6})[3\lambda _{1} n(t)\frac{\partial^{3} u}{\partial x_{1}^{2} \partial t_{1} } +\lambda _{2} r(t)\frac{\partial^{3} u}{\partial x_{1}^{3} }]\\&&+\frac{1}{\epsilon }\Delta t^{2}(\tau-2\tau^{2} ) \frac{\partial^{2} }{\partial x_{1} \partial t_{1}}[\lambda _{3} l (t)u^{2} +\lambda _{4} m (t)u^{p+1}] = 0. \end{eqnarray} $

采用LBM文献中所使用的符号, $ DnQm $格子Boltzmann离散速度模型是指在$ n $维空间中具有$ m $个离散速度的模型. 本文采用$ D1Q5 $格子Boltzmann离散速度模型, 其离散速度的集合为：$ c_{i} = \left \{ 0, \pm 1, \pm 2 \right \} $. 由方程$ (2.23)\times \epsilon +(2.27)\times \epsilon^{2} +(2.31)\times \epsilon^{3} $可以导出具有四阶截断误差的宏观方程

(2.32) $ \begin{eqnarray} &&u_{t}+2 \lambda _{3}\tau \Delta t l(t) u u_{x} +(p+1)\lambda _{4}\tau \Delta t m (t)u^{p}u_{x}+\Delta t(\frac{1}{2} -\tau )\lambda _{1} n(t)u_{xx}\\&&+\Delta t^{2}(\tau^{2}-\tau+\frac{1}{6})[3\lambda _{1}\epsilon u_{xx} \partial t_{1} n(t) +\lambda _{2} r(t)u_{xxx} ]\\& &+\epsilon \Delta t^{2}(\tau-2\tau^{2} ) [2\lambda _{3}uu_{x} \partial t_{1} l (t) +(p+1)\lambda _{4} u^{p}u_{x} \partial t_{1} m (t)]+O(\epsilon ^{4} ) = 0, \end{eqnarray} $

为恢复宏观方程(1.3), 令

(2.33) $ \begin{equation} \left\{\begin{array}{ll} { }\lambda _{1} = \frac{1}{\Delta t (\tau-0.5)}, \\ { }\lambda _{2} = \frac{1}{\Delta t^{2} (\tau^{2} -\tau+\frac{1}{6} )}, \\ { }\lambda _{3} = \frac{1}{2\Delta t \tau }, \\{ }\lambda _{4} = \frac{1}{(p+1)\Delta t \tau }. \end{array} \right. \end{equation} $

此时方程(2.32)为

(2.34) $ \begin{equation} u_{t} + l(t) u u_{x} + m (t)u^{p}u_{x}- n(t)u_{xx}+r(t)u_{xxx} = {\cal R}, \end{equation} $

$ {\cal R} = -\epsilon ^{2} \left \{ \frac{6\tau^{2}-6\tau +1 }{2\tau-1}u_{xx}\partial _{t_{1} }n(t)+(1-2\tau)[uu_{x}\partial _{t_{1} }l (t) +u^{p} u_{x} \partial _{t_{1} }m (t) ] \right \} = O(\epsilon ^{2} ). $

为了推导出局部平衡态分布函数的具体表达式, 我们应该引入一个新的约束条件, 考虑到$ f $的矩都是$ u $的多项式函数, 我们假设平衡分布函数满足以下约束

(2.35) $ \begin{equation} f_{4}^{\left ( 0 \right ) } = \sigma u, \end{equation} $

其中引入参数$ \sigma $来调节平衡态分布函数

(2.36) $ \begin{eqnarray} f_{i}^{(0)} = \left\{\begin{array}{lll}{ } {} [1+6\sigma -\frac{n(t)}{c^{2}\Delta t\left ( \tau -0.5 \right ) } +\frac{r(t)}{2c^{3}\Delta t^{2} \left (\tau ^{2} -\tau +\frac{1}{6} \right ) } ]u, \ \ &i = 0, \\{ } {} [-4\sigma +\frac{n(t)}{2 c^{2}\Delta t\left ( \tau -0.5 \right ) } -\frac{r(t)}{2 c^{3}\Delta t^{2} \left (\tau ^{2} -\tau +\frac{1}{6} \right )} ]u, \ \ &i = 1, \\ { } {}[-4\sigma +\frac{n(t)}{2 c^{2}\Delta t\left ( \tau -0.5 \right ) } -\frac{r(t)}{6 c^{3}\Delta t^{2} \left (\tau ^{2} -\tau +\frac{1}{6} \right )} ]u, \ \ &i = 2, \\{ }{} [\sigma +\frac{r(t)}{6 c^{3}\Delta t^{2} \left (\tau ^{2} -\tau +\frac{1}{6} \right )} ]u, \ \ &i = 3, \\ \sigma u, \ \ &i = 4.\end{array}\right. \end{eqnarray} $

耦合方程(2.14)–(2.16). 假设速度1和速度2方向上的修正函数相同, 我们可以推导出补偿函数为

(2.37) $ \begin{equation} h_{i } = \left\{\begin{array}{ll} { }\frac{1}{2} {\cal H}, &i = 0, \\{ } -\frac{1}{3} {\cal H}, &i = 1, \\{ } -\frac{1}{3} {\cal H}, &i = 2, \\{ } \frac{1}{3} {\cal H}, &i = 3, \\{ }-\frac{1}{6} {\cal H}, &i = 4.\end{array}\right. \end{equation} $

(2.38) $ \begin{equation} {\cal H} = [\frac{l (t)}{2} u^{2} +\frac{m(t)}{p+1} u^{p+1}]/\left (\Delta t \tau c \right ). \end{equation} $

考虑到方程的系数为时变系数. 因此在具体的数值模拟时, 松弛时间$ \tau $不应该是一个固定常数, 由方程(2.36), 假设$ \tau _{1} $满足

(2.39) $ \begin{equation} \tau _{1} = \tau ^{2} -\tau +\frac{1}{6}, \end{equation} $

由文献[11], 把参数$ k $设为一个给定的值

(2.40) $ \begin{equation} k = \frac{r(t)}{c^{3}\Delta t^{2} \tau _{1} }. \end{equation} $

为了确保本文建立的格子Boltzmann模型的稳定性, 因此必须保证松弛时间满足约束条件$ \tau _{1}>0.5 $ (参见文献[14]), 由此可得

(2.41) $ \begin{equation} \tau (t) = \frac{1}{2} +\sqrt{\frac{1}{12} +\frac{\Delta t}{\Delta x^{3}k } r(t) }, \end{equation} $

(2.42) $ \begin{equation} \frac{1}{12} +\frac{\Delta t}{\Delta x^{3}k } r(t)\ge 0\Rightarrow k\ge \frac{-12\Delta t r(t)}{\Delta x^{3} }. \end{equation} $

对参数取值总结如下.

(1) 参数$ k $是单松弛时间$ \tau $的调节因子, 一般取值为$ 0.0<k<1.5 $ (参见文献[11]);

(2) 常数$ \sigma $是平衡函数的一个调整因子, 它可以在演化过程中调整模拟结果, 一般$ \sigma $是一个负的微小值, $ -\frac{1}{12} <\sigma <-0.01 $ (参见文献[12]).

本文采用均方根误差$ E_{2} $, 最大误差$ E_{\infty } $, 全局相对误差$ GRE $来验证格子Boltzmann模型的有效性, 并采用最小二乘拟合来计算模型的精度.

(3.1) $ \begin{equation} E_{2} = \frac{1}{N} \sqrt{\sum\limits_{j = 1}^{N}\left [ u\left ( x_{j} , t \right )-u^{*} \left ( x_{j}, t \right ) \right ]^{2} }, \end{equation} $

(3.2) $ \begin{equation} E_{\infty } = \max\limits_{j = 0, 1, 2\dots, N }\left | u\left ( x_{j} , t \right )-u^{*} \left ( x_{j}, t \right ) \right |, \end{equation} $

(3.3) $ \begin{equation} GRE = \frac{ { \sum\limits_{j = 1}^{N}}\left | u\left ( x_{j} , t \right )-u^{*} \left ( x_{j}, t \right ) \right | }{ { \sum\limits_{j = 1}^{N} \left | u^{*} \left ( x_{j} , t \right ) \right | } }, \end{equation} $

其中, $ u\left ( x, t \right ) $和$ u^{\ast } \left ( x, t \right ) $分别表示数值解和解析解, $ N $是网格节点数. 并通过以下方法计算离散化方法收敛速度.

(3.4) $ \begin{equation} p \approx\left|\frac{\log E_{\mbox{new}}-\log E_{\mbox{old}}}{\log H_{\mbox{new}}-\log H_{\mbox{old}}}\right|, \end{equation} $

其中$ p $是收敛阶, $ E $代表误差, $ H $代表序列, 下标new和old分别代表新旧状态.

算例3.1   取$ l(t)\! = \!\alpha , m(t)\! = \!n(t)\! = \!0, r(t)\! = \!-\beta $时, 此时方程可以转化成KdV方程^[12]

(3.5) $ \begin{equation} u_{t} +\alpha uu_{x} = -\beta u_{xxx}. \end{equation} $

它被用来描述流体中的波传播和分层内波传播现象, 以及等离子体中的离子声波. 令$ \alpha = 6, $$ \beta = 1 $, KdV方程的双孤子解可以表示为

$ u(x, t) = \frac{\partial^{2}}{\partial x^{2}}\{2 \log [F(x, t)]\}, $

$ F(x, t) = 1+a_{1} e^{\kappa_{1} x-\kappa_{1}^{3} t}+a_{2} e^{\kappa_{2} x-\kappa_{2}^{3} t}+\frac{a_{1} a_{2}\left(\kappa_{1}-\kappa_{2}\right)^{2} e^{\left(\kappa_{1}+\kappa_{2}\right) x-\left(\kappa_{1}^{3}+\kappa_{2}^{3}\right) t}}{\left(\kappa_{1}+\kappa_{2}\right)^{2}}, $

其中, $ a_{1}, a_{2}, \kappa_{1}, \kappa_{2} $为常数, $ a_{1} = 2, a_{2} = 1, \kappa_{1} = 1.5, \kappa_{2} = 1 $, $ k = 0.11, \sigma = -0.0095 $, $ \Delta x = 0.05, \Delta t = 0.0001, c = \Delta x/\Delta t = 500 $, 计算区域固定在$ I = [-10, 20] $, 图 1给出了数值解和解析解在二维空间中的视觉对比图. 表 1给出了不同时刻数值结果的$ E_{2}, E_{\infty }, GRE $. 图 2给出了数值结果的时空演化图. 图 3左给出了$ \Delta x = \left \{0.01\quad0.02\quad0.03\quad0.04 \right \} $时, $ \Delta t $取0.001时, 不同时刻的空间精度图, $ t = 1 $时数值精度为1.9001, $ t = 2 $时数值精度为2.2003. 图 3右给出了$ \Delta t = \left \{ 0.0001\quad 0.0002\quad 0.0003\quad 0.0004 \right \} $, $ \Delta x $取0.01时, 不同位置的时间精度图, $ x = 0 $时斜率为2.1833, $ x = 5 $时斜率为2.0456. 这与方程(2.34)的精度是相同的.

注3.1   数值模拟时, 为区分数值解与解析解, 用matlab中的colormap default表示数值解, colormap jet表示精确解, 下同.

算例3.2   取$ m(t) = \alpha, n(t) = \beta \times t^{\frac{2}{n} -1} $, 方程转化为广义变系数Burgers方程$ B\left ( n, 1 \right ) $ (参见文献[15])

(3.6) $ \begin{equation} u_{t} +5 \alpha u^{p} u_{x} = \beta t^{\frac{2}{n}-1 } u_{xx}. \end{equation} $

初始条件和边界条件为

$ \left\{\begin{array}{ll}u(x, 0) = 0, &-\infty <x<\infty, \\ u(0, t) = q(t), &t>0, \\ { }\lim\limits_{x \to \infty} u(x, t) = 0, &t>0, \end{array}\right. $

其中, 取$ q(t) = \gamma \times t^{-\frac{1}{n} } $, 此时方程的解析解为

$ u(x, t) = \frac{1}{t^{\frac{1}{n}}}\left\{\frac{{\rm e}^{-\frac{x^{2}}{2 n \beta t^{\frac{2}{n}}}}}{\left[\gamma^{1-n}-\alpha \sqrt{\frac{\pi n(n-1)}{2 \beta}} {\rm erf}\left(\sqrt{\frac{n-1}{2 n \beta}} \frac{x}{t^{\frac{1}{n}}}\right)\right]^{\frac{1}{n-1}}}\right\}. $

上式中, 取$ n = 5, \alpha = 0.5, \beta = 2, \gamma = 0.5 $, $ k = 1.05, \sigma = -0.19, \Delta x = 0.01, \Delta t = 0.0001, $$ c = \Delta x/\Delta t = 100 $, 计算区域为$ I = [-10, 10] $. 图 4给出了数值结果的时空演化图, 图 5给出了数值解和解析解在二维空间中的视觉对比图, 从图中可以看出, 数值结果与解析解吻合良好, 图 6左给出了$ \Delta x = \left \{ 0.01 \quad 0.02 \quad 0.03\quad0.04 \right \} $, $ \Delta t = 0.001 $时, 不同时刻的空间精度图, $ t = 1 $时数值精度为1.9940, $ t = 2 $时数值精度为2.1730. 图 6右给出了$ \Delta t = \left \{ 0.0001\quad 0.0002\quad 0.0003\quad 0.0004 \right \} $, $ \Delta x = 0.01 $, 不同位置的时间精度图, $ x = 0 $时数值精度为2.0288, $ x = 10 $时数值精度为1.9942, 这与方程(2.34)的精度是相同的. 表 2给出了不同时刻数值结果的$ E_{2}, E_{\infty }, GRE $.

算例3.3   考虑如下形式的广义变系数KdV-mKdV方程^[15]

(3.7) $ \begin{equation} u_{t} -15\cos\left ( \pi t/2 \right ) u_{x} -6uu_{x} +0.06u^{p} u_{x} = 0.01u_{xxx}. \end{equation} $

$ p = 2 $时, 方程的解析解可以表示为

$ u(x, t) = -\frac{202}{\Phi} \exp \left[\frac{t}{100}+\frac{30}{\pi} \sin \left(\frac{\pi t}{2}\right)+x+\frac{1}{100}\right], $

$ \begin{eqnarray*} \Phi& = & 20, 200 \exp \left[\frac{t}{100}+\frac{30}{\pi} \sin \left(\frac{\pi t}{2}\right)+x+\frac{1}{100}\right]+10, 001 \exp \left[\frac{60}{\pi} \sin \left(\frac{\pi t}{2}\right)+2 x+\frac{1}{50}\right] \\&&+10, 201 \exp \left(\frac{t}{50}\right). \end{eqnarray*} $

上式中, 取$ k = 0.985, \sigma = -0.09, \Delta x = 0.01, \Delta t = 0.0005, c = \Delta x/\Delta t = 500 $, 计算区域固定在$ I = [-20, 10] $. 表 3给出了不同时刻数值结果的$ E_{2}, E_{\infty }, GRE $, 图 7给出了给出了数值结果的时空演化图, 图 8数值解和解析解在二维空间中的视觉对比图, 从图中可以看出, 数值结果与解析解吻合良好. 图 9左给出了$ \Delta x = \left \{ 0.01 \quad 0.02 \quad 0.03\quad0.04 \right \} $时, $ \Delta t $为0.001, 不同时刻的空间精度图, $ t = 1 $时数值精度为1.9532, $ t = 2 $时数值精度为1.9942, 图 9右给出了$ \Delta t = \left \{ 0.0001\quad 0.0002\quad 0.0003\quad 0.0004 \right \} $时, $ \Delta x $取0.01, 不同位置的时间精度图, $ x = 0 $时数值精度为2.0456, $ x = 5 $时数值精度为2.0503, 这与方程(2.34)的数值精度是相同的, 时间和空间均可达到2阶精度.

本文研究了一类变系数广义KdV-Burgers方程的格子Boltzmann模型. 首先, 通过选择平衡分布函数和适当地增加修正函数建立了方程(1.3). 利用D1Q5格子Boltzmann模型对不同形式的变系数KdV-Burgers方程进行数值模拟试验.其次, 证明了格子Boltzmann模型与宏观方程是相容的. 然后，引入自由参数$ \sigma $便于构造出平衡分布函数$ f_{i}^{eq} $, 引入参数$ \mathit{k} $调整松弛时间, 很好的避免了松弛时间$ \tau $的取值问题. 最后, 对数值算例的时空变化趋势和时间空间精度进行分析, 证明了模型可以从时间上和空间上达到二阶精度. 该模型推广到更高维的变系数Burgers方程, 变系数对流扩散方程, 变系数模糊类波动方程数值求解中.