ABSTRACT

This paper presents a CFD numerical simulation method for ship E/R ventilation. Through the discretization of the fluid control equation, the basic physical model of the engine room is established by using GAMBIT software, and the air intake of the engine room is optimized according to the air supply quantity. The physical model is established into FLUENT, and the boundary condition, the solution method and the precision of the model are selected reasonably. Finally, the air flow distribution in different ventilation modes is taken as an example, and the characteristics of each scheme are compared and analyzed, and the optimal ventilation scheme is obtained. The method proposed in this paper has important theoretical and practical significance for optimizing the ventilation of the ship room.

Keywords:Ship Engine Room, Air Ventilation, Numerical Simulation, Scheme Optimization

基于CFD技术的船舶机舱空气速度场数值模拟

陈建平¹，徐洁²，龚幼¹，黄鹏举¹，陈俊发¹

¹广州航海学院船舶工程学院，广东 广州

²广东工业大学自动化学院，广东 广州

收稿日期：2017年11月8日；录用日期：2017年11月22日；发布日期：2017年11月29日

摘 要

文章提出了船舶机舱通风的CFD数值模拟方法。首先建立起流场的基本控制方程，然后通过流体控制方程的离散化，使用GAMBIT软件建立了船舶机舱的基本物理模型，并根据供气量优化了机舱进气口分布。将物理模型建立导入FLUENT，合理选择模型的边界条件，获得模型的仿真数值解。最后，以不同通风方式下的空气流量分布为例，对各种方案的特点进行比较分析，得出最优通风方案。本文提出的方法对于优化船舶机房通风具有重要的理论和现实意义。

关键词 :船舶机舱，机舱通风，数值模拟，方案优化

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

船舶机舱为船舶主要动力装置提供空间，包括主机、锅炉和增压装置等。船上机舱通风换气系统与陆上相应系统比较具有明显的特点，例如船舶机舱空间狭窄就是其最大特点，单位容积余热量为陆上的5~10倍 [1] ，加热设备局部区域分布不均匀，温度梯度大等。同时随着科技的发展，节能环保也成为船舶通风技术改进的重要研究方向 [2] 。

在船舶机舱通风领域，国外文献资料很少，国内研究主要集中在大学和船舶研究机构。船舶机舱的通风属于大型空间通风，与火力电站和锅炉房通风系统相似，但机舱通风的特性决定了其通风设计的难度。郝维坚 [3] 对集装箱船机舱气流进行了仿真计算，通过CFD软件PHOENICS计算，分析了机舱的速度场和温度场，并通过调整空气通风口的位置，研究机舱通风设计数值模拟的可行性。周山 [4] 使用FLUENT对客滚轮机舱通风系统进行了模拟，通过对内部温度场和速度场分布的分析表明，机舱内温度分布不均匀，但改变入口位置和出口角度可以降低泵送温度。陈宁，张东 [5] 利用CFD软件AIRPAK对油轮机舱通风系统进行了模拟分析，详细介绍了该机舱通风系统的设计和计算，具有很大的参考价值。AIRPAK作为HVAC气流在特殊条件下的数值模拟软件设计简单，后期结果表明FLUENT更加方便，但在复杂舱室模型建模方面存在局限性，网格划分困难。刘凤荣，孟宪松等 [6] 使用FLUENT模拟船舶机舱的热环境，并分析了机舱内的三维温度场和速度场分布。索文超、王宪成 [7] 通过数值模拟，将机舱三维温度场和速度场与实验值进行比较，认为机舱空气流动的数值模拟结果是可行的。

本文通过船舶机舱通风设计，采用CFD技术对机舱通风系统进行数值模拟，比较各工作条件以获得最佳通风方案，为通风设计提供依据。研究舱室气流速度场的分布，为改善机舱气流死区和涡流，合理分布空气，减少热量累积，消除有害气体提供理论依据。GAMBIT是为了帮助分析者和设计者建立并网格化计算流体力学(CFD)模型和其它科学应用而设计的一个软件包。GAMBIT GUI可以直接实现建立模型、网格化模型、指定模型区域大小等基本步骤。

2. 流场基本控制方程

