求解對流擴散問題的積分方程法

魏濤，許明田，汪引

（山東大學土建與水利學院工程力學系，山東 濟南 250061）

求解對流擴散問題的積分方程法

魏濤，許明田，汪引

（山東大學土建與水利學院工程力學系，山東 濟南 250061）

提出了一種求解對流擴散問題的積分方程法。在這一方法中，首先利用Laplace方程的級數形式的格林函數將對流擴散方程轉化為積分方程，然后利用級數的正交性質，把積分方程進一步簡化為代數方程組，求解該方程組即可得到對流擴散方程的級數形式的近似解。最后，分別利用Chebyshev多項式和Fourier級數求解了3個典型的一維和二維對流擴散問題。該方法和有限體積法、有限元法和迎風差分法相比，展現出非常高的精度并且避免了由解的不連續性造成的虛假振蕩。

計算流體力學；積分方程法；對流擴散方程；有限體積法; 傳熱

引 言

對流擴散方程在化學、物理和工程中有著廣泛的應用，如流體力學、傳熱學、化學反應動力學以及環境科學中的許多現象都可以歸結為求解對流擴散方程問題[1-5]。有限元法、有限差分法和有限體積法是3種常用的求解對流擴散方程的數值方法[6-12]。其中，由于有限體積法能保證控制體積內的通量守恒，所以成為求解對流擴散方程的一種重要方法。有限體積法是20世紀60年代提出的一種數值計算方法，它對偏微分方程在每一個小的子區域上進行積分，并且保持物理量在控制容積內的守恒。該方法概念簡單，物理守恒原理和數值算法之間的關系易于理解，因此被廣泛應用于求解流體力學和傳熱學問題。但是結構化網格的高階離散降低了它的計算效率[13-14]，并且解的不連續性經常會引起高階有限體積法的數值振蕩，這就需要在變量重構過程中引入限制器來抑制這個振蕩。van Leer[15]通過在經典MUSCL（monotonic upstream-centered scheme for conservation laws）格式中添加斜率限制器（slope limiter）來處理不連續性，Venkatakrishnan[16-17]通過對限制器的修正得到了求解穩定解的高效收斂方法。使用限制器可以有效消除虛假振蕩，但同時它也有一些不足，如使數值解收斂到穩定解的速度變慢，降低了光滑極值點附近的精度等。

后來，研究者提出了一種ENO（essentially non-oscillatory）格式用于替代MUSCL格式。該方法通過在每一個迭代步選取光滑通量模板，避免斜率限制器的使用，它可以得到高階精度的解。但是，由于這個過程內在的不可微性，導致數值解不能收斂到穩定解，于是，又發展出了一種WENO（weighted essentially non-oscillatory）格式來處理這個問題[18-20]。然而，這些方法收斂到穩定解的效率卻不如MUSCL格式高。

近來，研究者們基于積分方程理論也提出了一些新的數值方法[21-26]。Wen等[21]給出了一個比傳統有限差分法精度高很多的有限積分法。而由積分公式發展出來的非標準有限差分格式（nonstandard finite difference schemes）在求解反應擴散對流方程時，不僅提高了計算精度，而且減少了計算時間[22]。Xu等[23-26]建立了數值模擬磁流體自激發電和磁重聯現象的積分方程法，與其他方法相比較，該方法展現出很高的精度和可靠性。并且積分方程法的計算結果對法國VKS（von Kármán-sodium）磁流體自激發電實驗的成功起到了關鍵作用[24,26]。在本文中，將建立求解對流擴散問題穩態解的積分方程法。

1 數值方法

在本文中，考慮下面的穩態對流擴散問題

其中，φ1n、φ2n、φ3m、φ4m和式（10）、式（11）、式（12）、式（13）中給出的表達式相同。與求解格林函數時類似，式（22）、式（23）中也有4個與其他線性相關。從式（21）中選出(M?1)× (N?1)個

從式（22）和式（23）中選出2M+2N個線性無關的，然后它們組成一個線性代數方程組。求解該方程組即可得到未知系數Tmn，再將這些系數代入式（18）就得到了式（1）滿足邊界條件式（2）的近似解。

2 數值實驗

2.1 一維對流擴散問題

選取正交多項式為Chebyshev多項式Pm(x)=cos(marccosx)。

例1：考慮如下一維穩態對流擴散問題

邊界條件

式中，T是溫度場，℃；u是流體的流動速度，m·s?1；ρ是密度，kg·m?3；Γ是擴散系數，kg·m?1·s?1。取下列量綱1參數

則式（25）、式(26)可以寫成下面的形式

邊界條件

這個問題的解析解為

考慮下面3種情形。

