面向無(wú)網(wǎng)格法的質(zhì)點(diǎn)生成算法

周詞揚(yáng) 鄭澎 馬智博 劉敏娟

摘要：

介紹一種面向無(wú)網(wǎng)格數(shù)值模擬方法的質(zhì)點(diǎn)生成算法。將四面體、三角形網(wǎng)格生成算法分別用于空間平面、曲面和實(shí)體模型，將網(wǎng)格單元的屬性賦予單元內(nèi)某一點(diǎn)作為質(zhì)點(diǎn)，并生成對(duì)應(yīng)質(zhì)點(diǎn)集。為研究不同網(wǎng)格生成算法和質(zhì)點(diǎn)生成算法對(duì)質(zhì)點(diǎn)集的影響，提出一種質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)，開(kāi)展對(duì)不同算法組合的質(zhì)量分析，得到網(wǎng)格生成算法和質(zhì)點(diǎn)生成算法中的最佳組合。

關(guān)鍵詞：

無(wú)網(wǎng)格法； 質(zhì)點(diǎn)生成； 質(zhì)點(diǎn)集； 網(wǎng)格生成； 評(píng)價(jià)標(biāo)準(zhǔn)

中圖分類號(hào)： TP391.9

文獻(xiàn)標(biāo)志碼： B

Particle generation algorithm for meshless method

ZHOU Ciyang1a， ZHENG Peng1， MA Zhibo2， LIU Minjuan1a

（1. a. Computer Application Institute， Mianyang 621900， Sichuan， China； b. Software Center for High Performance

Numerical Simulation， Beijing 100088， China， China Academy of Engineering Physics； 2. Beijing Institute of Applied

Physics and Computational Mathematics， Beijing 100094， China）

Abstract：

A particle generation method is introduced for meshless numerical simulation method. Etrahedron and triangle mesh generation algorithm is used for plane， surface and entity model； the properties of the mesh unit are given to a point inside the unit as particle and the corresponding particle sets are generated. In order to study the influence of different mesh generation and particle generation algorithms on particle sets， a quality evaluation standard is proposed. The quality analysis of different algorithm combinations is carried out， and the optimal combination of mesh generation algorithm and particle generation algorithm is obtained.

Key words：

meshless method； particle generation； particles set； mesh generation； evaluation standard

收稿日期： 2017-12-19

修回日期： 2018-01-23

基金項(xiàng)目： 國(guó)防基礎(chǔ)科研計(jì)劃（C1520110002）；中國(guó)工程物理研究院雙百基金（HT1609-01-KY124-019）

作者簡(jiǎn)介：

周詞揚(yáng)（1993—），男，湖北仙桃人，碩士研究生，研究方向?yàn)楦咝阅軘?shù)值模擬前處理，（E-mail）1095798605@qq.com

0 引 言

隨著高性能計(jì)算技術(shù)的發(fā)展，有限元法在核物理學(xué)、爆轟、材料學(xué)等領(lǐng)域發(fā)揮著越來(lái)越大的作用。但是，有限元法在處理材料變形、高速碰撞等問(wèn)題時(shí)，由于網(wǎng)格發(fā)生形變，重新劃分困難，導(dǎo)致計(jì)算速度下降，計(jì)算精度降低，嚴(yán)重時(shí)甚至出現(xiàn)計(jì)算崩潰的情況。采用無(wú)網(wǎng)格法解決這類問(wèn)題時(shí)，由于不依賴網(wǎng)格或只使用部分網(wǎng)格，可避免網(wǎng)格重新劃分等問(wèn)題，能有效解決有限元法所遇到的困難。然而，目前無(wú)網(wǎng)格法仍處于研究初級(jí)階段，其相應(yīng)的工程軟件較少，部分功能尚未實(shí)現(xiàn)。[1]無(wú)網(wǎng)格法的基本思想是使用一系列無(wú)網(wǎng)格節(jié)點(diǎn)排列，采用一種與權(quán)函數(shù)相關(guān)的近似，使某個(gè)域上的節(jié)點(diǎn)可以影響研究對(duì)象任意一點(diǎn)的力學(xué)特性，擺脫不連續(xù)性的束縛，保證求解的精度。[2]因此，在前處理過(guò)程中如何布置和劃分用于表達(dá)連續(xù)體信息的質(zhì)點(diǎn)成為重要研究方向。本文設(shè)計(jì)一個(gè)面向無(wú)網(wǎng)格法、用于實(shí)現(xiàn)幾何模型質(zhì)點(diǎn)集生成的算法，生成模型的四面體/三角形網(wǎng)格，由網(wǎng)格得到對(duì)應(yīng)的質(zhì)點(diǎn)集。討論質(zhì)點(diǎn)生成算法的不同實(shí)現(xiàn)方法并比較其對(duì)結(jié)果的影響，通過(guò)對(duì)同一個(gè)模型進(jìn)行測(cè)試，比較不同實(shí)現(xiàn)方案結(jié)果的差異。