空气作为管道和机舱中的流体在其中流动，必须遵循相应流体运动物理定律，这些定律包括质量守恒、动量守恒和能量守恒等。围绕这些规律，建立流体运动的基本方程式根据流体的实际物理状态，求解这些方程将能够获得流体的流动状态和规律，这里的实际物理状态指的是层流或湍流，组份转移或化学反应 [8] 。

2.1. 质量守恒方程

质量守恒定律是自然界的普遍规律，机舱气体流动同样必须遵循该规律。质量守恒定律可以简单描述为任意瞬间流体控制体的质量增加等于该瞬间流入该控制体的流体的净质量。根据这个规律，可以获得机舱气体质量守恒方程

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$ (1)

式中， $i$ 分别表示坐标轴 $x$ ， $y$ 和 $z$ 方向， $u$ 为气体流动速度， $\rho$ 为气体密度， $t$ 为时间。

在通风和空调工程领域，空气通常是认为不可压缩的流体，即密度 $\rho$ 不随时间变化，则其连续性方程可以简化为

$\frac{\partial u_i}{\partial x_i}=0$ (2)

2.2. 动量守恒方程

在惯性坐标中，动量方程为如下描述

$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}+\rho g_i+F_i$ (3)

式中， $\rho g_i$ 为 $i$ 方向体积力， $p$ 为流体静压， ${\tau }_{ij}$ 为单元流体中法向为 $i$ 方向的表面上沿 $j$ 方向作用的切应力， $F_i$ 包含其它力源项。 ${\tau }_{ij}$ 由式(4)确定

${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_{ij}}{\partial x_i}\right)-\frac{2\mu }{3}\frac{\partial u_i}{\partial x_i}{\delta }_{ij}$ (4)

式中， $\mu$ 为动力粘性系数。等式右边第二项是考虑了气体体积膨胀的作用， ${\delta }_{ij}$ 为克罗内克(Kronecker)函数， $i=j$ 时， ${\delta }_{ij}=1$ ， $i\ne j$ 时， ${\delta }_{ij}=0$ 。

2.3. 能量守恒方程

在通风和空调工程领域，流体之间的能量传递形式通常用焓 $h$ ( $h={\displaystyle {\int }_{298.15}^T{c}_p\text{d}T}$ ，其中298.15是参考温度)来表示

$\frac{\partial \left(\rho h\right)}{\partial t}+\frac{\partial \left(\rho u_{i}h\right)}{\partial x_i}=\frac{\partial \left(k+k_t\right)}{\partial x_i}\frac{\partial T}{\partial x_i}+S_h$ (5)

式中， $t$ 为温度， $k$ 为流体分子运动产生的导热系数， $k_t$ 为湍流附加的导热系数( $k_t=c_p{\mu }_t/P{r}_t$ ， $P{r}_t$ 为湍流普朗克数， $P{r}_t=0.85$ )， $S_h$ 为源项，可以是任何体积热源。

对于固体区域的导热，可用导热方程表示

$\frac{\partial \left(\rho h\right)}{\partial t}=\frac{\partial }{\partial x_i}\left(k\frac{\partial T}{\partial x_i}\right)+q_{\nu }$ (6)

其中， $q_{\nu }$ 是内热源。

如果舱室内空气流动导热和对流并存时，必须联立式(6)与式(5)进行求解。

2.4. 组分输运方程

求解组份传输方程式的思路是通过第 $i$ 种物质的对流扩散方程来预测每种物质的质量分数 $Y_i$ 。 组分输运方程用以下通式表示

$\frac{\partial \left(\rho Y_i\right)}{\partial t}+\nabla \cdot \left(\rho u{Y}_i\right)=-\nabla \cdot J_i+R_i+S_i$ (7)

其中， $S_i$ 是由扩散项及用户定义的源项所导致的额外产生速率， $R_i$ 是由化学反应引起的该组分净生产速率。如果系统中存在着 $N$ 种物质相互作用，则需要求解 $\left(N-1\right)$ 个该形式方程。由于质量分数和必须为1，所以第 $N$ 种物质的分数等于1减去已解得的 $\left(N-1\right)$ 个质量分数之和。

1) 层流中的质量扩散

式(7)中， $J_i$ 为物质 $i$ 的扩散通量，由浓度梯度产生，可以采用稀释近似表示

