网络交通流模型的高分辨率迎风格式

A High Resolution Upwind Scheme for Traffic Flow on Networks

摘要: 本文利用双曲守恒律方程描述车流密度的变化规律，选择一种同时满足CBC (Convection Boundedness Criterion)准则和TVD (Total Variation Diminishing)限制条件的HPUS高阶有界格式，并利用该格式对交通流问题中的一些典型算例进行数值模拟，对比结果表明该格式在间断处可有效抑制非物理振荡，并具有良好数值逼近效果。

Abstract: In this paper, the hyperbolic conservation law is used to describe the change law of the traffic flow density. Some typical examples of traffic flow problems are numerically simulated with a HPUS high-order bounded scheme that satisfies both CBC condition and TVD constraint. The results show that the scheme can effectively suppress non-physical oscillations near the discontinuities, and has good performance of numerical approximation.

1. 引言

交通流问题一直是计算流体力学研究的主要方向。现代交通流理论研究始于上世纪30年代初期，1935年Greenshields提出用速度-密度线性模型描述道路车流量与车速间的关系 [1]。1936年Adams在论文中提出将道路车流量看作随机独立变量 [2]，利用概率论知识来讨论车流量与车速之间的关系。1950年左右，人均汽车使用率显著提高，城市道路交通系统建设加快，James Lighthill、Gerald Whitham [3] 和Richards [4] 提出将流体力学的概念应用于交通问题中，即宏观交通流连续介质理论(经典LWR理论)，该理论旨在将车流比拟为某种连续流体，利用相关双曲守恒律方程来描述车流密度的基本变化规律，此外还提出了相关的激波理论。

交通流问题的关键在于求解上述双曲守恒律方程，一般很难求得该类方程的精确解，故需要构造一些数值格式来求其数值解。但主要的求解困难在于，即使初始条件充分光滑，其数值解也会随时间的变化在某一时刻产生间断。当前求解该类方程的数值计算方法有很多，如有限差分法、有限元法、有限体积法等。有限差分法在处理简单问题时效率较高，但不适于复杂区域问题；有限元方法在处理复杂边界问题时具有较大优势，但面对大梯度问题时较为困难，有限体积法是基于无结构网格下，处理强间断问题或大梯度问题的一种高效方法。

现有的经典数值格式有一阶迎风格式FOU (First-Order Upwind)、Lax-Friedrichs格式 [5]，这些格式精度较低，数值耗散较大。而中心差分格式CD (Central Difference)、二阶迎风格式SOU (Second-Order Upwind)、QUICK格式 [6] (Quadratic Upstream Interpolation for Convective Kinematic)、CUI格式 [7] (Cubic Upwind Interpolation)数值耗散较低，在光滑初值下可达到二阶及以上的精度，但在间断初值上讨论时格式不满足对流有界性，易在间断解处产生非物理振荡。为了在间断处抑制振荡，Gaskell和Lau提出CBC准则 [8] (Convection Boundedness Criterion)，可以保证数值格式的对流有界性，但并不能保证数值格式的精度。1983年，Harten提出高分辨率TVD [9] (Total Variation Diminishing)全变差格式，在满足有界性的同时保证了格式的高分辨率。本文利用一种既满足TVD限制条件又满足BAIR准则的HPUS [10] (Hermite Polynomial Upwind Scheme)高阶有界格式对以下交通流密度方程进行数值模拟得到流体的运动规律，进而解决流体运动中遇到的实际问题。

本文主题内容安排如下：第二节介绍网络交通流模型的相关背景，并简要分析两种情况下交叉点处流通量的判别情况；第三节简要分析HPUS格式的构造过程；第四节利用HPUS格式对网络交通流模型中的两种情况进行数值模拟，并对结果进行分析；第五节对上述工作进行简要总结。

2. 网络交通流模型

2.1. LWR模型

在单一路径中，非线性交通流问题可转化为如下双曲守恒律方程，并利用此方程来描述车流密度的变化规律 [3] [4]