采用無(wú)網(wǎng)格法中的光滑粒子流體動(dòng)力學(xué)（smoothed partical hydrodynamics， SPH）方法，通過(guò)四面體網(wǎng)格得到模型質(zhì)點(diǎn)。SPH方法的主要思想是任何一個(gè)連續(xù)系統(tǒng)可離散為一系列任意分布的質(zhì)點(diǎn)，所有有關(guān)這一系統(tǒng)的量都被認(rèn)為集中在這些質(zhì)點(diǎn)上。[3-4]生成四面體網(wǎng)格后，用網(wǎng)格單元內(nèi)的某點(diǎn)代替網(wǎng)格，網(wǎng)格單元的屬性值集中于該點(diǎn)，使生成的質(zhì)點(diǎn)集能用于SPH方法及其衍生方法的計(jì)算。

1 三維實(shí)體內(nèi)部質(zhì)點(diǎn)集生成算法

1.1 算法概述

用無(wú)網(wǎng)格數(shù)值模擬方法對(duì)模型進(jìn)行仿真計(jì)算時(shí)，需要將模型離散為若干質(zhì)點(diǎn)，用質(zhì)點(diǎn)代表周?chē)鷧^(qū)域，質(zhì)點(diǎn)的屬性即為周?chē)鷧^(qū)域的屬性。

生成模型的四面體網(wǎng)格，每個(gè)網(wǎng)格單元對(duì)應(yīng)一個(gè)質(zhì)點(diǎn)，質(zhì)點(diǎn)質(zhì)量為網(wǎng)格單元質(zhì)量。此外，算法需要一個(gè)標(biāo)準(zhǔn)量化網(wǎng)格質(zhì)點(diǎn)分布狀況，以便于比較不同生成算法之間的效果。

質(zhì)點(diǎn)生成算法流程為：

（1）生成模型的初始網(wǎng)格。

（2）細(xì)化初始網(wǎng)格到用戶指定程度。

（3）采用質(zhì)點(diǎn)生成算法求出細(xì)化網(wǎng)格每個(gè)網(wǎng)格單元的質(zhì)心，得到模型對(duì)應(yīng)的質(zhì)點(diǎn)集。

（4）對(duì)質(zhì)點(diǎn)集采用相應(yīng)評(píng)價(jià)標(biāo)準(zhǔn)得到具體數(shù)值，與設(shè)置的閾值比較，如果不符合要求則重復(fù)步驟（2）和（3），直到滿足要求為止。

在上述步驟中，網(wǎng)格生成算法和質(zhì)點(diǎn)生成算法均會(huì)對(duì)最終結(jié)果產(chǎn)生影響，不同的算法組合會(huì)得到不同的結(jié)果。本文討論每個(gè)算法的多種實(shí)現(xiàn)方式，比較不同方法所生成質(zhì)點(diǎn)集的質(zhì)量并確定一種最佳方法。

1.2 網(wǎng)格生成算法

本算法對(duì)四面體網(wǎng)格生成算法的要求是網(wǎng)格質(zhì)量良好、網(wǎng)格密度可控：網(wǎng)格質(zhì)量良好則四面體單元更接近于正四面體，相鄰質(zhì)點(diǎn)之間的距離方差更小，質(zhì)點(diǎn)分布更加均勻；網(wǎng)格密度可控是為了保證能細(xì)化網(wǎng)格到指定程度。

