固定網格中基于邊界描述的CAD模型快速分析方法

李心耀，張衛紅，陳亮

西北工業大學 工程仿真與宇航計算技術實驗室，西安 710072

固定網格與傳統貼體網格的最大區別在于其采用獨立于結構模型的規則網格，將傳統有限元分析中的網格劃分問題轉換為單元識別問題［1－2］。實現單元識別的前提就是能夠快速地判斷固定網格中任意一個節點是否位于CAD模型區域之內，然后根據一個單元所包含所有節點的位置確定單元的屬性。如果一個單元或者被細化而成的子單元所包含的節點分別位于CAD模型之內和之外，則該單元或該子單元需要采用四叉樹或八叉樹技術進行細化。

目前，用來實現固定網格中節點屬性判斷的方法主要有2種：一種是采用隱式函數描述CAD模型，通過函數在固定網格中每一個節點處的正負值來判斷該點的位置［3－6］；另一種是采用像素模型描述CAD模型，通過密集分布的像素點來識別固定網格中任意一個節點的位置［7－9］。然而，這2種方法都不能直接利用通過設計軟件建立的結構CAD模型，均需要對CAD模型進行重構，且對于復雜結構重構困難、費時。

本文研究的基于邊界描述的CAD模型在固定網格中的快速單元識別方法不同于現有的單元識別方法，它直接利用設計軟件建立的CAD模型信息，避免了對CAD模型進行重構的操作，降低了模型轉換過程中可能造成的精度損失，同時也降低了固定網格中單元識別的難度。結合加權B樣條有限胞元法，本文建立了一種固定網格中由結構CAD模型轉換成有限元分析模型，進而完成有限元仿真計算的新的設計框架。

1 CAD模型至分析模型的快速轉換

1.1 CAD模型的邊界描述法

在計算機輔助設計軟件中，CAD模型的常用描述方法有邊界表示法，構造實體幾何法，單元分解表示法和特征表示法等［10］。其中采用邊界表示法描述CAD模型時，模型的相關數據都可以顯式地給出，有利于后續的數控加工和3D打印等操作。STL格式是CAD模型邊界表示的一種重要方法，它采用若干個三角面片來近似CAD模型的邊界，同時這些三角面片的法向量和頂點坐標可以直接獲得。以一個帶孔立方體為例，在給定近似精度的前提下，它可以利用包含96個三角面片的STL格式進行表示，如圖1所示。

圖1 帶孔立方體的CAD模型及其STL格式Fig.1 CAD model of cube with hole and its STL format

1.2 基于射線交點法的模型轉換

1.2.1 射線交點法原理

空間任意一點P（xP，yP，zP）與前述帶孔立方體（STL格式）的位置關系如圖2所示。在笛卡爾坐標系中，以點P為起點，沿x軸正向作一條射線LP，統計射線LP與帶孔立方體的交點的個數，并記為nc。

圖2 任意一點P與帶孔立方體的位置關系Fig.2 Relative position between point Pand cube with hole

若nc為偶數，則點P位于帶孔立方體的外部；若nc為奇數，則點P位于帶孔立方體的內部。另外，點P位于帶孔立方體邊界上的情況可以通過坐標比較等方式確定。

1.2.2 交點統計原則

點與多邊形的位置關系判斷是計算機圖形學的基本問題之一，國內外相關學者對此進行了深入的研究［11－13］，相關技術發展較為成熟。比較而言，點與多面體的位置關系判斷是點與多邊形位置關系判斷的延伸和拓展。然而，增加的一維空間使得點與多面體的位置關系判斷變得更為困難。目前，點與多面體位置關系判斷的方法主要有以下幾種：Kalay［14］、Lane［15］和 Horn［16］等分別提出通過分解或降維的算法實現點與多面體位置關系的判斷，但是這些方法應用較為復雜，并且容易出現奇異現象。像素空間法［17］通過密集分布的像素點來識別待判斷點的位置。八叉樹法［18］借助對空間的逐層劃分進行點的位置判斷。二元空間劃分方法［19］通過對空間信息進行分類來判斷點的位置。基于Jordan曲線定理的方法［20］是應用最普遍的一種方法，它通過統計從待判斷點沿任意方向所引射線與多面體交點個數的策略實現點與多面體相對位置的識別。Feito-Torres方法［21－22］借助待判斷點，將多面體拆分為多個四面體，然后基于這些四面體的符號體積實現點的相對位置信息判斷。文獻［23］對上述方法進行了詳細的比較，并對每一種方法的優缺點和適用范圍進行了闡述。