${\partial }_t\rho +{\partial }_{x}f\left(\rho \right)=0$

其中 $\rho =\rho \left(t,x\right)\in \left[0,{\rho }_{\max }\right]$ 表示t时刻时不同位置的车流密度， ${\rho }_{\max }$ 表示该路段可达到的最大车流密度。 $f\left(\rho \right)=\rho v$ 表示流通量，指单位时间单位面积内的车流通量，该函数为严格凹函数，v是给定的只与密度有关的光滑单调递减函数，通常我们假设 $f\left(0\right)=f\left({\rho }_{\max }\right)=0$。由于函数 $f\left(\rho \right)$ 为严格凹函数故必存在唯一的极大值点 $\sigma$，称其为临界密度，且满足 $f^{\prime }\left(\rho \right)=0$。当 $0\le \rho \le \sigma$ 时，上述流量被认为是自由流动的，函数 $f\left(\rho \right)$ 为密度的增加函数；而当 $\sigma \le \rho \le {\rho }_{\max }$ 时，函数 $f\left(\rho \right)$ 为密度的递减函数，表明此时路上出现堵塞情况。

本文所研究的路网模型包含有限条路 [11]，这些路交汇于某些交叉点处。例如在交叉路口J处，包含n条入路记作 $I_i\left(i=1,2,\cdots ,n\right)$，m条出路记作 $I_j\left(j=n+1,n+2,\cdots ,n+m\right)$，记入路车流密度为 ${\rho }_i\left(i=1,2,\cdots ,n\right)$，出路车流密度为 ${\rho }_j\left(j=n+1,n+2,\cdots ,n+m\right)$，交叉口J处满足入路车流总流量等于出路车流总流量：

${\displaystyle \underset{i=1}{\overset{n}{\sum }}f\left({\rho }_i\left(t,b_i\right)\right)=}{\displaystyle \underset{j=n+1}{\overset{n+m}{\sum }}f\left({\rho }_j\left(t,a_j\right)\right)}$

每条出路的车流量分布可由如下交通流分布矩阵A表示 [12] ：

$A=\left(\begin{array}{cccc}{\alpha }_{n+1,1}& {\alpha }_{n+1,2}& \cdots & {\alpha }_{n+1,n}\\ {\alpha }_{n+2,1}& {\alpha }_{n+2,2}& \cdots & {\alpha }_{n+2,n}\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{n+m,1}& {\alpha }_{n+m,2}& \cdots & {\alpha }_{n+m,n}\end{array}\right)$

其中 ${\alpha }_{ji}\left(0\le {\alpha }_{ji}\le 1\right)$ 为流量分配系数，表示车量从入路i驶入并从出路j驶出的百分比。

网络交通流模型中道路交叉点的处理起到至关重要的作用，一个固定交叉点J处的Riemann问题即可表示为一个带初值的Cauchy问题，且其初值在交叉点处的每条道路上都是常数。上述双曲守恒方程在整个路网结构下的Riemann解是由每一交叉点处的Riemann解演化而来，为了在交叉点处取到Riemann问题的唯一解，需满足以下条件 [13] [14] ：

1) 矩阵A每一列的百分比相加和为1，即：

${\displaystyle \underset{j=n+1}{\overset{n+m}{\sum }}{\alpha }_{ji}}=1,\left(j=n+1,\cdots ,n+m;i=1,2,\cdots ,n\right)$

2) 当满足条件1时，流通量函数 $f\left(\rho \right)$ 的最大值在交叉点处取到。若入路数比出路数多时，即当 $n>m$ 时，此结论不成立；

3) 当入路数比出路数多时，即当 $n>m$ 时。现假设交叉点入口处有a，b两条路，出口处只有c一条路，从a，b两条路驶入的车不一定都能从c驶出，可能有一部分会滞留在交叉口J处，故在出路c上要设定一个最大输出量Q，设 $q\left(0\le q\le 1\right)$ 为路径参数，记从入路a进来的车流量为 $q\cdot Q$，从入路b进来的车流量为 $\left(1-q\right)\cdot Q$。