網(wǎng)格生成算法采用Delaunay算法保證網(wǎng)格質(zhì)量。四面體Delaunay網(wǎng)格的特性是對(duì)任意網(wǎng)格單元，其外接球內(nèi)不存在除了構(gòu)成網(wǎng)格單元頂點(diǎn)之外的點(diǎn)[5]，網(wǎng)格整體質(zhì)量良好，畸形網(wǎng)格數(shù)量較少，并且能保證限定條件在網(wǎng)格中的一致性，網(wǎng)格尺度和質(zhì)量可控。采用自適應(yīng)網(wǎng)格生成算法，通過(guò)限制網(wǎng)格單元的體積或相鄰點(diǎn)距離實(shí)現(xiàn)功能。

限制期望距離的Delaunay細(xì)化算法需要先指定一個(gè)距離L，使生成的網(wǎng)格邊長(zhǎng)在該值附近，然后為算法設(shè)置一個(gè)點(diǎn)生成規(guī)則，檢查每個(gè)網(wǎng)格單元，如果滿足點(diǎn)生成規(guī)則時(shí)則生成新的點(diǎn)并細(xì)化網(wǎng)格，直到所有單元均無(wú)法再細(xì)化為止。對(duì)應(yīng)的點(diǎn)生成規(guī)則[6]為：

（1）當(dāng)某條限定邊不存在于網(wǎng)格中或網(wǎng)格中某條邊被干涉，則將該邊的中點(diǎn)加入到四面體網(wǎng)格中并重構(gòu)網(wǎng)格。

（2）如果某個(gè)限定面不存在于網(wǎng)格中或網(wǎng)格中的三角面片被干涉，則將該面的外心加入到四面體網(wǎng)格中并重構(gòu)網(wǎng)格。如果出現(xiàn)規(guī)則（1）中的情況，則用R1的方法處理干涉邊。

（3）如果四面體的Radius-edge比大于閾值，或者外接球半徑大于L，則將四面體單元的外心加入到網(wǎng)格中并重構(gòu)網(wǎng)格，如果網(wǎng)格出現(xiàn)規(guī)則（1）和（2）中的情況，則分別按第1.1節(jié)中的（1）和（2）處理對(duì)應(yīng)線面。

在點(diǎn)生成規(guī)則中，三角面片被干涉是指以該三角面片外心為球心，外接圓半徑為球半徑的球含有除三角形頂點(diǎn)外的其他點(diǎn)。線段被干涉是指以線段中點(diǎn)為球心，線段長(zhǎng)度的1/2為半徑的球含有除線段端點(diǎn)以外的點(diǎn)。

基于上述點(diǎn)生成規(guī)則，網(wǎng)格細(xì)化算法流程如下：

（1）將初始四面體網(wǎng)格的網(wǎng)格單元存儲(chǔ)在隊(duì)列中；

（2）從隊(duì)列取出一個(gè)四面體單元，按照點(diǎn)生成規(guī)則生成新的點(diǎn)并重構(gòu)網(wǎng)格，刪除隊(duì)列中已不存在于網(wǎng)格中的四面體單元，將新生成的四面體單元加入隊(duì)列；（3）重復(fù)步驟（2），直到隊(duì)列為空。

限制最大體積的自適應(yīng)Delaunay算法需要指定一個(gè)體積V，生成的網(wǎng)格中所有網(wǎng)格單元體積均小于V。對(duì)于限制最大體積的網(wǎng)格生成算法，對(duì)應(yīng)的點(diǎn)生成規(guī)則與限制距離的規(guī)則大致相同，唯一的不同之處是點(diǎn)生成規(guī)則中第（3）條加入外心的條件是四面體的Radius-edge比大于閾值或者體積大于V。算法流程與限制距離的完全一樣。

1.3 質(zhì)點(diǎn)生成算法

質(zhì)點(diǎn)生成算法要求確保生成的質(zhì)點(diǎn)在網(wǎng)格內(nèi)。若某個(gè)網(wǎng)格單元對(duì)應(yīng)的質(zhì)點(diǎn)有在單元外的可能性，則對(duì)于位于邊界的單元會(huì)出現(xiàn)對(duì)應(yīng)質(zhì)點(diǎn)在幾何模型外的情況，同時(shí)內(nèi)部質(zhì)點(diǎn)分布不均勻。設(shè)計(jì)2種質(zhì)點(diǎn)生成方法，取四面體內(nèi)心和取四面體重心。