本文基于結構CAD模型的STL格式中三角面片的頂點和法向量信息，在Jordan曲線定理方法的基礎上，提出了向量運算和坐標參數相融合的算法，解決了Jordan曲線定理方法中射線過頂點、邊或面（共面）時，需要重新確定射線方向的問題，實現了初始射線與三角面片任意相交情況的判定。

1）一般情形

假定△ABC是任意一個結構CAD模型的STL格式中的任意一個三角形，A（xa，ya，za）、B（xb，yb，zb）、C（xc，yc，zc） 和 n△ABC分 別 為△ABC的頂點和法向量。在笛卡爾坐標系中，一般情形下射線LP與△ABC的相交情況如圖3所示。

圖3 一般情形下射線LP與△ABC的相交情況Fig.3 Intersection between LPand△ABC in common situation

如圖3所示，P0為射線LP與 △ABC的交點。如果P0位于 △ABC內部，則向量可以由向量表示，其數學表達式為

式中：0≤m≤1，0≤n≤1，0≤m＋n≤1。

另外，射線LP的參數化方程可以表示為

式中：k≥0，nx＝ （1，0，0）。

根據式（1）和式（2）可得

求解式（3），可得m，n和k的值如表1所示。表中，類別Ⅰ中P0位于△ABC的內部，因此記射線LP與△ABC之間有1個交點。類別Ⅳ中，可以直接確定P0位于結構CAD模型的邊界上。類別Ⅱ和Ⅲ中，P0位于△ABC的一條邊或一個頂點上，需要進行深入的討論。

表1 m、n和k的取值與P0的位置Table 1 Values of m，n，kand positions of P0

2）交點位于三角形的一條邊上

由CAD模型的STL格式可知，三角形的任意一條邊都同時屬于相鄰的2個三角形，也就是說當P0位于△ABC的一條邊上時，它同時也位于與△ABC相鄰三角形的邊上。如圖4所示，2個相鄰的三角形T1和T2，其法向量分別為nT1和nT2。nx＝ （1，0，0）為沿x 軸正向的單位法向量。

計算得到向量nT1和nx之間的夾角θ1，向量nT2和nx之間的夾角θ2。根據θ1和θ2之間的關系，可以得到此種情況下射線與三角形之間的交點統計方法如下：①θ1和θ2均為銳角或鈍角時，記射線與CAD模型在此處有1個交點；②θ1和θ2分別為銳角和鈍角時，無交點。

圖4 交點位于三角形的一條邊上Fig.4 Intersection lies on one edge of triangle

3）交點位于三角形的一個頂點上

由CAD模型的STL格式可知，三角形的任意一個頂點至少屬于相鄰的3個三角形，也就是說當P0位于△ABC的一個頂點上時，它同時也位于與△ABC共頂點的相鄰三角形的頂點上。如圖5所示，共頂點的4個三角形T1、T2、T3和T4，其法向量分別為nT1、nT2、nT3和nT4。

計算得到向量nT1和nx之間的夾角θ1，向量nT2和nx之間的夾角θ2，向量nT3和nx之間的夾角θ3，向量nT4和nx之間的夾角θ4。根據θ1、θ2、θ3和θ4之間的關系，可以得到此種情況下射線與三角形之間的交點統計方法如下：①θ1，θ2，θ3和θ4均為銳角或鈍角時，記射線與CAD模型在此處有1個交點；② 其他情況，無交點。

4）射線與三角形位于同一個平面上

連接點P與△ABC的任意一個頂點（例如：頂點A），得到向量，如果向量與△ABC的法向量n△ABC之間是正交關系，則可以確定點P與△ABC位于同一個平面上。此時，可以分以下幾種情況進行討論：

當點P在△ABC內部時，式（1）可以寫成如下矩陣形式：

