盛天爽*

?

一種適用于任意復(fù)雜結(jié)構(gòu)的曲六面體網(wǎng)格生成算法

盛天爽*

（南京師范大學(xué)附屬中學(xué)高二（4）班，江蘇南京，210094）

曲六面體網(wǎng)格自動生成是時域譜元法發(fā)展的瓶頸問題。本文采用基于體基函數(shù)的協(xié)變投影算法，將一個10節(jié)點(diǎn)曲四面體單元分解為四個20節(jié)點(diǎn)曲六面體單元，從而形成一種適用于任意復(fù)雜結(jié)構(gòu)的曲六面體網(wǎng)格生成算法。使用該算法能夠?qū)⑶拿骟w單元均勻分割，且確保新形成的曲六面體單元的翹曲最小。典型算例表明，使用該算法建立的曲六面體離散模型，完全滿足時域譜元法的計算精度要求。

網(wǎng)格自動生成；曲四面體；曲六面體；協(xié)變投影

引言

在計算電磁學(xué)領(lǐng)域，時域譜元法（SETD）是一種較為新型的數(shù)值仿真算法。這種算法不僅具有類似于FDTD算法的顯式求解格式，而且隨著算法中所使用的正交基函數(shù)階數(shù)的增加，其計算精度呈指數(shù)增長，即具有譜精度特性[1]。在時域譜元法中，通常使用具有20個節(jié)點(diǎn)的曲六面體單元作為基本離散單元，以確保仿真計算的譜精度性，但是受到六面體拓?fù)浣Y(jié)構(gòu)的限制，實(shí)現(xiàn)可靠、高質(zhì)量的曲六面體網(wǎng)格自動生成一直是時域譜元法發(fā)展的瓶頸問題。為了繞開這一問題，人們只能使用分區(qū)域離散的方式建模[2]，對復(fù)雜結(jié)構(gòu)使用曲四面體建模，其余部分采用曲六面體建模。但這樣的操作不僅增加了算法的復(fù)雜性，而且也破壞了算法所具有的全局顯式求解格式，降低了算法的效率。

事實(shí)上，人們很早就認(rèn)識到使用曲六面體建模具有空間逼近度高、計算精度高的特點(diǎn)，也設(shè)計了很多種方法來實(shí)現(xiàn)曲六面體網(wǎng)格的自動生成，但相比能適用于任意復(fù)雜結(jié)構(gòu)的曲四面體來說，實(shí)現(xiàn)曲六面體網(wǎng)格的自動生成還沒有通用的方法[3]。目前通常使用的方法主要有單元映射法[4]、柵格分解法[5]、四面體轉(zhuǎn)換法[6]、逐層推進(jìn)的波前法[7][8]等等。其中，單元映射法雖然最易實(shí)現(xiàn)，卻不適用于復(fù)雜三維結(jié)構(gòu)的離散。柵格分解法只能生成規(guī)則化網(wǎng)格，且邊界單元的質(zhì)量較差。逐層推進(jìn)的波前法生成的六面體網(wǎng)格質(zhì)量最優(yōu)，但其實(shí)現(xiàn)難度最大，易產(chǎn)生單元間的裂縫以及形狀尺寸不合理的單元。

四面體轉(zhuǎn)換法是目前最具有通用性的六面體自動生成算法，該算法首先使用直四面體建模，在此基礎(chǔ)上，再將每個直四面體單元分離為4個直六面體單元。但使用直四面體建模時，往往需要很大的剖分密度才能保證對復(fù)雜結(jié)構(gòu)的擬合精度，從而導(dǎo)致生成更多的直六面體單元，極大增加了未知量的個數(shù)。事實(shí)上，使用具有曲邊和曲面的高階離散單元可以更好的擬合復(fù)雜結(jié)構(gòu)，并且能避免過大的剖分密度。于是國內(nèi)有學(xué)者提出了采用從曲四面體單元中分離出六面體單元的算法[9]，該算法在區(qū)分曲面和直面的基礎(chǔ)上，使用面基函數(shù)，通過空間映射得到曲四面體的曲面形心坐標(biāo)，再根據(jù)常規(guī)的四面體轉(zhuǎn)換法，獲得4個準(zhǔn)曲六面體單元。該算法操作起來較為繁瑣，且形成的準(zhǔn)曲六面體單元中即包含曲面，又包含直面，并非常規(guī)的包含20個幾何節(jié)點(diǎn)的曲六面體單元，因此無法直接用于時域譜元法。