$J_i=-\rho D_{i,m}\nabla Y_i$ (8)

式中， $D_{i,m}$ 为混合物中第 $i$ 种物质的扩散系数。

2) 湍流中的质量扩散

湍流中的质量扩散表示为

$J_i=-\left(\rho D_{i,m}+\frac{u_t}{S{c}_t}\right)\nabla Y$ (9)

式中， $S{c}_t$ 为湍流施密特数(一般默认为0.7)。

3) 能量方程中的物质输送处理

在多组分流动的过程中，组分的扩散伴随着能量的传递。对于焓场，这种扩散具有重要意义，因此不能简单地忽略。特别是所有物质的Lewis数

$L{e}_i=\frac{k}{\rho c_p{D}_{i,m}}$ (10)

当 $L{e}_i$ 远大于1时，忽略此项会产生严重误差。

综上，各守恒方程可统一描述为

$\frac{\partial \left(\rho \varphi \right)}{\partial t}+div\left(\rho u\varphi \right)=div\left(\Gamma grad\varphi \right)+S$ (11)

将其展开为

$\frac{\partial \left(\rho \varphi \right)}{\partial t}+\frac{\partial \left(\rho u\varphi \right)}{\partial x}+\frac{\partial \left(\rho v\varphi \right)}{\partial y}+\frac{\partial \left(\rho w\varphi \right)}{\partial z}=\frac{\partial }{\partial x}\left(\Gamma \frac{\partial \varphi }{\partial x}\right)+\frac{\partial }{\partial y}\left(\Gamma \frac{\partial \varphi }{\partial y}\right)+\frac{\partial }{\partial z}\left(\Gamma \frac{\partial \varphi }{\partial z}\right)+S$ (12)

式中， $\varphi$ 为代表 $u$ ， $v$ ， $w$ 和 $T$ 等求解变量， $S$ 为广义源项， $\Gamma$ 为广义扩散系数。

对于未确定的方程式， $\varphi$ 、 $\Gamma$ 和 $S$ 以特定形式表达。在通风和空调工程中，空气通常被视为不可压缩气体，而在微观上，相同分子之间的距离保持恒定，即气体的密度 $\rho$ 在宏观上不会改变，即式(11)或式(12)中 $\frac{\partial \rho }{\partial t}=0$ 。

3. 数值求解方法

流体运动控制方程是由一系列偏微分方程组成，求解方法可分为分析、实验和数值方法。一般来说，通过分析方法获得方程的解析解很难。对于流体流动和热传递的一些问题，必须采用实验和数值方法。数值解是空间中原始物理坐标的连续场量(如速度场，温度场和浓度场)之间的一个有限的集合值，通过使用某些规则，将这些节点变量与这些规则相联系的一组代数方程来替代原始代数方程组而获得的近似解 [9] 。CFD就是使用数值分析方法来获取变量的近似解。根据流体控制方程的离散模式和代数方程的近似求解，数值分析方法包括有限差分法，有限元法，有限元分析法和有限体积法。有限体积法将计算区域划分为一系列连续体积，并用基点代替每个控制体积。从有限体积法得到的离散方程可以保证守恒性，离散方程的系数具有明确的物理意义，被广泛应用于CFD计算。文章使用有限体积法。

3.1. 计算区域的离散化

计算域的离散化是数值方法的基础。区域离散化是将计算区域划分为几个子区域，并确定每个区域的节点，这个过程也称为网格生成。选择离散计算区域，有限数量的离散点表示原始连续空间。在计算区域离散化完成后，获得以下4个几何元素 [10] ：

1) 节点。待求解的未知物理量的几何位置；

2) 控制容积。适用于控制方程或守恒定律的最小几何单位；

3) 界面。它指定对应于每个节点的控制容积的分界面位置；

4) 网线。通过在坐标轴方向连接两个相邻节点形成的曲线簇。

在模拟机舱通风问题时，通常将模型分为若干控制体，每个控制体选择一个节点进行更换，并将控制体温度、速度，浓度等物理量存储在该节点上。

3.2. 控制方程的离散化

