求解交通流LWR模型的高分辨率熵相容格式

High-Resolution Entropy Consistent Scheme for the LWR Model of Traffic Flow

摘要: 针对交通流LWR模型，构造了求解交通流LWR模型的熵相容格式，并将一种基于MUSCL型重构方法的新型斜率限制器应用于该格式中，得到了求解交通流LWR模型的高分辨率熵相容格式。将新构造的格式应用于多个实际交通流问题的求解中，数值结果表明，该格式对激波有良好的捕捉效果，没有非物理振荡，且在稀疏波区域能够平滑地过渡。

Abstract: In this paper, a high-resolution entropy consistent scheme of the LWR model in traffic flow is pro-posed. Based on MUSCL-type reconstruction method, a new kind of slope limiter is constructed. The high-resolution entropy consistent scheme for solving LWR model can be obtained by adding a slope limiter to the entropy consistent scheme. Applying the newly constructed scheme to the solu-tion of multiple practical traffic flow problems, the numerical results show that the scheme can capture shock waves well, and there are no non-physical oscillations in the discontinuous region and smooth transition in the rarefaction region.

1. 引言

连续交通流模型在交通模拟、交通控制等许多应用中具有重要的现实意义。由Lighthill和Whitham [1] 和Richards [2] 分别独立提出的LWR模型是所有其他交通流模型的先驱。该模型是一个宏观模型，描述了均匀单向公路交通流量的总体特征。虽然LWR模型只是提供了交通流的粗略表示，但它能够定性地再现大量的真实交通现象，例如交通堵塞、车流疏散、信号灯放行等。

LWR模型背后有两个主要关系，一个是连续性方程，表明道路上的车辆数量是守恒的，另一个是流量和流量密度之间的基本关系，LWR模型将这两种关系相结合，得到了一个反映车流密度的时间变化和流量的空间变化的偏微分方程。该方程在数学上属于非线性的双曲守恒律方程，因此，在模拟LWR模型方程时，可以采用非线性双曲守恒律方程的数值计算方法，如有限差分法、有限体积法和有限元法等。非线性双曲守恒律方程的一个重要的特征是：即使初始条件连续，其解也可能会出现间断的情况。为了解决间断问题，1954年，Lax提出了弱解 [3] 的概念，这一概念允许间断解的存在，但弱解并不唯一，必须增加限制条件，才能够求出与物理相关的唯一弱解。在数学上，满足“粘性消失”的Cauchy问题具有唯一、有物理意义的解，这个解叫做熵解。“熵”是求解双曲守恒律方程的关键，只要所选择的数值格式能够保持总熵不变，该格式被称为熵守恒格式。然而，在存在间断解的情况下，熵守恒格式缺少熵的耗散机制，使得它的数值解具有较大不稳定性。所以Tadmor研究指出 [4]：一个三点格式只需含有比熵守恒格式更多的粘性则是熵稳定的。2006年，Roe在熵守恒格式的基础上加入Roe格式的数值粘性项，得到了熵稳定的ERoe格式 [5]。随后在2009年，Roe和Ismail提出了熵相容格式 [6]，该格式对数值粘性项的量化更精确。2015年，刘友琼、封建湖、任炯 [7] 提出了一种复合型通量限制器，并在此基础上给出了一种高分辨率熵相容格式。2021年，任璇、封建湖基于MUSCL型重构方法，提出了一种新的斜率限制器 [8]，有效地提高了熵相容格式的精度和分辨率。随后沈亚玲等 [9] 将该斜率限制器应用到理想磁流体动力学(MHD)方程的数值求解中，并通过数值算例验证了该格式是一种较为理想的数值方法。