可以发现在交叉口处求解Riemann问题相等于在每条路上求解方程的初边值问题，由于方程的非线性性，初边值问题的解可能达不到边界值，所以必须对每条道路的最大流量进行限制 [15]。

现假设每条路的初始车流密度为： ${\rho }_{1,0},{\rho }_{2,0},\cdots ,{\rho }_{n,0},{\rho }_{n+1,0},\cdots ,{\rho }_{n+m,0}$，密度严格最大值即临界密度为 $\sigma$，此时入路的最大流通量 ${\gamma }_i^{\max }$ 满足：

${\gamma }_i^{\max }=\{\begin{array}{l}f\left({\rho }_{i,0}\right),{\rho }_{i,0}\in \left[0,\sigma \right]\\ f\left(\sigma \right),{\rho }_{i,0}\in \left[\sigma ,{\rho }_{\max }\right]\end{array}i=1,\cdots ,n$

出路的最大流通量 ${\gamma }_j^{\max }$ 满足：

${\gamma }_j^{\max }=\{\begin{array}{l}f\left({\rho }_{j,0}\right),{\rho }_{j,0}\in \left[0,\sigma \right]\\ f\left(\sigma \right),{\rho }_{j,0}\in \left[\sigma ,{\rho }_{\max }\right]\end{array}j=n+1,\cdots ,n+m$

其中 $\sigma$ 为密度最大值点，有 $f^{\prime }\left(\sigma \right)=0$ 成立，由上式可知在交叉口处入路流通量小于等于需求流通量，即 $f\left({\rho }_{i,0}\right)\le f\left(\sigma \right)$，且交叉点处的出路流通量小于等于供应流通量，即 $f\left({\rho }_{j,0}\right)\le f\left(\sigma \right)$，一些车流密度会在交叉点处根据上述式子重新确定。

2.2. 交叉点处流通量

交通流的问题关键在于如何确定交叉点J的流通量，本文主要讨论以下两种路网结构的黎曼问题：(A)一条入路和一条出路(Bottelneck问题)；(B)两条入路和一条出路。

2.2.1. 一条入路a和一条出路b

考虑最简单的由一条入路a和一条出路b构成的交通路网问题。首先根据入路、出路最大流通量的定义式分别确定入路a和出入b在交叉点J处最大流通量 ${\gamma }_a^{\max }$ 、 ${\gamma }_b^{\max }$，之后取二者的最小值记为交叉点J处的流通量，记作：

${\hat{\gamma }}_a={\hat{\gamma }}_b=\hat{\gamma },\hat{\gamma }=\min \left\{{\gamma }_a^{\max },{\gamma }_b^{\max }\right\}$

则J处入流道和出流道的车流量取为：

${\hat{f}}_{N_a+\frac{1}{2}}^a={\hat{\gamma }}_a,{\hat{f}}_{\frac{1}{2}}^b={\hat{\gamma }}_b$

2.2.2. 两条入路a、b和一条出路c

考虑由两条入路a、b和一条出路c构成的交通路网模型，首先同样根据入路、出路最大流通量的定义分别确定入路a、入路b和出路c在交叉点J处的最大流通量 ${\gamma }_a^{\max }$ 、 ${\gamma }_b^{\max }$ 、 ${\gamma }_c^{\max }$，此时交叉点处的流通量需分两种情况：

Case1：两条入路在J处的最大流通量之和小于出入在J处的最大流通量，即当 ${\hat{\gamma }}_a^{\max }+{\hat{\gamma }}_b^{\max }\le {\hat{\gamma }}_c^{\max }$ 时，两条入路的车流均可顺利通过交叉点到达出路c，则记

${\hat{\gamma }}_a={\gamma }_a^{\max },{\hat{\gamma }}_b={\gamma }_b^{\max },{\hat{\gamma }}_c={\hat{\gamma }}_a+{\hat{\gamma }}_b$

Case2：两条入路在J处的最大流通量之和大于出入在J处的最大流通量，即当时，两条入路的车流量之和超出出路c可承受车流量，故两条入路的车流需按2.1中条件3中所给出的比例进入c车道。现取路径参数为 $q\left(0\le q\le 1\right)$，此时c车道上的最大流通量为 ${\gamma }_c^{\max }$，则由路a驶入路c的车流量为 $q\cdot {\gamma }_c^{\max }$，由路b驶入路c的车流量为 $\left(1-q\right)\cdot {\gamma }_c^{\max }$，则交叉点J处三条路的实际流通量应满足以下三种情况：

1) 当 ${\gamma }_a^{\max }\ge q{\gamma }_c^{\max },{\gamma }_b^{\max }\ge \left(1-q\right){\gamma }_c^{\max }$ 时，则有

${\hat{\gamma }}_a=q{\gamma }_c^{\max },{\hat{\gamma }}_b=\left(1-q\right){\gamma }_c^{\max },{\hat{\gamma }}_c={\gamma }_c^{\max }$

2) 当 ${\gamma }_a^{\max }<q{\gamma }_c^{\max },{\gamma }_b^{\max }\ge \left(1-q\right){\gamma }_c^{\max }$ 时，则有

