復(fù)雜幾何模型接觸關(guān)系并行識(shí)別方法及應(yīng)用

劉田田，劉偉杰，楊 洋，鄭 澎

(1. 中物院高性能數(shù)值模擬軟件中心，北京 100088；2. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100094；3. 中國(guó)工程物理研究院計(jì)算機(jī)應(yīng)用研究所，四川 綿陽(yáng) 621900)

1 引言

高性能計(jì)算機(jī)能力的不斷提升使得結(jié)構(gòu)分析計(jì)算能針對(duì)數(shù)千部件的復(fù)雜系統(tǒng)進(jìn)行模擬[1]。復(fù)雜系統(tǒng)的多部件之間通過(guò)粘接、鉸接、鉚接等方式裝配在一起， 存在著接觸、摩擦等相互作用， 獲得零部件之間的相互作用是模擬計(jì)算過(guò)程中的重要一環(huán)[2]， 也是決定結(jié)構(gòu)分析結(jié)果準(zhǔn)確可靠的關(guān)鍵之一[3，4]。

在進(jìn)行接觸分析計(jì)算時(shí)， 需要先確定哪些部件間會(huì)存在接觸關(guān)系。在工業(yè)生產(chǎn)部門中， 數(shù)值模擬分析人員得到的來(lái)自于設(shè)計(jì)部門的幾何模型往往缺少裝配信息， 對(duì)于多部件、裝配關(guān)系復(fù)雜的系統(tǒng)級(jí)模型， 不得不花費(fèi)大量的時(shí)間人工尋找接觸關(guān)系。因此， 發(fā)展自動(dòng)識(shí)別接觸關(guān)系的方法顯得尤為重要。近年來(lái)， 基于幾何信息的接觸關(guān)系識(shí)別方法得到各界學(xué)者的廣泛研究[5-14]。但是， 目前已有的方法大都只適用于簡(jiǎn)單模型， 局限于效率或精度的原因， 對(duì)復(fù)雜模型和高階曲面不適用。

面向復(fù)雜模型， 接觸關(guān)系識(shí)別算法達(dá)到實(shí)用化標(biāo)準(zhǔn)需要解決的兩個(gè)關(guān)鍵問(wèn)題是： 提高識(shí)別精度和提高識(shí)別效率。模型中存在幾何體自身的誤差和幾何體之間的誤差等多種誤差來(lái)源， 選取合適的計(jì)算容差是提高算法識(shí)別精度的關(guān)鍵。已有的方法大都采用全局容差進(jìn)行接觸識(shí)別[10]， 但全局容差的識(shí)別精度不高。已有的方法大都通過(guò)空間篩選技術(shù)提高算法效率[15]。

本文面向復(fù)雜B-Rep幾何模型， 提出了一種適用于任意曲面類型的接觸關(guān)系識(shí)別算法， 并將該算法應(yīng)用到真實(shí)工程模型的結(jié)構(gòu)分析實(shí)例中。算法基于幾何模型的離散數(shù)據(jù)進(jìn)行接觸識(shí)別， 根據(jù)模型的幾何特征自動(dòng)計(jì)算適用于相鄰幾何體和幾何面之間的局部容差， 并采用基于局部容差的碰撞檢測(cè)技術(shù)提高識(shí)別精度。采用多線程并行技術(shù)， 組合使用空間劃分技術(shù)和包圍盒篩選技術(shù)提高識(shí)別效率。

2 接觸關(guān)系識(shí)別算法

接觸表示兩個(gè)幾何體之間有某些面、線或點(diǎn)重合， 但是沒(méi)有體積上的重合[16]。幾何體間的接觸關(guān)系有六類： 面-面接觸、面-線接觸、面-點(diǎn)接觸、線-線接觸、線-點(diǎn)接觸和點(diǎn)-點(diǎn)接觸。不同類型的接觸識(shí)別方法都是類似的， 由于計(jì)算時(shí)大部分接觸約束是通過(guò)面聯(lián)接的， 因此， 本文只介紹面-面接觸的識(shí)別方法。

2.1 算法概述

B＿Rep幾何模型通常由大量的三維幾何體組成， 包含多種類型的曲面。為了適用于包含任意曲面類型的模型， 算法基于幾何面的離散數(shù)據(jù) (三角面片) 進(jìn)行接觸探測(cè)。目前， 商業(yè)和開源幾何引擎都提供B＿Rep模型的離散化接口， 獲取離散數(shù)據(jù)非常簡(jiǎn)單。