取四面體內(nèi)心法：已知四面體4頂點(diǎn)坐標(biāo)，求內(nèi)心坐標(biāo)。設(shè)4頂點(diǎn)坐標(biāo)分別為A1（x1，y1，z1），A2（x2，y2，z2），A3（x3，y3，z3），A4（x4，y4，z4），內(nèi)心坐標(biāo)設(shè)為（xa，ya，za），則求解公式[7]為

xa=4i=1（sixi）/4i=1xi

ya=4i=1（siyi）/4i=1yi

za=4i=1（sizi）/4i=1zi

（1）

式中：si為點(diǎn)Ai所對(duì)側(cè)面的面積。

取四面體重心法：已知四面體4頂點(diǎn)坐標(biāo)，求重心坐標(biāo)。設(shè)重心坐標(biāo)為（xg，yg，zg），則求解公式為

xg=4i=1xi/4

yg=4i=1yi/4

zg=4i=1zi/4

（2）

1.4 質(zhì)點(diǎn)集質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)

對(duì)于生成的質(zhì)點(diǎn)集，質(zhì)點(diǎn)分布越均勻，后續(xù)計(jì)算中計(jì)算精度和計(jì)算效率越高，計(jì)算出錯(cuò)的可能性越低。同時(shí)，無(wú)網(wǎng)格法用質(zhì)點(diǎn)代表對(duì)應(yīng)的四面體，每個(gè)質(zhì)點(diǎn)的質(zhì)量為對(duì)應(yīng)四面體的質(zhì)量，質(zhì)點(diǎn)的質(zhì)量值越集中，說(shuō)明模型劃分得越均勻，有助于提高后續(xù)計(jì)算精度。假設(shè)模型密度恒定，則可以用體積值代替質(zhì)量。因此，評(píng)價(jià)質(zhì)點(diǎn)集的質(zhì)量可從質(zhì)點(diǎn)分布質(zhì)量和質(zhì)點(diǎn)所在單元體積兩方面進(jìn)行。

為了量化質(zhì)點(diǎn)集分布質(zhì)量，生成質(zhì)點(diǎn)集的Delaunay網(wǎng)格，將Delaunay網(wǎng)格質(zhì)量作為質(zhì)點(diǎn)集質(zhì)量。由于Delaunay網(wǎng)格與Voronoi圖為對(duì)偶關(guān)系，Voronoi圖中每個(gè)空間單元均為對(duì)應(yīng)點(diǎn)的鄰域，而鄰域共面的兩點(diǎn)在Delaunay網(wǎng)格中是相連的，因此Delaunay網(wǎng)格中線段兩端點(diǎn)為空間相鄰點(diǎn)，因此該網(wǎng)格的質(zhì)量可表示質(zhì)點(diǎn)集的分布質(zhì)量。

網(wǎng)格的質(zhì)量可以用能量函數(shù)[8]、平均曲率[9]等方式表示，本文用網(wǎng)格單元的平均質(zhì)量表示網(wǎng)格整體質(zhì)量。每個(gè)網(wǎng)格單元的質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)如下。

Radius-edge ratio標(biāo)準(zhǔn)公式[10]為

d=r/l （3）

式中：r為四面體外接球半徑；l為四面體最短邊長(zhǎng)。

ParaQ標(biāo)準(zhǔn)公式[11]為

Q=CdV/（1i

式中：V為四面體體積；lij為端點(diǎn)是頂點(diǎn)i、j的線段長(zhǎng)度；Cd為使正四面體的質(zhì)量度量值取1而采用的比例因數(shù)，取Cd=1 832.820 8。

ParaV標(biāo)準(zhǔn)公式[11]為

Q=723V/（1i

l2ij）1.5

（5）

式中：V為四面體體積；lij為端點(diǎn)是頂點(diǎn)i、j的線段長(zhǎng)度。

每個(gè)四面體單元對(duì)應(yīng)質(zhì)點(diǎn)的屬性為該四面體單元體積，單元體積越集中，屬性值越集中，后續(xù)計(jì)算的精度和效率越高。用最大體積與最小體積之比衡量質(zhì)點(diǎn)集屬性的分布狀況。體積比計(jì)算公式為

