基于非结构RKDG实现多介质界面追踪的数值模拟方法

摘要

本发明提出一种基于非结构RKDG实现多介质界面追踪的数值模拟方法。在非结构网格上利用一系列标志点对多介质界面进行离散，在标志点法向构造黎曼问题并采用双激波近似方法求解。黎曼问题的解不仅用于推进界面，而且可以直接更新界面附近真实流体一侧的流体状态并外推得到对应的虚拟流体状态。同时，采用空间紧致性较好的非结构RKDG方法求解各个单介质流场。整个求解方案编程简单，计算高效。通过与基于有限体积格式得到的结果进行对比，本发明得到的守恒误差较小，数值稳定性更好。

说明书

技术领域

本发明属于数值模拟技术领域，具体涉及一种基于非结构RKDG实现多介质界面追 踪的数值模拟方法。

背景技术

可压缩多介质问题的求解主要包括两部分：一是对于各个单介质流场进行高精度 数值模拟，二是精确地处理多介质界面及其附近区域。在近二十年里，间断伽辽金(DG)方法 在求解单介质流场上逐渐成为研究热点。在对DG方法进行的众多研究中，值得一提的是 Cockburn等人提出的龙格-库塔间断伽辽金(RKDG)方法：时间离散上采用总变差减小的 (TVD)高阶龙格-库塔(RK)方法，在空间离散上通过求解黎曼问题获得网格界面通量，应用 总变差有界(TVB)限制器抑制数值解的虚假振荡。RKDG方法的主要优点：在流场的光滑区域 易于构造高精度格式；在流场任意位置都可以得到多项式形式的数值解；空间紧致性好，网 格的边界通量只与边界的左右网格单元有关。基于这些优势，RKDG方法在求解单介质流场 问题上取得了较大的成功。在多介质界面处理上，早期基于γ模型、质量分数模型、level set函数的方法在多介质界面处往往会产生较大的数值耗散。然而，对于界面追踪方法，界 面是显式追踪的而且界面位置在计算过程中能够始终保持非常清晰。虚拟流体方法主要用 于定义界面边界条件从而使得每种介质可以单独进行求解。但是当压力或者速度在界面附 近存在较大的梯度时，该方法的表现效果较差。事实上，虚拟流体状态的定义应当恰当地考 虑界面处波系的相互作用以及不同介质的流体属性。后续相关学者提出了一些改进版本的 虚拟流体方法，例如真实虚拟流体方法(RGFM)。该方法通过求解界面黎曼问题来更新界面 真实流体一侧的流体状态并外推得到虚拟流体状态。因为界面处的波系可以向上游和下游 同时传播，所以真实虚拟流体方法定义的界面边界条件具有较小的守恒误差。通过这些方 法的综合应用能对可压缩多介质问题进行较好的数值模拟。

到目前为止，国内外对于DG方法与虚拟流体方法耦合实现多介质界面追踪的研究 仍然较少，而且多介质界面问题多是在均匀结构网格上进行求解的。为了实现对单介质流 场的精确求解，早期的研究往往采用高阶精度的有限差分格式或者有限体积格式。然而，高 阶格式由于空间紧致性较差往往需要较多的模板点，在界面附近需要定义更多的虚拟流体 状态。这在一些具有复杂界面的多介质问题中往往是不利的。因为界面追踪方法在很大程 度上依赖于界面上的几何信息，在远离界面处求解法向量等几何参数存在较大的误差，从 而进一步导致远离界面处的虚拟流体状态存在较大的不精确性。

发明目的

由于现有技术在多介质界面追踪和虚拟流体状态求解上多是单独进行的，这样就 增加了编程和计算的复杂性。另外，在求解单介质流场时，高阶精度数值格式需要较多的模 板点，这样就不可避免地需要定义远离界面处的虚拟流体状态。然而在界面追踪方法中远 离界面的虚拟流体状态往往存在较大的误差，这就导致了现有技术在实现多介质界面追踪 上具有较大的不精确性。因此，本发明提出了一种更简单、更精确的多介质界面追踪方法。 在多介质界面构造黎曼问题，黎曼问题的解不仅用于推进界面，而且可以直接用于更新真 实流体状态和外推得到虚拟流体状态。另外，采用空间紧致性较好的非结构RKDG方法求解 单介质流场，大大减少了所需虚拟流体状态的信息。