流体运动控制方程的离散化是将其偏微分方程的控制方程转化为每个节点上的代数方程组。 当使用有限体积法时，控制体的物理量及其导数需要通过节点的物理量进行插值运算。在机舱通风区域，设有任意的网格节点P，其在每个垂直交叉方向的相邻网格点称为E，W，S，N，T和B，如图1所示。

在如图1所示三维定常坐标系中，控制方程的一般形式(11)的表示为典型形式

$a_P{\varphi }_P=a_E{\varphi }_E+a_W{\varphi }_W+a_S{\varphi }_S+a_N{\varphi }_N+a_T{\varphi }_T+a_B{\varphi }_B+S_c$ (13)

简写形式

$a_P{\varphi }_P={\displaystyle \sum a_{nb}{\varphi }_{nb}}+S_c$ (14)

式中， $a$ 为差分方程的系数， $S_c$ 为源项，为建立差分方程所要求解的独立变量。

通过式(13)，控制体 $P$ 在计算域上进行积分，可以获得由控制体积法建立的一般方程

${\displaystyle \underset{e,w,s,n,t,b}{\sum }\left[\left(\rho u\varphi -{\Gamma }_{\varphi }grad\varphi \right)\times A\right]}={\left(S_{\varphi }V\right)}_P$ (15)

式中， $e$ 、 $w$ 、 $s$ 、 $n$ 、 $t$ 和 $b$ 为控制体东西南北上下各界面的中心， $V$ 为控制体的体积， $A$ 为控制体某

图1. 控制体单元

界面的面积，方向为该界面的外法线方向。

式(15)中对流项和扩散项中已经通过物理节点的插值转换为控制体界面上的物理量。有限体积法中的离散方程必须保持控制体界面的连续性，源项的斜率为负，中心系数等于相邻点系数和，符合上述原则，离散方程则可以实现数学和物理的正确性和合理性。

4. 数值模拟实例

4.1. 机舱通风模型

机舱的主要设备是汽轮机、燃油锅炉、涡轮增压器、高温除氧器、油水分离过滤器、部分烟管通道等。机舱一般可分为三层，底部通常为各舱室的顶板，顶板为机舱顶甲板，在底层和顶甲板之间为两层格栅。风机通风出口布置在格栅下面。由于舱内有许多设备和管道，结构复杂，必须简化模型，使模型适合于数值模拟。这种模型简化的基本原理是：

1) 复杂结构形状设备的正则化和简化减小对仿真结果的影响。机舱主机的形状非常不规则，所有的绘制都为网格划分，计算难度大大增加，所以简化成规则形状的组合，同样道理锅炉和涡轮增压机组和辅助设备都被正规体所取代。

2) 忽略对流场的热影响和空间流场产生的作用较小的设备。机舱里有管道密布，包括蒸汽管道、烟管和各种水电管道，列出这么多管道是非常困难的。本文模型只列出部分主机入口蒸汽管道和部分烟气管道。

3) 设备壁面的光滑和无厚度处理。通风数值模拟一般认为空气和固体壁面粘性力忽略，空气流经壁无摩擦损失，所以将机舱舱壁和设备等壁面进行无厚度光滑壁面简化。

4) 管道处理。在机舱顶部的主风道，虽然尺寸较大，但对位置较低的流场没有影响。同时支风管尺寸较小，本文不考虑内外热交换管。本文模型管道简化为在风口段为一个立方体，出风口简化为底面。

5) 考虑排、抽风围阱位于机舱甲板顶部，于流场影响不大，以风口面替代。

4.2. 数学模型选择

船舶机舱结构的复杂性是其主要特征，这也导致了机舱空间狭小，空气流动不规则，气流复杂。为了提高计算效率，需要对计算对象进行必要假设：

1) 恒定流场

在实际的空调和通风过程，当空气流以均匀的速度经出风口送到各个舱室，一段时间后，送风空气和舱内空气充分混合达到一个稳定的状态，于是舱内整个流场可视为定常流。

2) 机舱内空气不可压缩。

机舱送风速度较小，密度变化可以忽略，将机舱空气认为是不可压缩的。

3) 符合Boussinesq假设

舱内空气符合Boussinesq假设，即不考虑压力对密度的影响。

4) 机舱为密闭空间

为使机舱通风保持正压，舱内送排空气应保持平衡，因此应视整个机舱为一个密闭系统。