${\hat{\gamma }}_a={\gamma }_a^{\max },{\hat{\gamma }}_b={\gamma }_c^{\max }-{\gamma }_a^{\max },{\hat{\gamma }}_c={\gamma }_c^{\max }$

3) 当 ${\gamma }_a^{\max }\ge q{\gamma }_c^{\max },{\gamma }_b^{\max }<\left(1-q\right){\gamma }_c^{\max }$ 时，则有

${\hat{\gamma }}_a={\gamma }_c^{\max }-{\gamma }_b^{\max },{\hat{\gamma }}_b={\gamma }_b^{\max },{\hat{\gamma }}_c={\gamma }_c^{\max }$

则交叉点J处入路a、b和出路c的车流量分别取为：

3. HPUS格式

由前文可知一维非线性交通流双曲守恒方程为

${\partial }_t\rho +{\partial }_{x}f\left(\rho \right)=0$ (1)

其中 $\rho =\rho \left(x,t\right)$ 表示t时刻时不同位置的车流密度， $f\left(\rho \right)=\rho v$ 表示流通量，指单位时间单位面积的车流通量，现利用有限体积法对上述方程进行空间离散。

3.1. 空间离散

首先将区间 $\left[a,b\right]$ 平均剖分成N个控制单元，记： $a=x_{\frac{1}{2}}\le x_{\frac{3}{2}}\le \cdots \le x_{\left(N-1\right)+\frac{1}{2}}\le x_{N+\frac{1}{2}}=b$，每个控制单元记为 $I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right]\left(j=1,2,\cdots ,N\right)$，现定义控制单元的积分平均值为： $\bar{\rho }\left(x_j,t\right)=\frac{1}{\Delta x}{\displaystyle {\int }_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\rho \left(x,t\right)\text{d}x}$。

在区间上对方程(1)进行积分化简得：

$\frac{\text{d}{\bar{\rho }}_j\left(t\right)}{\text{d}t}=-\frac{1}{\Delta x}\left({\hat{f}}_{j+\frac{1}{2}}-{\hat{f}}_{j-\frac{1}{2}}\right)$

其中 ${\bar{\rho }}_j\left(t\right)$ 、 ${\hat{f}}_j$ 分别为 $\bar{\rho }\left(x_j,t\right)$ 、 $f\left({\rho }_j,t\right)$ 的近似值。