发明内容

为了克服现有技术中的不足，本发明提供了一种基于非结构RKDG实现多介质界面 追踪的数值模拟解决方案，具体如下：

步骤1，在计算域上生成非结构网格；

步骤2，根据界面的位置确定虚拟流体区域；

步骤3，利用一系列标志点离散多介质界面，在标志点法向构造黎曼问题，采用双 激波近似方法求解。根据黎曼问题的解得到标志点的运动速度；

步骤4，直接利用标志点上黎曼问题的解更新界面附近真实流体一侧的流体状态， 并外推得到对应的虚拟流体状态；

步骤5，采用非结构RKDG方法对每种介质单独求解；

步骤6，由标志点的速度推进多介质界面，重构界面，得到新的界面位置和标志点 并根据新的界面确定整个多介质流场的数值解；

步骤7，判断是否满足计算终止条件，如计算还需继续进行，返回步骤2进行下一时 刻的计算。

本发明采用上述技术方案，具有以下有益效果：

(1)非结构网格对于计算域边界具有较强的适应性；

(2)在界面构造黎曼问题不仅用于推进界面，而且可直接用于定义虚拟流体状态， 编程简单，计算高效。

(3)采用空间紧致性较好的RKDG方法进行计算使得所需虚拟流体状态信息大为减 少，提高了流场计算的精确性。

附图说明

图1是斜率限制器的构造。

图2是在标志点处构造黎曼问题。

图3是数值结果(算例5.1)。

图4是计算域(算例5.2)。

图5是不同时刻的流场密度(算例5.2)。

图6是特征点运动图(算例5.2)。

图7是相对质量误差(算例5.2)。

图8是计算域(算例5.3)。

图9是不同时刻的密度云图(算例5.3)。

图10是特征点运动图(算例5.3)。

图11是SF₆介质相对质量误差(算例5.3)。

图12是不同时刻的密度云图(算例5.4)。

图13是不同时刻密度(左)和压力(右)对比结果(算例5.4)。

图14是空气介质相对质量误差(算例5.4)。

具体实施方式

本发明在非结构网格上结合RKDG方法与界面追踪方法并且通过RGFM定义界面边 界条件来求解多介质问题。除了非结构网格对于计算区域边界具有较强的适应能力外，由 于RKDG方法具有较好的空间紧致性，远离界面处的虚拟流体状态没有参与计算。通过与基 于有限体积格式得到的结果进行对比，本发明得到守恒误差较小，数值稳定性更好。

下面结合附图和算例对发明内容作进一步说明：

(1)非结构RKDG方法

控制方程为：

$\frac{\partial \overrightarrow{U}}{\partial t}+▿·\overrightarrow{F}\left(\overrightarrow{U}\right)=0---\left(1\right)$

其中$\overrightarrow{U}={[\rho ,\rho u,\rho v,E]}^T,\overrightarrow{F}\left(\overrightarrow{U}\right)=[{\overrightarrow{F}}_1\left(\overrightarrow{U}\right),\overrightarrow{F_2}\left(\overrightarrow{U}\right)],{\overrightarrow{F}}_1\left(\overrightarrow{U}\right)={[\rho u,{\mathrm{\rho u}}^2+p,\rho uv,(E+p)u]}^T,$${\overrightarrow{F}}_2\left(\overrightarrow{U}\right)={[\rho v,\rho uv,{\mathrm{\rho v}}^2+p,(E+p)v]}^T.$ρ是密度,u和v是速度,p是压力，E是单位体积的总 能量，其定义为：

E＝ρe+ρ(u²+v²)/2(2)

其中e是单位质量的内能。状态方程为：

p＝(γ-1)ρe-γB(3)

γ和B是流体常量。对于理想气体γ是比热比，B取为0。

对于非结构网格K₀，数值解和基函数空间由给定(其 中P^k(K₀)是网格K₀上多项式次数小于或等于k的多项式集合)。数值解可以表示为：

${\overrightarrow{U}}_{K_0}^h(x,y,t)=\underset{i=1}{\overset{N}{\Sigma }}{\overrightarrow{U}}_{K_0}^{\left(i\right)}\left(t\right){\phi }_i(x,y),\forall (x,y)\in K_0---\left(4\right)$