2.1.1 情形1:η=1 對于該種情形，由本文的積分方程法、解析解和有限體積法得到的結果都列在表1中。從該表可以看出當η=1，M=5時，與解析解相比由有限體積法得到的結果絕對誤差為10?3量級，而由積分方程法得到的結果絕對誤差為10?6量級。

2.1.2 情形2:η=1/25 對于該種情形，由積分方程法、解析解和有限體積法得到的結果列在表2中。當ξ趨向于1時，梯度的絕對值越來越大，從該表可以看出，有限體積法的解在梯度大的點處誤差增大，而積分方程法則表現出了較好的穩定性。

表1η=1、M=5時由解析解、有限體積法和積分方程法得到的結果Table 1 Results obtained by analytical method, FVM and integral equation approach （IEA） forη=1, M=5

表2η=1/25、M=20時由解析解、有限體積法和積分方程法得到的結果Table 2 Results obtained by analytical method, FVM and IEA forη=1/25, M=20

2.1.3 情形3:η=1/60 和前兩種情形相比較，本情形考慮的量綱1擴散系數η較小。為了更清晰地展示對應于不同η的對流擴散問題的解的不同變化趨勢，將3種情形下的解析解都畫在圖1中。由圖可見，對于較大的η，場變量Θ的變化比較平穩，但是當η較小時解在邊界ξ=1處的變化非常劇烈，這就是導致有限體積法不穩定的原因[27]。當η=1/60時，由積分方程法、解析解、迎風差分格式及有限元法求得的結果都列在表3和表4中。從這兩個表格中可以看出當M=20時，由積分方程法所得解的誤差比由采用有理基函數（rational basis function）的有限元法大，而當M=40時積分方程法的精度要比采用有理基函數的有限元法高。

圖1 例1中對應不同η的解析解Fig.1 Analytic solution of Example 1 for differentη

表3 當η=1/60、M=20時由各種方法得到的結果Table 3 Results obtained by various methods forη=1/60, M=20

表4 當η=1/60、M=40時由各種方法得到的結果Table 4 Results obtained by various methods forη=1/60, M=40

2.2 二維對流擴散問題

選取的正交多項式為Pm(x)=sin(mπx/Lx)，Pn(y)=sin(nπy/Ly)。

例2：在式（1）中取ν=1 m2·s?1，u1=u2= ?2 m·s?1，S(x,y)= ?4x(y?y2) ?4y(x?x2)，以及邊界條件φ(x,y)=0，則該對流擴散問題的解析解為

其中，(x,y)∈[0,1]×[0,1]。

在式（24）中取M=N=20，相當于采用20×20的網格。用下面的誤差公式計算絕對誤差和相對誤差

表5 例2中由積分方程法所得數值解的收斂性Table 5 Convergence of numerical solution obtained by IEA for Example 2

對流占優的對流擴散方程在工程問題中有著廣泛的應用，有限體積法、有限差分法和有限元法是求解此類問題常用的方法。然而，當熱擴散率比較小時，由這些方法得到的數值結果經常會產生數值擴散或者數值振蕩。在此，通過比較由積分方程法得到的數值解和由式（30）給出的解析解之間的誤差來檢驗該方法的穩定性。取M=N=100，u1=u2=?2 m·s?1。對于不同的熱擴散率數值解和解析解之間的誤差列于表6中。從該表可以看出當熱擴散率變小時誤差有所增加，但數值振蕩并沒有發生。這說明積分方程法在求解對流占優的對流擴散問題時具有很好的穩定性。

表6 積分方程法所得對流占優對流擴散問題數值解的穩定性Table 6 Stability of numerical solution obtained by IEA for convection-dominated diffusion problems

為了進一步驗證該方法，給出了下面的例子。

例3：在本例中，取ν=1 m2·s?1，u1=u2= ?2 m·s?1，源項S(x,y)=ex+y[(6-5π2)sin(πx)sin(2πy)+ 4πcos(πx) sin(2πy)+8πsin(πx)cos(2πy)]，以及邊界條件φ(x,y)=0，則式（1）的解析解為

其中，(x,y)∈[0,1]×[0,1]。

由積分方程法和有限體積法[27]所得的結果列于表7中，所有計算用Fortran語言編寫程序，并在DELL OPTIPLEX 790臺式計算機上完成運算。從該表中可以看出在同等網格數量下用積分方程法進行計算所需的時間略多于有限體積法，但是積分方程法的精度要明顯好于有限體積法。圖3(a)、(b)分別展示了解析解和由積分方程法得到的數值解。通過比較表7和表5發現例3中的絕對誤差εmax比例2中的大，由圖3和圖2可以看到例3中解的結構復雜并且很多點處的梯度也比較大，這應是表7中的絕對誤差比表5大的原因。但是表7和表5中的相對誤差εL2在數量級上并沒有明顯的差別。