关于LWR模型的求解，郑素佩等 [10] 提出了求解LWR模型的高分辨率熵相容格式，并在空间方向采用三阶CWENO重构，避免了膨胀激波现象的产生。冯娟娟等 [11] 采用WENO-Z+重构的熵稳定格式对LWR模型进行了数值模拟，结果表明该算法对激波有良好的捕捉效果。本文为了给出一种求解交通流LWR模型的更好的数值格式，在熵相容格式的基础上，加入了MUSCL型重构方法的新型斜率限制器，进一步提高了格式的计算精度。本文数值算例由熵相容格式、四阶CWENO熵相容格式以及加入新型斜率限制器的熵相容格式计算得到，模拟了LWR交通流模型随着时间的推移距离和密度之间的变化关系，证明了新算法的有效性。

2. 交通流LWR模型

LWR模型的质量守恒方程及速度表达式如下 [12]

$\frac{\partial k}{\partial t}+\frac{\partial f\left(k\right)}{\partial x}=0,f=ku,u=u_f\left(1-\frac{k}{k_m}\right).$

k为密度；t为时间；x为距离； $u_f$ 为最大速度； $k_m$ 为阻塞密度。显然，该模型方程属于非线性双曲守恒律方程。

2.1. 方程离散

空间上采用均匀网格上的半离散格式

$\frac{\text{d}{k}_j\left(t\right)}{\text{d}t}=-\frac{1}{\Delta x}\left(f_{j+1/2}-f_{j-1/2}\right),$ (1)

其中 $k_j\left(t\right)=\frac{1}{\Delta x}{\displaystyle {\int }_{x_{j-1/2}}^{x_{j+1/2}}k\left(x,t\right)\text{d}x}$。

时间方向上采用三阶Runge-Kutta方法 [13]：

$\{\begin{array}{l}{k}_j^{\ast }=k_j^k+\Delta tL\left(k_j^k\right),\hfill \\ k_j^{\ast \ast }=\frac{3}{4}{k}_j^k+\frac{1}{4}{k}_j^{\ast }+\frac{\Delta t}{4}L\left(k_j^{\ast }\right),\hfill \\ k_j^{k+1}=\frac{1}{3}{k}_j^k+\frac{2}{3}{k}_j^{\ast \ast }+\frac{2\Delta t}{3}L\left(k_j^{\ast \ast }\right).\hfill \end{array}$ (2)

其中 $L\left(k\right)=-\frac{1}{\Delta x}\left(f_{j+1/2}-f_{j-1/2}\right)$。

2.2. 熵相容格式

对于交通流LWR方程的熵守恒格式，要求满足离散熵等式同时保持总熵不变：

$\frac{\text{d}}{\text{d}t}E\left(k_j\left(t\right)\right)+\frac{1}{\Delta x}\left(F_{j+1/2}-F_{j-1/2}\right)=0,$

其中，取熵函数 $E\left(k\right)=k^2$，计算出熵通量函数 $F=u_f\left(k^2-\frac{4k^3}{3k_m}\right)$，熵变量 $v=E\prime \left(k\right)=2k$ 和熵势

$\psi \left(\nu \left(k\right)\right)=\nu \left(k\right)f\left(k\right)-F\left(k\right)=u_f\left(k^2-\frac{2k^3}{3k_m}\right).$

则熵守恒格式的数值通量为 [4]

$f_{j+1/2}^C=\frac{\psi \left({\nu }_{j+1}\right)-\psi \left({\nu }_j\right)}{{\nu }_{j+1}-{\nu }_j}=\frac{u_f\left(k_{j+1}^2-k_j^2\right)-\frac{2u_f}{3k_m}\left(k_{j+1}^3-k_j^3\right)}{2\left(k_{j+1}-k_j\right)}.$ (3)

然而，熵守恒格式在解的间断处会表现出严重的不稳定性，而要使格式达到熵稳定，总熵必须有所耗散，Ismail [6] 提出在熵守恒格式的数值通量中加入一个二阶的迎风项，得到了如下熵稳定格式的数值通量

$f_{j+1/2}^{ES}=f_{j+1/2}^C-\frac{1}{2}\left|\bar{a}\right|\left[k\right]=f_{j+1/2}^C-\frac{1}{2}{u}_f\left|1-\frac{k_{j+1}+k_j}{k_m}\right|\left(k_{j+1}-k_j\right),$ (4)

其中 $\left|\bar{a}\right|=\frac{a_{j+1}+a_j}{2}$， $a_j=f\prime \left(k_j\right)$， $\left[k\right]=k_{j+1}-k_j$。

熵稳定格式包含的数值粘性大于熵守恒格式，但在实际的数值实验中，仍不能避免伪振荡的产生。Ismail等在文献 [6] 中给出了一维双曲守恒律方程的熵相容通量的通用表达式，并将其推广到一维Euler方程组。本文将Ismail等人的熵相容通量的通用表达式推广到交通流LWR模型，并给出了求解交通流LWR模型的熵相容格式的数值通量

$\begin{array}{c}{f}_{j+1/2}^{EC}=f_{j+1/2}^C-\frac{1}{2}\left(\left|\bar{a}\right|+\alpha \left|\left[a\right]\right|\right)\left[k\right]\\ =\frac{u_f\left(k_{j+1}^2-k_j^2\right)-\frac{2u_f}{3k_m}\left(k_{j+1}^3-k_j^3\right)}{2\left(k_{j+1}-k_j\right)}-\frac{1}{2}{u}_f\left|1-\frac{k_{j+1}+k_j}{k_m}\right|\left(k_{j+1}-k_j\right)\\ -\frac{1}{6}{u}_f\left|\frac{k_j-k_{j+1}}{k_m}\right|\left(k_{j+1}-k_j\right),\end{array}$ (5)

取 $\alpha =\frac{1}{6}$， $\left[k\right]=k_{j+1}-k_j$。

2.3. 高分辨率熵相容格式

2.3.1. 斜率限制器

在对LWR模型方程进行数值模拟时，为了使熵相容格式的精度和分辨率得到进一步的提高，重构各控制单元 $I_j=\left[x_{j-1/2},x_{j+1/2}\right]$ 上的参考解 $k_j\left(x\right)$，并将重构后的函数表示为 $R_j\left(x\right)$，若满足以下条件

$\left|k_j\left(x\right)-R_j\left(x\right)\right|\le O\left({\left(\Delta x\right)}^{k+1}\right),$

则称该重构函数 $R_j\left(x\right)$ 具有k阶精度。假设 $k_j^n$ 表示单元 $I_j$ 在时间层 $t^n$ 上的平均值，即

$k_j^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{R}_j}\left(x\right)\text{d}x=\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{k}_j}\left(x\right)\text{d}x.$

