CN107784141B - 一种二维有限控制体积计算的加速方法 - Google Patents

Abstract

本发明公开了一种二维有限控制体积计算的加速方法，具有方法简单，易于移植，稳定性好、运行效率提高等显著优势。本发明人通过近几年来模拟实践，在改进无结构网格有限体积显式算法方面取得了一些进展，着重对优化网格设置、双时间轴变时步长法等进行讨论。本发明的加速方法步骤如下：第一步：优化网格设置，第二步：双时间步长推进法。

Description

一种二维有限控制体积计算的加速方法

技术领域

本发明涉及一种二维有限控制体积计算的加速方法，属于无结构网格有限体积计算领域。

背景技术

近年来，泥沙数学模型研究已取得了长足的进展，广泛应用于水库、河道、湖泊以及河口的河床冲淤变形预测。为了提高预测精度，除了致力于河床变形机理和计算模式的研究外，水沙运动方程组的数值离散算法的改进和创新亦是不可缺少的重要环节。有限差分法因其简单实用，已成为水沙数值模拟的主要算法，但在实际模拟工作中总感到不够满意，主要表现在：(1)有限差分法除非矩形网格很细，否则难以准确描述计算域边界和地形，从而导致局部流场和泥沙浓度场模拟精度欠佳；(2)由于有限差分法属于节点算法，难以准确实现计算单元以及整个计算域上水量、动量以及泥沙量的平衡，影响了水沙运动数值模拟的计算精度和数值稳定性；(3)边界条件处理繁复或不准确，内点与边界点格式不一致。为了克服上述困难，近年来重点开发研究无结构网格的有限体积法水沙数值模拟算法。

但由于构建无结构网格有限体积法的隐式算法的复杂性和庞大稀疏矩阵的迭代不稳定性，目前大多数采用无结构网格有限体积显式算法，由于显式算法的时间步长受限制，在进行长历时计算时效率不高，已成为制约该算法广泛运用的主要瓶颈。发明人通过近几年来模拟实践，在改进无结构网格有限体积显式算法方面取得了一些进展，提出优化网格设置、双时间轴变时步长法等新方法。

本发明旨在提供一种二维有限控制体积算法的新改进，具有方法简单，易于移植，稳定性好、运行效率提高等显著优势。

发明内容

基于背景技术存在的技术问题，本发明针对背景技术存在的问题，提供一种二维有限控制体积算法的新改进，具有方法简单，易于移植，稳定性好、运行效率提高等显著优势。

在二维浅水流动计算中，常遇到如何处理复杂边界形状以及计算区域内有堤防、公路、铁路等天然分界这类问题。采用任意三角形或多边形网格剖分是非常合适的，既可以克服矩形网格锯齿形边界所造成的流动失真，也可以避免生成有结构贴体曲线网格的复杂计算和其它困难。因此，为了更好地拟合湖区的复杂岸形和地形，研制了由三角形和四边形单元混合组成的网格单元体系。在无结构三角形或四边形单元中，采用有限体积法对水流微分方程组进行数值离散，其实质是逐单元进行水量和动量平衡，物理意义清晰，准确满足积分形式的守恒律，成果无守恒性误差，且能够处理含间断或陡梯度的流动。所以，基于无结构网格有限体积法的算法的优点是显而易见的。

但由于构建无结构网格有限体积法的隐式算法的复杂性和庞大稀疏矩阵的迭代不稳定性，目前大多数采用无结构网格有限体积显式算法，由于显式算法的时间步长受限制，在进行长历时计算时效率不高，已成为制约该算法广泛运用的主要瓶颈。本发明人通过近几年来模拟实践，在改进无结构网格有限体积显式算法方面取得了一些进展，着重对优化网格设置、双时间轴变时步长法等进行讨论。

为了提高非定常流场的计算效率，20世纪90年代中期，国外在ADI、LU分解等定常高效隐式算法的基础上，通过“亚迭代”技术，在空气动力学领域发展了一种所谓“双时间推进法”的非定常流场计算方法。该算法具有结构简单、易于推广和计算效率高等特点。

