付 宏，呂 游，徐 靜，黃 山，于建群

（吉林大學(xué)a.計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，長春 130012；b.生物與農(nóng)業(yè)工程學(xué)院，長春 130025）

非規(guī)則曲面的離散元法分析模型建模軟件

付 宏a，呂 游a，徐 靜a，黃 山a，于建群b

（吉林大學(xué)a.計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，長春 130012；b.生物與農(nóng)業(yè)工程學(xué)院，長春 130025）

在采用離散元法分析機(jī)械部件與顆粒材料接觸作用時，需要建立機(jī)械部件（邊界）的離散元法分析模型。分析可知，機(jī)械部件中與顆粒材料接觸作用的零件表面，存在不能用初等解析函數(shù)表達(dá)的非規(guī)則曲面。為此，采用推進(jìn)波前法（AFT：Advancing Front Technique）進(jìn)行非規(guī)則曲面網(wǎng)格劃分，把非規(guī)則曲面離散成三角形平面片的組合，同時添加運(yùn)動屬性和材料特性參數(shù)，由此建立非規(guī)則曲面邊界的離散元法分析模型。在對PRO／E軟件進(jìn)行二次開發(fā)的基礎(chǔ)上，研制了非規(guī)則曲面邊界建模軟件。通過實(shí)例驗(yàn)證，初步證明了基于AFT邊界建模方法和軟件的可行性，為復(fù)雜結(jié)構(gòu)機(jī)械部件工作過程的仿真分析奠定了基礎(chǔ)。

離散元法；邊界建模；非規(guī)則曲面；推進(jìn)波前法

0 引 言

在采用離散元法分析顆粒材料與機(jī)械部件的接觸作用時，需要建立機(jī)械部件（邊界）的離散元法分析模型［1］。分析可知，機(jī)械部件中與顆粒材料接觸作用的零件表面，一般可分為兩類：一類是能用初等解析函數(shù)表達(dá)的規(guī)則曲面，如平面、球面、柱面和錐面等；另一類是不能用初等解析函數(shù)表達(dá)的非規(guī)則曲面。

對于規(guī)則曲面，可采用基于圖元的邊界建模方法和實(shí)邊界與虛邊界的方法［2］，建立邊界的離散元法分析模型。對于非規(guī)則曲面，筆者采用推進(jìn)波前法（AFT：Advancing Front Technique）［3－8］進(jìn)行曲面網(wǎng)格劃分，把非規(guī)則曲面離散成三角形平面片的組合，同時添加運(yùn)動屬性和材料特性參數(shù)，由此建立非規(guī)則曲面邊界的離散元法分析模型［9］。

在對PRO／E軟件進(jìn)行二次開發(fā)的基礎(chǔ)上，研制了非規(guī)則曲面邊界的建模軟件，并與規(guī)則曲面邊界建模軟件集成，從而開發(fā)一種通用的離散元法邊界建模軟件［9］。通過實(shí)例驗(yàn)證，初步證明了基于AFT邊界建模方法和邊界建模軟件的可行性，為進(jìn)一步完善筆者提出的基于CAD模型的離散元法邊界建模方法［10］，以及采用離散元法分析復(fù)雜結(jié)構(gòu)機(jī)械部件的工作奠定了基礎(chǔ)。

1 AFT法的基本原理和實(shí)現(xiàn)

1.1 基本原理

如圖1所示，AFT法的基本原理為：首先由離散非規(guī)則曲面的邊界環(huán)，形成前沿線段的集合，然后由每條前沿線段確定一個新生成點(diǎn)，在新生成點(diǎn)和已有前沿線段點(diǎn)中搜索一個最優(yōu)點(diǎn)。將最優(yōu)點(diǎn)投影到曲面上，計(jì)算出該點(diǎn)在曲面上的投影點(diǎn)并作為最優(yōu)頂點(diǎn)。由該前沿線段和最優(yōu)頂點(diǎn)形成一個三角形劃分，更新前沿線段集合，去除已生成三角形的邊并加入新邊，如此反復(fù)，直到未劃分區(qū)域只剩下一個三角形，則網(wǎng)格劃分結(jié)束。

圖1 AFT法的基本原理Fig.1 The basic principle of AFT method

1.2 實(shí)現(xiàn)步驟

如圖2所示，基于PRO／E軟件的AFT法實(shí)現(xiàn)步驟如下：

1）應(yīng)用PRO／Toolkit選取對象的方法，選擇待劃分的曲面，并遍歷該曲面的邊界環(huán)；

2）逐一離散邊界環(huán)中的各條邊，形成前沿線段集合Ω；