σ=Vmax/Vmin （6）

式中：σ為體積比；Vmax為最大四面體單元體積；Vmin為最小四面體單元體積。

得到質(zhì)點(diǎn)集后，綜合對(duì)比3個(gè)網(wǎng)格評(píng)價(jià)標(biāo)準(zhǔn)，比較不同方法生成的質(zhì)點(diǎn)集的分布質(zhì)量，并討論質(zhì)點(diǎn)集的體積比。

2 空間面片質(zhì)點(diǎn)生成

在無(wú)網(wǎng)格數(shù)值模擬方法中，常常需要離散空間平面或曲面。得到的質(zhì)點(diǎn)集一般用于作為邊界條件求解碰撞發(fā)生時(shí)間、狀態(tài)等；除此之外，還可能用于求解該平面或曲面的物理屬性，例如求解材料變形時(shí)表面材料的狀態(tài)和形狀。對(duì)于前者，離散后的質(zhì)點(diǎn)必須在平面或曲面上，質(zhì)點(diǎn)沒(méi)有屬性；對(duì)于后者，質(zhì)點(diǎn)與平面或曲面之間允許有一定距離，質(zhì)點(diǎn)可看作將所在微元面集中于該點(diǎn)，該點(diǎn)屬性代表所在微元面屬性。

空間面片質(zhì)點(diǎn)生成算法輸入為一組空間平面或曲面，輸出為一組質(zhì)點(diǎn)集?？臻g平面的質(zhì)點(diǎn)生成算法的流程為：

（1）生成空間平面初始網(wǎng)格，將這些初始網(wǎng)格單元按所屬平面分類，以隊(duì)列的形式存儲(chǔ)。

（2）遍歷初始網(wǎng)格隊(duì)列，細(xì)化每一個(gè)外表面網(wǎng)格，采用質(zhì)點(diǎn)生成方法生成每個(gè)網(wǎng)格單元質(zhì)點(diǎn)并合并為質(zhì)點(diǎn)集。

（3）重復(fù)步驟（2）直到遍歷完隊(duì)列所有單元，整合所有質(zhì)點(diǎn)集為一個(gè)總質(zhì)點(diǎn)集并輸出。

空間曲面的質(zhì)點(diǎn)生成算法流程為：

（1）生成曲面的初始三角形網(wǎng)格。

（2）采用曲面逼近算法生成新點(diǎn)并細(xì)化網(wǎng)格

（3）根據(jù)實(shí)際應(yīng)用的不同取三角形網(wǎng)格的頂點(diǎn)、所有網(wǎng)格單元內(nèi)心或重心的集合作為質(zhì)點(diǎn)集。

曲面質(zhì)點(diǎn)生成算法根據(jù)選擇頂點(diǎn)方式不同又分為映射法和直接法2類，本文采用直接法中的網(wǎng)格前沿法。

2.1 平面三角形網(wǎng)格生成算法

由于在二維網(wǎng)格中Delaunay三角形網(wǎng)格具有空外接圓性和最優(yōu)性，因此采用Delaunay算法作為三角形網(wǎng)格生成算法。細(xì)化網(wǎng)格采用限制網(wǎng)格邊長(zhǎng)的方法確保能細(xì)化到指定程度。

細(xì)化算法質(zhì)點(diǎn)生成規(guī)則為：

（1）如果邊被干涉或限定邊不存在于網(wǎng)格中，則將邊的中點(diǎn)加入網(wǎng)格并重新劃分。

（2）如果三角形被干涉或三角形外接球半徑大于指定值L，則將三角形的外心加入網(wǎng)格并重新劃分；如果外心與某一條邊干涉，則用規(guī)則（1）的方式處理干涉邊。

細(xì)化算法檢查每個(gè)三角形單元，如果滿足上述規(guī)則，那么細(xì)化該單元。整個(gè)生成算法流程與四面體的相同，只是處理的網(wǎng)格單元為三角形單元。

2.2 曲面三角形網(wǎng)格生成算法

