多尺度有限元法結(jié)合分層網(wǎng)格模擬二維奇異攝動(dòng)的兩端邊界層問(wèn)題

孫美玲，江山

（1.南通職業(yè)大學(xué) 數(shù)學(xué)教研室，江蘇 南通 226007; 2.南通大學(xué)理學(xué)院，江蘇 南通 226019）

多尺度有限元法結(jié)合分層網(wǎng)格模擬二維奇異攝動(dòng)的兩端邊界層問(wèn)題

孫美玲1，2，江山2*

（1.南通職業(yè)大學(xué) 數(shù)學(xué)教研室，江蘇 南通 226007; 2.南通大學(xué)理學(xué)院，江蘇 南通 226019）

通過(guò)攝動(dòng)系數(shù)建立分層網(wǎng)格，用多尺度有限元法捕捉對(duì)流擴(kuò)散方程的兩端邊界層，研究二維奇異攝動(dòng)模型?；诜謱泳W(wǎng)格并利用多尺度基函數(shù)刻畫(huà)了邊界層的微觀信息，用有限的計(jì)算資源、較短的計(jì)算時(shí)間，得到了不依賴于攝動(dòng)系數(shù)、一致穩(wěn)定的模擬結(jié)果。

奇異攝動(dòng)；自適應(yīng)網(wǎng)格；兩端邊界層；多尺度有限元；一致穩(wěn)定

0 引言

文獻(xiàn)［3-9］主要考慮一維奇異攝動(dòng)的有效求解方法。KADALBAJOO等［3］采用三次B-樣條配點(diǎn)法在Shishkin網(wǎng)格上得到了關(guān)于最大模的收斂結(jié)果。但Shishkin分片等距網(wǎng)格存在一定局限性，粗略估計(jì)過(guò)渡點(diǎn)位置可能造成方法精度不高。GENG等［4］用再生核方法模擬兩端邊界層現(xiàn)象；楊宇博等［5］研究非對(duì)稱內(nèi)罰間斷有限元法在分層網(wǎng)格中的一致收斂性，在一維情形下對(duì)左端附近的分層網(wǎng)格進(jìn)行加密構(gòu)建。受其啟發(fā)，本文將其拓展為左右兩端附近皆可用自適應(yīng)迭代的網(wǎng)格，適用于二維情形下的方向。ZHENG等［6］用混合有限差分格式處理擬線性邊值問(wèn)題；基于Bakhvalov-Shishkin網(wǎng)格，江山等［7］、鄭權(quán)等［8］分別得到多尺度有限元法、混合差分法的奇異攝動(dòng)高精度結(jié)果；CHENG［9］利用局部間斷有限元法有效模擬了雙參數(shù)模型，并基于各范數(shù)給出了穩(wěn)定性估計(jì)。

對(duì)二維奇異攝動(dòng)的研究已取得一定成果，如FRANZ等［10］利用單元邊界穩(wěn)定化技術(shù)，證明了高階有限元格式總能得到理想的一致超收斂；BRDAR等［11］面向二維對(duì)流反應(yīng)擴(kuò)散的雙參數(shù)方程，基于Duran-Lombardi與Duran-Shishkin網(wǎng)格，用傳統(tǒng)雙線性有限元法得到了比Shishkin網(wǎng)格精度更高的能量范數(shù)誤差；JIANG等［12］提出了獨(dú)立構(gòu)造試探函數(shù)空間、檢驗(yàn)函數(shù)空間的Petrov-Galerkin多尺度有限元基函數(shù)，在等級(jí)網(wǎng)格上自動(dòng)消除共振誤差，得到了精確、高效、穩(wěn)定的一致收斂模擬；LI等［13］針對(duì)反應(yīng)擴(kuò)散方程，基于分層網(wǎng)格證明了高階Galerkin有限元的優(yōu)化理論；XU等［14］、CAI等［15］分別給出了組合有限元、高階有限元的魯棒性分析與數(shù)值驗(yàn)證，但在二維情形下這類有限元模擬的計(jì)算代價(jià)較大。

本文針對(duì)二維奇異攝動(dòng)中小參數(shù)導(dǎo)致的兩端邊界層問(wèn)題，利用基于分層網(wǎng)格的多尺度有限元計(jì)算格式，實(shí)現(xiàn)較傳統(tǒng)有限元計(jì)算格式數(shù)值精度更高、計(jì)算代價(jià)更小、運(yùn)算時(shí)間更短的一致穩(wěn)定結(jié)果。

