三維泡泡布點方法及網(wǎng)格生成?

沙鑫池,聶玉峰,張偉偉

(西北工業(yè)大學理學院應用數(shù)學系，西安 710129)

1 引言

有限元方法作為一種在科學研究和工程計算領域重要的數(shù)值模擬方法，需要節(jié)點布置作為其前處理的工作．目前已經(jīng)有一些高效的節(jié)點布置方法[1]，其中泡泡布點方法由于其人工干預少、計算精度高等特點在近幾年得到了深入的研究．

泡泡布點方法[2]是借鑒泡泡網(wǎng)格化方法[3]和分子動力學方法的思想和處理技術提出的一種單純布點方法．該方法是將區(qū)域中的點看做有相互作用力的泡泡，通過泡泡之間的相互作用力和阻尼力移動點的位置，同時結合分子動力學中更新泡泡鄰近鏈表的方法，獲得新時刻的臨近泡泡．泡泡布點法與泡泡網(wǎng)格化方法相比，可以避免不斷的網(wǎng)格重劃分，具有較高的效率，與一般的節(jié)點布置方法相比生成的節(jié)點質(zhì)量較高．

目前，泡泡布點方法在二維區(qū)域和雙參數(shù)曲面區(qū)域[4]已經(jīng)得到了很好的應用．二維區(qū)域和雙參數(shù)曲面區(qū)域的泡泡布點方法對區(qū)域有很好的適應性，生成的網(wǎng)格質(zhì)量高并且具有潛在的并行性，但是計算速度慢．針對這一缺點，通過泡泡運動方程中粘性系數(shù)的變化，改變傳統(tǒng)泡泡布點的數(shù)值求解方法和對交疊率設置閥值等一系列操作得到快速泡泡布點方法[5]．

對于三維區(qū)域，本文首先在邊界上運用雙參數(shù)曲面的泡泡布點方法進行節(jié)點布置，然后在區(qū)域內(nèi)部借鑒泡泡布點方法的思想進行節(jié)點布置，最后對生成的節(jié)點運用基于節(jié)點的局部網(wǎng)格生成算法[6-8]進行網(wǎng)格剖分．

2 三維區(qū)域泡泡布點方法

為了在三維區(qū)域得到高質(zhì)量的節(jié)點集，運用泡泡布點方法分別在三維區(qū)域的邊界面和內(nèi)部先后布點．對于區(qū)域的邊界部分，運用基于黎曼度量的雙參數(shù)曲面泡泡布點方法[4]進行布點．對于區(qū)域內(nèi)部，可以借鑒二維平面區(qū)域的泡泡布點方法[2]，將區(qū)域中的泡泡看做位于一個力場中．首先進行初始布點，然后將泡泡間的相互作用力和粘性阻尼力代入經(jīng)典運動方程，并運用數(shù)值方法進行迭代計算，得到泡泡的平衡位置，最后通過計算鄰近泡泡交疊率的大小來進行增點和刪點操作，得到符合間隔控制函數(shù)要求的高質(zhì)量節(jié)點集．下面給出具體實現(xiàn)過程．

2.1 三維泡泡布點方法的輸入信息

泡泡布點方法的輸入信息包括求解區(qū)域的幾何邊界和節(jié)點的理想間隔控制函數(shù)．幾何邊界利用雙參數(shù)曲面函數(shù)表示．

節(jié)點的理想間隔控制函數(shù)由用戶根據(jù)誤差估計定義，比如在結構分析時應力集中區(qū)域的間隔函數(shù)要適當小，其它部分可以適當?shù)姆糯螅@樣可以減少總體誤差．

2.2 區(qū)域邊界面的節(jié)點布置

運用雙參數(shù)曲面泡泡布點方法在區(qū)域邊界面上進行節(jié)點布置[4]．該方法利用三維曲面上的黎曼度量矩陣和雙參數(shù)映射計算出對應參數(shù)域上的尺寸控制矩陣，然后在二維參數(shù)域上運用泡泡布點方法進行節(jié)點布置，根據(jù)雙參數(shù)映射將參數(shù)域上的節(jié)點和拓撲結構映射回空間三維曲面．

2.3 區(qū)域內(nèi)部泡泡的初始布置