对于P¹的情形，自由度是网格各条边中点的值，基函数φ_i(x,y)是线性函 数，其在第i条边的中点上取值为1，在其他边的中点取值为0。对于P²的情形，自由度是网格顶点和各条边中点的值，基函数φ_i(x,y)是二次函数，其在上述第i个点上取值为1， 在其他点取值为0。把式(4)代入式(1),对式(1)乘以基函数φ_i(x,y)并在网格K₀上分部积 分：

$M\frac{{\mathrm{dU}}_{K_0}\left(t\right)}{dt}=-\underset{e\in \partial K_0}{\Sigma }{\int }_e\overrightarrow{F}\left({\overrightarrow{U}}^h\right(x,y,t\left)\right)·{\overrightarrow{n}}_{e,K_0}\phi (x,y)d\Gamma +{\int }_{K_0}\overrightarrow{F}\left({\overrightarrow{U}}^h\right(x,y,t\left)\right)·grad\phi (x,y)dxdy---\left(5\right)$

其中$M_{N\times N}={\left({\int }_{K_0}{\phi }_i\right(x,y\left){\phi }_j\right(x,y\left)dxdy\right)}_{ij}$是质量矩阵，$U_{K_0}\left(t\right)={[{\overrightarrow{U}}_{K_0}^{\left(1\right)},{\overrightarrow{U}}_{K_0}^{\left(2\right)},...,{\overrightarrow{U}}_{K_0}^{\left(N\right)}]}^T,$φ(x,y)＝[φ₁(x,y),φ₂(x,y),…,φ_N(x,y)]^T，是网格K₀上边e的单位法向量。 采用Lax-Friedrichs通量近似。右端积分项采用高斯求积法求解。 式(5)的半离散形式表示为：

${\overrightarrow{U}}_t=L\left(\overrightarrow{U}\right)---\left(6\right)$

采用3阶TVDRK方法进行时间离散：

${\overrightarrow{U}}^{\left(1\right)}={\overrightarrow{U}}^n+\Delta tL\left({\overrightarrow{U}}^n\right)$

${\overrightarrow{U}}^{\left(2\right)}=\frac{3}{4}{\overrightarrow{U}}^n+\frac{1}{4}{\overrightarrow{U}}^{\left(1\right)}+\frac{1}{4}\Delta tL\left({\overrightarrow{U}}^{\left(1\right)}\right)---\left(7\right)$

${\overrightarrow{U}}^{n+1}=\frac{1}{3}{\overrightarrow{U}}^n+\frac{2}{3}{\overrightarrow{U}}^{\left(2\right)}+\frac{2}{3}\Delta tL\left({\overrightarrow{U}}^{\left(2\right)}\right)$

当流场出现间断时应用斜率限制器克服可能出现的虚假振荡问题。如图1所示， K₁，K₂和K₃是与网格K₀相邻的三个网格，m₁，m₂和m₃是网格K₀上三条边的中点,b₀，b₁，b₂和b₃分 别表示网格的形心。由几何分解得到：

$x_{m_1}-x_{b_0}={\alpha }_1(x_{b_1}-x_{b_0})+{\alpha }_2(x_{b_2}-x_{b_0})$

$y_{m_1}-y_{b_0}={\alpha }_1(y_{b_1}-y_{b_0})+{\alpha }_2(y_{b_0}-y_{h_0})$

其中α₁和α₂为非负实数。网格的平均值为：

${\bar{U}}_{K_i}=\frac{1}{K_i}{\int }_{K_i}{\overrightarrow{U}}^{h}dxdy,i=0,1,2,3$

记：

${\overline{\stackrel{‾}{U}}}^h(m_1,K_0)\equiv {\overrightarrow{U}}_{K_0}^h\left(m_1\right)-{\bar{U}}_{K_0}---\left(8\right)$

$\Delta \bar{U}(m_1,K_0)\equiv {\alpha }_1({\bar{U}}_{K_1}-{\bar{U}}_{K_0})+{\alpha }_2({\bar{U}}_{K_2}-{\bar{U}}_{K_0})---\left(9\right)$

其他两个边中点上的和可以类似给出。对于分段线性函数可以得到：

