李洪建，韓立國，鞏向博

(吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春130026)

隨著油氣勘探開發(fā)的不斷深入，油氣藏勘探逐漸從常規(guī)構(gòu)造油氣藏向復(fù)雜隱伏構(gòu)造油氣藏和隱蔽油氣藏發(fā)展，復(fù)雜的地下地質(zhì)結(jié)構(gòu)和地層條件對地震勘探方法技術(shù)提出了更高的要求。地震波場數(shù)值模擬是研究復(fù)雜構(gòu)造地震波傳播問題的重要工具，人們充分利用現(xiàn)有的地質(zhì)和地球物理資料，不斷提高數(shù)值模擬的精度，為油氣藏勘探開發(fā)提供更準(zhǔn)確的依據(jù)[1]。

正演數(shù)值模擬在油氣勘探中應(yīng)用廣泛，李軍鋒[2]分析了高精度有限差分波場正演模擬方法，并對張性斷層和逆掩推覆斷層進行了模擬計算；王秀明等[3]利用高階交錯網(wǎng)格有限差分法模擬了地震波在烴類儲集層中的傳播；劉春園等[4]利用變網(wǎng)格有限差分法對孔洞模型等碳酸鹽巖儲層進行了研究。有限差分法計算速度快，但網(wǎng)格剖分方法無法適應(yīng)不規(guī)則異常地質(zhì)界面的情況，且存在數(shù)值頻散。劉炯等[5]提出了錯格偽譜法并模擬了不均勻縫洞型儲層介質(zhì)中彈性波的傳播問題；劉淵博等[6]提出了分段點最佳化的網(wǎng)格優(yōu)化算法，提高了偽譜法對非光滑最優(yōu)控制問題的求解效率。偽譜法模擬精度高，但偽譜法網(wǎng)格剖分同樣無法很好地適應(yīng)邊界的變化。劉廣峰等[7]利用有限元方法開展了油氣儲層地應(yīng)力模擬研究；王月英[8]采用MPI并行計算方法，提升了有限元正演模擬數(shù)據(jù)容量和計算速度。有限元法可以很好地模擬復(fù)雜邊界構(gòu)造，但大規(guī)模的計算量使它難以發(fā)揮優(yōu)勢。

譜元法數(shù)值模擬方法結(jié)合了偽譜法高精度和有限元法靈活處理復(fù)雜地質(zhì)結(jié)構(gòu)的優(yōu)點，在模擬復(fù)雜構(gòu)造方面具有很大的優(yōu)勢，在應(yīng)用到地震勘探領(lǐng)域的20多年里，國內(nèi)外專家、學(xué)者對譜元法數(shù)值模擬技術(shù)及其應(yīng)用做了一系列的改進，使其更具有實用性。Patera[9]最早將譜元法應(yīng)用于流體動力學(xué)計算；Komatitsch等[10-17]編寫了從一維到三維的譜元法計算程序，在地震波場計算中取得了很好的應(yīng)用效果；王秀明等[18]提出了逐元算法，大大提高了譜元法計算效率。但是，常規(guī)的譜元網(wǎng)格剖分方法仍受到地層厚度限制，導(dǎo)致部分區(qū)域網(wǎng)格剖分過于致密或稀疏，模擬精度受到很大影響。對于復(fù)雜構(gòu)造的地震正演模擬，發(fā)展一種既快速方便，又能靈活適應(yīng)地層邊界變化的正演算法尤為重要。

在譜元法數(shù)值模擬中，靈活使用四邊形網(wǎng)格對模型進行空間離散，能更好地逼近復(fù)雜構(gòu)造。在三角形網(wǎng)格內(nèi)部插入一個節(jié)點，連接這一節(jié)點和三角形質(zhì)心，網(wǎng)格可以被剖分為3個四邊形，但這種方法剖分出的四邊形往往形狀不規(guī)則，單元質(zhì)量不高；對于合適的幾何域，四邊形映射網(wǎng)格剖分(quadrilateral mapped meshing)[19]或子映射(sub-mapping)[20]方法可以將整個區(qū)域離散為一組結(jié)構(gòu)化的四邊形，每個內(nèi)部節(jié)點都連接有4個單元，這種方法快速且有效，但是僅僅適用于某些特定的問題。為此，更為普遍的非結(jié)構(gòu)化算法被提出。非結(jié)構(gòu)化四邊形網(wǎng)格算法可以分為直接法和間接法兩類。直接法四邊形網(wǎng)格剖分方法可以借助于域分解[21-22]把區(qū)域直接剖分為簡單的可變換凸四邊形單元，但是其剖分結(jié)果存在大量不穩(wěn)定單元。區(qū)別于直接法，間接法四邊形網(wǎng)格剖分方法首先將域剖分為三角形網(wǎng)格，進而合成四邊形網(wǎng)格。Johnston等[23]使用拆分和交換方法，三角形單元被盡可能全部重組為四邊形單元，但仍有一部分三角形網(wǎng)格余留；Lee等[24]和Lo[25]提出了首先在邊界處將三角形重組形成四邊形單元，進而利用波前法在區(qū)域內(nèi)部依次將剩余所有的三角形重組為四邊形單元的方法，這種重組方法得到的四邊形網(wǎng)格質(zhì)量較高，剖分效果較好。