在显式有限控制体积算法的基础上，提高计算效率的方法主要有以下几个方面：①优化网格布置，按照各单元的CFL接近1来划分网格，尽可能采用均匀四边形网格；②双时间步长推进法，在显式算法的基础上，引入空气动力学领域一种所谓“双时间推进法”。采用“亚迭代”技术，即在选择合适的物理时间步长的基础上，运用虚拟迭代时间步长，加速方程组数值解的收敛过程，达到提高计算效率目的。

二维有线控制体积算法有：

1二维浅水方程组

为了保证格式的守恒性，以及适用于含间断或陡梯度的流动，采用二维不恒定浅水方程组的守恒形式[4]：

其中守恒物理量W、x向和y向通量向量F和G，以及源项向量D分别为：

式中h为水深、u和v分别是x和y方向垂线平均的水平流速分量，g为重力加速度；

和

分别为x和y方向的水底底坡；

和

分别为x和y方向的摩阻坡度，采用二维形式的曼宁阻力公式计算；q为单元旁侧入流。

2有限体积离散

采用有限体积法对水沙微分方程组进行数值离散，其实质是逐单元进行水量、动量和沙量平衡，物理意义清晰，准确满足积分形式的守恒律，成果无守恒性误差，且能够处理含间断或陡梯度的流动。

为了适应复杂几何形状流场和泥沙浓度场的数值计算，采用任意三角形或四边形网格剖分和网格中心格式，即将流动变量定义在单元形心，且控制体与单元本身重合。记Ωi为单元域边界，利用格林公式，可得方程组(1)的有限体积近似[4]：

式中Ai为第i单元Ωi的面积；(cosφ，sinφ)为Ωi的外法向单元向量；d1为线积分微元；

为非齐次项在单元Ωi上的某种平均；若记F_n＝F·cosφ+G·sinφ为跨单元界面的法向水流数值通量，时间积分采用显向前格式，那么(4)式为

式中求和号下的指标j表示某一单元i的第j边，lj为单元界面边长，上标n表示时间，算法的核心是如何计算法向数值通量Fn。

3法向数值通量计算

水沙法向数值通量计算的关键是在无结构网格上进行跨单元界面法向数值通量沿特征的逆风分解。对于二维水流数值计算，大部分单元采用Osher格式计算水量和动量的法向通量；对于洪水漫堤和分洪的特殊单元，则采用水力衔接方程计算法向通量。

本发明的目的通过以下技术方案实现：

一种二维有限控制体积计算的加速方法，其步骤如下：

第一步：优化网格设置

尽量少用三角型网格，按照CFL接近1时的Δt的计算公式：

单元为三角形时，

单元为四边形时，

其中，Δr为三角形内接园半径；ΔX、ΔY分别是四边形的边长，显而易见，相同面积的四边形的时间步长几乎是三角形的时间步长的1倍；

第二步：双时间步长推进法

在空气动力学领域发展了一种所谓“双时间推进法”的非定常流场计算方法。即求解非恒定方程组时，采用物理时间和虚拟时间双步长方法。物理时间步长的选取取决于精度要求，虚拟时间步长参与方程求解的迭代过程，也是提高计算效率的关键环节。将求解空气动力学方程的双时间推进法，移植到浅水动力学方程数值计算中。双时间步长推进法的核心技术是基于虚拟时间步长的亚迭代技术，核心问题是尽可能加速获得的收敛解，本发明提出局部时间步长法和残差消除加速收敛法。

(1)局部时间步长加速收敛法