式中：0≤m≤1，0≤n≤1，0≤m＋n≤1。

為了方便對式（4）的解的討論，定義

圖5 交點位于三角形的一個頂點上Fig.5 Intersection lies at one vertex of triangle

若矩陣G的秩為2，并且矩陣N的秩也為2，則式（4）有唯一解。同時，若此唯一解滿足0≤m≤1，0≤n≤1和0≤m＋n≤1的條件，則可以確定點P位于結構CAD模型的邊界上。

當點P位于△ABC外部時，若點P在△ABC的右側，此時射線LP與 △ABC無交點。因此，只需考慮點P位于△ABC左側的情況，如圖6所示。

實際上在結構CAD模型的STL格式中，射線LP可能會同時與多個三角形共面，如圖7所示。

圖6 點P位于△ABC的左側Fig.6 Point Plies on the left side of△ABC

圖7 射線與三角形共面情況的兩種分類（除了共面曲面外，其余曲面上的離散三角形未畫出）Fig.7 Two different categories of coplanar situation of ray and triangle （Besides on coplanar surface，discrete triangles on other surfaces are not plotted）

圖7 （a）中，S2是與點P共面的曲面，S1和S3分別是與S2左相鄰和右相鄰的曲面，nS1和nS3分別為面S1和S3的法向量。JL1是射線LP與面S2在最左側的交點，JR1是射線LP與面S2在最右側的交點。圖7（b）中，S5是與點P共面的曲面，S4和S6分別是與S5左相鄰和右相鄰的曲面，nS4和nS6分別為面S4和S6的法向量。JL2是射線LP與面S5在最左側的交點，JR2是射線LP與面S5在最右側的交點。

通過計算可以分別得到nS1和nx之間的夾角θ1，nS3和nx之間的夾角θ2，nS4和nx之間的夾角θ3，nS6和nx之間的夾角θ4。根據θ1和θ2，θ3和θ4之間的關系，可以得到此種情況下射線與三角形之間的交點統計方法為：①θ1和θ2同時為銳角或鈍角，記射線與CAD模型在此處有1個交點；②θ3和θ4分別為銳角和鈍角時，無交點。

特別地，CAD模型的STL格式中可能存在幾個面與射線同時共面的情況，在交點統計時，只關注最左和最右端的兩個面的法向量與射線的夾角之間的關系即可，上述射線與三角形共面時的交點統計方法仍適用。

1.3 單元識別與局部細化

以帶孔立方體結構為例，由CAD模型向固定網格中有限元分析模型的轉換過程如圖8所示。

圖8 CAD模型向固定網格中分析模型的轉換流程Fig.8 Transform flow from CAD model to analysis model at fixed grid

在笛卡爾坐標系中，首先建立一個尺寸略大于結構CAD模型的六面體，并采用各向等比的方法將其劃分成規則的網格結構；然后將STL格式的CAD模型嵌入網格結構中；接著根據1.2節闡述的射線交點法判斷固定網格中每一個節點與CAD模型的位置關系，進而確定單元的屬性，完成CAD模型向固定網格中有限元分析模型的轉換。相對于當前固定網格中主要采用的2種模型轉換方法而言，本文所提方法避免了采用隱式函數描述CAD模型方法中將CAD模型人為地分解成若干基元，再通過布爾運算進行組合的過程；避免了采用像素模型描述CAD模型方法中需借助專門工具生成像素模型的過程，不需要任何的預先處理，直接根據CAD模型的信息得到其在固定網格中對應的分析模型。

為提高固定網格中有限元分析模型對CAD模型的近似精度，采用八叉樹技術對邊界單元進行細化處理，并對細化后的子單元進行新的單元識別，如圖9所示。

一般而言，固定網格越密，對邊界單元的細化層級越多，有限元分析模型越接近CAD模型，但所耗費的時間也越長。以帶孔立方體為例，固定網格中單元數量與CPU耗時之間的關系如圖10所示。由圖10可知，網格越密，轉換模型所耗費的時間也越多。因此需要在模型近似精度和計算效率上進行合理的平衡。

圖9 基于八叉樹技術的邊界單元細化Fig.9 Refinement of boundary element based on octree technique