之后引入全局Lax-Friedrichs数值流通量来确定 ${\hat{f}}_{j+\frac{1}{2}},{\hat{f}}_{j-\frac{1}{2}}$ 的值 [16]

$\begin{array}{l}{\hat{f}}_{j+\frac{1}{2}}=\frac{1}{2}\left[\hat{f}\left({\rho }_{j+\frac{1}{2}}^{+}\right)+\hat{f}\left({\rho }_{j+\frac{1}{2}}^{-}\right)-\alpha \left({\rho }_{j+\frac{1}{2}}^{-}-{\rho }_{j+\frac{1}{2}}^{+}\right)\right]\\ {\hat{f}}_{j-\frac{1}{2}}=\frac{1}{2}\left[\hat{f}\left({\rho }_{j-\frac{1}{2}}^{+}\right)+\hat{f}\left({\rho }_{j-\frac{1}{2}}^{-}\right)-\alpha \left({\rho }_{j-\frac{1}{2}}^{-}-{\rho }_{j-\frac{1}{2}}^{+}\right)\right]\end{array}$

其中 $\alpha =\underset{u}{\max }\left|f^{\prime }\left(\rho \right)\right|$，之后利用HPUS格式重构 ${\rho }_{j+\frac{1}{2}}^{+}\rho {,}_{j+\frac{1}{2}}^{-},{\rho }_{j-\frac{1}{2}}^{+},{\rho }_{j-\frac{1}{2}}^{-}$。

3.2. 对流有界性

现取三个相邻的网格节点U、C、D分别表示上游、中心、下游网格节点如图1所示，这三个节点所

Figure 1. Three neighboring mesh points and the mesh face

图1. 三个相邻节点与单元边界

对应的函数值分别为 ${\varphi }_U,{\varphi }_C,{\varphi }_D$，取f表示均匀网格上的界面，则 ${\varphi }_f$ 表示均匀网格的界面值，通过相邻的三个节点对其进行重构，可将一些常见的线性格式统一为如下表达式：

${\varphi }_f={\varphi }_C+\left[\frac{1+\kappa }{4}\left({\varphi }_D-{\varphi }_C\right)+\frac{1-\kappa }{4}\left({\varphi }_C-{\varphi }_U\right)\right]$

其中 $k\in \left[-1,1\right]$，此时 ${\varphi }_f$ 的取值与 ${\varphi }_U,{\varphi }_C,{\varphi }_D$ 有关，由如下Leonard正则化思想 [17]

$\hat{\varphi }=\frac{\varphi -{\varphi }_U}{{\varphi }_D-{\varphi }_U}$

将上述格式转化为如下正则化形式(Normalized Variable Formulation, NVF)为

${\hat{\varphi }}_f=\left(1-\frac{1}{2}k\right){\hat{\varphi }}_c+\frac{1}{4}\left(1+k\right)$

发现正则化后 ${\varphi }_f$ 的值只与 ${\varphi }_C$ 和参数k有关，表1给出几种常见线性对流格式正则化前后的形式：

Table 1. System resulting data of standard experiment

表1. 几种常见线性对流格式及其正则化格式

这些格式在光滑初值下都可达到二阶及以上的精度，但对于间断初值不满足对流有界性，在间断解处产生非物理振荡。为了在间断处有效抑制振荡，Gaskell和Lau提出了CBC准则 [8]，但CBC准则只能保证数值格式的有界性，并不能保证数值格式的精度，为了提高数值格式的精度，Wei和Hou提出BAIR准则，即改进的CBC准则 [18] [19] ：