局部时间步长加速收敛法就是采用局部单元时间步长法。在使用显格式时，一般采用计算单元最小时间步长，对于大多数计算网格来说，计算时步偏小、效率低，且单元CFL值小于1的计算精度小于CFL接近1时的精度。为此，根据有限控制体法的离散方程，即U^n＝1＝Uⁿ+ΔtR(Uⁿ)可知，各控制体单元允许采用各自单元的时步长。考虑到各单元时步长不一致会对时间过程积分产生误差，但对于收敛解没有影响。具体来说，具体来说，计算四边形网格的时步长Δt，计算三角形网格的时步长Δt计算；

按照柯兰德数，可令

假设β代替c_r，则通过试验确定β，即逐步缩小β直至计算稳定。

此方法的思路是使扰动波尽快离开各控制体，利用此方法可使得水深浅流速小的地方收敛加速较多，数值实验表明仅这一项改进，实际物理时间步长可以增加一倍左右；

(2)残差消除加速收敛法

按照隐式残差光滑法和多重网格插值消除残差方法的思路，在双时间步长推进法的运用中引入常规的消除残差的计算方法。在虚拟时间迭代过程中，按照收敛解的准则，自动判断当前数值解是否接近收敛解，通常水位计算精度较高，宜以流速误差作为指标，即前后时步流速最大误差满足小于某一设定值。将这判别准则延伸至进出单元法向通量之和的变化值为某一设定小值。将进出单元法向通量之和的变化值与设定小值之差作为残差。在空间上仿照多重网格插值消除残差方法的思路，可以在相邻网格单元进行平均残差，在时间上采用前后多个时步平均消除残差，具体残差消除加速收敛法如下：

按照多重网格插值消除残差方法的思路，采用简化空气动力学多重网格插值消除残差方法，即采用粗细两重网格单元消除残差，具体来说，基于生成的计算网格，以计算单元为中心，将与其连接的周边单元共同组成一个粗网格，将计算单元的残差，在该粗网格内进行平均，即残差局部光滑，依次类推，进行一次空间分布上的残差消除；在时间上，采用前后多个时步的残差修正。

数值试验表明，采用上述简化的残差消除方法，没有增加CPU的耗时，反而增加了数值计算的稳定性和加快了计算收敛，使得物理时步长增加60％。

本发明的有益之处在于：

1.本发明方法简单，易于移植；

2.本发明稳定性好、运行效率大幅提高；

3.法向数值通量的计算误差极低，误差率低于0.05％。

具体实施方式

实施例：

一种二维有限控制体积计算的加速方法，其步骤如下：

第一步：优化网格设置

尽量少用三角型网格，按照CFL接近1时的Δt的计算公式：

单元为三角形时，

单元为四边形时，

其中，Δr为三角形内接园半径；ΔX、ΔY分别是四边形的边长，显而易见，相同面积的四边形的时间步长几乎是三角形的时间步长的1倍；

第二步：双时间步长推进法

在空气动力学领域发展了一种所谓“双时间推进法”的非定常流场计算方法。即求解非恒定方程组时，采用物理时间和虚拟时间双步长方法。物理时间步长的选取取决于精度要求，虚拟时间步长参与方程求解的迭代过程，也是提高计算效率的关键环节。将求解空气动力学方程的双时间推进法，移植到浅水动力学方程数值计算中。双时间步长推进法的核心技术是基于虚拟时间步长的亚迭代技术，核心问题是尽可能加速获得的收敛解，本发明提出局部时间步长法和残差消除加速收敛法。

(1)局部时间步长加速收敛法

局部时间步长加速收敛法就是采用局部单元时间步长法。在使用显格式时，一般采用计算单元最小时间步长，对于大多数计算网格来说，计算时步偏小、效率低，且单元CFL值小于1的计算精度小于CFL接近1时的精度。为此，根据有限控制体法的离散方程，即U^n＝1＝Uⁿ+ΔtR(Uⁿ)可知，各控制体单元允许采用各自单元的时步长。考虑到各单元时步长不一致会对时间过程积分产生误差，但对于收敛解没有影响。具体来说，四边形网格的时步长Δt计算公式(4)式，三角形网格的时步长Δt计算、公式见(5)式。按照柯兰德数，可令