我們研究并提出了復(fù)雜構(gòu)造網(wǎng)格剖分及其高精度譜元數(shù)值模擬方法。首先基于Delaunay三角網(wǎng)格剖分方法將區(qū)域剖分為三角形網(wǎng)格；然后采用波前法將所有三角形網(wǎng)格重組，得到譜元法所需的四邊形網(wǎng)格。應(yīng)用基于復(fù)雜構(gòu)造網(wǎng)格剖分的高精度譜元法對Marmousi速度模型和二維隨機介質(zhì)模型進行數(shù)值模擬，并與常規(guī)譜元法的模擬結(jié)果進行對比，以驗證其有效性和應(yīng)用效果。

1 方法原理

根據(jù)Galerkin加權(quán)余量原理，與微分方程初、邊值問題對應(yīng)的波動方程弱形式可以表示為

(1)

將模型剖分為若干個離散的單元，在每個單元上做離散化運算。對于測試函數(shù)空間內(nèi)的位移值uh，有

(2)

式中：gh表示邊界位移；vh表示除邊界外測試函數(shù)空間的位移值。用截斷的Chebyshev正交多項式或者Legendre多項式確定單元基函數(shù)NA(x)，測試函數(shù)wh，vh，gh在單元內(nèi)各個節(jié)點上的節(jié)點值分別為ciA(t)，diB(t)，giB(t)，則有

(3)

式中：η為全局節(jié)點集合；ηgi為滿足邊界條件的節(jié)點集合。經(jīng)過計算，譜元方程矩陣形式可以簡寫為

(4)

(5)

式中：ei，ej為單位向量；fi為單元體積的體力;hi為邊界應(yīng)力；cijkl為彈性系數(shù)。

圖1 Delaunay三角形網(wǎng)格重組示意圖解

針對某一起伏地表層狀模型，圖2給出了常規(guī)譜元法網(wǎng)格剖分結(jié)果和本文復(fù)雜構(gòu)造網(wǎng)格剖分法的剖分結(jié)果。從圖2中可以看出，在起伏地表和起伏界面處，二者都能很好地適應(yīng)邊界的變化；但對比可見圖2a剖分效果相對較差，因為常規(guī)譜元法的網(wǎng)格剖分方法必須保證每層具有特定的網(wǎng)格數(shù)目，由于地層厚度的變化，難免會導(dǎo)致部分區(qū)域網(wǎng)格剖分過于致密或稀疏，從而導(dǎo)致網(wǎng)格大小不一，差別很大；而圖2b則有效地避免了這一問題的產(chǎn)生，四邊形網(wǎng)格大小均勻，質(zhì)量較好。

圖2 起伏地表層狀模型的常規(guī)譜元法(a)和復(fù)雜構(gòu)造網(wǎng)格剖分方法(b)的網(wǎng)格剖分結(jié)果

Delaunay三角網(wǎng)格剖分的任意性導(dǎo)致重組之后的網(wǎng)格形狀為不規(guī)則四邊形，無法直接用于運算，需要通過參數(shù)變換，才能將單元運算法則應(yīng)用到形狀不規(guī)則的網(wǎng)格單元中。將參考單元映射到實際物理區(qū)域的變換方法稱為等參變換，即單元幾何形狀變換，單元內(nèi)場函數(shù)采用相同的節(jié)點參數(shù)及相同的插值函數(shù)，參與變換的單元稱為等參單元。圖3為具有4個節(jié)點單元的等參變換的基本示意圖。