對(duì)模型中的任意兩個(gè)幾何面， 判斷是否存在接觸關(guān)系的方法是對(duì)兩個(gè)幾何面的離散數(shù)據(jù)進(jìn)行碰撞檢測(cè)。復(fù)雜工程模型的尺寸特征跨度大，為了提升算法精度，采用基于局部容差的碰撞檢測(cè)。由于復(fù)雜裝配模型中往往存在數(shù)千上萬(wàn)個(gè)幾何面， 離散數(shù)據(jù)可能包含數(shù)百萬(wàn)個(gè)三角面片， 直接進(jìn)行三角形接觸檢測(cè)計(jì)算量太大。本文采用兩種策略提高接觸關(guān)系的識(shí)別效率：

1)采用空間劃分和包圍盒篩選相結(jié)合的技術(shù)篩選幾何面對(duì)；

2)采用多線程并行技術(shù)加速篩選和接觸檢測(cè)過(guò)程。

為方便算法描述， 記{Bi,i=1,…，n}為模型中所有三維幾何體的集合，{Fij，j=1,…,m}為第i個(gè)幾何體包含的所有幾何面的集合。Box(Bi)表示幾何體Bi的包圍盒，Box(Fij)表示幾何體Bi中第j個(gè)面Fij的包圍盒。Box(Bi)±ε表示Bi的包圍盒向各方向外擴(kuò)ε。

2.2 局部容差計(jì)算

局部容差是根據(jù)模型的局部特征長(zhǎng)度計(jì)算得到的。 幾何體的特征長(zhǎng)度定義為包圍盒的最短邊長(zhǎng)

L(B)=min{lx，ly，lz}，

(1)

其中l(wèi)x，ly，lz分別表示包圍盒在x，y，z方向的邊長(zhǎng)。幾何面的特征長(zhǎng)度定義為組成該面的最短邊長(zhǎng)度或包圍盒的對(duì)角線長(zhǎng)度

(2)

其中k為該面包含的邊個(gè)數(shù)。

定義幾何體的特征容差εb如下

εb=max{L(Bi)/R，i=1，…，n}，

(3)

其中n表示模型包含的幾何體個(gè)數(shù)，R是容差計(jì)算因子， 取值大于0。R越大， 局部容差越小， 根據(jù)經(jīng)驗(yàn)，R默認(rèn)取值20。

定義幾何面的特征容差εf如下

εf=max{L(Fj)/R，j=1，…，m}

(4)

其中m表示模型包含的幾何面?zhèn)€數(shù)。

定義兩個(gè)相鄰幾何面之間的局部容差如下

tol(Fj1，F(xiàn)j2)=min{L(Fj1)/R，L(Fj2)/R}

(5)

2.3 面-面接觸檢測(cè)

記Fj1和Fj2是兩個(gè)待檢測(cè)的幾何面，Sj1和Sj2分別為兩幾何面的離散數(shù)據(jù)， 兩面之間的局部容差根據(jù)式(5)計(jì)算得到。判斷兩個(gè)幾何面之間是否存在接觸關(guān)系， 只需判斷兩個(gè)集合面的離散數(shù)據(jù)中是否存在三角形共面相交，即存在s1∈Sj1，s2∈Sj2滿足

①s1和s2的法向相反；

②s1和s2在容差范圍內(nèi)共面相交 (或包含)。

三角形之間的相交檢測(cè)是圖形學(xué)中的經(jīng)典問(wèn)題， 有多種求解方法， 如分離軸算法[17]、Moller算法[18]、Devillers算法[19]、Shen算法[20]、Tropp算法[21]等。采用的是分離軸算法。為了減少不必要的碰撞操作， 對(duì)每個(gè)幾何面的三角面片構(gòu)建k-d樹， 基于k-d樹進(jìn)行鄰域內(nèi)的三角形搜索。

2.4 空間劃分+包圍盒分層篩選技術(shù)

由于幾何模型的零件之間是緊密排列的， 且進(jìn)行的碰撞檢測(cè)是靜態(tài)的， 為了提高效率， 算法選用空間劃分和包圍盒分層篩選相結(jié)合的方法來(lái)減少面-面之間的接觸檢測(cè)。根據(jù)B-Rep模型的層級(jí)關(guān)系， 如果兩個(gè)幾何體的包圍盒不相交 (含相切)， 那么這兩個(gè)幾何體之間一定不存在接觸面對(duì)。如果兩個(gè)幾何面的包圍盒不相交， 那么這兩個(gè)幾何面一定不接觸。幾何面對(duì)篩選的算法如下：

Step 1：計(jì)算所有幾何體、幾何面的包圍盒。