3）應(yīng)用PRO／Toolkit繪制線段的方法，繪制出集合Ω中所有的前沿線段；

4）選擇最短的前沿線段，如圖3中AB，作為下一個生成三角形的前沿，計(jì)算前沿線段AB的中點(diǎn)M、推進(jìn)方向和單元尺寸函數(shù)h；

圖2 AFT法網(wǎng)格劃分算法流程Fig.2 The meshing procedure of AFT method

5）由前沿線段中點(diǎn)M和推進(jìn)方向，計(jì)算新生成點(diǎn)N，N點(diǎn)應(yīng)滿足：①到前沿線段的距離小于AB邊長的0.7倍，以保證三角形單元的質(zhì)量；②NA和NB與前沿線段集合Ω中任一前沿線段不相交；

如果以上兩個條件有一個不能滿足，則當(dāng)前節(jié)點(diǎn)N不可選，進(jìn)入6）繼續(xù)選擇；反之，選擇當(dāng)前節(jié)點(diǎn)為最優(yōu)頂點(diǎn)。

6）以N點(diǎn)為圓心，0.9倍的單元尺寸函數(shù)h（見圖3）為半徑R，建立搜索范圍，在已存在的前沿線段中尋找最優(yōu)頂點(diǎn)T，T點(diǎn)應(yīng)滿足：①三角形ABT不包含前沿線段集合Ω中任一前沿線段的端點(diǎn)；②TA和TB與前沿線段集合Ω中任一前沿線段不相交；③點(diǎn)T到前沿線段集合中其他前沿線段的距離大于AB邊長的一半，以避免狹長單元的產(chǎn)生。

對滿足條件的三角形計(jì)算其單元質(zhì)量評測函數(shù)，選擇使生成三角形質(zhì)量最高的點(diǎn)T，作為最優(yōu)頂點(diǎn)；如果以上3個條件有一個不能滿足，則當(dāng)前前沿線段不能生成新頂點(diǎn)，待后續(xù)過程中有合適的點(diǎn)再重新計(jì)算。

7）如果當(dāng)前前沿存在最優(yōu)頂點(diǎn)，應(yīng)用PRO／Toolkit繪制線段的方法在PRO／E環(huán)境中繪出新的三角形；更新前沿線段集合Ω，刪除舊的前沿，添加新的前沿；否則，執(zhí)行8）；

8）如果當(dāng)前集合Ω為空，網(wǎng)格劃分過程結(jié)束；否則，繼續(xù)執(zhí)行4），生成新單元。

圖3 最優(yōu)頂點(diǎn)的搜索范圍Fig.3 The search scope of the best point

2 關(guān)鍵問題的處理方法

2.1 曲面邊界離散方法

在Pro／E軟件中，曲面分為開曲面、半閉合曲面和閉合曲面3種。對于開曲面和半閉合曲面，筆者采用的離散方法如下：

1）確定邊界環(huán)中最短邊的離散密度，以防止最短邊不夠分的情況發(fā)生；

2）在參數(shù)域?qū)⒆疃踢厔澐譃橹付ǖ姆輸?shù)，然后將劃分的參數(shù)坐標(biāo)從參數(shù)域轉(zhuǎn)換成物理域三維坐標(biāo)，由此完成最短邊的劃分；

3）如果邊界環(huán)中存在下一條邊，轉(zhuǎn)步驟4）；否則，轉(zhuǎn)步驟5）；

4）根據(jù)最短邊離散的平均尺寸計(jì)算當(dāng)前邊的離散份數(shù)，然后轉(zhuǎn)步驟2）；

5）連接每條邊的離散結(jié)果，邊界環(huán)離散過程結(jié)束。

對于閉合曲面，根據(jù)旋轉(zhuǎn)參數(shù)方向?qū)⑶鎰澐譃閮蓚€開曲面，然后按照開曲面的離散方法離散。

2.2 新節(jié)點(diǎn)的計(jì)算方法

在新節(jié)點(diǎn)計(jì)算之前，要對單元尺寸h的確定方法進(jìn)行說明。選定前沿線段后，需要計(jì)算單元尺寸函數(shù)，以確定單元尺寸，筆者使用的單元尺寸函數(shù)是區(qū)域自主的方法。該方法的原理是，使新生成單元的大小與鄰近前沿的尺寸接近。在三維條件下，以當(dāng)前選定前沿線段的中點(diǎn)為中心構(gòu)造一個球形的區(qū)域，凡是在此區(qū)域內(nèi)的其他前沿線段都屬于新生成單元的尺寸確定參數(shù)。確定單元尺寸后，將按照下述步驟計(jì)算新節(jié)點(diǎn)的位置。