本文對該算法做了改進(jìn)，不必區(qū)分曲四面體中的曲面和直面，直接使用體基函數(shù)，通過諧變投影映射方式，在一個10節(jié)點(diǎn)曲四面體單元中獲得4個20節(jié)點(diǎn)曲六面體單元。

2 算法原理

笛卡兒坐標(biāo)系下的二次曲面四面體單元包含10個幾何特征節(jié)點(diǎn)，分別定義于四面體的4個頂點(diǎn)和6條棱邊的中點(diǎn)，通常稱為10節(jié)點(diǎn)曲四面體單元，如圖1(a)所示。

圖1 10節(jié)點(diǎn)曲四面體單元

圖2 參量空間的的直角四面體單元

根據(jù)協(xié)變投影原理，真實(shí)空間中的10節(jié)點(diǎn)曲四面體單元與參量空間的直角四面體單元之間的空間映射關(guān)系為

顯然，在參量空間中，直角四面體單元里任意一點(diǎn)的坐標(biāo)由相應(yīng)的4個參變量唯一確定。據(jù)此，可以將任意一個曲四面體單元投影為體坐標(biāo)系中的標(biāo)準(zhǔn)直角四面體單元，如圖2所示。其中，直角頂點(diǎn)位于坐標(biāo)系原點(diǎn)，三條直角棱邊分別位于三坐標(biāo)軸，邊長均為1。

上述形函數(shù)N的具體表達(dá)形式為

將式（3）帶入空間映射關(guān)系式（1），即可以獲得體坐標(biāo)系中直角四面體單元內(nèi)任意一點(diǎn)在相應(yīng)曲四面體單元中的笛卡爾坐標(biāo)?；谶@一原理，我們可以先在體坐標(biāo)系中將一個標(biāo)準(zhǔn)的直角四面體單元分裂為四個直六面體單元，如圖3所示。

圖4 20個幾何節(jié)點(diǎn)的曲六面體單元

其中，編號1~10的點(diǎn)對應(yīng)于曲四面體的10個幾何節(jié)點(diǎn)，編號11的點(diǎn)為體中心點(diǎn)，編號12~15的點(diǎn)分別為4個面的形心。將相應(yīng)各點(diǎn)連接起來，便可以在一個直四面體中獲得4個直六面體。表1列舉了將一個四面體單元分裂為4個六面體單元時，各六面體中8個頂點(diǎn)的序列分布。

表1 四個六面體的頂點(diǎn)編號

以上操作在體坐標(biāo)系中進(jìn)行，最后獲得4個直六面體單元。但曲六面體元的幾何節(jié)點(diǎn)有20個，分別位于8個頂點(diǎn)和12條棱邊的中點(diǎn)，如圖4所示。因此，我們還需要在體坐標(biāo)系中，針對新形成的4個直六面體單元，分別按照曲六面體20個幾何節(jié)點(diǎn)的分布規(guī)律，重新建立相應(yīng)幾何點(diǎn)的體坐標(biāo)。最后，再通過協(xié)變投影，在真實(shí)空間獲得4個曲六面體單元。