${\overrightarrow{U}}_{K_0}^h(x,y,t)=\underset{i=1}{\overset{3}{\Sigma }}{\overrightarrow{U}}_{K_0}^h\left(m_i\right){\phi }_i(x,y)={\bar{U}}_{K_0}+\underset{i=1}{\overset{3}{\Sigma }}{\overline{\stackrel{‾}{U}}}^h(m_i,K_0){\phi }_i(x,y),\forall (x,y)\in K_0---\left(10\right)$

限制器的具体构造如下，记：

${\Delta }_i=\bar{m}\left({\overline{\stackrel{‾}{U}}}^h\right(m_i,K_0),v\Delta \bar{U}(m_i,K_0\left)\right)---\left(11\right)$

其中ν＝1.5，是修正形式的minmod函数：

其中M是限制器常量。

如果$\underset{i=1}{\overset{3}{\Sigma }}{\Delta }_i=0,$则有：

${\overrightarrow{U}}_{K_0}^h(x,y,t)={\bar{U}}_{K_0}+\underset{i=1}{\overset{3}{\Sigma }}{\Delta }_i{\phi }_i(x,y)---\left(13\right)$

否则，令：

$pos=\underset{i=1}{\overset{3}{\Sigma }}max(0,{\Delta }_i),neg=\underset{i=1}{\overset{3}{\Sigma }}max(0,-{\Delta }_i)$

${\theta }^{+}=min(1,\frac{neg}{pos}),{\theta }^{-}=min(1,\frac{pos}{neg})$

取：

${\bar{\Delta }}_i={\theta }^{+}max(0,{\Delta }_i)-{\theta }^{-}max(0,-{\Delta }_i)$

则有：

${\overrightarrow{U}}_{K_0}^h(x,y,t)={\bar{U}}_{K_0}+\underset{i=1}{\overset{3}{\Sigma }}{\bar{\Delta }}_i{\phi }_i(x,y)---\left(14\right)$

(2)界面追踪方法

如图2所示，在tⁿ时刻介质1和介质2被界面分开。标记点是界面和网格线的交点。 和分别是标记点的法向量和切向量。是标记点P(x_P,y_P)的单位法向 量。点A(x_A,y_A)和点B(x_B,y_B)在不同介质中,与标记点P相距为h(h为网格外接圆的最大直 径)。这里：

x_A＝x_P+h·N_Px,y_A＝y_P+h·N_Py

x_B＝x_P-h·N_Px,y_B＝y_P-h·N_Py

点A和点B的状态由式(4)直接得到：

$\left(\begin{array}{c}{\overrightarrow{U}}_A=\underset{i=1}{\overset{N}{\Sigma }}{\overrightarrow{U}}_{K_A}^{\left(i\right)}\left(t\right){\phi }_i(x_A,y_A)\\ {\overrightarrow{U}}_B=\underset{i=1}{\overset{N}{\Sigma }}{\overrightarrow{U}}_{K_B}^{\left(i\right)}\left(t\right){\phi }_i(x_B,y_B)\end{array}\right)---\left(15\right)$

其中K_A和K_B分别是包含点A和点B的网格。计算点A和点B的密度、法向速度和压力， 分别表示为和沿着标志点P的法线方向构造黎 曼问题，其初始条件为：

${\overrightarrow{W}}_0=\left(\begin{array}{c}{\overrightarrow{W}}_B,\\ {\overrightarrow{W}}_A.\end{array}\right)---\left(16\right)$

采用双激波近似方法求解上述黎曼问题并记其解为其中 上标I表示界面，下标L和R表示界面左右两侧。标志点P的切向速度定义为：

$v_T^I=\left(\begin{array}{cc}{v}_T^B,& u_N^I\ge 0,\\ v_T^A,& u_N^I<0,\end{array}\right)---\left(17\right)$

其中和是标志点P的法向速度和切向速度，和是点A和点B的切向速度。

当求出所有标志点的速度后，采用三阶TVDRK方法获得标志点新的位置：

${\overrightarrow{x}}_f^{\left(1\right)}={\overrightarrow{x}}_f^n+\Delta t·{\overrightarrow{v}}_f\left({\overrightarrow{x}}_f^n\right),$

${\overrightarrow{x}}_f^{\left(2\right)}=\frac{3}{4}{\overrightarrow{x}}_f^n+\frac{1}{4}{\overrightarrow{x}}_f^{\left(1\right)}+\frac{1}{4}\Delta t·{\overrightarrow{v}}_f\left({\overrightarrow{x}}_f^{\left(1\right)}\right),---\left(18\right)$

