TVD格式求解对流扩散方程的最优限制器研究

2023-03-29戚园春侯庆志

计算机仿真 2023年2期

戚园春，刘昉，侯庆志，2

(1. 天津大学建筑工程学院，天津 300354；2. 天津大学智能与计算学部，天津 300354)

1 引言

对流扩散方程作为描述自然界流动、传热及物质输运过程的主要控制方程，其可充分反映流体在流动过程中所携带的物理量的变化规律，在流体力学、能源开发及环境科学等领域发挥着重要的作用[1]。

对流扩散方程的求解方法主要有有限差分法、有限元法和有限体积法[3-6]。有限体积法物理意义明确、计算效率较高，基于该方法求解对流扩散方程时，常用的离散格式有一阶迎风、二阶迎风、中心差分、Lax-Friedrichs、Lax-Wendroff、Beam-Warming、Fromm等，然而这些格式会导致计算结果出现不同程度的数值耗散或者数值振荡[7]，精度过低会抹平数值解中的间断，精度过高又会导致计算结果在间断附近出现伪振荡。近年来，高分辨率格式的发展有效缓和了这一问题，其可敏锐地捕捉间断，提高激波的分辨率并消除伪振荡。常用的高分辨率离散方法有FCT (Flux Corrected Transport)方法和TVD (Total Variation Diminishing)方法。

FCT方法是几乎所有保单调和非振荡流体输运算法的起源，该方法的本质是在输运格式里面施加一个扩散和反扩散项，使之能平滑激波区和流场陡变区，从而避免伪振荡。而TVD方法是通过保证变量的变化不随时间增加，从而消除伪振荡，其在对流扩散方程求解中表现出了良好的性能。早期，Harten[8]依据TVD性质构造了一个通量限制函数，并与高阶格式相结合，首次提出了“高分辨率方法”这一概念。Tamamidis等[9]采用3种高阶单调格式分别求解了4类二维对流扩散问题，并与无通量限制器格式的计算结果进行对比，表明带通量限制器的高阶格式可显著提高计算结果的准确性。Hutter等[10]采用多种数值格式求解一维对流扩散问题，表明带有Superbee限制器的MTVDLF方法是求解对流占优问题的最佳方法。郭建红等[11]基于有限体积通量修正，应用高阶迎风对流格式和带TVD限制器的高分辨率格式研究了自然空泡湍流流动的数值计算方法，发现后者能更好地捕捉空泡界面附近物理量的阶跃特性。冯定华等[12]应用8种差分格式和5种限制器计算分析了激波管Riemanna问题，并对比了各类限制器的压缩性和耗散性。Zhang等[13]对现有的TVD格式进行归类，表明合理的通量限制器是TVD格式的关键要素，并在此基础上提出了一种改进的用于稳态计算的SS-TVD(Steady-state TVD)限制器TCDF。

因此，为进一步优选求解不同类型对流扩散问题的最优算法，本文共采用了3种典型的TVD格式：MUSCL、TVDLF、MTVDLF格式，10种常用的通量限制器：Superbee[14]、Minmod[15]、Woodward[10]、Van Leer[9]、UMIST[16]、WACEB[17]、Albada[18]、OSPRE[19]、TCDF[13]、Koren[13]，求解了阶跃型纯对流问题、阶跃型对流扩散问题、线性高斯型对流扩散问题、拟线性高斯型对流扩散问题、Burgers方程，共5类典型的对流扩散问题。同时，为了促进该问题的数值模拟研究发展，本文所有格式的相关代码可在以下网站自行下载：http：∥114.55.218.152/wp-blog/，以期为国产自主软件开发提供一定的技术支持。

2 数学模型

对流扩散方程的一维形式和守恒形式为

(1)

(2)

3 TVD方法与通量限制器

有限体积法守恒性较好，物理意义明确，以下主要采用该方法离散方程(2)中的导数，通过在[xj-1/2，xj+1/2]×[tn，tn+1]上积分，可以得到该方程的离散形式为

(3)

(4)

(5)

(6)

(7)

(8)

(9)

φj=φ(θj)

(10)

(11)

式中，σj为斜率限制器，θj为连续梯度比，φ(θ)函数的选用见表1。

表1 TVD格式中的常用通量限制器

本文主要对比带有以上通量限制器的三种TVD格式：MUSCL、TVDLF和MTVDLF格式(见表2)，并通过数值实验找出适应于不同类型问题的最优格式。由于在对流扩散问题中，扩散项对不同的计算格式并不敏感，一般采用中心差分格式处理就可以满足实际需要。因此，以下仅采用TVD格式构造对流通量，三种格式中对流通量的构造方法如下

(12)

表2 三种TVD格式

可以看出，当对流速度为常数时，MUSCL格式和MTVDLF格式本质上是相同的。