圖10 固定網格數量與CPU耗時之間的關系Fig.10 Relationship between number of fixed grids and CPU times

2 加權B樣條有限胞元方法

2.1 B樣條有限胞元方法

B樣條有 限胞元方法 （B-Spline Finite Cell Method，BSFCM）［3－6，24-25］作為計算力學領域中新興的先進數值方法之一，其本質是一種采用高階連續的B樣條基函數作為待求物理場的插值函數，并使用四叉樹／八叉樹自適應積分策略處理被結構邊界所分割的網格單元的虛擬區域法。

為求解結構Ωphy上的線彈性力學問題，虛擬區域法先將結構Ωphy嵌入到規則的計算域Ω中（如圖11所示），然后采用固定網格離散Ω，并基于變分原理建立離散后的平衡方程組，最后求解此方程組得到結構的待求物理場。

圖11 結構域Ωphy擴展到簡單規則的計算域ΩFig.11 Embedded domainΩextendedfrom structural domainΩphy

定義在計算域Ω上的線彈性問題為

式中：σ為柯西應力張量；u為位移向量；g為邊界ΓD上的給定位移；f為體積力向量；t和n分別為邊界ΓN上的外界力和單位外法向量。標量因子λ定義為

相容方程及本構方程描述為

式中：ε為應變張量；D為Ωphy上的彈性本構張量。

定義以下2個可容位移場空間

式中：H1（Ω）代表一階希爾伯特空間。對于嵌入域Ω，式（5）的等效積分弱形式可以表示為

式中：

假定Shu是Su的有限維子空間，Shv是Sv的有限維子空間，因此上述問題的離散問題為尋找一個離散的解uh∈Shu，使其能夠滿足

類似于有限單元法，計算域Ω可以離散成獨立于結構域Ωphy的規則的網格單元（這些單元在BSFCM中也稱為“胞元”），根據文章第1節所述的方法，可以將這些單元區分為內部單元、邊界單元和外部單元。在BSFCM中，利用定義在固定網格上的高階B樣條形函數來離散插值位移向量uh和vh。對于三維問題，三變量B樣條基函數可以表示為

式中：Ni，p（ξ）、Nj，q（η）和 Nh，r（γ）分別是p、q和r次單變量B樣條基函數，分別定義在參數域坐標系下 的 Ξ ＝ ｛ξ1，ξ2，…，ξns＋p＋1｝，H ＝ ｛η1，η2，…，ηms＋q＋1｝和ψ ＝ ｛γ1，γ2，…，γks＋r＋1｝節點矢量上，他們的個數分別為ns、ms和ks。對于計算域Ω內的任意一點P ＝ ［x，y，z］T，可由以下的映射關系確定

式中：Ps＝ ［xs，ys，zs］T表示三維空間中的第s個控制點。計算域Ω內的點PΩ，可由以下簡單的線性映射確定

式中：x0和xend分別為計算域在x方向上的最小和最大值；y0和yend分別為計算域在y方向上的最小和最大值；z0和zend分別為計算域在z方向上的最小和最大值。

本文中采用統一開放的節點矢量來定義B樣條基函數［26］，根據考克斯-德布爾遞推公式，定義在Ξ上的第i個p次B樣條基函數Ni，p（ξ）的表達式為

由式（8）可知，對于vh而言，基函數Ms在邊界ΓD上的值需為0，則存在一組向量A，使得vh滿足線性插值形式

假設U是一組待求位移向量，表示控制點Ps的位移。對于給定的函數gh∈Sh，如果滿足gh｜ΓD＝g，則位移的離散形式可以表示為

特別地，對于施加齊次Dirichlet邊界條件的情況，g＝0。根據Riz-Galerkin方法，將式（16）和式（18）代入式（11），可得到BSFCM 中的剛度方程

式中：K是整體剛度矩陣；F是整體載荷矩陣；具體為

式中所表達的總體剛度矩陣K和總體載荷矩陣F在BSFCM中均由單元的剛度矩陣和載荷向量組裝而成。根據第1節所確定的單元屬性不同，其剛度矩陣計算的方法也不同。對于外部單元，由于其對K和F沒有貢獻，因此無需計算；對于內部單元，直接按照傳統有限元方法的單元積分方式進行計算；對于邊界單元，根據1.3節的細化策略，采用自適應積分方法沿著CAD模型邊界加密高斯積分點，對邊界單元中的有材料部分進行高精度的積分計算。