根據(jù)網(wǎng)格生成方法的不同，空間曲面網(wǎng)格生成算法可分為直接法和映射法。其中直接法又可以分為曲面分解法和網(wǎng)格前沿法2類。本文采用網(wǎng)格前沿法。網(wǎng)格前沿法[12]最早由LHNER提出，用于生成二維或三維非結(jié)構(gòu)網(wǎng)格，其中生成三角形網(wǎng)格時(shí)要求面為平面。

將網(wǎng)格前沿法用于空間曲面三角形網(wǎng)格生成時(shí)，由于網(wǎng)格與模型不重合，因此需要對(duì)算法進(jìn)行修改。LAU等[13]提出一種基于網(wǎng)格前沿法的空間曲面網(wǎng)格生成算法，該算法與平面網(wǎng)格不同之處是生成新邊界點(diǎn)的方式不同。

假設(shè)邊AB為下一個(gè)從前沿中刪除的邊，h為新網(wǎng)格單元預(yù)期高度，d為前沿中某條邊與曲面相切且與邊夾角為90°的向量，不同前沿點(diǎn)的d值不同，則新點(diǎn)坐標(biāo)為

Cz=M+hd （7）

式中：M為前沿邊隊(duì)列中某一個(gè)端點(diǎn)。選擇的端點(diǎn)不同，得到的新點(diǎn)不同。比較不同點(diǎn)與邊AB構(gòu)成的三角形的質(zhì)量，選擇質(zhì)量最好的網(wǎng)格單元作為新網(wǎng)格單元，同時(shí)更新網(wǎng)格前沿和網(wǎng)格單元的隊(duì)列。

2.3 質(zhì)點(diǎn)生成算法

質(zhì)點(diǎn)生成算法輸入一個(gè)三角形網(wǎng)格，輸出一組質(zhì)點(diǎn)。一般情況下，每個(gè)網(wǎng)格單元對(duì)應(yīng)一個(gè)質(zhì)點(diǎn)，計(jì)算時(shí)該質(zhì)點(diǎn)代表對(duì)應(yīng)單元。質(zhì)點(diǎn)生成方法與三維質(zhì)點(diǎn)生成算法相似：如果質(zhì)點(diǎn)有屬性，則取三角形的內(nèi)心或三角形的重心；如果質(zhì)點(diǎn)沒(méi)有屬性值，則可以取三角形網(wǎng)格的頂點(diǎn)。

已知三角形三頂點(diǎn)坐標(biāo)為A1（x1，y1，z1）、A2（x2，y2，z2）、A3（x3，y3，z3），設(shè)內(nèi)心的坐標(biāo)為Aa（xa，ya，za），重心的坐標(biāo)為Ag（xg，yg，zg），則坐標(biāo)值為

xg=3i=1xi/3

yg=3i=1yi/3

zg=3i=1zi/3

xa=（ax1+bx2+cx3）/（a+b+c）

ya=（ay1+by2+cy3）/（a+b+c）

za=（az1+bz2+cz3）/（a+b+c）

（8）

式中：a的值為|BC|；b的值為|AC|；c的值為|AB|。

2.4 質(zhì)點(diǎn)集質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)

二維質(zhì)點(diǎn)集的評(píng)價(jià)標(biāo)準(zhǔn)與三維類似，不同之處是用于質(zhì)量評(píng)價(jià)的網(wǎng)格是由模型表面生成的三角形網(wǎng)格，網(wǎng)格質(zhì)量為所有網(wǎng)格單元的平均質(zhì)量，將該值作為質(zhì)點(diǎn)集的質(zhì)量。該方法得到的質(zhì)量能用于比較不同網(wǎng)格生成方法得到質(zhì)點(diǎn)集之間分布質(zhì)量，無(wú)法用于比較不同質(zhì)點(diǎn)生成方式得到質(zhì)點(diǎn)集的質(zhì)量。網(wǎng)格單元的質(zhì)量有以下幾種標(biāo)準(zhǔn)。

Radius-edge ratio標(biāo)準(zhǔn)公式為

d=R/l （9）

式中：R為三角形外接圓半徑；l為三角形最短邊長(zhǎng)。

縱橫比標(biāo)準(zhǔn)公式為

d=R/r （10）

式中：r為三角形內(nèi)切圓半徑。