1 模型

考慮二維情形的對(duì)流擴(kuò)散方程

其中，雙線性形式為

2 分層網(wǎng)格與解的分解

2.1 兩端邊界層的分層網(wǎng)格

為得到式（2）變分形式的有效近似，用有限維逼近無(wú)限維的思想進(jìn)行區(qū)域離散。因常規(guī)的一致網(wǎng)格難以有效求解奇異攝動(dòng)小參數(shù)問(wèn)題，即使剖分?jǐn)?shù)很大，其等距步長(zhǎng)也無(wú)法滿足，故難以形成可靠的分辨率。對(duì)二維區(qū)域，先在方向形成適合左右兩端邊界層的優(yōu)化分層網(wǎng)格，再類似處理方向的上下兩端邊界層。

分層（graded）網(wǎng)格［5］由迭代格式生成，用攝動(dòng)系數(shù)和網(wǎng)格參數(shù)計(jì)算。方向的節(jié)點(diǎn)為

從而形成一端稠密、另一端稀疏的分層網(wǎng)格。為滿足兩端均有邊界層的情況，將式（3）改進(jìn)為

下文將驗(yàn)證分層網(wǎng)格是一種能自適應(yīng)地逼近邊界層位置及寬度的優(yōu)化網(wǎng)格，其剖分?jǐn)?shù)并非簡(jiǎn)單地成倍增加，從而突破了一致網(wǎng)格和Shishkin網(wǎng)格剖分?jǐn)?shù)偶數(shù)倍加密的局限，得到了更好的數(shù)值精度與穩(wěn)定結(jié)果。

2.2 解的多尺度分解

依據(jù)解的多尺度性質(zhì)，將其分解為若干部分之和，即

用式（4）對(duì)兩端邊界層的分層網(wǎng)格進(jìn)行區(qū)域離散，再采用多尺度有限元計(jì)算格式，使其更好地逼近子分塊上多尺度解式（5）的局部形態(tài)。

圖1 區(qū)域的子分塊Fig.1 Sub-domains of domain

3 有限元法與多尺度有限元法

3.1 有限元的變分原理

傳統(tǒng)有限元法（FEM）通過(guò)分片多項(xiàng)式構(gòu)造基底以形成有限維函數(shù)空間。如選定一組基，記有限元空間，其變分形式對(duì)應(yīng)為尋求，使得

3.2 多尺度有限元的變分原理

不同于傳統(tǒng)有限元，多尺度有限元法（MsFEM）在構(gòu)造有限維空間時(shí)，不采用顯式多項(xiàng)式函數(shù)，而采用基于與原問(wèn)題相同的微分算子在粗風(fēng)格單元中求解非顯式基函數(shù)。在每個(gè)粗網(wǎng)格單元中用有限元法求對(duì)應(yīng)的齊次子問(wèn)題：

3.3 多尺度有限元解的誤差估計(jì)

度量誤差，結(jié)合分層網(wǎng)格得到收斂結(jié)果。

定理1

為二維分層網(wǎng)格的生成函數(shù)。

4 數(shù)值驗(yàn)證

用已有文獻(xiàn)算例和程序結(jié)果驗(yàn)證相應(yīng)方法的精度和效率，本文僅討論當(dāng)很小時(shí)產(chǎn)生奇異攝動(dòng)邊界層求解困境的情況。將分層網(wǎng)格上與傳統(tǒng)有限元、多尺度有限元對(duì)應(yīng)的結(jié)果分別記作FEM（G）和MsFEM（G），用真解、近似解和誤差的三維圖示、范數(shù)值分析度量實(shí)際模擬效果。

圖2 當(dāng)=時(shí)的真解Fig.2 Exact solutions when and ，respectively

圖3 當(dāng)時(shí)傳統(tǒng)有限元法與多尺度有限元法的解Fig.3 The solutions of FEM（G）and MsFEM（G） when