式（20）中，B是應變位移矩陣（幾何矩陣），表達式為

2.2 Dirichlet邊界條件施加

不失一般性，本文僅考慮BSFCM框架下齊次Dirichlet邊界條件的精確施加策略，并選擇了簡單且能精確施加齊次Dirichlet邊界條件的加權B樣條方法［27－28］。該方法的本質就是將B樣條基函數Ms在齊次位移約束邊界ΓD上的值修正為0，修正后的式（18）變為

式中：權函數ω（x，y，z）是隱式函數，且滿足ω（x，y，z）｜ΓD＝0。對于實際CAD模型的約束邊界，可以構造其隱式函數表達式，并將其直接作為權函數。

在使用權函數ω（x，y，z）修正B樣條基函數M之后，式（21）定義的幾何矩陣B也相應地變化為

此外，在求解一些對稱問題時，常對結構的1／2、1／4或1／8結構模型進行力學分析。此時的位移邊界條件需施加在對稱邊界上，如：在邊界ΓDX上要求x方向的位移為0，在邊界ΓDY上要求y方向的位移為0，在邊界ΓDZ上要求z方向的位移為0，也即分別構造在ΓDX上的權函數ωx＝0，在ΓDY上的權函數ωy＝0和在ΓDZ上的權函數ωz＝0。因此，基函數矩陣變化為

3 CAD模型分析流程

綜合本文第2節和第3節的內容，可以建立一種由結構CAD模型至固定網格仿真分析的設計框架，如圖12所示。

圖12 固定網格中CAD模型分析流程Fig.12 Analytic flow of CAD model in fixed grid

首先通過設計軟件建立結構的CAD模型，并將其保存成采用邊界形式表示的STL格式；然后將CAD模型嵌入建好的固定網格中，并采用射線交點法進行單元識別，同時對識別出的邊界胞元采用八叉樹技術進行細化，以提高有限元分析模型的建模精度；再者對固定網格區域采用BSFCM，并通過加權B樣條方法施加Dirichlet邊界條件和載荷；最后求解獲得CAD模型的有限元仿真結果。

4 數值算例

4.1 平板帶孔結構

平板帶孔結構作為計算力學中的經典算例，相關學者對此進行了深入的研究［8］。不失一般性，本文選用線彈性平板帶孔結構算例來說明文中建立的基于邊界描述的CAD模型快速分析方法的有效性。假定其楊氏模量和泊松比分別為2.074×105MPa和0.3。如圖13所示，平板帶孔結構沿±z向承受的均布載荷為t1＝22.5MPa，沿±y向承受的均布載荷為t2＝45MPa。為提高計算效率，將平板帶孔結構簡化為1／8模型進行分析，簡化后模型的尺寸為R＝1mm，L＝4mm，H＝0.5mm。

簡化后的平板帶孔結構模型的STL格式如圖14所示。

圖13 均布載荷作用下的平板帶孔結構Fig.13 Plate with a circular hole under uniform loads

圖14 平板帶孔結構簡化模型的STL格式Fig.14 STL format of simplified model for plate with a circular hole

根據平板帶孔結構的特點，建立一個尺寸略大于平板帶孔結構且分辨率為8×8×2的固定網格，然后將平板帶孔結構嵌入固定網格中，并應用射線交點法進行單元屬性識別，建立固定網格中的有限元分析模型，如圖15所示。

為提高固定網格中有限元分析模型的精度，采用八叉樹技術對邊界單元進行分層細化，然后在固定網格中應用BSFCM，并通過加權B樣條方法施加齊次Dirichlet邊界條件和載荷，可以得到平板帶孔結構的應力云圖如圖16所示。由圖16可知，最大應力發生在孔的上邊緣處，且最大應力值約為133.04MPa。