圖3 4個節(jié)點單元的等參變換示意圖解

為將自然坐標(biāo)(ξ，η)中形狀規(guī)則的單元轉(zhuǎn)換成笛卡爾坐標(biāo)(x,y)中形狀扭曲的單元，需要建立一個坐標(biāo)變換：

(6)

并建立以下映射：

(7)

式中：Na表示基函數(shù)；nen表示單元中的節(jié)點數(shù)，此例中nen=4。由于

(8)

式中：Na,ξ，Na,η分別表示基函數(shù)Na在ξ，η方向的導(dǎo)數(shù)。通過推導(dǎo)，得到基函數(shù)Na關(guān)于x，y的導(dǎo)數(shù)：

(9)

Na,ξ和Na,η來自于幾何形狀規(guī)則的單元計算,容易求得。并且

(10)

其中，J為Jacobian矩陣。

(11)

即可得

(12)

通過等參變換，即可以利用規(guī)則坐標(biāo)來計算我們剖分所得的不規(guī)則網(wǎng)格。把所有單元按照一定的順序進行組裝，求解總體譜元方程，即能得到全局的近似解。

2 數(shù)值模擬算例

2.1 Marmousi速度模型

首先利用基于復(fù)雜構(gòu)造網(wǎng)格剖分方法的高精度譜元法和常規(guī)譜元法對Marmousi速度模型(圖4a) 進行模擬對比研究，模型大小為4170m×1810m。為了測試復(fù)雜構(gòu)造網(wǎng)格剖分方法的優(yōu)勢，按照大致的速度邊界將整個模型區(qū)域劃分為15個區(qū)域，在每個速度區(qū)域進行Delaunay三角網(wǎng)格剖分和間接法重組，得到譜元法計算所需要的四邊形網(wǎng)格(圖4b)，網(wǎng)格總數(shù)約為120000個。由于網(wǎng)格剖分較細(xì)，將圖4b中黑框區(qū)域放大如圖5a，再將圖5a中黑框區(qū)域放大如圖5b。對比圖6 所示的常規(guī)譜元法剖分結(jié)果(將模型區(qū)域剖分為400×300=120000個網(wǎng)格)，可以清楚地看出，復(fù)雜構(gòu)造網(wǎng)格剖分方法不但能很好地適應(yīng)復(fù)雜構(gòu)造速度邊界和傾角的變化，且網(wǎng)格均勻，質(zhì)量較高，對復(fù)雜地層適應(yīng)性更好。

圖4 Marmousi速度模型(a)及其復(fù)雜構(gòu)造網(wǎng)格剖分法的剖分結(jié)果(b)

圖5 Marmousi模型復(fù)雜構(gòu)造網(wǎng)格剖分法剖分結(jié)果(圖4a)按黑框區(qū)域的依次放大顯示a 圖4b中黑框區(qū)域的放大顯示; b 圖5a中黑框區(qū)域的放大顯示

圖6 Marmousi模型的常規(guī)譜元法網(wǎng)格剖分結(jié)果(a)及其黑框區(qū)域的放大顯示(b)

選取頻率為40Hz的雷克子波作為震源，激發(fā)點位置在模型表面(2085m,25m)處；在(100m,25m)到(4100m,25m)處均勻排列共251個接收道，道間距為16m；取采樣點數(shù)9000，采樣率0.2ms，總記錄長度為1.8s，分別采用基于復(fù)雜構(gòu)造網(wǎng)格剖分方法的高精度譜元法和常規(guī)譜元法對Marmousi速度模型做正演數(shù)值模擬。圖7為高精度譜元法數(shù)值模擬0.5，0.6，0.7，0.8s時刻的波場快照；圖8 為高精度譜元法和常規(guī)譜元法正演模擬的單炮記錄對比結(jié)果；圖9是對應(yīng)圖8中黑框區(qū)域(50～240道，4500～9000采樣點(0.9～1.8s))的放大顯示對比。

從高精度譜元法數(shù)值模擬的波場快照(圖7)可以看出，彈性波遇到Marmousi模型的速度邊界時，在界面處發(fā)生較清晰的反射和透射，波形清晰連續(xù)。但由于該模型比較復(fù)雜，其反射、透射也很復(fù)雜。