比較不同質(zhì)點(diǎn)集質(zhì)量時(shí)，使用上述標(biāo)準(zhǔn)得到質(zhì)點(diǎn)集的多個(gè)質(zhì)量值，綜合比較這些數(shù)值，判斷質(zhì)點(diǎn)集質(zhì)量的優(yōu)劣。

3 實(shí)驗(yàn)驗(yàn)證

3.1 三維實(shí)體質(zhì)點(diǎn)生成算法

爆炸裝置模型見(jiàn)圖1，該模型共5個(gè)部件，每個(gè)部件的內(nèi)容物不同。無(wú)網(wǎng)格算法模擬裝置的爆炸過(guò)程，判斷裝置從最上方起爆，爆轟波能否到達(dá)最下面。設(shè)半圓柱最上方部件的頂面和側(cè)面為固壁。

采用限制最大距離的自適應(yīng)Delaunay算法時(shí)，期望距離設(shè)為0.000 23 m，模型質(zhì)點(diǎn)生成結(jié)果見(jiàn)圖2。采用限制最大體積的自適應(yīng)Delaunay算法時(shí)，設(shè)置最大體積為1 mm3，得到的質(zhì)點(diǎn)集見(jiàn)圖3。

圖 1 爆炸裝置模型

圖 2 限制最大距離的自適應(yīng)算法模型質(zhì)點(diǎn)集

圖 3

限制最大體積的自適應(yīng)算法模型質(zhì)點(diǎn)集

2×2×12長(zhǎng)方體質(zhì)點(diǎn)集生成質(zhì)量見(jiàn)表1，按照上文介紹的3種質(zhì)量評(píng)價(jià)標(biāo)準(zhǔn)進(jìn)行質(zhì)量分析。由此可知：內(nèi)心生成的質(zhì)點(diǎn)集分布質(zhì)量?jī)?yōu)于重心的；質(zhì)點(diǎn)集生成的網(wǎng)格中，四面體單元整體質(zhì)量?jī)?yōu)于后者，但部分四面體單元質(zhì)量較差，從而導(dǎo)致在用Radius-edge標(biāo)準(zhǔn)判斷質(zhì)量時(shí)數(shù)值過(guò)大。

表 1 2×2×12長(zhǎng)方體質(zhì)點(diǎn)集生成質(zhì)量

整理半圓柱模型最上方部件網(wǎng)格的各網(wǎng)格單元體積分布，見(jiàn)圖4。由此可知，原四面體各網(wǎng)格單元體積取值范圍很大，最大值為1.750E-11，最小值為1.849E-13，最大值約為最小值的95倍。說(shuō)明四面體網(wǎng)格生成算法得到的網(wǎng)格單元的大小不均勻，需要改進(jìn)。

圖 4 四面體網(wǎng)格單元體積分布

3.2 表面質(zhì)點(diǎn)生成

生成固壁對(duì)應(yīng)的質(zhì)點(diǎn)集時(shí)，質(zhì)點(diǎn)取對(duì)應(yīng)三角形單元的重心，固壁質(zhì)點(diǎn)見(jiàn)圖5。

圖 5 固壁質(zhì)點(diǎn)

生成整個(gè)模型的表面質(zhì)點(diǎn)，質(zhì)點(diǎn)取網(wǎng)格單元的節(jié)點(diǎn)，見(jiàn)圖6。

圖 6 整個(gè)模型的表面質(zhì)點(diǎn)

對(duì)固壁的2個(gè)表面進(jìn)行質(zhì)量分析，結(jié)果見(jiàn)表2。由此可知：由于縱橫比接近于2，因此固壁上的三角形網(wǎng)格中每個(gè)網(wǎng)格單元形狀接近正三角形，網(wǎng)格質(zhì)量良好。

表 2 模型表面質(zhì)點(diǎn)生成質(zhì)量

固壁對(duì)應(yīng)的三角形網(wǎng)格單元面積見(jiàn)圖7。由此可知，三角形面積的最大值為7.22E-8，最小值為3.00E-8，最大值約為最小值的2.4倍。結(jié)果表明：對(duì)固壁進(jìn)行三角形劃分時(shí)得到的網(wǎng)格單元面積值集中，有助于后續(xù)計(jì)算保持高精度。