表7 例3中由積分方程法和有限體積法得到的結果Table 7 Results obtained by IEA and FVM for Example 3

圖2 例2中當ν=1 m2·s?1、u1=u2= ?2 m·s?1時解析解和取M=N=60時數值解的結構圖Fig. 2 Profiles of analytical and numerical solutions withM=N=60 for Example 2 withν=1 m2·s?1，u1=u2= ?2 m·s?1

圖3 例3中解析解和取M=N=60時數值解的結構圖Fig. 3 Profiles of analytical solution, and numerical solution withM=N=60 for Example 3

3 結 論

本文提出了一種積分方程法用于研究穩態對流擴散問題。該方法首先利用Laplace方程的級數形式的格林函數將對流擴散方程轉化為積分方程。并利用級數的正交性，進一步將積分方程變成一個代數方程組，最后求得對流擴散問題的近似解。借助Chebyshev多項式和Fourier級數，應用積分方程法求解了非齊次邊界條件的一維對流擴散問題和齊次邊界條件的二維對流擴散問題。通過和有限體積法、有限元法和迎風差分法的比較，發現積分方程法具有很好的收斂性和很高的精度。而且即使選取較少的有限和項數用于逼近級數形式的解，積分方程法也取得了較好的精度，并且在處理對流占優的對流擴散問題時該方法并沒有出現數值振蕩，展示出非常好的穩定性。

符 號 說 明

L——長度，m

Lx——x方向的長度，m

Ly——y方向的長度，m

S——源項，K·s?1

T——溫度場，℃

T0——0處的溫度，℃

u——速度向量

u——速度，m·s?1

u1——x方向的速度，m·s?1

u2——y方向的速度，m·s?1

Г——擴散系數，kg·m?1·s?1

η——量綱1擴散系數

Θ——量綱1溫度

ν——熱擴散率，m2·s?1

ξ——量綱1長度

ρ——密度，kg·m?3

[1] Fudym O, Pradère C, Batsale J C. An analytical two-temperature model for convection-diffusion in multilayered systems: application to the thermal characterization of microchannel reactors [J].Chem. Eng. Sci., 2007, 62(15): 4054-4064.

[2] Li Huanyong, Gu Zhi, Zhang Haiyang, Li Wenwei. Thermodynamic analysis and growth of ZnSe single crystals in Zn-Se-I2system [J].J. Cryst. Growth, 2015, 415: 158-165.

[3] Lü Hexiang (呂和祥), Qiu Kunyu (邱崑玉), Chen Jianfeng (陳建峰). Solution under convective coordinate for convection-diffusion chemical reaction kinetic equations [J].CIESC Journal(化工學報), 2005, 56(1): 41-46.

[4] Yang Xiaohua, Yang Zhifeng, Yin Xinan, Li Jianqiang. Chaos gray-coded genetic algorithm and its application for pollution source identifications in convection-diffusion equation [J].Commun. Nonlinear Sci. Numer. Simulat., 2008, 13(8): 1676-1688.

[5] Song Wei (宋偉), Kong Qingyuan (孔慶媛), Li Hongmei (李洪枚). Characteristics of VOC mass transfer from wood furniture surface [J].CIESC Journal(化工學報), 2013, 64(5): 1549-1560.

[6] Stynes M. Finite volume methods for convection- diffusion problems [J].J. Comput. Appl. Math., 1995, 63(1/2/3): 83-90.

[7] Sonar T. On the construction of essentially non-oscillatory finite volume approximations to hyperbolic conservation laws on general triangulations: polynomial recovery, accuracy and stencil selection [J].Comput. Methods Appl. Mech. Engrg., 1997, 140(1/2): 157-181.

[8] Farhloul M, Serghini Mounim A. A mixed-hybrid finite element method for convection-diffusion problems [J].Appl. Math. Comput., 2005, 171(2): 1037-1047.

[9] Stynes M. Steady-state convection-diffusion problems [J].Acta Numerica, 2005, 14: 445-508.

[10] Tian Z F, Ge Y B. A fourth-order compact ADI method for solving two-dimensional unsteady convection-diffusion problems [J].J. Comput. Appl. Math., 2007, 198(1): 268-286.

[11] Shih Y, Cheng J Y, Chen K T. An exponential-fitting finite element method for convection-diffusion problems [J].Appl. Math. Comput., 2011, 217(12): 5798-5809.