如圖4所示，首先定義3個向量n，e，d。其中n是前沿線段AB的中點(diǎn)M 的法向量，這個法向量取兩個端點(diǎn)A和B法向量的平均值，在實(shí)際的應(yīng)用過程中，這種假定能滿足網(wǎng)格劃分的要求。e是前沿線段AB的單位方向向量。d是前沿線段AB推進(jìn)方向的單位向量。根據(jù)推進(jìn)方向的計(jì)算方法，3個向量的關(guān)系滿足

在d方向上能與AB組成正三角形的點(diǎn)N′的位置，可計(jì)算如下

圖4 曲面上新節(jié)點(diǎn)的計(jì)算方法簡圖Fig.4 The diagram of computing new point on the surface

由于點(diǎn)N′不在曲面上，因此新節(jié)點(diǎn)N是N′在曲面上的投影。

2.3 線段與線段的重疊判定

在空間對曲面進(jìn)行網(wǎng)格劃分時，對于前沿線段重疊的判定可以映射到參數(shù)平面上，若空間內(nèi)兩條線段在參數(shù)平面上相交，則必在空間內(nèi)重疊。

圖5 參數(shù)域線段相交判斷示意圖Fig.5 The diagram of lines'intersection judgement on parametric space

設(shè)空間兩條線段A′B′和C′D′，采用幾何特征計(jì)算方法直接映射，依次求出三維空間點(diǎn)A′、B′、C′、D′的二維參數(shù)點(diǎn)A、B、C、D。由此得到兩條參數(shù)平面上的線段AB和線段CD。對AB和CD相交的判定原理如圖5所示。如果AB和CD相交，則必互相穿過，因此兩條線段相交的充要條件是A、B點(diǎn)位于線段CD的兩側(cè)且C、D點(diǎn)位于線段AB的兩側(cè)。若A、B點(diǎn)位于線段CD的兩側(cè)，等價(jià)于位于線段CD的兩側(cè)，即叉乘得到的平面法向量n1與叉乘得到的向量n2方向相反且夾角為π，即等價(jià)于n1與n2的點(diǎn)乘小于0，即

滿足式（1）即可保證A、B點(diǎn)位于線段CD的兩側(cè)，同理，為了保證C、D點(diǎn)位于線段AB的兩側(cè)，類似有

式（1），式（2）均成立時，線段AB和線段CD相交，參數(shù)域的相交同時也表明在物理域的兩條線段重疊，即線段A′B′和線段C′D′在空間內(nèi)重疊。

2.4 單元質(zhì)量評測

單元質(zhì)量評測用來比較不同的三角形單元的質(zhì)量，對于給定的三角形ABC，單元質(zhì)量評測函數(shù)為

采用這種方式計(jì)算的λ值，最大值為1，即三角形單元為正三角形時，其質(zhì)量最優(yōu)。當(dāng)三角形單元為鈍角或銳角三角形時，其質(zhì)量評測函數(shù)值會小于1。

2.5 運(yùn)動屬性和材料特性設(shè)定

為了仿真機(jī)械部件的工作過程，在提取了離散后得到的三角形平面片的幾何參數(shù)后，還應(yīng)添加運(yùn)動屬性和材料特性，以便建立邊界的物理模型。

目前開發(fā)的軟件，邊界運(yùn)動包括平動（含振動）、轉(zhuǎn)動和平轉(zhuǎn)動組合運(yùn)動。邊界運(yùn)動形式包括勻速運(yùn)動、變速運(yùn)動、伺服運(yùn)動和往復(fù)運(yùn)動等。伺服運(yùn)動是一種特殊的運(yùn)動形式，即當(dāng)外力達(dá)到給定閾值時，邊界才發(fā)生運(yùn)動，運(yùn)動可以是平動或轉(zhuǎn)動。

平動的屬性設(shè)置包含運(yùn)動方向向量、運(yùn)動速度和運(yùn)動的最大距離等。轉(zhuǎn)動屬性設(shè)置包含轉(zhuǎn)軸方向向量和角速度。

在一個機(jī)械部件中，不同零件可由不同材料制成。由于不同材料對應(yīng)的力學(xué)性質(zhì)不同，因此在提取不同零件上的曲面時，應(yīng)予以標(biāo)識。

上述運(yùn)動屬性和材料特性的添加，均是通過人機(jī)交互方式設(shè)置不同的特征值，以便在離散元法計(jì)算和計(jì)算結(jié)果仿真顯示時，采用不同的處理方式。