${\overrightarrow{x}}_f^{n+1}=\frac{1}{3}{\overrightarrow{x}}_f^n+\frac{2}{3}{\overrightarrow{x}}_f^{\left(2\right)}+\frac{2}{3}\Delta t·{\overrightarrow{v}}_f\left({\overrightarrow{x}}_f^{\left(2\right)}\right),$

其中和是标志点在tⁿ和tⁿ⁺¹时刻的位置，是标志点的速度，Δt是时间步 长。

(3)确定虚拟流体区域

由于各个介质单独进行求解，因此有必要确定每种介质的虚拟流体区域。如图2所 示，规定所有标志点的法向矢量均由介质1指向介质2。以网格点C为例，首先计算其与周围 标志点的距离。假定标志点P与网格点C的距离最小。计算数量积其中是由点P 指向点C的位置矢量。如果数量积大于0，则网格点C在介质1的虚拟流体区域，否则在介质2 的虚拟流体区域。对于其他网格点可以采用类似的方法。虚拟流体区域上的法向矢量采用 面积加权方法求得。由于RKDG方法的空间紧致性较好，远离界面处的法向矢量没有参与计 算从而不必进行求解。

(4)真实虚拟流体方法

在真实虚拟流体方法中，在界面附近通过构造黎曼问题更新真实流体状态并外推 得到虚拟流体状态。由于在界面追踪部分已经构造标志点处的黎曼问题，其黎曼问题的解 可以直接用于定义虚拟流体状态。如图2所示，点P，Q，R是格点C附近的三个标志点，是网 格C的法向矢量。网格C的状态可由与其法向夹角最小的标志点上的黎曼解进行更新。如果 标志点P是与网格C法向夹角最小的标志点，那么把标记点P处的黎曼解投影到特征空间可 以更新网格C的平均值。其它界面附近真实流体状态更新方法类似。

虚拟流体状态通过求解如下方程得到：

$\frac{\partial \phi }{\partial t}\pm \overrightarrow{N}·▿\phi =0,---\left(19\right)$

其中φ表示密度，法向速度，切向速度和压力。如图2所示，在求解介质1的边界条 件时采用“+”号，否则采用“-”号。采用迭代方法求解该方程，空间上采用一阶迎风格式离 散，时间上采用三阶TVDRK方法进行离散。由于RKDG方法空间紧致性很好，本文仅采用2～3 个网格作为虚拟流体区域。通过把虚拟流体状态投影到基函数空间，得到虚拟网格状态的 平均值，这样即可确定界面边界条件。

(5)新方法验证

算例5.1.激波管问题。本算例用于研究高压空气在水中的膨胀问题。计算域设计 为圆形，其圆心位于(1,1)、直径为2。初始时刻，直径为0.8的圆形气泡放置在计算域的中 心。气泡内外的初始条件为:

外边界均采用无反射边界条件。对计算域进行网格划分，如图3(a)所示，区域内部有4746个 网格，格点数为2443。采用P²RKDG方法，M取为0.1。图3(b)给出了在t＝0.002时刻流场的密 度云图。其中可以清楚地看到界面，激波和膨胀波。膨胀波向着计算域中心移动而激波在向 计算域的边界移动。为了和具有相同空间离散精度的有限体积方法进行对比，这里把RKDG 方法替换为三阶有限体积WENO方法并保持程序的其他部分不变。RKDG方法、WENO方法与界 面追踪方法的结合分别简记为RKDG-FT方法和WENO-FT方法。图3(c)比较了基于RKDG-FT方 法和WENO-FT方法密度沿着y＝1的分布情况。整体看来，密度分布基本一致。然而，在激波和 稀疏波的捕捉上，RKDG-FT方法更具优势。可以明显地看出基于RKDG-FT方法得到的激波区 域更窄。为了进一步对两种方法作出比较，图3(d)展示了空气的相对质量误差Me。基于 RKDG-FT方法得到的误差明显较小，这充分展示了RKDG-FT方法在界面处理上的优势。