從模擬單炮記錄(圖8)及其放大結(jié)果(圖9)可明顯看出，與常規(guī)譜元法相比，基于復(fù)雜構(gòu)造網(wǎng)格剖分的高精度譜元法由于網(wǎng)格剖分更適應(yīng)實際速度界面，同相軸更清晰、連續(xù)，且頻散很小，模擬精度高，具有很高的信噪比和分辨率。

圖7 Marmousi模型高精度譜元法數(shù)值模擬不同時刻的波場快照a 0.5s; b 0.6s; c 0.7s; d 0.8s

圖8 Marmousi模型的高精度譜元法(a)和常規(guī)譜元法(b)正演模擬單炮記錄

圖9 Marmousi模型的高精度譜元法(a)和常規(guī)譜元法(b)正演模擬單炮記錄中黑框區(qū)域的放大顯示

2.2 隨機介質(zhì)模型

在石油勘探中，實際儲層(如碳酸鹽巖儲層)往往是非均勻的，大量微小異常極不規(guī)則地分布在非均勻介質(zhì)中，正演模擬中如果采用波長遠(yuǎn)大于異常尺度的地震波，非均勻異常對地震波的傳播影響較小。但在如今高精度的地球物理探測中，對于高、寬頻地震信號，這些小尺度異常卻無法忽略，因為它們對地震波的傳播產(chǎn)生十分明顯的影響。為了驗證復(fù)雜構(gòu)造高精度譜元法的應(yīng)用效果，建立如圖10 所示的二維水平地表層狀隨機介質(zhì)模型。模型大小為1200m×1200m，其中隨機介質(zhì)的物理參數(shù)如表1所示。

圖10 層狀隨機介質(zhì)速度模型

采用基于復(fù)雜構(gòu)造網(wǎng)格剖分的高精度譜元法將層狀隨機介質(zhì)模型區(qū)域剖分為大約90000個網(wǎng)格，結(jié)果示于圖11a，由于網(wǎng)格剖分較細(xì)，圖11b給出了圖11a中心區(qū)域的局部放大顯示。對比圖12所示的常規(guī)譜元法網(wǎng)格剖分結(jié)果(將模型區(qū)域剖分為300×300=90000個網(wǎng)格)，可以看出，復(fù)雜構(gòu)造網(wǎng)格剖分法的剖分結(jié)果大小均勻，質(zhì)量較高，對復(fù)雜界面適應(yīng)性更好。

為考慮非均勻性對地震波的影響，地震波長要盡可能接近異常的尺度，隨機介質(zhì)模型中非均勻異常不規(guī)則且尺度相對微小，因此需要選用高頻地震信號進行譜元法數(shù)值模擬。在模型中心(600m,480m)位置處選取頻率為150Hz的雷克子波作為震源，并在(250m,480m)到(950m,480m)均勻排列201個檢波器，取采樣點數(shù)3360，采樣率0.1ms，得到每個時刻的波場快照和單炮記錄。圖13為高精度譜元法數(shù)值模擬0.08，0.12，0.16，0.20s時刻的波場快照；圖14為二維隨機介質(zhì)速度模型高精度譜元法和常規(guī)譜元法正演模擬單炮記錄對比；圖15為對應(yīng)圖14中黑框區(qū)域(25～201道，2400～3360采樣點(0.240～0.336s))的放大顯示。

從高精度譜元法數(shù)值模擬的波場快照(圖13)

中可以看出，直達(dá)波的波前面不再是在均勻介質(zhì)中傳播時呈現(xiàn)的圓形，可見隨機介質(zhì)中存在速度方面的微小差異，具有明顯的不均勻性；在速度模型內(nèi)部明顯存在各種微小的散射波，它們由非均勻地質(zhì)體的速度差異和密度差異引起。

表1 隨機介質(zhì)物性參數(shù)

圖11 隨機介質(zhì)模型的復(fù)雜構(gòu)造網(wǎng)格剖分法剖分結(jié)果(a)及其中心區(qū)域的局部放大顯示(b)

圖12 隨機介質(zhì)模型的常規(guī)譜元法網(wǎng)格剖分結(jié)果(a)及其中心區(qū)域的局部放大顯示(b)

圖13 隨機介質(zhì)模型高精度譜元法數(shù)值模擬不同時刻的波場快照a 0.08s; b 0.12s; c 0.16s; d 0.20s

圖14 隨機介質(zhì)模型的高精度譜元法(a)和常規(guī)譜元法(b)正演模擬單炮記錄