3 軟件設(shè)計(jì)和實(shí)例驗(yàn)證

在對PRO／E軟件進(jìn)行二次開發(fā)的基礎(chǔ)上，設(shè)計(jì)了基于AFT法的非規(guī)則曲面邊界建模軟件，并實(shí)現(xiàn)了與規(guī)則曲面邊界建模軟件的集成，從而開發(fā)出一種通用的離散元法邊界建模軟件，其結(jié)構(gòu)和流程如圖6所示。

圖6 基于PRO／E軟件的離散元法邊界建模軟件結(jié)構(gòu)和流程圖Fig.6 The structure and flow chart of the DEM boundary modeling software based on PRO／E

以一種開曲面——非完整直圓柱面、一種半閉合曲面——圓臺面和一種閉合曲面——球面為例，采用基于AFT的網(wǎng)格劃分方法，由曲面的CAD模型實(shí)現(xiàn)的曲面網(wǎng)格劃分，把曲面離散成三角形面片的組合。其流程如圖7所示。

圖7 由CAD模型實(shí)現(xiàn)的3種類型曲面網(wǎng)格劃分Fig.7 Three type surface's meshing created by CAD model

采用AFT的網(wǎng)格劃分方法，由一種非規(guī)則曲面——開溝器工作面的三維CAD模型，建立的該邊界的三維離散元法分析模型如圖8所示，圖9為開溝器工作過程的三維離散元法仿真結(jié)果。圖10～圖12分別為混合采用規(guī)則曲面和非規(guī)則曲面的建模方法，由一種排肥器、一種排種器和一種螺旋輸送機(jī)的三維CAD模型，建立了該邊界的三維離散元法分析模型，然后采用三維離散元法分析排肥器、排種器和輸送機(jī)的工作過程。

圖8 由非規(guī)則曲面——開溝器工作面的三維CAD模型建立的其三維離散元法分析模型Fig.8 The 3DDEM simulation analysis model of non-analytical surface，furrower's work surface created by 3DCAD model

圖9 非規(guī)則曲面——開溝器工作面開溝過程的三維離散元法仿真分析Fig.9 The 3DDEM simulation analysis of work process of non-analytical surface，furrower's work surface

圖10 由CAD模型實(shí)現(xiàn)的排肥器工作過程的三維離散元法仿真分析Fig.10 The 3DDEM simulation analysis of fertilizer apparatus's work process created by CAD model

圖11 由CAD模型實(shí)現(xiàn)的排種器工作過程的三維離散元法仿真分析Fig.11 The 3DDEM simulation analysis of metering device work process created by CAD model

圖12 用CAD模型實(shí)現(xiàn)的螺旋輸送機(jī)工作過程的三維離散元法仿真分析Fig.12 The 3DDEM simulation analysis of screw conveyer's work process created by CAD model

以上實(shí)例初步證明了筆者建立的離散元法邊界建模方法及其軟件的可行性和有效性，為采用離散元法分析復(fù)雜結(jié)構(gòu)機(jī)械部件的工作過程奠定了基礎(chǔ)。

4 結(jié) 語

筆者采用AFT法進(jìn)行非規(guī)則曲面網(wǎng)格劃分，把非規(guī)則曲面離散成三角形平面片的組合，同時添加運(yùn)動屬性和材料特性參數(shù)，由此建立了非規(guī)則曲面邊界的離散元法分析模型。在對PRO／E軟件進(jìn)行二次開發(fā)的基礎(chǔ)上，研制了非規(guī)則曲面邊界建模軟件，并實(shí)現(xiàn)了與規(guī)則曲面邊界建模軟件集成，從而開發(fā)出一種通用的離散元法邊界建模軟件。通過實(shí)例驗(yàn)證，初步證明了基于AFT邊界建模方法和軟件的可行性，為完善筆者提出的基于CAD模型的離散元法邊界建模方法和采用離散元法分析復(fù)雜結(jié)構(gòu)機(jī)械部件的工作過程奠定了基礎(chǔ)。

［1］KREMMER M，F(xiàn)AVIER J F.A Method for Representing Boundaries in Discrete Element Modelling-PartⅡ：Kinematics［J］.International Journal for Numerical Methods in Engineering，2001，51（12）：1423-1436.

［2］付宏，烏蘭，黃萬風(fēng)，等.基于圖元的三維離散元法邊界建模方法［J］.計(jì)算機(jī)集成制造系統(tǒng)，2008，14（12）：2328-2333.