本文運用三維均勻“桶”結構來進行區(qū)域內(nèi)部泡泡的初始布置．主要思想是首先將計算區(qū)域用一個正六面體包圍起來，然后將包含計算區(qū)域的正六面體劃分為一系列規(guī)則的子區(qū)域，這些規(guī)則的子區(qū)域被稱為“桶”，“桶”的中心位置即為泡泡的初始位置．對于均勻布點即理想間隔控制函數(shù)為常數(shù)的情況，桶的邊長取為理想間隔控制函數(shù)；對于非均勻布點，桶的邊長取為理想間隔控制函數(shù)的最大值和最小值和的一半．對于已知的計算區(qū)域，只要找出區(qū)域中三個坐標最大的點和三個坐標最小的點就能得到每一個“桶”的中心節(jié)點．針對“桶”的中心節(jié)點可分為如下三種情況：

1)IN，桶的中心節(jié)點在計算區(qū)域內(nèi)；

2)OUT，桶的中心節(jié)點在計算區(qū)域外；

3)BOUNDARY，桶的中心節(jié)點在計算區(qū)域的邊界上．

對于第一種情況，將桶的中心節(jié)點保留，作為計算區(qū)域內(nèi)的初始泡泡；對于后兩種情況，則將桶的中心節(jié)點刪除，不予考慮．這樣得到的均勻“桶”中心節(jié)點即為初始布置的泡泡節(jié)點．

2.4 計算泡泡的新位置和速度

三區(qū)域中泡泡間相互作用力[9]選取為

其中w為兩個泡泡中心的實際距離l與理想距離l0的比值，這里理想距離l0為兩個泡泡中心理想尺寸和的一半．

泡泡的運動遵循牛頓第二定律，運動方程如下

其中mi是第i個泡泡的質(zhì)量，ci是阻尼系數(shù)，xi表示第i個泡泡中心點的位置，fi(t)是系統(tǒng)中所有鄰近泡泡作用在第i個泡泡上的合力，其中ij表示第j個泡泡對第i個泡泡的作用力．對于此二階常微分方程組采用歐拉預估矯正方法求解．在迭代求解過程中，借鑒快速泡泡布點方法的策略，對粘性系數(shù)C做適當?shù)恼{(diào)整，這樣不僅可以使泡泡盡快運動到一個合適的位置，而且有利于進行下一輪模擬時確保泡泡分布質(zhì)量．這里特別指出，對于那些被移動到區(qū)域外面的泡泡要被移回到區(qū)域內(nèi)，并與邊界保持一定的距離，這里將距離取為從泡泡中心向區(qū)域邊界所做法線與邊界交點處理想間隔的0.6倍．

2.5 泡泡鄰接鏈表的建立和更新

在動態(tài)模擬時，需要計算每個泡泡所受的合力．由于泡泡之間的相互作用力是短程力，只有鄰近的泡泡之間才有相互的作用．因此，為每個泡泡定義一個鄰接鏈表來儲存鄰近泡泡的信息是十分必要的，這樣可以大量節(jié)約搜索時間，提高計算效率．當計算每個泡泡所受合力時，只需計算來自鄰接鏈表中泡泡的相互作用力．

不同于泡泡網(wǎng)格化方法利用網(wǎng)格拓撲確定鄰近泡泡，三維區(qū)域泡泡布點方法中節(jié)點的鄰接泡泡鏈表使用節(jié)點初始布置中的三維均勻“桶”結構來建立，“桶”的邊長取為理想間隔函數(shù)最大值的1.7倍．每個桶都有一個節(jié)點列表，包含這個桶中的所有節(jié)點，在當前節(jié)點進行搜索時，只需要搜索該節(jié)點所在的桶和相鄰桶，這樣就可以大大減少搜索節(jié)點時所用的時間．之后，需要判斷搜索到的泡泡是否就是當前泡泡的鄰近泡泡．若是，則加入到鄰接鏈表中，反之則繼續(xù)搜索，直到節(jié)點所在的桶和相鄰桶搜索完為止．

具體方法如下：

1)找到區(qū)域中三個坐標取最小的點：A=(m in x,m in y,m in z)和區(qū)域中三個坐標取最大的點：B=(max x,max y,max z)，設

2)桶的邊長取為bl=1.7D，其中D為理想間隔距離的最大值，則三維桶的個數(shù)為

平面一層桶的個數(shù)為

包含區(qū)域的正方體的邊長為

3)計算每個泡泡所屬桶的編號，并建立每個桶所包含的泡泡信息．對任意的泡泡，設其中心坐標為(x,y,z)，則其所在桶的序號為