$\{\begin{array}{l}\frac{3}{2}{\hat{\varphi }}_c\le f\left({\hat{\varphi }}_c\right)<\frac{1}{2}\left({\hat{\varphi }}_c+1\right),0<{\hat{\varphi }}_c<\frac{1}{2}\\ \frac{1}{2}\left({\hat{\varphi }}_c+1\right)\le f\left({\hat{\varphi }}_c\right)\le \frac{3}{2}{\hat{\varphi }}_c,f\left({\hat{\varphi }}_c\right)\le 1,\frac{1}{2}\le {\hat{\varphi }}_c<1\\ {\hat{\varphi }}_f={\hat{\varphi }}_c,{\hat{\varphi }}_c\ge 1\\ {\hat{\varphi }}_f={\hat{\varphi }}_c,{\hat{\varphi }}_c\le 0\end{array}$

此外，TVD约束条件在归一化变量公式可定义为

$\{\begin{array}{l}{\hat{\varphi }}_f\le 1,{\hat{\varphi }}_f\in \left[{\hat{\varphi }}_c,2{\hat{\varphi }}_c\right],0<{\hat{\varphi }}_c<1\\ {\hat{\varphi }}_f={\hat{\varphi }}_c,{\hat{\varphi }}_c\le 0\\ {\hat{\varphi }}_f={\hat{\varphi }}_c,{\hat{\varphi }}_c\ge 1\end{array}$

现以 ${\hat{\varphi }}_c$ 图为横坐标， ${\hat{\varphi }}_f$ 为纵坐标将BAIR准则和TVD约束条件表示在图2中，当格式的NV线落在图中阴影部分时表示该格式既满足BAIR准则又满足TVD约束条件。

Figure 2. The regions of the TVD (shaded) and BAIR (hatched)

图2. BAIR区域(虚线部分)和TVD区域(阴影部分)

3.3. 对流项离散

现已知高分辨率格式在归一化变量公式下可统一表示为： ${\hat{\varphi }}_f=f\left({\hat{\varphi }}_c\right)$，由Zijlema和Wesseling所提出的结论可知 [20]，下列性质适用于这类格式：

1) 若对流格式的NV线过点 $\left(\frac{1}{2},\frac{3}{4}\right)$，则该式具有二阶局部截断误差。

2) 对流格式的NV线需满足 $f\left(\frac{1}{2}\right)=\frac{3}{4},f^{\prime }\left(\frac{1}{2}\right)=1-\frac{1}{2}k$。

3) 为使对流格式同时满足BAIR准则和TVD约束条件，格式的图像应落在图中阴影部分且 ${\hat{\varphi }}_c\in \left[0,1\right]$。

故所取对流格式需满足：当 ${\hat{\varphi }}_c\in \left[0,1\right]$ 时，方程需过点 $\left(0,0\right)$，点 $\left(\frac{1}{2},\frac{3}{4}\right)$ 与点 $\left(1,1\right)$，且满足 $f^{\prime }\left(0\right)={\theta }_1$， $f^{\prime }\left(\frac{1}{2}\right)=1-\frac{1}{2}k$， $f^{\prime }\left(1\right)={\theta }_2$，由图可知当 ${\theta }_1\in \left[\frac{3}{2},2\right]$， ${\theta }_2\in \left[0,\frac{1}{2}\right]$ 时，格式的图像落在阴影部分。

现取 ${\theta }_1=2$， ${\theta }_2=0$ 得满足上述条件的HPUS格式为：