Step 2：對(duì)所有幾何體的包圍盒構(gòu)建k-d樹， 記為BTree。對(duì)第每個(gè)幾何體包含的幾何面的包圍盒構(gòu)建k-d樹， 記為FTree(i)， i取值1～n。

Step 3：對(duì)任一幾何體Bi， 在BTree中查找與Box(Bi)±εb相交的包圍盒列表{Box(Bik),k=1,…,K}。遍歷Bik， 為避免重復(fù)計(jì)算， 若Bik的實(shí)體編號(hào)小于等于Bi， 則跳過(guò)； 否則， 執(zhí)行Step 4。

Step 4：對(duì)OBi的任一幾何面Fij， 在FTree(ik)中查找與Box(Fij)±εf相交的包圍盒列表{Box(Fikl),l=1,…,L}。若Fikl的實(shí)體編號(hào)小于等于Fij， 則跳過(guò)； 否則，F(xiàn)ij與Fikl時(shí)可能存在接觸關(guān)系的面對(duì)，需進(jìn)行面-面接觸檢測(cè)。

圖1 接觸算法并行設(shè)計(jì)流程

2.5 多線程并行

串行算法難以滿足復(fù)雜模型在計(jì)算效率方面的應(yīng)用需求。目前的處理器都是多核系統(tǒng)， 采用多線程分配任務(wù)的策略是提高算法性能的有效手段。本文采用OpenMP實(shí)現(xiàn)多線程并行算法。

算法中可以被并行分解的任務(wù)內(nèi)容包括： 包圍盒計(jì)算、幾何特征長(zhǎng)度計(jì)算、多個(gè)k-d樹的構(gòu)建、包圍盒篩選和面對(duì)接觸檢測(cè)。因此， 并行設(shè)計(jì)流程如圖1所示， 其中I、II、III、IV表示循環(huán)。由于幾何體包圍盒篩選、幾何面包圍盒篩選和接觸檢測(cè)之間存在層級(jí)依賴， 為了充分利用并行資源， 將并行粒度放在幾何體包圍盒篩選的循環(huán)上。

3 算法分析

利用圖2中的裝配模型對(duì)算法效率和算法精度進(jìn)行了測(cè)試， 該模型包含76個(gè)體， 14467個(gè)面，含多種類型的曲面， 模型內(nèi)包含505對(duì)面-面接觸關(guān)系。如果對(duì)所有的幾何面進(jìn)行兩兩接觸檢測(cè)， 則需要檢測(cè)209，294，089次。通過(guò)包圍盒篩選后， 只需要檢測(cè)557，367次， 檢測(cè)次數(shù)降到了0.27%。幾何面對(duì)篩選和面-面接觸檢測(cè)過(guò)程中， 算法采用k-d樹進(jìn)行鄰近搜索， k-d樹的構(gòu)建復(fù)雜度是O(NlogN)， 搜索的復(fù)雜度是O(logN)， 大大提高了計(jì)算效率。

圖2 公開發(fā)動(dòng)機(jī)模型

利用該模型對(duì)算法的并行效率進(jìn)行分析， 分別采用1/2/4/8線程進(jìn)行測(cè)試， 結(jié)果如表1所示。隨著線程數(shù)增加， 計(jì)算時(shí)間逐漸下降。8核并行的運(yùn)行時(shí)間是串行時(shí)間的1/3左右， 表明多線程并行的有效性。

表1 不同線程數(shù)目算法運(yùn)行時(shí)間

由于識(shí)別方法主要基于幾何信息， 并未考慮物理模型， 因此對(duì)于識(shí)別結(jié)果， 還需要分析人員進(jìn)行手工剔除，把需要被剔除的接觸關(guān)系稱為“冗余接觸關(guān)系”。分別測(cè)試了本文的接觸識(shí)別方法和商業(yè)軟件的自動(dòng)識(shí)別方法， 結(jié)果如表2所示。測(cè)試環(huán)境CPU型號(hào)為Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.4GHz， 8核。從表中看出， 采用本文的方法識(shí)別結(jié)果更精確， 只識(shí)別出少量的冗余接觸， 效果好于商業(yè)軟件ABAQUS和ANSYS。分析原因， 主要是因?yàn)楸疚牡淖R(shí)別方法是基于局部容差進(jìn)行的， 而商業(yè)軟件采用一致化的全局容差， 為識(shí)別最大特征， 不得不設(shè)置最大容差值， 從而識(shí)別出大量的冗余接觸。對(duì)于模型尺寸跨度比較大的復(fù)雜模型， 選取合適的全局容差非常困難， 使用基于幾何特征的局部容差計(jì)算效果明顯優(yōu)于全局容差。 從表中的運(yùn)行時(shí)間來(lái)看， 本文中的算法效率略低于ANSYS，但明顯優(yōu)于ABAQUS。