4)對于每一個計算泡泡，對其所屬桶及鄰近桶內(nèi)的其他泡泡進行判斷，將落在以計算泡泡中心(x,y,z)為圓心，以1.7σ為半徑的球域內(nèi)的泡泡加入其鄰接Verlet鏈表中，其中σ為計算泡泡與其他泡泡中心理想尺寸之和的一半．

2.6 三維區(qū)域泡泡數(shù)目的自適應調(diào)整

動態(tài)模擬時，為了找到一個合適的三維泡泡集合來覆蓋整個三維區(qū)域，使得泡泡彼此之間具有最小的裂縫和重疊，需要對泡泡的數(shù)目進行自適應調(diào)整．其數(shù)目的調(diào)整需要使用交疊率，其計算公式[2]為

其中ri和rj表示泡泡i和泡泡j的半徑，lij表示泡泡i和泡泡j的實際距離，N是鄰近鏈表中的節(jié)點數(shù)．

在理想情況下，曲面上的標準交疊率為6，三維實體內(nèi)部點的標準交疊率是12．在交疊率較小的泡泡附近添加新泡泡并且在交疊率較大的泡泡附近刪除泡泡，從而達到自適應調(diào)整的目的．為避免重復添加和刪除，對已做過操作的泡泡的鄰近泡泡做標記，不再進行添加和刪除的判斷．圖1表示了三維情況下中心泡泡周圍的臨近泡泡情況，對于其中的第一種和第三種情況應該分別進行適當?shù)脑黾优菖莺蛣h除泡泡操作．

圖1:泡泡的交疊率

3 節(jié)點集質(zhì)量的評價

本文分別針對區(qū)域邊界和區(qū)域內(nèi)部的點集進行網(wǎng)格生成．對于區(qū)域邊界上的節(jié)點集，運用高質(zhì)量點集的快速局部網(wǎng)格生成算法(BLMG)[7]進行網(wǎng)格生成．該方法利用泡泡布點過程中生成的泡泡鄰接鏈表，并結合Delaunay三角剖分的外接圓準則，從中心節(jié)點的鄰接鏈表中剔除非衛(wèi)星點，可以快速的生成局部網(wǎng)格．因為三維區(qū)域邊界是雙參數(shù)曲面，所以本文首先在參數(shù)域內(nèi)對曲面上生成的點集運用BLMG算法進行網(wǎng)格生成然后再將參數(shù)域中的網(wǎng)格映射回計算區(qū)域．對于區(qū)域內(nèi)部的點集，直接運用Delaunay三角剖分[8]進行網(wǎng)格生成．

在網(wǎng)格生成結束之后，通過計算生成的Delaunay網(wǎng)格質(zhì)量來評價節(jié)點集的質(zhì)量．近20年來，人們從不同角度提出了各種各樣的四面體網(wǎng)格單元質(zhì)量度量標準，本文采用內(nèi)切球–外接球半徑度量法[10]，本方法通過計算內(nèi)切球–外接球半徑比來判斷四面體網(wǎng)格單元質(zhì)量，其計算公式為ρ=3r/R，其中r,R分別為四面體的內(nèi)切圓和外接圓半徑．

半徑比ρ越接近1說明四面體網(wǎng)格質(zhì)量越好，將區(qū)域內(nèi)所有的四面體網(wǎng)格質(zhì)量求出后采用統(tǒng)計方法來評估整個區(qū)域的節(jié)點集質(zhì)量．

4 數(shù)值算例

運用上述方法，在半徑為1.0的三維球體區(qū)域內(nèi)，分別進行均勻和非均勻的泡泡節(jié)點布置．

算例1對一個半徑是1.0的三維球體進行泡泡布點．間隔函數(shù)d(x,y,z)取為常數(shù)0.22，這里取粘性系數(shù)的初值c0=1，變化斜率[5]k1=0.04,k2=0.004．節(jié)點布置的效果見圖2(a)，經(jīng)過Delaunay網(wǎng)格剖分后表面網(wǎng)格的生成情況見圖2(b)．最終三維區(qū)域生成的節(jié)點數(shù)是715，網(wǎng)格單元數(shù)是3367，網(wǎng)格質(zhì)量采用內(nèi)切球–外接球半徑度量法計算，網(wǎng)格質(zhì)量分布的統(tǒng)計結果見表1，網(wǎng)格的平均質(zhì)量是0.9474．

圖2:區(qū)域均勻布點及網(wǎng)格化結果

表1:網(wǎng)格單元質(zhì)量分布統(tǒng)計信息