${\hat{\varphi }}_f=\{\begin{array}{l}{\hat{\varphi }}_c\left(-4{\hat{\varphi }}_c^4+10{\hat{\varphi }}_c^3-8{\hat{\varphi }}_c^2+{\hat{\varphi }}_c+2\right),0<{\hat{\varphi }}_c<1\\ {\hat{\varphi }}_c,其他\end{array}$

3.4. 时间离散

空间离散之后方程变为只与时间有关的常微分方程，为了有效的抑制振荡，利用三阶Runge-Kutta方法对时间进行离散，三阶R-K格式如下：

$\{\begin{array}{l}{\rho }^{\left(1\right)}={\rho }^n+\Delta tL\left({\rho }^n\right)\\ {\rho }^{\left(2\right)}=\frac{3}{4}{\rho }^n+\frac{1}{4}\left[{\rho }^{\left(1\right)}+\Delta tL\left({\rho }^{\left(1\right)}\right)\right]\\ {\rho }^{n+1}=\frac{1}{3}{\rho }^n+\frac{2}{3}\left[{\rho }^{\left(2\right)}+\Delta tL\left({\rho }^{\left(2\right)}\right)\right]\end{array}$

4. 数值算例

本节我们将运用上述HPUS格式对文献 [13] [14] [16] 中两个较为复杂的路网模型进行数值模拟，并将Jiang、Shu所提出的五阶WENO格式 [21] 在4000个网格点上对方程进行数值模拟的结果记为参考解，图像上表示为“Reference”。以下数值算例中，三阶(SSP)-RK方法均在CFL条件下使用。

4.1. Bottleneck问题

Bottleneck问题为交通流的不同宽度的单一路径模型，由一条较宽的入路Road a和一条较窄的出路Road b构成，J表示两条道路的相交点即节点。由于两条路具有不同的宽度，从而两条道路具有不同的流通量函数。Bottleneck问题的实质即讨论不同宽度单一路径模型中从宽的入路a到窄的出路b的车流密度变化规律。

首先在计算区间[0,2]上考虑交通流双曲守恒方程 ${\partial }_t\rho +{\partial }_{x}f\left(\rho \right)=0$，取 $J=1$，[0,1]区间表示较宽的入路车道a，所对应的流通量函数为

$f_1\left(\rho \right)=\rho \left(1-\rho \right),\rho \in \left[0,1\right]$

当 ${\sigma }_1=\frac{1}{2}$ 时，可取到入路a上的最大流通量： $f_1\left({\sigma }_1\right)=\underset{\rho \in \left[0,1\right]}{\max }{f}_1\left(\rho \right)=\frac{1}{4}$，

则[1,2]区间表示较窄的出路车道b，所对应的流通量函数为

$f_2\left(\rho \right)=\rho \left(1-\frac{3}{2}\rho \right),\rho \in \left[0,\frac{2}{3}\right]$

当 ${\sigma }_2=\frac{1}{3}$ 时，可取到出路b上的最大流通量： $f_2\left({\sigma }_2\right)=\underset{\rho \in \left[0,\frac{2}{3}\right]}{\max }{f}_2\left(\rho \right)=\frac{1}{6}$。

现假设 ${\rho }_1\left(t,x\right)、{\rho }_2\left(t,x\right)$ 分别表示入路a和出路b的t时刻不同位置的车流密度，取N = 200，CFL = 0.05，下面讨论不同初值条件及边界条件下的瓶颈问题。

情形1：初值条件及边界条件为：

$\begin{array}{l}\{\begin{array}{l}{\rho }_1\left(0,x\right)=0.66,x\in \left[0,1\right]\\ {\rho }_2\left(0,x\right)=0.66,x\in \left[1,2\right]\end{array}\\ {\rho }_{1,b}\left(t,0\right)=0.25\end{array}$

其中初始条件所给车流密度趋于出路b的最大车流密度，随着时间推移我们可以观察到交通堵塞的形成，详见 [13]。现利用HPUS格式求解时间T分别为0.5，1.0，4.0时的数值解，并在图3中与五阶WENO格式所给出的参考解进行对比分析。

情形2：初值条件及边界条件为：

$\begin{array}{l}\{\begin{array}{l}{\rho }_1\left(0,x\right)=0,x\in \left[0,1\right]\\ {\rho }_2\left(0,x\right)=0,x\in \left[1,2\right]\end{array}\\ {\rho }_{1,b}\left(t\right)={\rho }_1\left(t,0\right)=0.4\end{array}$

利用HPUS格式计算时间T为2.0，4.0，10.0时的数值解，并在图4中与五阶WENO格式所给出的参考解进行对比分析。

情形3：初值条件及边界条件为：

$\begin{array}{l}\{\begin{array}{l}{\rho }_1\left(0,x\right)=0.4+0.2\sin \left(5\pi x\right),x\in \left[0,1\right]\\ {\rho }_2\left(0,x\right)=0.66,x\in \left[1,2\right]\end{array}\\ {\rho }_{1,b}\left(t\right)={\rho }_1\left(t,0\right)=0.25\end{array}$

利用HPUS格式计算时间时间T为0.2，0.5，0.7时的数值解，并在图5中与五阶WENO格式所给出的参考解进行对比分析。

由图3~5可知HPUS格式可以很好的近似参考解，在间断处及大梯度处都可以有效地抑制非物理震荡。

Figure 3. Bottleneck problem with initial and boundary data—Case 1

图3. 情形1初边值条件下Bottleneck问题的数值解

Figure 4. Bottleneck problem with initial and boundary data—Case 2

图4. 情形2初边值条件下Bottleneck问题的数值解

Figure 5. Bottleneck problem with initial and boundary data—Case 3

图5. 情形3初边值条件下Bottleneck问题的数值解

4.2. 两条入路一条出路模型

现考虑图6所示的包含两条入路车道a，b及一条出路车道c的交叉口情况，每一条路的计算区间均取[0,1]，每条路的流通量函数皆满足 $f\left(\rho \right)=\rho (\left(1-\rho \right),\rho \in \left[0,1\right]$。

Figure 6. A junction with two incoming and one outgoing roads

图6. 两条入路一条出路的交叉口示意图

现给出如下初值条件：

$\begin{array}{l}{\rho }_1\left(0,x\right)=\{\begin{array}{l}0.1+0.1{\sin }^3\left(5\pi x\right),x\in \left[0,0.2\right]\cup \left[0.4,0.6\right]\cup \left[0.8,1\right]\\ 0.2,其他\end{array}\\ {\rho }_2\left(0,x\right)=\{\begin{array}{l}0,x\in \left[0,0.2\right]\cup \left[0.4,0.6\right]\cup \left[0.8,1\right]\\ 0.1+0.1\sin \left(5\pi x\right),其他\end{array}\\ {\rho }_3\left(0,x\right)=0.1\end{array}$

及边界条件：

$\begin{array}{l}{\rho }_{1,b}\left(t,0\right)=0.1\\ {\rho }_{2,b}\left(t,0\right)=0.1\end{array}$

图7~9分别表示入路a，入路b，出路c在N = 500，CFL = 0.01，时间T分别取0.5，1.0时的HPUS格式数值解及五阶WENO格式参考解。由所给数值结果可知，HPUS格式数值解可以很好的近似五阶WENO格式参考解，且相较于参考解HPUS格式可以在较粗网格下达到良好的逼近效果，从而有效的提高了计算速度。

Figure 7. Numerical results for Road a at T = 0.5, 1.0.

图7. T = 0.5，1.0时Road a的数值解

Figure 8. Numerical results for Road b at T = 0.5, 1.0.

图8. T = 0.5，1.0时Road b的数值解

Figure 9. Numerical results for Road c at T = 0.5, 1.0.

图9. T = 0.5，1.0时Road c的数值解

5. 结论

本文基于交通流理论，将交通流基本规律转化为相关的非线性双曲守恒律方程，分析讨论两种简单路网模型交叉点处的流通量，在有限体积法下对上述双曲守恒方程进行数值离散，利用高分辨率Hermite多项式迎风格式逼近对流项，通过相关交通流算例分析，HPUS格式在间断处及大梯度处可以有效的抑制非物理振荡，且其数值解可在较粗尺寸的网格下达到良好的逼近效果，从而可以正确描述交通流特性，准确实时的模拟交通过程，以便于对交通堵塞等问题提前做出预判。

致谢

作者感谢审稿专家给予本文宝贵的修改意见和建议，感谢期刊编辑对于本文的付出。

基金项目

本文由内蒙古大学科研发展基金项目(21100-5187133)支持。