周詞揚 鄭澎 馬智博 劉敏娟

摘要：

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

關鍵詞：

無網(wǎng)格法； 質點生成； 質點集； 網(wǎng)格生成； 評價標準

中圖分類號： TP391.9

文獻標志碼： 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

基金項目： 國防基礎科研計劃（C1520110002）；中國工程物理研究院雙百基金（HT1609-01-KY124-019）

作者簡介：

周詞揚（1993—），男，湖北仙桃人，碩士研究生，研究方向為高性能數(shù)值模擬前處理，（E-mail）1095798605@qq.com

0 引 言

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

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

1 三維實體內(nèi)部質點集生成算法

1.1 算法概述

用無網(wǎng)格數(shù)值模擬方法對模型進行仿真計算時，需要將模型離散為若干質點，用質點代表周圍區(qū)域，質點的屬性即為周圍區(qū)域的屬性。

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

質點生成算法流程為：

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

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

（3）采用質點生成算法求出細化網(wǎng)格每個網(wǎng)格單元的質心，得到模型對應的質點集。

（4）對質點集采用相應評價標準得到具體數(shù)值，與設置的閾值比較，如果不符合要求則重復步驟（2）和（3），直到滿足要求為止。

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 質點生成算法

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

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

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

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

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

（1）

式中：si為點Ai所對側面的面積。

取四面體重心法：已知四面體4頂點坐標，求重心坐標。設重心坐標為（xg，yg，zg），則求解公式為

xg=4i=1xi/4

yg=4i=1yi/4

zg=4i=1zi/4

（2）

1.4 質點集質量評價標準

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

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

網(wǎng)格的質量可以用能量函數(shù)[8]、平均曲率[9]等方式表示，本文用網(wǎng)格單元的平均質量表示網(wǎng)格整體質量。每個網(wǎng)格單元的質量評價標準如下。

Radius-edge ratio標準公式[10]為

d=r/l （3）

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

ParaQ標準公式[11]為

Q=CdV/（1i

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

ParaV標準公式[11]為

Q=723V/（1i

l2ij）1.5

（5）

式中：V為四面體體積；lij為端點是頂點i、j的線段長度。

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

σ=Vmax/Vmin （6）

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

得到質點集后，綜合對比3個網(wǎng)格評價標準，比較不同方法生成的質點集的分布質量，并討論質點集的體積比。

2 空間面片質點生成

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

空間面片質點生成算法輸入為一組空間平面或曲面，輸出為一組質點集。空間平面的質點生成算法的流程為：

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

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

（3）重復步驟（2）直到遍歷完隊列所有單元，整合所有質點集為一個總質點集并輸出。

空間曲面的質點生成算法流程為：

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

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

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

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

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

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

細化算法質點生成規(guī)則為：

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

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

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

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

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

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

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

Cz=M+hd （7）

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

2.3 質點生成算法

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

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

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 質點集質量評價標準

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

Radius-edge ratio標準公式為

d=R/l （9）

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

縱橫比標準公式為

d=R/r （10）

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

比較不同質點集質量時，使用上述標準得到質點集的多個質量值，綜合比較這些數(shù)值，判斷質點集質量的優(yōu)劣。

3 實驗驗證

3.1 三維實體質點生成算法

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

采用限制最大距離的自適應Delaunay算法時，期望距離設為0.000 23 m，模型質點生成結果見圖2。采用限制最大體積的自適應Delaunay算法時，設置最大體積為1 mm3，得到的質點集見圖3。

圖 1 爆炸裝置模型

圖 2 限制最大距離的自適應算法模型質點集

圖 3

限制最大體積的自適應算法模型質點集

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

表 1 2×2×12長方體質點集生成質量

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

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

3.2 表面質點生成

生成固壁對應的質點集時，質點取對應三角形單元的重心，固壁質點見圖5。

圖 5 固壁質點

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

圖 6 整個模型的表面質點

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

表 2 模型表面質點生成質量

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

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

4 結束語

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

參考文獻：

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

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

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

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

進展， 2005， 35（3）： 323-337.

[5] 武曉波，王世新，肖春生.Delaunay三角網(wǎng)的生成算法研究[J]. 測繪學報， 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)心坐標公式及其應用[J].平頂山工學院學報， 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三角剖分理論及可視化應用研究[M]. 哈爾濱： 哈爾濱工業(yè)大學出版社， 2010： 66-67.

[11] 聶春戈，劉劍飛，孫樹立.四面體網(wǎng)格質量度量準則的研究[J]. 計算力學學報， 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.