基于上述假设和流体流动理论及其模型的研究，本文选择标准 $k-\epsilon$ 双方程模型作为通风系统的数学模型。在近壁问题中，采用壁面函数法进行处理。由有关数理理论，文章选定模型在湍流和近壁层流区的计算都能满足计算要求，是通风空调领域常用的数学模型之一。

4.3. 风口位置对气流组织的影响

通过对气流组织各种形式的比较分析，本算例采用下送上回的送风方式，选择在两层格栅的底部均匀布置送风口，这样可以使下部区域的空气均匀，顺利地消除舱内下部的热量。根据风口布置的原理，通常将风口布置散热量较大的设备附近，所以主机和锅炉附近区域应增大送风量。根据空气量的计算结果，与上述安排，本算例初步安排如下风口位置、数量以及尺寸为在风口的尺寸是0.4 × 0.2 m，方向为垂直向下，其中第一层格栅个数25，第二层为28个，具体分布见图2。

据送风口和回风风口的分布情况，文章选择了六种风口布置方案进行模拟比较。方案1设置单回程回风，回风口位于舱顶中央。方案2为双回风布置，回风围阱距离为6 m。方案3为四回出口，并均匀设置在锅炉的顶部位置，侧向送风口通过风阀调节开闭状态。具体方案如表1所示。

4.4. 速度场分析

经过模拟运算后，可以从热负荷分布看， $y=0\text{m}$ 横截面之间的两锅炉主机的右侧部位，是最大的散热面积，该区域与回风口正交，故流场也最为复杂。下面选取 $x=-8\text{m}$ 和 $z=-3.5\text{m}$ 截面来分析机舱下部区域的气流组织，选定截面位置见图3。

速度矢量图可以显示空气的流向和速度，因此本文用速度矢量图反映机舱内气流情况。通过对仿真结果的分析，在机舱控制室下部增加送风口和回风口位置的变化之间的影响非常小，其速度矢量图与气流结构相似，因此，文章只选取在 $y=0\text{m}$ 的截面上的方案1至4的速度矢量图进行分析，如图4所示。

如图4所示，我们可以看到在 $y=0\text{m}$ 截面，每个方案的送风量是相同的，但气流组织形式存在较大的差异的。

1) 在方案1中，机舱上部的空气流向顶部中间区域，而在回风口附近的气流速度较快。同时在锅炉和主机右侧有大量回流现象。

2) 方案2中由于设置了双回风口，增大了回风面积，但总回风量保持不变。对主机两侧的气体可以有效地抽除，而锅炉侧仍有大量回流，虽然回风面积有了加大，但单个回风口的风量却减少了一半，表

图2. 机舱风口分布图；(a) 第一层格栅；(b) 第二层格栅

表1. 风口布置方案

图3. 机舱分析截面示意图

图4. 在y = 0 m截面上的速度场云图；(a) 方案1；(b) 方案2；(c) 方案1；(d) 方案2

明回风速度下降了，所以上升到顶部的气体不能及时排除，同时阻碍下面气体的有效传输。

3) 方案3是从整个舱室的余热分配考虑的回风方案。锅炉侧热约为主侧的两倍，风口位置设计时侧部布置较多，风量也大。但主机侧没有安装回风口，使主机侧气流向锅炉侧整体运动，一部分从机舱中部排出，而一部分则继续测流与锅炉侧气流合流，从而在机舱中部形成一个较大的漩涡场。

4) 方案4的流态与方案1基本相同。类似地，方案5和方案6类似于方案2方案3。

通过上述比较和分析，如将回风口设置在锅炉侧的顶部，有利于锅炉侧气体的顺利消除，而机舱气流的整体结构也较为合理。根据分析，可以根据舱内空气量的分布确定回风口的位置，因此方案3和6的回风形式是上述方案中最好的。

选取 $x=-8\text{m}$ 截面及 $z=-3.5\text{m}$ 截面来分析增加的侧向风口对机舱气流组织的影响。数值模拟结果如图5和图6所示。