為更清晰地展現(xiàn)相應(yīng)方法的精確性與穩(wěn)定性，通過(guò)網(wǎng)格加密的方法觀察數(shù)值變化。表格格式與文獻(xiàn)［5］的表1與表2一致，區(qū)別在于文獻(xiàn)［5］處理的是一維問(wèn)題，行數(shù)較多、單方向剖分?jǐn)?shù)較小，本文研究的是二維問(wèn)題，因受算力限制，行數(shù)較少、單方向兩端邊界層的剖分?jǐn)?shù)較大。由表1知，無(wú)論攝動(dòng)參數(shù)如何選取，依據(jù)網(wǎng)格參數(shù)的遞減，由迭代式（4）自適應(yīng)生成兩端疏密不同的分層網(wǎng)格，用于離散化計(jì)算。橫向看，表1中單方向剖分?jǐn)?shù)、范數(shù)誤差均微增，縱向看，其范數(shù)誤差隨遞減呈穩(wěn)定收斂。表1為用傳統(tǒng)有限元法求解二維問(wèn)題，若網(wǎng)格剖分?jǐn)?shù)較大，計(jì)算消耗很大，繼續(xù)剖分將超出單機(jī)的運(yùn)行限定，具有一定的局限性。表2采用的是多尺度有限元法，僅在較粗的分層網(wǎng)格上計(jì)算誤差的能量范數(shù)，求式（7）時(shí)其子單元剖分?jǐn)?shù)，對(duì)應(yīng)行的精度略遜于表1，但其計(jì)算消耗小、時(shí)間短，繼續(xù)剖分能得到更精確的結(jié)果，其收斂結(jié)果與理論估計(jì)式（12）一致。另需指出，多尺度有限元法的精度在剖分?jǐn)?shù)=222，262，300時(shí)較傳統(tǒng)有限元法在剖分?jǐn)?shù)=240，288，328時(shí)更高。當(dāng)然此優(yōu)化結(jié)果也有計(jì)算消耗，主要用以刻畫(huà)奇異攝動(dòng)問(wèn)題的邊界層微觀屬性。

表1 不同下傳統(tǒng)有限元法在分層網(wǎng)格上的誤差能量范數(shù)Table 1 FEM（G）with different parameters for errors of energy norm

表1 不同下傳統(tǒng)有限元法在分層網(wǎng)格上的誤差能量范數(shù)Table 1 FEM（G）with different parameters for errors of energy norm

h誤差誤差誤差1364.510×10-31604.664×10-31924.748×10-32401.310×10-32881.362×10-33281.402×10-34483.870×10-45364.054×10-46164.189×10-48881.083×10-41 0481.142×10-41 2001.186×10-4

表2 不同下多尺度有限元法在分層網(wǎng)格上的誤差能量范數(shù)Table 2 MsFEM（G）with different parameters for errors of energy norm

表2 不同下多尺度有限元法在分層網(wǎng)格上的誤差能量范數(shù)Table 2 MsFEM（G）with different parameters for errors of energy norm

h誤差誤差誤差348.866×10-2409.084×10-2489.083×10-2602.270×10-2722.303×10-2822.329×10-21124.686×10-31344.938×10-31544.955×10-32229.228×10-42629.565×10-43009.627×10-4

（a） 傳統(tǒng)有限元法（b） 多尺度有限元法

圖4 當(dāng)時(shí)傳統(tǒng)有限元法于和多尺度有限元法于的誤差
Fig.4 Errors of FEM（G） onand MsFEM（G）onwhen

圖5 當(dāng)時(shí)傳統(tǒng)有限元法于和多尺度有限元法于的誤差Fig.5 Errors of FEM（G） on and MsFEM（G）on when

在上述數(shù)值精度與穩(wěn)定分析的基礎(chǔ)上，考慮2種數(shù)值方法所需的運(yùn)行時(shí)間和效率，表3給出了當(dāng)時(shí)臺(tái)式機(jī)Intel Core i9 CPU 3.7 GHz運(yùn)行相應(yīng)程序所需的CPU時(shí)間，可見(jiàn)在較密二維網(wǎng)格上用傳統(tǒng)有限元法所需的CPU時(shí)間是同一行較粗二維網(wǎng)格用多尺度有限元法的近10倍，顯然多尺度有限元法的計(jì)算效率更高。進(jìn)一步，圖6為攝動(dòng)參數(shù)取更?。ㄅc）時(shí)相應(yīng)方法的CPU時(shí)間與剖分?jǐn)?shù)的對(duì)數(shù)比例關(guān)系，再次證實(shí)了多尺度有限元法的計(jì)算代價(jià)更小、計(jì)算效率更高。