在本文介紹的算法中，將體坐標(biāo)系中的直六面體單元映射到真實(shí)空間時采用了協(xié)變投影，因此得到的曲六面體單元的翹曲最小。并且，在體坐標(biāo)系中定義4個六面體的分界點(diǎn)時采用的是體積變量，因此可以準(zhǔn)確控制曲四面體的均勻分割。例如，四面體的中心點(diǎn)（圖3中編號為11的點(diǎn)），可以直接定義為（1/4, 1/4, 1/4, 1/4）。通過協(xié)變投影，該點(diǎn)在相應(yīng)曲六面體中也必然位于體積均分的中心點(diǎn)位置。根據(jù)這一原則，我們給出圖3中15個幾何點(diǎn)的體坐標(biāo)定義，如表2所示。

表2 圖3中15個幾何點(diǎn)的體坐標(biāo)定義

根據(jù)表2可以確定體坐標(biāo)系中每個直六面體的8個頂點(diǎn)坐標(biāo)。在此基礎(chǔ)上，計算每條棱邊的中點(diǎn)坐標(biāo)，從而獲得所有20個節(jié)點(diǎn)的體坐標(biāo)。帶入式（1），即可求得相應(yīng)曲六面體單元的20個幾何節(jié)點(diǎn)的真實(shí)坐標(biāo)。

2 算例

使用協(xié)變投影獲得的曲六面體單元具有最小翹曲，因此可以確保使用該算法所獲得的曲六面體單元在用于時域譜元法時計算的準(zhǔn)確性，下面通過一個簡單的算例加以驗(yàn)證。一個半徑為1米的球形金屬諧振腔，利用商業(yè)電磁仿真軟件HFSS，得到三個最低諧振頻率分別為132.089MHz、186.439MHz、216.551MHz。在ANSYS軟件中建立模型，并且分別使用10節(jié)點(diǎn)的曲四面體單元和20節(jié)點(diǎn)的曲六面體單元對模型進(jìn)行離散。使用曲六面體離散時，需要先對球體進(jìn)行合理切割，以確保各子區(qū)域的空間拓?fù)浣Y(jié)構(gòu)滿足曲六面體離散的要求，剖分后獲得448個曲六面體單元。使用曲四面體離散時，對模型無需特殊處理，直接離散后獲得72個曲四面體單元。再根據(jù)本文算法，將這72個曲四面體離散為288個曲六面體，計算出每個曲六面體單元中20個幾何節(jié)點(diǎn)的坐標(biāo)，并對其進(jìn)行全局編號。依據(jù)兩種不同離散方式建立的模型，使用時域譜元法分別計算該球形諧振腔的諧振頻率，仿真結(jié)果如圖5所示。

圖5 半徑1米球形諧振腔的三個最低諧振頻率

圖中，灰色豎線為HFSS仿真得到的三個最低諧振頻率值；黑實(shí)線為使用曲四面體轉(zhuǎn)換曲六面體方法建模時，仿真所獲得的諧振頻率值；虛線為直接使用曲六面體建模時，仿真所獲得的諧振頻率值。顯然，使用本文算法所獲得的曲六面體離散模型，完全滿足時域譜元法的計算精度要求。

3 結(jié)束語

使用時域譜元法進(jìn)行電磁仿真時，實(shí)現(xiàn)復(fù)雜計算模型的曲六面體網(wǎng)格離散是個難題。本文對常規(guī)的“四面體轉(zhuǎn)換法”做了改進(jìn)，根據(jù)協(xié)變投影原理，采用體基函數(shù)映射的方式，在一個10節(jié)點(diǎn)曲四面體單元中獲得4個20節(jié)點(diǎn)曲六面體單元。實(shí)際操作中，由于將體坐標(biāo)系中的直六面體單元映射到真實(shí)空間時采用了協(xié)變投影，因此所獲得的曲六面體單元具有最小翹曲。同時，在體坐標(biāo)系中使用體積變量定義4個六面體的分界點(diǎn)，可以更準(zhǔn)確地控制曲四面體的均勻分割。從給出的典型算例也可以看出，由于曲四面體適用于曲面擬合，因此在不降低計算精度的前提下，可以極大減少建模時所需要的離散單元數(shù)量，并最終減少曲六面體單元的數(shù)量，從而減少未知量，加快仿真速度。該曲六面體網(wǎng)格自動生成算法具有很強(qiáng)的通用性，適于任意復(fù)雜結(jié)構(gòu)的建模。