[12] Ge Yongbin, Cao Fujun. Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems [J].J. Comput. Phys., 2011, 230(10): 4051-4070.

[13] Zingg D W, De Rango S, Nemec M, Pulliam T H. Comparison of several spatial discretizations for the Navier-Stokes equations [J].J. Comput. Phys., 2000, 160(2): 683-704.

[14] de Rango S, Zingg D W. Higher-order spatial discretization for turbulent aerodynamic computations [J].AIAA Journal, 2001, 39(7): 1296-1304.

[15] van Leer B. Towards the ultimate conservative difference scheme (Ⅴ): A second-order sequel to Godunov?s method [J].J. Comput. Phys., 1979, 32(1): 101-136.

[16] Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions [R]. Reno, NV: AIAA, 1993.

[17] Venkatakrishnan V. Convergence to steady state solutions of the Euler equations on unstructured grids with limiters [J].J. Comput. Phys., 1995, 118(1): 120-130.

[18] Jiang Guangshan, Shu Chiwang. Efficient implementation of weighted ENO schemes [J].J. Comput. Phys., 1996, 126(1): 202-228.

[19] Friedrich O. Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids [J].J. Comput. Phys., 1998, 144(1): 194-212.

[20] Hu Changqing, Shu Chiwang. Weighted essentially non-oscillatory schemes on triangular meshes [J].J. Comput. Phys., 1999, 150(1): 97-127.

[21] Wen P H, Hon Y C, Li M, Korakianitis T. Finite integration method for partial differential equations [J].Appl. Math. Model., 2013, 37(24): 10092-10106.

[22] Hernandez-Martinez E, Puebla H, Valdes-Parada F, Alvarez-Ramirez J. Nonstandard finite difference schemes based on Green?s function formulations for reaction-diffusion-convection systems [J].Chem. Eng. Sci., 2013, 94: 245-255.

[23] Xu Mingtian, Stefani F, Gerbeth G. Integral equation approach to time-dependent kinematic dynamos in finite domains [J].Phys. Rev. E, 2004, 70(5): 056305.

[24] Xu Mingtian, Stefani F, Gerbeth G. The integral equation approach to kinematic dynamo theory and its application to dynamo experiments in cylindrical geometry [J].J. Comput. Phys., 2008, 227(17): 8130-8144.

[25] Xu Mingtian, Wei Tao, Zhang Longzhou, Hao Aiqin. An integral equation approach to two-dimensional incompressible resistive magnetohydrodynamics [J].Geophys. Astrophys. Fluid Dyn., 2014, 108(4): 463-478.

[26] Xu Mingtian. Effect of soft-iron impellers on the von Kármán-sodium dynamo [J].Phys. Rev. E, 2014, 89(1): 013012.

[27] Versteeg H K, Malalasekera W. An Introduction to Computational Fluid Dynamics: The Finite Volume Method[M]. 2nd ed. London: Pearson Prentice Hall, 2007: 115-155.

[28] van Niekerk F D, van Niekerk A. An optimal rational basis function in the finite element method for convection-diffusion problems in one space variables [J].Mathl. Comput. Modelling, 1994, 19(5): 43-50.

Integral equation approach to convection-diffusion problems

WEI Tao, XU Mingtian, WANG Yin
(Department of Engineering Mechanics,Shandong University,Jinan250061,Shandong,China)

In the present work, an integral equation approach is developed to solve the convection-diffusion equations. In this approach, Green?s function of the Laplace equation in the form of series is employed to transform the convection-diffusion equation into an integral equation. With the help of orthogonal polynomials, the integral equation is reduced to an algebraic equation system with a finite number of unknown variables. Finally, this integral equation approach is examined by three examples. The Chebyshev polynomial is used to approximate the one-dimensional convection-diffusion problem with nonhomogeneous boundary conditions and the Fourier series is for the two-dimensional convection-diffusion problem with homogeneous boundary conditions. The comparisons with the finite volume method, finite element method and upwind difference method show that the integral equation approach is more accurate and stable. The stability is also proved by the convection-dominated diffusion problems. Furthermore, it can achieve a satisfactory accuracy even with a small number of grid points.

computational fluid dynamics; integral equation approach; convection-diffusion equation; finite volume method; heat transfer

Prof. XU Mingtian, mingtian@sdu.edu.cn

10.11949/j.issn.0438-1157.20150364

O 351.2

：A

：0438—1157（2015）10—3888—07

2015-03-20收到初稿，2015-06-03收到修改稿。

聯系人：許明田。

：魏濤（1981—），男，博士研究生，講師。

國家自然科學基金項目（11272187）。

Received date: 2015-03-20.

Foundation item: supported by the National Natural Science Foundation of China (11272187).