圖 7 固壁對(duì)應(yīng)的三角形網(wǎng)格單元面積

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

對(duì)多個(gè)模型生成對(duì)應(yīng)的內(nèi)部質(zhì)點(diǎn)集和表面質(zhì)點(diǎn)集。由分析結(jié)果可知：網(wǎng)格生成算法為限制距離的自適應(yīng)Delaunay算法時(shí)，生成的質(zhì)點(diǎn)分布質(zhì)量?jī)?yōu)于限制最大體積的自適應(yīng)Delaunay算法；質(zhì)點(diǎn)生成方式中，取內(nèi)心的方式得到的質(zhì)點(diǎn)集整體分布比取重心方式的更加均勻，但局部地區(qū)分布狀況不如后者。目前的算法得到的質(zhì)點(diǎn)集仍存在問(wèn)題，當(dāng)給每個(gè)質(zhì)點(diǎn)賦予對(duì)應(yīng)四面體體積值時(shí)，數(shù)值分布范圍過(guò)大，數(shù)值分布比較發(fā)散，采用適當(dāng)?shù)膬?yōu)化算法優(yōu)化由模型得到的四面體網(wǎng)格，使網(wǎng)格單元體積集中于某個(gè)值附近。表面質(zhì)點(diǎn)集的質(zhì)點(diǎn)分布均勻，屬性值集中，能用于實(shí)際SPH無(wú)網(wǎng)格法的計(jì)算中。

參考文獻(xiàn)：

[1] 王宇新，孫明，張建臣.無(wú)網(wǎng)格MPM法三維前處理系統(tǒng)設(shè)計(jì)[J]. 計(jì)算力學(xué)學(xué)報(bào)， 2008， 25（3）： 392-396.

[2] 劉天祥，劉更，朱均，等.無(wú)網(wǎng)格法的研究進(jìn)展[J]. 機(jī)械工程學(xué)報(bào)， 2002， 38（5）： 7-12.

[3] 王宇新，陳震，張洪武，等.多層抗暴結(jié)構(gòu)沖擊響應(yīng)無(wú)網(wǎng)格MPM法分析[J]. 工程力學(xué)， 2007， 24（12）： 186-192.

[4] 顧元通，丁樺.無(wú)網(wǎng)格法及其最新進(jìn)展[J]. 力學(xué)

進(jìn)展， 2005， 35（3）： 323-337.

[5] 武曉波，王世新，肖春生.Delaunay三角網(wǎng)的生成算法研究[J]. 測(cè)繪學(xué)報(bào)， 1999， 28（1）： 28-35.

[6] SI H. On refinement of constrained Delaunay Tetrahedralizations[C]//Proceedings of the 15th International Meshing Roundtable Conference. Birmingham， USA， 2006： 509-528.

[7] 王曉鳳，李永利.四面體的內(nèi)心坐標(biāo)公式及其應(yīng)用[J].平頂山工學(xué)院學(xué)報(bào)， 2004， 13（4）： 75-78.

[8] HOPPE H， DEROSE T， DUCHAMP T， et al. Mesh optimization[J]. Conference on Computer Graphics & Interactive Techniques， 1993， 27： 19-26.

[9] KARBACHER S， HAEUSLER G. New approach for the modeling and smoothing of scattered 3D data[R]. San Jose， USA： International Society for Optical Engineering， 1998.

[10] 李海生.Delaunay三角剖分理論及可視化應(yīng)用研究[M]. 哈爾濱： 哈爾濱工業(yè)大學(xué)出版社， 2010： 66-67.

[11] 聶春戈，劉劍飛，孫樹(shù)立.四面體網(wǎng)格質(zhì)量度量準(zhǔn)則的研究[J]. 計(jì)算力學(xué)學(xué)報(bào)， 2003， 20（5）： 579-582.

[12] LHNER R， PARIKH P. Generation of three-dimensional unstructured grids by the advancing-front method[J]. International Journal for Numerical Methods in Fluids， 1988， 8（10）： 1135-1149.

[13] LAU T S， LO S H. Finite element mesh generation over analytical curved surfaces[J].Computer & Structures， 1996， 59（2）： 301-309.