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

沙鑫池,聶玉峰,張偉偉

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

1 引言

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

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

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

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

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

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

2.1 三維泡泡布點(diǎn)方法的輸入信息

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

節(jié)點(diǎn)的理想間隔控制函數(shù)由用戶根據(jù)誤差估計(jì)定義，比如在結(jié)構(gòu)分析時(shí)應(yīng)力集中區(qū)域的間隔函數(shù)要適當(dāng)小，其它部分可以適當(dāng)?shù)姆糯?，這樣可以減少總體誤差．

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

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

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

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

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

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

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

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

2.4 計(jì)算泡泡的新位置和速度

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

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

泡泡的運(yùn)動遵循牛頓第二定律，運(yùn)動方程如下

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

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

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

不同于泡泡網(wǎng)格化方法利用網(wǎng)格拓?fù)浯_定鄰近泡泡，三維區(qū)域泡泡布點(diǎn)方法中節(jié)點(diǎn)的鄰接泡泡鏈表使用節(jié)點(diǎn)初始布置中的三維均勻“桶”結(jié)構(gòu)來建立，“桶”的邊長取為理想間隔函數(shù)最大值的1.7倍．每個桶都有一個節(jié)點(diǎn)列表，包含這個桶中的所有節(jié)點(diǎn)，在當(dāng)前節(jié)點(diǎn)進(jìn)行搜索時(shí)，只需要搜索該節(jié)點(diǎn)所在的桶和相鄰?fù)埃@樣就可以大大減少搜索節(jié)點(diǎn)時(shí)所用的時(shí)間．之后，需要判斷搜索到的泡泡是否就是當(dāng)前泡泡的鄰近泡泡．若是，則加入到鄰接鏈表中，反之則繼續(xù)搜索，直到節(jié)點(diǎn)所在的桶和相鄰?fù)八阉魍隇橹梗?/p>

具體方法如下：

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

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

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

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

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

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

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

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

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

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

圖1:泡泡的交疊率

3 節(jié)點(diǎn)集質(zhì)量的評價(jià)

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

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

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

4 數(shù)值算例

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

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

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

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

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

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

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

5 結(jié)論

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

參考文獻(xiàn)：

[1]聶玉峰，劉瑩.非結(jié)構(gòu)網(wǎng)格布點(diǎn)方法研究進(jìn)展[J].計(jì)算機(jī)工程與應(yīng)用，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]劉瑩，聶玉峰.泡泡布點(diǎn)方法及其并行性[J].計(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].計(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]齊楠，聶玉峰，張偉偉.快速泡泡布點(diǎn)方法[J].計(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é)點(diǎn)的局部網(wǎng)格生成并行算法[J].西北工業(yè)大學(xué)學(xué)報(bào)，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ì)量點(diǎn)集的快速節(jié)點(diǎn)生成算法[J].計(jì)算力學(xué)學(xué)報(bào)，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]王建華，徐強(qiáng)勛，張銳.任意形狀三維物體的Delaunay網(wǎng)格生成算法[J].巖石力學(xué)與工程學(xué)報(bào)，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ì)量度量標(biāo)準(zhǔn)的研究[J].計(jì)算力學(xué)學(xué)報(bào)，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