當網格單元數量較少時，網格單元數量不同所得到的最大應力也不同，但當網格單元數量滿足一定的值時，最大應力趨于穩定，如圖17所示。圖17中網格單元數量從16（4×4×1）～4 096（32×32×4）的變化過程中，當網格數量超過128（8×8×2）時，平板帶孔結構的最大應力逐漸收斂。

圖15 平板帶孔結構在固定網格中的分析模型Fig.15 Analytic model for plate with a circular hole in fixed grid

圖16 固定網格中平板帶孔結構的應力云圖Fig.16 Contour of Von-Mises stress of plate with a circular hole in fixed grid

作為參考，當采用傳統有限元方法求解上述平板帶孔結構的應力分布時，使用線彈性四面體網格或六面體網格離散結構區域，獲得的單元數量為2 637個，如圖18所示。平板帶孔結構的應力云圖如圖19所示，最大應力發生在孔的上邊緣處，最大應力值約為130.49MPa。在傳統有限元方法中，最大應力隨網格數量的變化曲線如圖20所示。

圖17 固定網格中不同網格數量下的最大應力Fig.17 Maximum Von-Mises stress with different elements in fixed grid

圖18 平板帶孔結構在傳統有限元方法中的分析模型Fig.18 Analytic model for plate with a circular hole in classical finite element method

圖19 傳統有限元方法中平板帶孔結構的應力云圖Fig.19 Contour of Von-Mises stress of plate with a circalarhole obtained by classical finite element method

圖20 傳統有限元方法中不同網格數量下的最大應力Fig.20 Maximum Von-Mises stress with different elements in classical finite element method

對比本文方法和傳統有限元方法可知，本文方法可以實現平板帶孔結構CAD模型至固定網格中分析模型的直接快速轉換，并且通過較少的單元即可獲得與傳統有限元方法同等精度的計算結果。由圖16和圖19可知，2種方法獲得的應力云圖基本一致，最大應力的相對誤差約為2.2%。

4.2 連接結構

為了進一步說明本文所提方法的有效性，下面通過一個更復雜的例子——航空連接結構進行驗證。連接結構的CAD模型及其載荷和約束情況如圖21所示。

連接結構的主要尺寸為：L11＝80mm，H11＝80mm，H12＝20mm，H13＝20mm，R11＝5mm，R12＝10mm。連接結構所承受的載荷為集中力，位于連接結構左端孔中心處，其值為Fzz＝5000N。連接結構的約束條件為其右端兩個孔的位移約束。同時假定連接結構的楊氏模量和泊松比分別為2.074×105MPa和0.3。

連接結構CAD模型的STL格式如圖22所示。建立一個尺寸略大于連接結構且分辨率為25×60×50的固定網格，然后將連接結構嵌入固定網格中，并應用射線交點法進行單元屬性識別，建立固定網格中的分析模型，如圖23所示。

圖22 連接結構的STL格式Fig.22 STL format of connector structure

對邊界單元采用八叉樹進行分層細化，提高分析模型對CAD模型的近似精度。然后在固定網格中應用BSFCM，并通過加權B樣條方法施加Dirichlet邊界條件和載荷，可以得到連接結構的應力云圖如圖24（a）所示。由圖24（a）可知，最大應力發生在連接結構左端過渡處，且最大應力值約為222.16MPa。采用傳統有限元方法計算得到連接結構的應力云圖如圖24（b）所示。由圖24（b）可知，最大應力同樣發生在連接結構的左端過渡處，且最大應力值約為219.02MPa。由圖24可知，2種方法所得的應力云圖基本一致，最大應力相對誤差約為1.5%。

圖23 連接結構在固定網格中的分析模型Fig.23 Analytic model for connector in fixed grid

圖24 連接結構應力云圖Fig.24 Contour of Von-Mises stress of connector

5 結 論

本文研究了復雜結構在固定網格中的力學分析問題，提出了一種基于CAD模型邊界表示的結構模型至固定網格中分析模型的快速轉換方法，在保證模型轉換精度的同時，提高了模型轉換的效率、降低了模型轉換的復雜度。同時結合計算力學領域新興的加權B樣條有限胞元方法，建立一種固定網格中結構CAD模型的快速分析框架，實現了結構CAD設計與分析的有效集成，為結構設計與分析的一體化建設提供了一種思路。