求解對(duì)流擴(kuò)散問(wèn)題的積分方程法

魏濤，許明田，汪引

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

求解對(duì)流擴(kuò)散問(wèn)題的積分方程法

魏濤，許明田，汪引

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

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

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

引 言

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

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

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

1 數(shù)值方法

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

其中，φ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è)線性無(wú)關(guān)的，然后它們組成一個(gè)線性代數(shù)方程組。求解該方程組即可得到未知系數(shù)Tmn，再將這些系數(shù)代入式（18）就得到了式（1）滿足邊界條件式（2）的近似解。

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

2.1 一維對(duì)流擴(kuò)散問(wèn)題

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

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

邊界條件

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

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

邊界條件

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

考慮下面3種情形。

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

2.1.2 情形2:η=1/25 對(duì)于該種情形，由積分方程法、解析解和有限體積法得到的結(jié)果列在表2中。當(dāng)ξ趨向于1時(shí)，梯度的絕對(duì)值越來(lái)越大，從該表可以看出，有限體積法的解在梯度大的點(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ù)η較小。為了更清晰地展示對(duì)應(yīng)于不同η的對(duì)流擴(kuò)散問(wèn)題的解的不同變化趨勢(shì)，將3種情形下的解析解都畫在圖1中。由圖可見(jiàn)，對(duì)于較大的η，場(chǎng)變量Θ的變化比較平穩(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中對(duì)應(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 二維對(duì)流擴(kuò)散問(wèn)題

選取的正交多項(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，則該對(duì)流擴(kuò)散問(wèn)題的解析解為

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

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

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

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

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

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

符 號(hào) 說(shuō) 明

L——長(zhǎng)度，m

Lx——x方向的長(zhǎng)度，m

Ly——y方向的長(zhǎng)度，m

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

T——溫度場(chǎng)，℃

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長(zhǎng)度

ρ——密度，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—），男，博士研究生，講師。

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

Received date: 2015-03-20.

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