表3 當(dāng)時(shí)傳統(tǒng)有限元法與多尺度有限元法的CPU時(shí)間Table 3 FEM（G）and MsFEM（G）apos;s CPU time when

表3 當(dāng)時(shí)傳統(tǒng)有限元法與多尺度有限元法的CPU時(shí)間Table 3 FEM（G）and MsFEM（G）apos;s CPU time when

FEM（G）CPU時(shí)間/sMsFEM（G）CPU時(shí)間/s4.13.349126598410 2771 106

圖6 =和時(shí)2種方法的剖分?jǐn)?shù)與CPU時(shí)間的log-log圖示Fig.6 Two methodsapos; log-log on partition and CPU time when and

綜上所述，多尺度有限元法只需在較粗分層網(wǎng)格上進(jìn)行計(jì)算，消耗的計(jì)算資源較少，且能保證穩(wěn)定收斂的有效精度，因此，多尺度有限元法在高維奇異攝動(dòng)問(wèn)題求解中具有廣闊的應(yīng)用前景。

5 結(jié)束語(yǔ)

基于自適應(yīng)的分層網(wǎng)格生成機(jī)制，主要利用多尺度有限元法處理奇異攝動(dòng)的二維對(duì)流擴(kuò)散變系數(shù)方程。用分層迭代精確逼近邊界層位置與寬度，結(jié)合多尺度計(jì)算格式有效捕捉了兩端邊界層的微觀信息，實(shí)現(xiàn)了不依賴攝動(dòng)系數(shù)的精確高效模擬結(jié)果，充分展現(xiàn)了多尺度有限元法結(jié)合分層網(wǎng)格求解高維奇異攝動(dòng)問(wèn)題的一致穩(wěn)定性和優(yōu)勢(shì)。

［1］蘇煜城，吳啟光. 奇異攝動(dòng)問(wèn)題數(shù)值方法引論［M］. 重慶：重慶出版社， 1991.

SU Y C， WU Q G. An Introduction to Numerical Methods for the Singular Perturbation Problems［M］. Chongqing： Chongqing Publishing House，1991.

［2］MILLER J J H， O’RIORDAN E，SHISHKIN G I. Fitted Numerical Methods for Singular Perturbation Problems （Revised Edition）［M］. Singapore： World Scientific，2012.

［3］KADALBAJOO M K， GUPTA V. A parameter uniform B-spline collocation method for solving singularly perturbed turning point problem having twin boundary layers［J］. International Journal of Computer Mathematics， 2010，87（14）：3218-3235.

［4］GENG F Z， QIAN S P. Reproducing kernel method for singularly perturbed turning point problems having twin boundary layers［J］. Applied Mathematics Letters， 2013，26（10）：998-1004. DOI：10.1016/j.aml.2013.05.006

［5］楊宇博，祝鵬，尹云輝. 分層網(wǎng)格上奇異攝動(dòng)問(wèn)題的一致NIPG分析［J］. 計(jì)算數(shù)學(xué)，2014， 36（4）：437-448.

YANG Y B， ZHU P，YIN Y H. Uniform analysis of the NIPG method on graded meshes for singularly perturbed convection-diffusion problems［J］. Mathematica Numerica Sinica， 2014，36（4）： 437-448.

［6］ZHENG Q， LI X Z，GAO Y. Uniformly convergent hybrid schemes for solutions and derivatives in quasilinear singularly perturbed BVPs［J］. Applied Numerical Mathematics， 2015，91： 46-59. DOI：10. 1016/j.apnum.2014.12.010

［7］江山，孫美玲. 多尺度有限元結(jié)合Bakhvalov-Shishkin網(wǎng)格法高效處理邊界層問(wèn)題［J］. 浙江大學(xué)學(xué)報(bào)（理學(xué)版），2015， 42（2）：142-146. DOI：10. 3785/j.issn.1008-9497.2015.02.004

JIANG S， SUN M L. Combining the multiscale finite element and Bakhvalov-Shishkin grid to solve the boundary layer problems［J］. Journal of Zhejiang University （Science Edition）， 2015，42（2）： 142-146. DOI：10.3785/j.issn.1008-9497.2015.02.004