4 数值实验

针对以下五个典型的对流扩散问题，分别开展数值实验。为了定量地比较这些数值格式的数值精度，进行了误差分析，误差定义为

(13)

4.1 线性对流扩散问题

(14)

(15)

(16)

将以上三式代入方程(1)，可得线性对流扩散方程的无量纲形式

(17)

4.1.1 阶跃型纯对流问题

该问题的初始浓度分布为

(18)

式中x0=0.1，为初始间断位置，精确解为

(19)

模拟时间T=0.5，Δx=0.005，Δt=0.000005。分别采用MUSCL格式和TVDLF格式进行计算，数值结果如表3所示。

表3 算例1计算误差

由上表可以看出，在网格数相同的情况下，选用同种通量限制器时，TVDLF格式的计算误差远大于MUSCL格式。且在两种TVD格式中，Superbee限制器的误差最小，而Minmod误差最大。两种格式的计算结果如图1和图2。

图1 带有Superbee和Minmod限制器的MUSCL格式计算结果

图2 带有Superbee和Woodward限制器的TVDLF格式

可以看出，基于TVD格式的数值解中没有出现明显的震荡，只存在小部分激波被抹掉的现象，而基于TVDLF格式的计算结果出现了严重的数值耗散。因此，对于阶跃型纯对流问题，应优先选用带有Superbee限制器的MUSCL格式进行求解。

4.1.2 阶跃型对流扩散问题

该问题的初始浓度分布及精确解为

c(x，0)=0， 0

(20)

(21)

该问题中模拟时间T=1，Δx=0.01，Δt=0.0001，Pe=100。分别采用MUSCL和TVDLF格式计算，发现TVDLF格式的计算结果仍会出现不可接受的数值耗散，故不再考虑该格式。MUSCL格式的计算结果见表4。

表4 算例2计算误差

可以看出，基于MUSCL格式的计算结果与精确解高度吻合，计算误差均在0.8%以下，该问题中限制器的类型基本不影响计算结果，Minmod限制器的计算结果如图3。

图3 带有Minmod限制器的MUSCL格式计算结果

4.1.3 高斯型对流扩散问题

该问题中高斯分布的初始位置为x0=1/4，σ=1/30，Pe=200，初始浓度分布及精确解为：

(22)

(23)

该问题中模拟时间T=0.4，Δx=0.01，Δt=0.00001。分别采用MUSCL格式和TVDLF格式计算，TVDLF格式的计算结果仍会出现不可接受的数值耗散，故不再考虑。MUSCL格式的计算结果见表5。

表5 算例3计算误差

可以看出，基于MUSCL格式的计算误差介于1.2%～3.6%之间，并且 Woodward限制器的计算误差最小，Minmod限制器的误差最大，基于该限制器的计算结果如图4。

图4 带有Woodward限制器的MUSCL格式计算结果

4.2 拟线性对流扩散问题

为了比选出TVD方法求解拟线性对流扩散问题时性能最优的通量限制器，求解如下方程：

(24)

该问题中对流速度a=0.02，扩散系数Γ=0.0001，Pe=200，高斯分布的初始位置x0=1/4，σ=1/30，初始浓度分布及精确解为

(25)

(26)

模拟时间T=10s，Δx=0.01m，Δt=0.0005s。控制体界面两侧对应同一个a(1+x)，故MUSCL和MTVDLF格式相同。基于TVDLF格式的计算结果仍表现出了严重的耗散，故不再考虑。MUSCL格式的数值计算结果见表6。

表6 算例4计算误差

可以看出，选用MUSCL格式计算时，Minmod限制器误差最大，而Woodward限制器误差最小。因此，对于拟线性高斯型对流扩散问题，应优先选用带有Woodward限制器的MUSCL或MTVDLF格式求解。选用Woodward和Minmod限制器的计算结果如图5。

图5 带有Woodward和Minmod限制器的MUSCL格式计算结果

4.3 非线性对流扩散问题-Burgers 方程

Burgers方程是模拟冲击波传播和反射的非线性偏微分方程，为找出TVD方法求解该方程时的最优限制器，设置了本算例。该方程的守恒形式为

(27)

(28)

该问题中计算域为[-1，2]，初始流速分布及精确解为

(29)

(30)

模拟时间T=2s，Δx=0.01m，Δt=0.001s，三种格式的数值计算结果见表7。

表7 算例5计算误差

可以看出，MUSCL格式和MTVDLF格式的计算误差均很小，而TVDLF格式的数值计算误差较大，尤其是采用Albada和OSPRE限制器的情况。因此，对于Burgers方程，应优先选用MUSCL或MTVDLF格式求解，且对于该问题限制器的类型基本不影响计算结果。带有Superbee限制器的三种格式的计算结果如图6。