圖15 隨機介質(zhì)模型高精度譜元法(a)和常規(guī)譜元法(b)模擬單炮記錄中黑框區(qū)域的放大顯示

從圖14和圖15的模擬單炮記錄來看，與常規(guī)譜元法相比，基于復(fù)雜構(gòu)造網(wǎng)格剖分的高精度譜元法模擬結(jié)果具有更高的信噪比和分辨率，同相軸連續(xù)性好，精度更高。高精度譜元法能將微觀層次的有效地震信息更清晰地反映在地震記錄剖面上，有助于在高精度地震勘探中發(fā)現(xiàn)小尺度的儲層非均質(zhì)體。

3 結(jié)論和討論

在空間上基于Delaunay三角網(wǎng)格剖分方法對模型區(qū)域進行剖分，再使用間接法首先在邊界處將三角形重組，進而利用波前法在區(qū)域內(nèi)部依次將剩余所有的三角形重組為譜元法計算所需的四邊形單元，形成了適應(yīng)復(fù)雜構(gòu)造模擬網(wǎng)格剖分的新方法。與常規(guī)網(wǎng)格剖分方法相比，復(fù)雜構(gòu)造網(wǎng)格剖分方法不但能靈活地適應(yīng)速度邊界的變化，且網(wǎng)格大小均勻、質(zhì)量較高，有效地解決了常規(guī)譜元法數(shù)值模擬在復(fù)雜構(gòu)造區(qū)域存在網(wǎng)格剖分過密或過疏的不足。

對Marmousi模型進行譜元法數(shù)值模擬的研究結(jié)果表明，在復(fù)雜網(wǎng)格空間離散條件下，基于復(fù)雜構(gòu)造網(wǎng)格剖分的高精度譜元法能夠更好地模擬復(fù)雜構(gòu)造，有效信號清晰且連續(xù)，數(shù)值頻散小，精度很高。隨機介質(zhì)的物理性質(zhì)(如速度、密度等)在微觀層次的微小差異導(dǎo)致介質(zhì)內(nèi)部存在波阻抗差，并影響其散射波場的強弱和分布，針對二維層狀隨機介質(zhì)模型的高精度譜元法數(shù)值模擬可以將這些有效信息更清晰地記錄下來。數(shù)值模擬結(jié)果展示了基于復(fù)雜構(gòu)造網(wǎng)格剖分的高精度譜元法在高精度地震勘探中具有較好的應(yīng)用前景，提高了勘探精度，有助于發(fā)現(xiàn)小尺度的儲層非均質(zhì)體。

參 考 文 獻(xiàn)

[1] 王錫文,秦廣勝,趙衛(wèi)鋒,等.正演模擬技術(shù)在地震采集設(shè)計中的應(yīng)用[J].地球物理學(xué)進展,2012,27(2)：642-650

Wang X W,Qing G S,Zhao W F,et al.Application of forward modeling in seismic acquisition design[J].Progress in Geophysics,2012,27(2)：642-650

[2] 李軍峰.高精度有限差分波場正演模擬方法[D].成都：成都理工大學(xué),2012

Li J F.High precision wave field forward modeling of finite difference method[D].Chengdu：Chengdu University of Technology,2012

[3] 王秀明,張海瀾,王東,等.利用高階交錯網(wǎng)格有限差分法模擬地震波在非均勻孔隙介質(zhì)中的傳播[J].地球物理學(xué)報,2003,46(6)：842-849

Wang X M,Zhang H L,Wang D,et al.Simulations of seismic wave propagation in inhomogeneous porous media using high-order staggered finite difference method[J].Chinese Journal of Geophysics,2003,46(6)：842-849

[4] 劉春園,朱生旺,魏修成,等.隨機介質(zhì)地震波正演模擬在碳酸鹽巖儲層預(yù)測中的應(yīng)用[J].石油物探,2010,49(2):133-139

Liu C Y,Zhu S W,Wei X C,et al.Random medium seismic forward modeling in carbonate reservoir prediction[J].Geophysical Prospecting for Petroleum,2010,49(2):133-139

[5] 劉炯,馬堅偉,楊慧珠,等.縫洞型儲層錯格偽譜法模擬[J].石油地球物理勘探,2008,43(6):723-727

Liu J,Ma J W,Yang H Z,et al.Vuggy reservoir simulation staggered grid pseudo-spectral method[J].Oil Geophysical Prospecting,2008,43(6):723-727