[1] Joon-Ho Lee, and Qing H. Liu, A 3-D Spectral-Element Time-Domain Method for Electromagnetic Simulation, IEEE Trans. on M.T.T., 2007, 55(5): 983-991.

[2] Qiang Ren, Luis E. Tobón, Qingtao Sun, Qing Huo Liu, A New 3-D Nonspurious Discontinuous Galerkin Spectral Element Time-Domain (DG-SETD) Method for Maxwell’s Equations, IEEE Trans. on A.P., 2015, 63(6): 2585-2594.

[3] 呂軍, 王忠金, 王仲仁, 有限元六面體網(wǎng)格的典型生成方法及發(fā)展趨勢, 哈爾濱工業(yè)大學(xué)學(xué)報，2001, 33(4): 485-490.

[4] 蔣浩民, 劉潤廣, 王忠金, 等, 基于映射法的三維有限元網(wǎng)格自動劃分, 塑性工程學(xué)報, 1998 , 5(3): 27-31.

[5] TEKKAYAAE. Fully automatic simulation of bulk metal forming processes. Proc NUMIFORM' 98. Rotterdam: Netherlands, 1998.

[6] G Xie, JAH Ramaekers, Graded mesh generation and transformation. Finite Elements in Analysis and Design, 1994, 17(1): 41-55.

[7] TD Blacker, RJ Meyers, A 3-D hexahedral mesh generation algorithm. Engineering with Computers, 1993, 9(2): 83-93.

[8] 尚菲菲, 甘洋科, 郭宇飛, 劉劍飛, 基于片分割的六面體網(wǎng)格生成方法, 計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報, 2017, 29(4): 662-669.

[9] 左旭, 衛(wèi)原平, 陳軍, 阮雪榆, 三維六面體有限元網(wǎng)格自動劃分中的一種單元轉(zhuǎn)換優(yōu)化算法, 計算力學(xué)學(xué)報, 1999, 16(3): 343-348.

An Algorithm for Generating Curved Hexahedron Mesh for Arbitrary Complex Structures

SHENG Tianshuang*

(High School Affiliated to Nanjing Normal University, Grade 2, Class 4, Jiangsu Nanjing, 210094, China)

Automatic generation of curved hexahedron mesh is a bottleneck in the development of spectral-element time-domain method. In this paper, the covariant projection algorithm based on volume basis function is adopted to decompose a 10-node curved tetrahedron element into four 20-node curved hexahedron elements. Thus, a curved hexahedron mesh generation algorithm suitable for arbitrary complex structures is proposed. Using this algorithm, the curved tetrahedron element can be partitioned evenly and the minimal warping can be ensured for these new hexahedron elements. The typical example shows that the curved hexahedron model established by this algorithm fully meets the precision requirement of the spectral-element time-domain method.

Automatic mesh generation; Curved tetrahedron; Curved hexahedron; Covariant projection

TP39

A

10.19551/j.cnki.issn1672-9129.2018.01.010

1672-9129(2018)01-0022-03

盛天爽. 一種適用于任意復(fù)雜結(jié)構(gòu)的曲六面體網(wǎng)格生成算法[J]. 數(shù)碼設(shè)計, 2018, 7(1): 22-24.

SHENG Tianshuang. An Algorithm for Generating Curved Hexahedron Mesh for Arbitrary Complex Structures[J]. Peak Data Science, 2018, 7(1): 22-24.

2017-11-01；

2018-01-05。

盛天爽(2001-)，男，江蘇南京，南京師范大學(xué)附屬中學(xué)高二（4）班，學(xué)生。E-mail: shengyij@.edu.cn