表2 不同方法識(shí)別的接觸關(guān)系結(jié)果對(duì)比

4 結(jié)構(gòu)分析應(yīng)用實(shí)例

本節(jié)以某土工離心機(jī)模型的靜強(qiáng)度分析為例驗(yàn)證接觸關(guān)系算法的有效性。離心機(jī)是一種利用旋轉(zhuǎn)來(lái)模擬特定高重力場(chǎng)加速度環(huán)境的試驗(yàn)設(shè)備， 廣泛地應(yīng)用于科研和國(guó)民經(jīng)濟(jì)建設(shè)的多個(gè)領(lǐng)域。土工離心機(jī)為巖土離心模型試驗(yàn)提供了重要的研究手段， 廣泛用于堤壩、邊坡地基、擋土結(jié)構(gòu)、路基、隧道、重要建筑物基礎(chǔ)等方面的工程應(yīng)用研究， 是評(píng)估巖土工程設(shè)計(jì)結(jié)果的合理性、安全性和經(jīng)濟(jì)性等重要設(shè)備。

圖3中的離心機(jī)模型， 包含126個(gè)幾何實(shí)體， 821個(gè)幾何面。采用本文的接觸關(guān)系識(shí)別算法， 共檢測(cè)到393對(duì)面-面接觸關(guān)系， 如圖5所示， 圖中紅色表示面-面接觸。

圖3 某土工離心機(jī)結(jié)構(gòu)

考慮離心機(jī)結(jié)構(gòu)在旋轉(zhuǎn)離心力載荷作用下的剛度與變形分析，待求解的力學(xué)平衡方程可表示為如下形式

(6)

其中，Ω表示域內(nèi)，Γ表示力邊界，σ是二階應(yīng)力張量，f是體載荷，n是邊界外法向，T是邊界上的外載荷。

此外，如圖4示意，由于該離心機(jī)模型的幾何實(shí)體間劃分為非協(xié)調(diào)網(wǎng)格，為求解式(6)的方程，還需在每對(duì)交界面處引入約束方程

(7)

圖4 非協(xié)調(diào)網(wǎng)格綁定約束條件

通過(guò)采用罰函數(shù)法構(gòu)造約束變分修正泛函，可將該含附加約束問(wèn)題轉(zhuǎn)化為無(wú)約束修正泛函的變分問(wèn)題，進(jìn)而得到待求解問(wèn)題的線性方程組

(K+K2)·u=P-Q，

(8)

其中，K和P分別對(duì)應(yīng)式(6)原問(wèn)題的結(jié)構(gòu)剛度矩陣和外力向量，K2和Q分別對(duì)應(yīng)由式(7)的附加約束方程所引入的額外結(jié)構(gòu)約束剛度矩陣和約束外力向量，形式如下

(9)

通過(guò)自研前處理引擎SuperMesh[22]生成網(wǎng)格， 自研結(jié)構(gòu)強(qiáng)度數(shù)值仿真軟件JSolid開展了旋轉(zhuǎn)離心力載荷作用下的靜強(qiáng)度計(jì)算， 其變形位移云圖如圖6所示?？梢钥吹剑?離心機(jī)整體結(jié)構(gòu)的位移分布光滑、連續(xù)、無(wú)間斷， 驗(yàn)證了該結(jié)構(gòu)接觸關(guān)系識(shí)別的完整性。

圖5 算法識(shí)別出的面-面接觸關(guān)系

圖6 離心機(jī)結(jié)構(gòu)的位移變形圖

5 結(jié)論

本文面向復(fù)雜幾何模型的結(jié)構(gòu)分析應(yīng)用， 提出了一種基于模型離散數(shù)據(jù)的接觸關(guān)系識(shí)別算法， 該算法適用于包含任意曲面類型的復(fù)雜模型。通過(guò)公開的發(fā)動(dòng)機(jī)模型驗(yàn)證了算法的高效和高精度。最后將接觸識(shí)別的結(jié)果應(yīng)用到實(shí)際工程模型的結(jié)構(gòu)分析模擬中， 進(jìn)一步驗(yàn)證了算法的有效性。

目前該算法已集成在并行前處理引擎SuperMesh中。為進(jìn)一步提升識(shí)別效率， 下一步考慮發(fā)展CPU+GPU異構(gòu)并行識(shí)別算法。