［8］鄭權(quán)，劉穎，劉忠禮. 奇異攝動(dòng)問(wèn)題在修正的Bakhvalov-Shishkin網(wǎng)格上的混合差分格式［J］. 浙江大學(xué)學(xué)報(bào)（理學(xué)版）， 2020，47（4）： 460-468. DOI：10. 3785/j.issn.1008-9497.2020.04.009

ZHENG Q， LIU Y，LIU Z L. The hybrid finite difference schemes on the modified Bakhvalov-Shishkin mesh for the singularly perturbed problem［J］. Journal of Zhejiang University （Science Edition）， 2020，47（4）： 460-468. DOI：10.3785/j.issn.1008-9497.2020. 04.009

［9］CHENG Y. On the local discontinuous Galerkin method for singularly perturbed problem with two parameters［J］. Journal of Computational and Applied Mathematics， 2021，392： 113485. DOI：10. 1016/j.cam.2021.113485

［10］FRANZ S， LIN? T，ROOS H G， et al. Uniform superconvergence of a finite element method with edge stabilization for convection-diffusion problems［J］. Journal of Computational Mathematics， 2010，28（1）： 32-44. DOI：10.4208/jcm.2009.09-m1005

［11］BRDAR M， ZARIN H，TEOFANOV L. A singularly perturbed problem with two parameters in two dimensions on graded meshes［J］. Computers and Mathematics with Applications， 2016，72（10）：2582-2603. DOI：10.1016/j.camwa.2016.09.021

［12］JIANG S， PRESHO M，HUANG Y Q. An adapted Petrov-Galerkin multi-scale finite element method for singularly perturbed reaction-diffusion problems［J］. International Journal of Computer Mathematics， 2016，93（7）： 1200-1211. DOI：10.1080/00207160. 2015.1041935

［13］LI Z W， WU B，XU Y S. High order Galerkin methods with graded meshes for two-dimensional reaction-diffusion problems［J］. International Journal of Numerical Analysis and Modeling， 2016，13（3）： 319-343.

［14］XU S P， DENG W B，WU H J. A combined finite element method for elliptic problems posted in domains with rough boundaries［J］. Journal of Computational and Applied Mathematics， 2018，336： 235-248. DOI：10.1016/j.cam.2017.12.049

［15］CAI D F， CAI Z Q，ZHANG S. Robust equilibrated a posteriori error estimator for higher order finite element approximations to diffusion problems［J］. Numerische Mathematik， 2020，144： 1-21. DOI：10.1007/s00211-019-01075-1

Simulation of multiscale finite element method on graded meshes for two-dimensional singularly perturbed twin boundary layers problems

SUN Meiling1,2, JIANG Shan2

（1. Department of Mathematics，Nantong Vocational University，Nantong226007，Jiangsu Province，China；2. School of Science，Nantong University，Nantong226019，Jiangsu Province，China）

To solve a two-dimensional singularly perturbed model, a multiscale finite element method on graded meshes built from the perturbed parameter is presented for capturing the twin boundary layers of convection-diffusion equations effectively. Based on the graded meshes, the multiscale basis functions are capable of subtly describing the microscopic information in the boundary layers. No wonder, it just costs a handful of computing resource and short time to achieve the accurate and efficient results, and the results are independent of the perturbed parameter with uniform stability.

singular perturbation; adaptive meshes; twin boundary layers; multiscale finite element; uniform stability

O 241.82

A

1008?9497（2022）05?564?06

2021?08?24.

南通市基礎(chǔ)科學(xué)研究指令性項(xiàng)目（JC2021123）；國(guó)家自然科學(xué)基金面上項(xiàng)目（11771224）；江蘇省高校青藍(lán)工程優(yōu)秀骨干教師資助項(xiàng)目.

孫美玲（1981—），ORCID：https：//orcid.org/0000-0003-0061-5155，女，博士，副教授，主要從事偏微分方程數(shù)值解及其應(yīng)用研究，E-mail：sunmeiling81@163.com.

通信作者，ORCID：https：//orcid.org/0000-0001-7983-0012，E-mail：jiangshan@ntu.edu.cn.

10.3785/j.issn.1008-9497.2022.05.007