算例5.2.激波打击氦气泡问题。流场计算域如图4所示，其中参数a＝1,b＝0.5,c ＝2,d＝6.5,e＝0.89。空气中马赫数为1.22的激波打击氦气泡。由于流场关于中心轴线对 称，此处仅仅研究上面半个流场的流动问题。在中轴线上取对称边界条件，左右边界均取为 无反射边界条件，上边界为壁面边界条件。初始条件为：波前空气ρ＝1,u＝0,v＝0,p＝1/ 1.4,γ＝1.4,B＝0,波后空气ρ＝1.3764,u＝-0.3336,v＝0,p＝1.5698/1.4,γ＝1.4,B＝ 0,氦气泡ρ＝0.1819,u＝0,v＝0,p＝1/1.4,γ＝1.648,B＝0。对计算域进行网格划分，区域 内部生成58364个网格，格点数为29604。采用P²RKDG方法，M＝0.1。图5给出了不同时刻流场 密度的变化情况。在初始时刻，当界面被激波打击后，部分界面发生移动。在氦气泡的压缩 过程中，界面的长度变小并且一些标志点逐渐消失。氦气泡以逆时针方向卷起来，导致气泡 后面形成较大的压力从而推动气泡向左移动。激波在打击氦气泡之后分成两部分，其中氦 气泡中的激波具有更高的强度而且运动速度也比入射激波更快。为了和前人的研究结果进 行对比，图6展示了三个特征点(Jet,Downstream,Upstream)的运动图以及Terashima等人 的计算结果。可以看出这些结果基本一致。为了进一步和WENO-FT方法做出对比，此处分别 计算了氦气泡的相对质量误差，如图7所示。可以看出，在氦气泡破裂之前两种方法的相对 质量误差均小于7％。总体上基于RKDG-FT方法的误差较小。

算例5.3.气-气Richtmyer-Meshkov不稳定性问题。如图8所示，流场计算域为[0, 4]×[0,0.5]，位于x＝3.2、马赫数为1.24的激波向左运动打击气气界面。界面初始形状为： x＝2.9-0.1sin(2π(y+0.25)),0＜y＜0.5。流场初始条件为：SF₆ρ＝5.04,u＝0,v＝0,p＝1, γ＝1.093,B＝0,波前空气ρ＝1,u＝0,v＝0,p＝1,γ＝1.4,B＝0,波后空气ρ＝1.411,u＝- 0.39,v＝0,p＝1.628,γ＝1.4,B＝0。左右边界取无反射边界条件，上下边界为对称边界条 件。对计算域进行网格划分，区域内部生成105542个网格，格点数为53357。采用P²RKDG方 法，M＝0.1。不同时刻的密度云图如图9所示，可以清晰地看到激波在SF₆介质中的运动过程 以及界面的变化。作为验证，图10展示了特征点(Spike,Bubble)的运动图，同时也给出了 Terashima等人的计算结果作为对比。可以看出，这些运动曲线基本一致。和算例5.2类似， 图11展示基于RKDG-FT方法和WENO-FT方法的SF₆介质的相对质量误差。可以看出，在初始时 刻，误差相差不大。然而随着激波在SF₆介质中不断运动，界面形状变得复杂，基于WENO-FT 方法的误差增加较快，而基于RKDG-FT方法的误差变化平缓。

算例5.4.气-液Richtmyer-Meshkov不稳定性问题。采用和算例5.3相同的流场计 算域、界面初始形状和网格划分。液体中位于x＝3.025、马赫数为1.95的激波向左运动打击 气液界面。流场初始条件为：空气ρ＝1,u＝0,v＝0,p＝1,γ＝1.4,B＝0,波前液体ρ＝5,u＝ 0,v＝0,p＝1,γ＝4,B＝1,波后液体ρ＝7.093,u＝-0.7288,v＝0,p＝10,γ＝4,B＝1。采用 P²RKDG方法，M＝0.1。图12给出了流场中密度的变化情况。从图中可以清晰地看出激波打击 界面后的波系结构。为了验证程序的正确性，图13给出了基于γ模型、沿着y＝0.5的密度、 压力分布进行对比。可以看出结果基本一致，证明了RKDG-FT方法的准确性。为了进一步展 示RKDG-FT方法的优越性，图14给出了空气介质的相对质量误差随时间变化情况。可以看出 在初始阶段，误差基本一致。在激波进入空气介质后，随着时间的增大，基于WENO-FT方法的 误差增长迅速，而基于RKDG-FT方法得到的误差较小。