假设β代替c_r，则通过试验确定β，即逐步缩小β直至计算稳定。

此方法的思路是使扰动波尽快离开各控制体，利用此方法可使得水深浅流速小的地方收敛加速较多，数值实验表明仅这一项改进，实际物理时间步长可以增加一倍左右；

(2)残差消除加速收敛法

按照隐式残差光滑法和多重网格插值消除残差方法的思路，在双时间步长推进法的运用中引入常规的消除残差的计算方法。在虚拟时间迭代过程中，按照收敛解的准则，自动判断当前数值解是否接近收敛解，通常水位计算精度较高，宜以流速误差作为指标，即前后时步流速最大误差满足小于某一设定值。将这判别准则延伸至进出单元法向通量之和的变化值为某一设定小值。将进出单元法向通量之和的变化值与设定小值之差作为残差。在空间上仿照多重网格插值消除残差方法的思路，可以在相邻网格单元进行平均残差，在时间上采用前后多个时步平均消除残差，具体残差消除加速收敛法如下：

按照多重网格插值消除残差方法的思路，采用简化空气动力学多重网格插值消除残差方法，即采用粗细两重网格单元消除残差，具体来说，基于生成的计算网格，以该计算单元为中心，将与其连接的周边单元共同组成一个粗网格，将计算单元的残差，在该粗网格内进行平均，即残差局部光滑，依次类推，进行一次空间分布上的残差消除；在时间上，采用前后多个时步的残差修正。

数值试验表明，采用上述简化的残差消除方法，没有增加CPU的耗时，反而增加了数值计算的稳定性和加快了计算收敛，使得物理时步长增加60％以上，而且法向数值通量的计算误差极低，误差率低于0.05％。

由此可见，本发明的算法简单，易于移植；而且稳定性好、运行效率大幅提高。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种用于水沙数值模拟的二维有限控制体积计算的加速方法，其步骤如下：

第一步：优化网格设置

少用三角形网格，按照CFL接近1时的Δt的计算公式：

单元为三角形时，

单元为四边形时，

其中，Δr为三角形内接圆半径；ΔX、ΔY分别是四边形的边长，相同面积的四边形的时间步长是三角形的时间步长的1倍；

第二步：双时间步长推进法

(1)局部时间步长加速收敛法

根据有限控制体积计算方法的离散方程，即Uⁿ⁺¹＝Uⁿ+ΔtR(Uⁿ)可知，各控制体单元允许采用各自单元的时间步长，计算四边形网格的时间步长Δt，计算三角形网格的时间步长Δt；

按照柯兰德数，可令

假设β代替cr，则通过试验确定β，即逐步缩小β直至计算稳定；

(2)残差消除加速收敛法：

按照隐式残差光滑法和多重网格插值消除残差方法的思路，在双时间步长推进法的运用中引入常规的消除残差的计算方法；

在虚拟时间迭代过程中，按照收敛解的准则，自动判断当前数值解是否接近收敛解，以流速误差作为指标；

将这判别准则延伸至进出单元法向通量之和的变化值为某一设定小值；将进出单元法向通量之和的变化值与设定小值之差作为残差；

在空间上仿照多重网格插值消除残差方法的思路，可以在相邻网格单元进行平均残差，在时间上采用前后多个时间步长平均消除残差；

第二步，步骤(2)中具体残差消除加速收敛法如下：

按照多重网格插值消除残差方法的思路，采用简化空气动力学多重网格插值消除残差方法，基于生成的计算网格，以该计算网格为中心，将与其连接的周边单元共同组成一个粗网格，将计算单元的残差，在该粗网格内进行平均，即残差局部光滑，依次类推，进行一次空间分布上的残差消除；在时间上，采用前后多个时间步长的残差修正；

没有增加CPU的耗时，反而增加了数值计算的稳定性和加快了计算收敛，使得物理时间步长增加60％。