[6] 劉淵博,朱恒偉.偽譜法求解非光滑最優(yōu)控制問題的網(wǎng)格優(yōu)化[J].系統(tǒng)工程與電子技術(shù),2013,35(11)：2396-2399

Liu Y B,Zhu H W.Pseudo-spectral method for solving non-smooth mesh optimization optimal control problems[J].Systems Engineering and Electronics,2013,35(11)：2396-2399

[7] 劉廣峰,陸紅軍,何順利.有限元法開展油氣儲層地應(yīng)力研究綜述[J].科學(xué)技術(shù)與工程,2009,9(24):7430-7435

Liu G F,Lu H J,He S L.A Summary of oil and gas reservoirs to stress finite element method[J].Science Technology and Engineering,2009,9(24):7430-7435

[8] 王月英.基于MPI的三維波動方程有限元法并行正演模擬[J].石油物探,2009,48(3)：221-243

Wang Y Y.3-D Finite dimensional wave equation forward simulation based on MPI parallelmodeling[J].Geophysical Prospecting for Petroleum,2009,48(3)：221-243

[9] Patera A T.A spectral element method for fluid dynamics:laminar flow in a channel expansion[J].Journal of Computational Acoustics,1994,2(4)：371-422

[10] Komatitsch D.Spectral and spectral-element methods for the 2D and 3D elasto-dynamics equations in heterogeneous media[D].Paris:Insti-tut de Physique du Globe,1997

[11] Komatitsch D,Vilotte J P.The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures[J].Bullein of the Seismolgical Society of America,1998,88(2):368-392

[12] Tromp J,Komatitsch D.The spectral element method:spectral-element and adjoint methods in seismology[J].Communications in Computational Physics,2008,3(1):1-32

[13] Lee S J,Komatitsch D.Effects of topography on seismic-wave propagation:an example from Northern Taiwan[J].Bulletin of the Seismological Society of America,2009,99(1):314-325

[14] Lee S J,Chan Y C,Komatitsch D.Effects of realistic surface topography on seismic ground motion in the Yangminshan region of Taiwan based upon the spectral-element method and LiDAR DTM[J].Bulletin of the Seismological Society of America,2009,99(2):681-693

[15] Komatitsch D,Erlebacher G,G?ddeke G,et al.High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster[J].Journal of Computational Physics,2010,229(20):7692-7714

[16] Komatitsch D,Tromp J.Introduction to the spectral element method for three-dimensional seismic wave propagation[J].Geophysical Journal International,1999,139(3):806-822

[17] Komatitsch D,Tromp J.Spectral-element simulations of wave propagation in porous media[J].Geophysical Journal International,2008,175(1):301-345

[18] 王秀明,Serianig G,林偉軍.利用譜元法計算彈性波場的若干理論問題[J].中國科學(xué)(G輯:物理學(xué) 力學(xué) 天文學(xué)),2007,37(1)：41-59

Wang X M,Serianig G,Lin W J.Using spectral element method to calculate some theory problems about elastic wave field[J].Science in China(Physics,Mechanics & Astronomy),2007,37(1)：41-59

[19] Cook W A,Oakes W R.Mapping methods for generating three-dimensional meshes[J].Computers in Mechanical Engineering,1982,67-72

[20] White D R.Automated hexahedral mesh generation by virtual decomposition[J].Meshing Roundtable,1985,165-176

[21] Baehmann P L,Wittchen S L,Shephard M S,et al.

Robust geometrically-based,automatic two-dimensional mesh generation[J].International Journal for Numerical Methods in Engineering,1987,24:1043-1078

[22] Talbert J A,Parkinson A R.Development of an automatic,two dimensional nite element mesh generator using quadrilateral elements and Bezier curve boundary dentition[J].International Journal for Numerical Methods in Engineering,1991,29:1551-1567

[23] Johnston B P,Sullivan J M,Kwasnik A.Automatic

conversion of triangularnite element meshes to quadrilateral elements[J].International Journal for Numerical Methods in Engineering,1991,31:67-84

[24] Lee C K,Lo S H.A new scheme for the generation of a graded quadrilateral mesh[J].Computers and Structures,1994,52(5):847-857

[25] Lo S H.A new mesh generation scheme for arbitrary planar domains[J].International Journal for Numerical Methods in Engineering,1985,21:1403-1426