算例2對一個半徑是1.0的三維球體進行泡泡布點．間隔函數(shù)如下理想間隔距離隨位置發(fā)生變化，在曲面x2+y2+z2=0.09附近為加密區(qū)域，粘性系數(shù)的選取與算例1相同．運行后布置效果見圖3(a)，經(jīng)過Delaunay網(wǎng)格剖分后的表面網(wǎng)格生成情況見3(b)．最終三維區(qū)域內(nèi)生成的節(jié)點數(shù)是883，網(wǎng)格單元數(shù)是4415，四面體網(wǎng)格平均質(zhì)量是0.9423．

圖3:區(qū)域非均勻布點及網(wǎng)格化結果

可見，在給定三維區(qū)域內(nèi)，對于均勻和非均勻布點情況，利用泡泡布點方法生成的點集具有很好的結構，并且根據(jù)點集生成的Delaunay四面體網(wǎng)格質(zhì)量高，同時對于邊界附近的點也有非常好的效果．因此在三維區(qū)域內(nèi)泡泡布點方法可以得到很好的應用．

5 結論

泡泡布點方法在三維區(qū)域內(nèi)可以取得較好的布點效果．在對凸區(qū)域進行均勻或非均勻的節(jié)點布置時，生成的點集不僅滿足理想間隔控制函數(shù)的要求，并且具有很高的質(zhì)量．三維區(qū)域的泡泡布點方法是一種單純的布點方法，實施起來簡單，而且生成的點集既可直接用于無網(wǎng)格方法，也可以作為有限元網(wǎng)格生成的給定點集為網(wǎng)格生成做準備．

參考文獻：

[1]聶玉峰，劉瑩.非結構網(wǎng)格布點方法研究進展[J].計算機工程與應用，2008,44(32):35-40 Nie Y F,Liu Y.Survey of point placement for unstructured mesh generation[J].Computer Engineer and Application,2008,44(32):35-40

[2]劉瑩，聶玉峰.泡泡布點方法及其并行性[J].計算物理，2009,26(6):813-820 Liu Y,Nie Y F.Node placement method with bubble simulation and parallelism[J].Chinese Journal of Computational Physics,2009,26(6):813-820

[3]Shimada K,Gossard D C.Bubble mesh:automated triangular meshing of non-manifold geometry by sphere packing[C]//Proceedings of Solid Modeling Applications,Salt Lake City,1995:409-419

[4]張偉偉，聶玉峰，王磊.雙參數(shù)曲面的泡泡網(wǎng)格化方法[J].計算物理，2012,29(1):43-50 Zhang W W,Nie Y F,Wang L.Bubble meshing method for two-parametric surface[J].Chinese Journal of Computational Physics,2012,29(1):43-50

[5]齊楠，聶玉峰，張偉偉.快速泡泡布點方法[J].計算物理，2012,29(3):333-338 Qi N,Nie Y F,Zhang W W.A fast node placement method with bubble simulation[J].Chinese Journal of Computational Physics,2012,29(3):333-338

[6]聶玉峰，樊祥闊，常升，等.基于節(jié)點的局部網(wǎng)格生成并行算法[J].西北工業(yè)大學學報，2006,24(6):731-735 Nie Y F,Fan X K,Chang S,et al.A new and efficient node-based local mesh generation parallel algorithm[J].Journal of Northwestern Polytechnical University,2006,24(6):731-735

[7]陳蔚蔚，聶玉峰，張偉偉，等.高質(zhì)量點集的快速節(jié)點生成算法[J].計算力學學報，2012,29(5):704-709 Chen W W,Nie Y F,Zhang W W,et al.A fast local mesh generation method about high-quality node set[J].Chinese Journal of Computational Mechanics,2012,29(5):704-709

[8]王建華，徐強勛，張銳.任意形狀三維物體的Delaunay網(wǎng)格生成算法[J].巖石力學與工程學報，2003,22(5):717-722 Wang J H,Xu Q X,Zhang R.Delaunay algorithm and related procedure to generate the tetrahedron mesh for an object with arbitrary boundary[J].Chinese Journal of Rock Mechanics and Engineering,2003,22(5):717-722

[9]Dibben D C.3D mesh generation using bubble meshing for microwave applicators[J].IEEE Transactions on Magnetics,2000,36(4):1514-1518

[10]聶春戈，劉劍飛，孫樹立.四面體網(wǎng)格質(zhì)量度量標準的研究[J].計算力學學報，2003,20(5):579-582 Nie C G,Liu J F,Sun S L.Study on quality measures for tetrahedral mesh[J].Chinese Journal of Computational Mechanics,2003,20(5):579-582