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

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

1 引言

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

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

面向復雜模型， 接觸關系識別算法達到實用化標準需要解決的兩個關鍵問題是： 提高識別精度和提高識別效率。模型中存在幾何體自身的誤差和幾何體之間的誤差等多種誤差來源， 選取合適的計算容差是提高算法識別精度的關鍵。已有的方法大都采用全局容差進行接觸識別[10]， 但全局容差的識別精度不高。已有的方法大都通過空間篩選技術提高算法效率[15]。

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

2 接觸關系識別算法

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

2.1 算法概述

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

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

1)采用空間劃分和包圍盒篩選相結合的技術篩選幾何面對；

2)采用多線程并行技術加速篩選和接觸檢測過程。

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

2.2 局部容差計算

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

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

(1)

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

(2)

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

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

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

(3)

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

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

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

(4)

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

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

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

(5)

2.3 面-面接觸檢測

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

①s1和s2的法向相反；

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

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

2.4 空間劃分+包圍盒分層篩選技術

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

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

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

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

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

圖1 接觸算法并行設計流程

2.5 多線程并行

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

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

3 算法分析

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

圖2 公開發(fā)動機模型

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

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

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

表2 不同方法識別的接觸關系結果對比

4 結構分析應用實例

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

圖3中的離心機模型， 包含126個幾何實體， 821個幾何面。采用本文的接觸關系識別算法， 共檢測到393對面-面接觸關系， 如圖5所示， 圖中紅色表示面-面接觸。

圖3 某土工離心機結構

考慮離心機結構在旋轉離心力載荷作用下的剛度與變形分析，待求解的力學平衡方程可表示為如下形式

(6)

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

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

(7)

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

通過采用罰函數(shù)法構造約束變分修正泛函，可將該含附加約束問題轉化為無約束修正泛函的變分問題，進而得到待求解問題的線性方程組

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

(8)

其中，K和P分別對應式(6)原問題的結構剛度矩陣和外力向量，K2和Q分別對應由式(7)的附加約束方程所引入的額外結構約束剛度矩陣和約束外力向量，形式如下

(9)

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

圖5 算法識別出的面-面接觸關系

圖6 離心機結構的位移變形圖

5 結論

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

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