FU Hong，WU Lan，HUANG Wan-feng，et al.Method for Modeling Boundaries Based on the Graphic Elements in the Three-Dimensional DEM［J］.Computer Integrated Manufacturing Systems，2008，14（12）：2328-2333.

［3］關(guān)振群，宋超，顧元憲，等.有限元網(wǎng)格生成方法研究的新進(jìn)展［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2003，15（1）：1-14.

GUAN Zhen-qun，SONG Chao，GU Yuan-xian，et al.Recent Advances of Research on Finite Element Mesh Generation Methods［J］.Journal of Computer-Aided Design and Computer Graphics，2003，15（1）：1-14.

［4］杜群貴，劉勝，黃曉東.閉曲面有限元網(wǎng)格生成的邊界預(yù)調(diào)整方法［J］.華南理工大學(xué)學(xué)報(bào)：自然科學(xué)版，2007，35（2）：27-32.

DU Qun-gui，LIU Sheng，HUANG Xiao-dong.Boundary Preadjustment Method for Finite Element Mesh Generation of Closed Surface［J］.Journal of South China University of Technology：Natural Science Edition，2007，35（2）：27-32.

［5］LAU T S，LO S H.Finite Element Mesh Generation over Analytical Curved Surfaces［J］.Computer ＆Structures，1996，59（2）：301-309.

［6］LOHNER R.Extension and Improvements of the Advancing-Front Grid Generation Technique［J］.Communications in Numerical Methods in Engineering，1996，12（10）：683-702.

［7］ITO Y，SHIH A M，ERUKALA A K，et al.Parallel Unstructured Mesh Generation by an Advancing Front Method［J］.Mathematics and Computers in Simulation，2007，75（5／6）：200-209.

［8］GUAN Zhen-qun，SHAN Ju-lin，ZHENG Yao，et al.An Extended Advancing Front Technique for Closed Surfaces Mesh Generation［J］.International Journal for Numerical Methods in Engineering，2008，74（4）：642-667.

［9］黃山.基于CAD模型的三維離散元法邊界建模方法研究［D］.長春：吉林大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，2011.

HUANG Shan.The Research of Boundary Modeling Method of 3DDEM Based on CAD Model［D］.Changchun：College of Computer Science and Technology，Jilin University，2011.

［10］于建群，付宏，劉振宇，等.基于CAD模型的離散元法邊界建模方法：中國，200510016835.7［P］.2006-07-28.

YU Jian-qun，F(xiàn)U Hong，LIU Zhen-yu，et al.Method Based on CAD Model for DEM Boundary Modeling：China，200510016835.7［P］.2006-07-28.

Development of Software for Establishing Analytical Model of Irregular Surfaces in Discrete Element Method

FU Honga，Lü Youa，XU Jinga，HUANG Shana，YU Jian-qunb
（a.College of Computer Science and Technology，Jilin University，Changchun 130012，China；b.School of Biological and Agricultural Engineering，Jilin University，Changchun 130025，China）

It needs to establish analysis models of machine parts（boundaries），when use DEM（Discrete Element Method）to analyze the contact action between machine parts and granular materials.There exist irregular surfaces which can not be expressed by the elementary analytic function in the parts'surfaces which contact with granular materials.The AFT（Advancing Front Technique）was used to mesh and discrete irregular surfaces into the triangle planar units，parameters of movement characters and material properties were added in the same time，so the DEM analysis models of irregular surfaces was created.Based on the redevelopment of PRO／E software，the boundary modeling software of irregular surfaces was developed.By application examples，the feasibility of boundary modeling method and the software which based on the AFT was validated，which lays foundations for simulation and analysis of working process for machine parts with complex structure.

discrete element method；boundary modeling；irregular surfaces；advancing front technique

TP391.72

A

1671-5896（2012）01-0023-07

2011-10-10

國家自然科學(xué)基金資助項(xiàng)目（60973090，51175219，11172112）；吉林省科技發(fā)展計(jì)劃基金資助項(xiàng)目（20100313）

付宏（1960—），女，浙江寧波人，吉林大學(xué)教授，碩士生導(dǎo)師，主要從事數(shù)值模擬與仿真研究和軟件工程，（Tel）86-13086875669（E-mail）fuhong＠jlu.edu.cn；通訊作者：于建群（1958—），男，長春人，吉林大學(xué)教授，博士生導(dǎo)師，主要從事數(shù)字化設(shè)計(jì)研究，（Tel）86-13074320858（E-mail）yujiangqun＠jlu.edu.cn。

（責(zé)任編輯：劉俏亮）