魏濤，許明田，汪引

（山東大學(xué)土建與水利學(xué)院工程力學(xué)系，山東 濟(jì)南 250061）

求解對流擴(kuò)散問題的積分方程法

魏濤，許明田，汪引

（山東大學(xué)土建與水利學(xué)院工程力學(xué)系，山東 濟(jì)南 250061）

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

計(jì)算流體力學(xué)；積分方程法；對流擴(kuò)散方程；有限體積法; 傳熱

引 言

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

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

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

1 數(shù)值方法

在本文中，考慮下面的穩(wěn)態(tài)對流擴(kuò)散問題

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

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

2 數(shù)值實(shí)驗(yàn)

2.1 一維對流擴(kuò)散問題

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

例1：考慮如下一維穩(wěn)態(tài)對流擴(kuò)散問題

邊界條件

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

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

邊界條件

這個(gè)問題的解析解為

考慮下面3種情形。

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

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

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

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

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

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

表3 當(dāng)η=1/60、M=20時(shí)由各種方法得到的結(jié)果Table 3 Results obtained by various methods forη=1/60, M=20

表4 當(dāng)η=1/60、M=40時(shí)由各種方法得到的結(jié)果Table 4 Results obtained by various methods forη=1/60, M=40

2.2 二維對流擴(kuò)散問題

選取的正交多項(xiàng)式為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，則該對流擴(kuò)散問題的解析解為

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

在式（24）中取M=N=20，相當(dāng)于采用20×20的網(wǎng)格。用下面的誤差公式計(jì)算絕對誤差和相對誤差

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

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

表6 積分方程法所得對流占優(yōu)對流擴(kuò)散問題數(shù)值解的穩(wěn)定性Table 6 Stability of numerical solution obtained by IEA for convection-dominated diffusion problems

為了進(jìn)一步驗(yàn)證該方法，給出了下面的例子。

例3：在本例中，取ν=1 m2·s?1，u1=u2= ?2 m·s?1，源項(xiàng)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]所得的結(jié)果列于表7中，所有計(jì)算用Fortran語言編寫程序，并在DELL OPTIPLEX 790臺式計(jì)算機(jī)上完成運(yùn)算。從該表中可以看出在同等網(wǎng)格數(shù)量下用積分方程法進(jìn)行計(jì)算所需的時(shí)間略多于有限體積法，但是積分方程法的精度要明顯好于有限體積法。圖3(a)、(b)分別展示了解析解和由積分方程法得到的數(shù)值解。通過比較表7和表5發(fā)現(xiàn)例3中的絕對誤差εmax比例2中的大，由圖3和圖2可以看到例3中解的結(jié)構(gòu)復(fù)雜并且很多點(diǎn)處的梯度也比較大，這應(yīng)是表7中的絕對誤差比表5大的原因。但是表7和表5中的相對誤差εL2在數(shù)量級上并沒有明顯的差別。

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

圖2 例2中當(dāng)ν=1 m2·s?1、u1=u2= ?2 m·s?1時(shí)解析解和取M=N=60時(shí)數(shù)值解的結(jié)構(gòu)圖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時(shí)數(shù)值解的結(jié)構(gòu)圖Fig. 3 Profiles of analytical solution, and numerical solution withM=N=60 for Example 3

3 結(jié) 論

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

符 號 說 明

L——長度，m

Lx——x方向的長度，m

Ly——y方向的長度，m

S——源項(xiàng)，K·s?1

T——溫度場，℃

T0——0處的溫度，℃

u——速度向量

u——速度，m·s?1

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

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

Г——擴(kuò)散系數(shù)，kg·m?1·s?1

η——量綱1擴(kuò)散系數(shù)

Θ——量綱1溫度

ν——熱擴(kuò)散率，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(化工學(xué)報(bào)), 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(化工學(xué)報(bào)), 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收到修改稿。

聯(lián)系人：許明田。

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

國家自然科學(xué)基金項(xiàng)目（11272187）。

Received date: 2015-03-20.

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