从图5(a)即方案3可以看出，在机舱控制室下部的气流形成一个大涡流。涡流从舱壁的左侧开始，沿着主机壁上升并回流到风口位置，同时气流从漩涡边界至中心呈现逐渐减少的趋势。该漩涡可以通过图6(a)进一步分析，图形显示控制室下部向主机侧运动气流较多，但向机舱中部区域的气流较少，这表明该涡流处于一个相对静止的状态，气体大部分在控制室下部流动，其结果是导致空气置换效果差。通过分析，该区域气流的合理组织形式是让控制室的下部气流绕过它，再上升到顶部。从图5(b)和图6(b)中可以看出，在增加斜向送气口之后，控制室下部的气流发生了明显的变化，大漩涡消失，因此，方案6比方案3的下空气置换效果要好。

4.5. 局部气流组织的改进

由于机舱内气流分布非常复杂，在通风过程中总会出现死角及漩涡，但必须避免存在安全隐患或影响较大的流场死区和漩涡。在假定通风方式不变的情况下，消除通风死角或者漩涡可参考以下两种方法。一是设置气流障碍物，虽然障碍物对舱室内气流阻力的影响是非常大的，但这种方法无疑凭空增加了通风空间的不必要的障碍，在实际工程中很少使用。二是采用主动诱导通风，在死区和漩涡区布置风口，该方法需要重新调整机舱通风系统，在小型船舶影响显著 [11] 。

图5. 在x = −8 m截面上的速度场云图；(a) 方案3；(b) 方案6

图6. 在z = −3.5 m截面上的速度场云图；(a) 方案3；(b) 方案6

算例中气流至底部后自下而上流动，各种障碍物迎风面不存在气流滞留区，算例采用在主机右侧安装隔屏来消除此漩涡，隔屏位置如图7所示。

设隔屏的宽度为4米，高5米，隔屏的宽度与主体右侧的宽度相同，底部第一层格栅的上方，高度超过第二层格栅，这个高度正好是生产的漩涡的顶部，气流会在这个高度与隔屏相碰撞。模拟时选取送风温度为37℃，送风气流速度为9.4 m/s，模拟计算完成后，选取 $y=-0.5\text{m}$ 的截面来定性分析效果。隔屏安装后气流的分布如图8所示。

图9(b)显示了在安装隔屏后主机面板右侧的速度矢量图。图9(b)没有按照隔屏的气流速度矢量图9(a)相比，由于隔屏的阻隔效果，主机上部的气流得到了提前提升，并且让贴附主机右侧壁的气体能够直接排放到机舱顶部的位置，从而从回风口排出，而位于隔屏位置处的漩涡消失。如图9(b)所示，锅炉和机舱中间的气流被引入到锅炉控制台的顶部，其结果是，来自右侧舱壁的上升气流产生少量回流，但不会产生大的漩涡，在短距离与其他上升气流汇合后，不会产生热量积聚的问题。图9显示只要条件允许的情况下，通过合理设置隔障可以有效提高气流死角和漩涡的通风效果。

图7. 隔屏安装位置示意图

图8. 机舱气流流线图

图9. 加装隔屏前后机舱速度场云图比较；(a) 加装之前；(b) 加装之后

5. 结论

本文研究了船舶机舱内空气流动的基本控制方程及其离散化，并用CFD方法进行了数值分析和仿真。在实例研究中，对六种回风口布置方案进行了模拟和比较。得出在锅炉侧顶部布置回风口的方案是较为合理的，同时，在机舱控制室下侧增加侧向送风口有利于局部气流组织改善。最后，通过研究设置隔屏对机舱气流的影响，得出合理设置隔屏可以减少回流和消除漩涡的结论。本文提出的CFD数值模拟方法对于设计机舱通风方案具有很好的理论指导和实用价值。

基金项目

广东省交通科技计划项目(2015-02-050; 2017-02-050)；广东省自然科学基金项目(2014A030313792; 2016A030313125)；广东省教育厅特色创新项目(2016-100)。

文章引用

陈建平,徐 洁,龚 幼,黄鹏举,陈俊发. 基于CFD技术的船舶机舱空气速度场数值模拟
A Numerical Simulation Method of Air Velocity Field in Ship Engine Room Based on CFD Technology[J]. 交通技术, 2017, 06(06): 299-309. http://dx.doi.org/10.12677/OJTT.2017.66038

参考文献 (References)