在构造中，首先将 $R_j\left(x\right)$ 在 $x_j$ 处泰勒展开

$R_j\left(x\right)=k_j\left(x_j\right)+{\frac{\partial k_j}{\partial x}|}_{x_j}\left(x-x_j\right)+{\frac{{\partial }^2{k}_j}{\partial x^2}|}_{x_j}\frac{{\left(x-x_j\right)}^2}{2}+\cdots ,$ (6)

之后将 $k_j\left(x\right)$ 在点 $x_j$ 处泰勒展开，可以得到

$k_j^n=\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{k}_j}\left(x\right)\text{d}x=k_j\left(x_j\right)+\frac{1}{\Delta x}{\frac{{\partial }^2{k}_j}{\partial x^2}|}_{x_j}{\displaystyle {\int }_{I_j}\frac{{\left(x-x_j\right)}^2}{2}\text{d}x}+\cdots ,$ (7)

合并(6)和(7)，得到

$R_j\left(x\right)=k_j^n+{\frac{\partial k_j}{\partial x}|}_{x_j}\left(x-x_j\right)+\frac{1}{2}{\frac{{\partial }^2{k}_j}{\partial x^2}|}_{x_j}\left({\left(x-x_j\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{\left(x-x_j\right)}^2\text{d}x}\right)+\cdots .$

我们取

$R_j\left(x\right)=k_j^n+{\frac{\partial k_j}{\partial x}|}_{x_j}\left(x-x_j\right)+\frac{1}{2}{\frac{{\partial }^2{k}_j}{\partial x^2}|}_{x_j}\left({\left(x-x_j\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{\left(x-x_j\right)}^2\text{d}x}\right).$

上式满足 $\left|k_j\left(x\right)-R_j\left(x\right)\right|\le O\left({\left(\Delta x\right)}^3\right)$，也就是在重构之后，得到了二阶精度，为了方便，将 $R_j\left(x\right)$ 写成

$R_j\left(x\right)=k_j^n+H\left(x-x_j\right),$ (8)

其中

$H\left(x-x_j\right)={\frac{\partial k_j}{\partial x}|}_{x_j}\left(x-x_j\right)+\frac{1}{2}{\frac{{\partial }^2{k}_j}{\partial x^2}|}_{x_j}\left({\left(x-x_j\right)}^2-\frac{1}{\Delta x}{\displaystyle {\int }_{I_j}{\left(x-x_j\right)}^2\text{d}x}\right).$

对上式中的一阶和二阶偏导数，用中心差商来近似，即

${\frac{\partial k_j}{\partial x}|}_{x_j}=\frac{k_{j+1}^n-k_{j-1}^n}{2\Delta x},{\frac{{\partial }^2{k}_j}{\partial x^2}|}_{x_j}=\frac{k_{j+1}^n-2k_j^n+k_{j-1}^n}{\Delta x^2}.$

进一步，将重构(8)改进为

$R_j\left(x\right)=k_j^n+{\phi }_{j}H\left(x-x_j\right),$ (9)

其中， ${\phi }_j\in \left[0,1\right]$ 是对类似于MUSCL方法 [14] 中斜率 ${\Delta }_j/\Delta x$ 的进一步限制，称之为斜率限制器。重构(9)在控制单元 $I_j=\left[x_{j-1/2},x_{j+1/2}\right]$ 的边界外推值分别为

$R_j^L=k_j^n+{\phi }_{j}H\left(x_{j-1/2}-x_j\right),R_j^R=k_j^n+{\phi }_{j}H\left(x_{j+1/2}-x_j\right).$

限制 $R_j^L$ 和 $R_j^R$

$\{\begin{array}{l}\min \left(k_{j-1}^n,k_j^n\right)\le R_j^L\le \max \left(k_{j-1}^n,k_j^n\right),\\ \min \left(k_{j+1}^n,k_j^n\right)\le R_j^R\le \max \left(k_{j+1}^n,k_j^n\right),\end{array}$ (10)

来确保不会产生新的极值。在上述不等式的各边都减去 $k_j^n$ 得到

$\{\begin{array}{l}\min \left(k_{j-1}^n-k_j^n,0\right)\le R_j^L-k_j^n\le \max \left(k_{j-1}^n-k_j^n,0\right),\\ \min \left(k_{j+1}^n-k_j^n,0\right)\le R_j^R-k_j^n\le \max \left(k_{j+1}^n-k_j^n,0\right).\end{array}$ (11)

然后，利用公式(11)中的第一个不等式，对左斜率限制器的构造过程进行分析。

当 $k_{j-1}^n>k_j^n$ 时

$0=\min \left(k_{j-1}^n-k_j^n,0\right)\le R_j^L-k_j^n={\phi }_{j}H\left(x_{j-1/2}-x_j\right)\le k_{j-1}^n-k_j^n,$

也就是

${\phi }_j\le \frac{k_{j-1}^n-k_j^n}{H\left(x_{j-1/2}-x_j\right)},$ (12)

由于 ${\phi }_j\in \left[0,1\right]$，所以取

${\phi }_L=\min \left(1,\frac{k_{j-1}^n-k_j^n}{H\left(x_{j-1/2}-x_j\right)}\right);$

同理，右斜率限制器取

${\phi }_R=\min \left(1,\frac{k_{j+1}^n-k_j^n}{H\left(x_{j+1/2}-x_j\right)}\right).$

为了确保不等式(12)的成立，斜率限制器取 ${\phi }_j=\max \left(\min \left({\phi }_L,{\phi }_R\right),0\right)$ [8]，对于 ${\phi }_L$ 和 ${\phi }_R$，一般情况下，在计算的过程中，为了避免出现计算错误，通常在分母上加上一个非常小的常数 $\epsilon ={10}^{-6}$。

2.3.2. 高分辨率熵相容格式

在本节中，基于2.2节所提出的熵相容格式，将2.3.1节构造的新的斜率限制器应用到该格式中，得到了求解交通流LWR模型的高分辨率熵相容格式(EC-SL)。

交通流LWR模型的高分辨率熵相容格式的数值通量为

$f_{j+1/2}^{EC-SL}=f_{j+1/2}^C-\frac{1}{2}{u}_f\left|1-\frac{k_{j+1}^{SL}+k_j^{SL}}{k_m}\right|\left(k_{j+1}^{SL}-k_j^{SL}\right)-\frac{1}{6}{u}_f\left|\frac{k_j^{SL}-k_{j+1}^{SL}}{k_m}\right|\left(k_{j+1}^{SL}-k_j^{SL}\right),$ (13)

$k_j^{SL}$ 和 $k_{j+1}^{SL}$ 是重构后得到的变量值。重构过程如下，利用加入斜率限制器后的式子 $R_j\left(x\right)=k_j^n+{\phi }_{j}H\left(x-x_j\right)$，计算出 $R_j\left(x\right)$ 和 $R_{j+1}\left(x\right)$，并令 $k_j^{SL}$ 和 $k_{j+1}^{SL}$ 分别等于 $R_j\left(x_{j-1/2}\right)$ 和 $R_{j+1}\left(x_{j+1/2}\right)$，在数据重构完成之后，将 $k_j^{SL}$ 和 $k_{j+1}^{SL}$ 代入数值通量中，就能得到求解交通流LWR模型的高分辨率熵相容格式。

3. 数值算例

本节中数值算例的参考解是由四阶CWENO熵相容格式在2000个网格单元计算所得。EC和EC-SL分别表示由熵相容格式的数值通量(公式(5))和高分辨率熵相容格式的数值通量(公式(13))得到的半离散守恒型差分格式(公式(1))计算所得，时间方向的离散采用三阶Runge-Kutta格式(公式(2))。数值模拟时，取100个网格单元，横坐标为距离的无量纲化，纵坐标为车流密度的无量纲化。边界条件均采用Neumann边界条件。

算例1 信号灯由红变绿的极端车流情形。长为 $b=10\text{km}$ 的公路， $u_f=150\text{m}/\text{s}$， $k_m=0.15\text{veh}/\text{m}$，模拟时间t设为120 s、300 s、450 s、600 s。

初始条件为：

$k_j^0=\{\begin{array}{l}0.8k_{m}0\le x\le 7000,\text{}\text{}\\ 0.4k_{m}7000\le x\le 10000.\end{array}\text{}$

边界条件是 $k_0^n=0.8k_m$， $k_m^n=0.4k_m$。

计算结果如图1所示。

t = 120 s t = 300 s t = 450 s t = 600 s

Figure 1. Results of example 1

图 1. 算例1数值结果

图1分别显示了在120 s、300 s、450 s、600 s时车流密度与距离之间的关系。由于信号灯由红灯变为绿灯，上游高密度的交通流逐渐向下游疏散，从而形成了稀疏波。如图所示，随着车辆行驶距离的增加，稀疏波的过渡区域逐渐扩大。通过与EC格式的对比，EC-SL格式的结果在稀疏波区域更加接近于参考解。在表1中，通过两种格式精度的比较，说明了斜率限制器的确提高了熵相容格式的精度。

Table 1. The EC/EC-Sl scheme computers the order of error and convergence up to t = 120 s

表1. EC/EC-SL格式计算到t = 120 s的误差和收敛阶

由表1可知，EC-SL格式在算例1中可以达到二阶精度，而EC格式只有一阶精度。

算例2 路长为 $b=1.1\text{km}$， $u_f=16.67\text{m}/\text{s}$， $k_m=0.168\text{veh}/\text{m}$，初始车流密度为0.075 veh/m，在 $x=1000\text{m}$ 处有一个红色信号灯。模拟时间t设为30 s、60 s、90 s、120 s。

初始条件为：

$k_j^0=\{\begin{array}{l}0.0750\le x\le 1000,\text{}\text{}\\ k_{m}1000\le x\le 1100.\end{array}\text{}$

边界条件是 $k_0^n=0.075$， $k_m^n=0.168$。

计算结果如图2所示。

图2分别显示了在30 s、60 s、90 s、120 s时车流密度与距离之间的关系。由于该路段在 $x=1000\text{m}$ 处是红色信号灯，行驶中的车流不断被压缩，造成了下游的堵塞，使得交通流向上游传播形成了激波。由图可以看出，EC格式和EC-SL格式在间断处均出现抹平现象，但EC-SL格式较EC格式有了良好的改善。EC-SL格式在激波处过渡带较窄，跨越间断的单元数目也少于EC格式，并且EC-SL格式解的分辨率比EC格式的高。

t = 30 s t = 60 s t = 90 s t = 120 s

Figure 2. Results of example 2

图2. 算例2数值结果

算例3 交通信号灯放行后车流密度的变化。路长为 $b=2.5\text{km}$， $u_f=14\text{m}/\text{s}$， $k_m=1.0\text{veh}/\text{m}$。

初始条件为：

$k_j^0=\{\begin{array}{l}0.5500\text{m}\le x\le 1000\text{m},\text{}\text{}\\ 01000\text{m}\le x\le 1500\text{m},\\ 0.25\text{otherwise}.\end{array}\text{}$

这个算例反映的交通流状况是：在500米位置有一个红色信号灯，信号灯左侧车流密度为0.25；500米~1000米路段车流密度为0.5；1000米处有一个绿色信号灯，信号灯前方道路上1000米~1500米范围内没有车辆；1500米位置右侧，车流密度为0.25，1500米处是一个向右的激波。模拟时间t设为15 s、30 s、50 s、60 s，计算结果如图3所示。

图3分别显示了在15 s、30 s、50 s、60 s时车流密度与距离之间的关系。根据初始条件可知，在500米处，信号灯刚好由红灯变成绿灯，大量的车辆从信号灯左边驶向右边，使得信号灯右边车流密度迅速增加，从而产生了激波，并且密度在短时间内保持不变；又因为绿灯放行的原因，被放行的车辆在1000米后快速向前行驶，使车流密度逐渐降低到0，导致出现了稀疏波现象；在绿灯放行后，在第一辆放行的车辆还没有追上前面最后一辆运行的车辆前，道路上仍然有一部分路段没有车辆；在零密度与没有被红灯阻断的最后一辆车之间，则会形成激波间断。从图中看出，EC-SL格式在激波间断处和稀疏波处算的结果都优于EC格式，即EC-SL格式更贴近参考解。

t = 15 s t = 30 s t = 50 s t = 60 s

Figure 3. Results of example 3

图3. 算例3数值结果

4. 结束语

由于交通流LWR模型在数学上可以表示为双曲守恒律方程，所以我们将求解一般双曲守恒律方程的熵相容格式进一步推广到交通流LWR模型，并在重构时将MUSCL型斜率限制器加入熵相容格式，得到了求解交通流LWR模型的高分辨率熵相容格式。计算结果表明，新格式在稀疏波和激波处的逼近效果优于其他格式，更好地模拟了LWR交通流模型随着时间的推移距离和密度之间的变化关系，是一种较为理想的方法。该算法也可推广应用于交通流其他类似属性方程的数值求解中，格式的推导和程序的实现过程类似于LWR模型方程。

基金项目

国家自然科学基金(11971075)。

NOTES

^*通讯作者。