面向CAE的簡化模型誤差評估與信息重用

唐建國 曹魏娟 高曙明

浙江大學CAD＆CG國家重點實驗室，杭州，310058

0 引言

在工程分析中，復雜CAD模型的網格生成和有限元分析非常困難，通過模型簡化得到合理的分析模型在CAD/CAE集成系統中十分重要。模型簡化一般基于幾何的標準，如縱橫比、曲率、重心的變化量、體積比等［1－4］。然而，小尺寸的特征因為邊界條件的變化或者所做的分析不同，對分析結果的影響也不一樣。因此，評估模型簡化引起的分析誤差非常必要。Gopalakrishnan等［5－6］提出了特征靈敏度的概念，根據伴隨理論，利用邊界元的方法，求解孔洞、槽等負特征被抑制后對感興趣區域的影響。該方法計算效率高，但無法處理正特征被抑制的情況。Ferrandes等［7］從應變能的角度，提出了判斷抑制特征對分析結果影響的方法。該方法的不足之處在于局部計算時切割邊界面的選擇依賴于分析者的經驗。龔志輝等［8］為提高模具成形精度，提出了一種仿真誤差補償模型。該方法通過迭代計算來獲取模具補償面，計算代價較大。

劉曉平等［9］根據不同的分析需求，提出了多態模型的概念，根據特征體積比來探索各簡化模型與計算精度之間的關系。此外，劉曉平等［10］還通過比較模型簡化前后的網格剖分及場函數之間的差異，依據仿射等價協調元和等參協調元的插值精度定理來估計模型簡化前后分析解的差異，從而判斷模型簡化策略的合理性［10］。

通過有限元模型的再分析技術也可以實現對模型簡化的評價。例如利用Sherman－Morrison公式對模型剛度矩陣的局部變化快速求逆［11］，該方法在模型簡化前后的剛度矩陣變化不大時效率較高。Wu等［12－13］對結構體再分析問題進行了研究，根據模型變化將模型的有限元網格分為節點自由度增加和減少兩種情況，將簡化模型剛度矩陣作為預處理矩陣，通過預條件共軛梯度法求解。該方法的分析效率較高，缺點是無法處理結構體變動同時引起自由度增加和減少的情況。此外，Kirsh等［14－15］提出了近似組合法，將二次項級數作為簡化基更高精度的表示，由于多次將初始模型的總剛度矩陣參與計算，因此計算代價較大。另外，ANSYS中的生死單元法［16］可對模型變動進行分析。上述的再分析技術均需要對初始模型網格化，并對整個模型求解剛度矩陣，計算代價大。

本文結合原模型與簡化模型的邊界特性，研究了抑制特征對簡化模型分析結果的影響。當簡化模型的分析結果不滿足分析要求時，根據特征對分析結果的影響，調整簡化模型的特征組成并得到合理的簡化結果。

目前，大型稀疏矩陣的分解軟件（GSS、UMFPACK、PARDISO、SUPERLU 等）使得百萬階的矩陣在PC上的分解成為可能。因此本文還研究了初始簡化模型采用直接法求解后，如何重用已有的有限元網格、單元剛度矩陣以及分解后的三角矩陣等分析信息，從而減小分析計算的代價。

1 方法概述

如圖1所示，簡化模型的誤差評估以及信息重用的流程分為3個階段。第一階段，對輸入的CAD模型進行預處理：通過特征識別，將模型分為簡化模型和抑制特征，并對它們進行粘合操作和網格剖分。對抑制特征與簡化模型粘合操作的目的是確保網格剖分時接觸面處網格的一致性。然后對初始簡化模型進行有限元分析。第二階段，計算抑制特征對分析的影響，根據影響大小確定最終簡化模型。第三階段，重用初始簡化模型中的網格單元、剛度矩陣以及總剛度矩陣分解后的三角陣。通過矩陣的局部分解方法，可快速得到最終簡化模型剛度矩陣的分解結果，并重用原有邊界條件和材料屬性，完成對最終簡化模型的分析。

圖1 簡化模型的誤差評估與信息重用方法的流程

由圖1可以看出，簡化模型的分析信息作為后續計算的輸入，一方面可用于評價簡化模型的誤差，另一方面，可以提高最終簡化模型的分析效率。

鑒于目前已有大量的特征識別與模型簡化方法的研究［4，17］，本文研究的重點不在模型簡化的具體方法，而是如何評估模型簡化對分析結果的影響以及如何重用已有的分析信息，從而提高最終簡化模型的分析精度和分析效率。

2 抑制特征對簡化模型分析結果影響的評估

定義模型簡化前后的共同區域為公共域、模型簡化時被抑制的特征為抑制特征、抑制特征與簡化模型的公共面為接觸面。將原模型與簡化模型的公共域作為子結構：抑制正特征后，簡化模型的子結構對應的特征接觸面上的外載荷為零；抑制負特征后，簡化模型的子結構對應的特征接觸面上產生了外力。綜上所述，模型簡化使得子結構的邊界條件發生變化，從而產生了分析誤差。

2.1 評價抑制特征對分析結果影響的理論指標

簡化模型與原模型的公共域是分析者關注的區域，因此將特征抑制前后公共域應變能的變化ηi作為評估抑制特征對分析結果影響的理論指標，即有

式中，ΔWi為特征i抑制前后公共域應變能的變化量；Wm為簡化模型在公共域的應變能；ΓF為原始模型的受力邊界；F、UF分別為原模型的外載荷向量和外載荷邊界上的位移向量；p為負特征的總數；Γsi為抑制負特征i與簡化模型的接觸面；R′s為簡化模型的子結構在負特征接觸面上的等效外力；U′s為簡化模型在負抑制特征接觸面的位移向量。

當i為某個正特征時，ΔWi標記為

當i為某個負特征時，ΔWi標記為

式中，Γa、Γs分別為正負特征與簡化模型的接觸面；以模型簡化前后的公共域為子結構，Ra為原模型的子結構在正特征接觸面上的等效外力，Ua、Us分別為原模型在正負特征接觸面上的位移向量；U′a為簡化模型在正抑制特征接觸面的位移向量。

ΔWi是根據簡化模型與原模型在公共域邊界條件的不同，結合貝蒂定理推導而得（推導原理參見文獻［18］）的。本文基于已有研究［18］，采用公共域應變能的變化量來評估抑制特征對分析結果的影響，更精確地定義簡化模型的誤差評估指標。

2.2 一種實用的估算方法

求解式（1）中的ηi需要計算ΔWi中的Ra、Ua、Us。這3個值通常采用子模型法求取。子模型法中分割區域的初選依賴于分析經驗，分割邊界位置的確定需要通過反復計算予以校驗，這對于擁有多個細節特征的模型并不可取。因此，本文對式（1）作了近似：

式中，ΔW′a、ΔW′s分別為根據簡化模型計算出的抑制正負特征的應變能；ηa、ηs分別為正負抑制特征對簡化模型公共域分析影響的估計指標；－R′a、－R′s分別為以U′a、U′s為邊界條件求解抑制特征得到的相應接觸面的外力向量。

如表1和表2所示，圖2a和2b模型中抑制特征的估算影響與理論值比較接近，具有較高的精度。

表1 圖2a模型中抑制特征對簡化模型的理論與估計影響值

表2 圖2b模型中抑制特征對簡化模型的理論與估計影響值

圖2 模型的特征抑制

3 信息重用方法

根據各個特征對分析結果的影響，在初始簡化模型不滿足分析精度的情況下重新確定出最終簡化模型。在對最終簡化模型的分析求解過程中，重用初始簡化模型的分析信息，實現快速求解。圖1所示的流程表明，最終簡化模型重用了簡化模型公共域的有限元網格以及相應的剛度矩陣和剛度矩陣分解后的三角矩陣。因此，最終簡化模型剛度矩陣的分解僅限于局部的有限單元。

3.1 剛度矩陣局部分解原理

矩陣分解的原理如圖3a所示，矩陣分解為下三角陣L和上三角陣U，分解的推進方向為對角線方向：即先分解對角線元素，然后是該元素對應的列和行剩余的未分解元素（圖3a中粗線表示正處于分解的行和列的元素）。剛度矩陣中尚未分解的元素值并不會影響已經分解的L、U元素的值。由于剛度矩陣為對稱正定稀疏矩陣，采用Cholesky分解法分解時，僅需要分解矩陣中的上三角或下三角矩陣元素。

如圖3b所示，根據剛度矩陣的稀疏性，如果簡化模型的剛度矩陣元素處于原模型總剛度矩陣的左上角，而抑制特征處于右下角位置，則在簡化模型剛度矩陣分解結果已知的前提下，求解含有抑制特征的模型只需分解抑制特征的剛度矩陣元素。對簡化模型和被抑制特征的編號預處理，可以調整被抑制特征和接觸面對應剛度矩陣元素到總剛度矩陣的右下側位置。

圖3 模型的矩陣信息重用與局部化分解

以圖3b為例，初始簡化模型分解得到三角矩陣L、U以及其他數值項Lm、Um，在最終簡化模型的剛度矩陣中分解結果不變。因此，初始簡化模型三角矩陣的可重用部分為L、Lm、U和Um。其中，L和U為初始模型公共域內部節點的分解矩陣，Lm和Um為總剛度矩陣中特征接觸面節點與公共域內部節點之間對應自由度元素的分解結果。

對稱正定的剛度矩陣的Cholesky分解法可以用下面的公式給出：設A∈Rn×n，有A＝LLT分解，其中，L＝（lij）n×n是下三角矩陣，分解過程中只需對剛度矩陣A的對稱部分（對應L下三角位置的元素）采用列優先分解，則有

從式（4）可以看出，矩陣中對角線位置的任一元素aii的分解結果lii僅僅依賴于第i行的已分解值，而第i列的其他值aki的分解結果lki依賴于第k行以及第i行的已有分解結果。簡而言之，下三角矩陣某一列元素的分解僅與該列首行的分解結果以及該元素對應行的已有分解值有關。

由剛度矩陣的性質可知，不鄰接節點的對應位置矩陣元素值為零。因此，將特征接觸面的對應節點作為連接簡化模型與抑制特征對應網格的鄰接點，其他非鄰接點在矩陣中相應位置處的矩陣元素值置零。如圖3b所示，Ln和Un為總剛度矩陣中特征接觸面上節點之間的自由度對應元素的分解值，只有Lm、Um參與了矩陣元素Ln和Un的計算。Ln和Un的局部計算減小了最終簡化模型剛度矩陣分解的計算代價。

圖3b中，矩陣元素的布局對簡化模型的節點自由度編號提出了特殊要求。如果編號與位置不一致，則需要行列變換。理想的自由度編號使得抑制特征的剛度矩陣元素總是位于總剛度矩陣的右下位置，這樣局部分解可從接觸面上的節點自由度開始，計算量最小。

3.2 信息重用與矩陣局部分解的算法描述

初始簡化模型的分析完成之后，需要評估抑制特征對分析結果的影響，然后生成最終簡化模型并完成最終簡化模型的分析計算。算法流程如下：

（1）計算抑制特征對分析結果的影響值ηa和ηs。首先對簡化模型進行預處理，使特征接觸面的節點自由度編號比公共域內部節點自由度編號大；然后對初始簡化模型求解，計算初始簡化模型的位移UF、U′a和U′s。以接觸面位移U′a、U′s為邊界條件，分別計算抑制正負特征接觸面上的力R′a和R′s。最后，根據前兩步的計算結果計算正負特征的影響值ηa和ηs。

（2）初始化剛度矩陣Kn。保留初始簡化模型中接觸面節點所在自由度的剛度，并將其他節點自由度的值置零。初始化后的Kn僅僅包含初始簡化模型公共域的接觸面處節點的剛度矩陣元素值。

（3）剛度矩陣局部分解。這一步的執行流程如圖4所示，首先讀入簡化模型剛度矩陣的分解值Lm（Lm是矩陣分解值中一小部分，即接觸面節點與簡化模型內部節點的關聯項對應值，參見圖3b的接觸面重用部分）；然后逐一讀入特征影響值，根據精度要求判斷該特征是否可抑制，對于不能抑制的特征，讀入其剛度矩陣并組裝到總剛度矩陣Kn中；最后根據讀入的Lm，對Kn矩陣進行分解。

圖4 剛度矩陣局部分解算法流程

（4）生成最終簡化模型的三角矩陣。將初始簡化模型公共域對應的剛度矩陣分解值與Kn矩陣的分解值合并，形成最終簡化模型總剛度矩陣的分解矩陣。

（5）讀入已有的邊界條件和材料屬性信息，完成最終簡化模型的有限元分析。

4 分析實例

本節將通過分析實例展示簡化模型的誤差評估與信息重用過程。實驗的硬件環境為：2.66GHz的Intel CPU，4GB內存；軟件平臺為Windows XP，ANSYS 11.0，以及 APDL （ANSYS parametric design language）。

如圖5a所示，假設彈性模量為2×1012Pa，泊松比為0.3，受力面的壓力值為6MPa，原模型的上表面受壓力、底面固定（圖5b）。抑制特征的標識如圖5a、圖5b所示，圖5c所示為初始簡化模型。

圖5 簡化模型誤差評估與信息重用工程實例1

表3展示了精確評價指標與估計指標的比較結果。其中，橫坐標為特征編號，縱坐標為抑制特征對簡化模型的影響值。可以看出，估計指標與精確指標相比，具有較強的魯棒性。其中，抑制特征3和4對簡化模型分析精度的影響較大，因而不能抑制，圖5d為調整后得到的最終簡化模型。

表3 實例1各特征的精確評價指標與估計指標

圖6、圖7是工程中的兩個CAD模型。圖6、圖7中，彈性模量和泊松比的設置同圖5的模型，邊界條件如下：圖6中，模型受到6MPa的壓力，位移約束為左側表面固定；圖7中，模型的局部表面受到6MPa壓力，底部表面為固定約束。計算出抑制特征對分析結果的影響值后，得到最終簡化模型，如圖6d和圖7d所示。

圖6 簡化模型誤差評估與信息重用工程實例2

圖7 簡化模型誤差評估與信息重用工程實例3

剛度矩陣局部分解的方法對簡化模型的預處理提出了特殊的要求。初始簡化模型和最終簡化模型必須保證公共域網格一致，確保域內單元剛度矩陣的值不變。另外，在分析初始簡化模型之前，需要對有限元節點的編號進行調整，使得特征相交面節點編號值介于簡化模型與被抑制特征的節點編號之間。這樣最終簡化模型總剛度矩陣元素的布局和圖3b一致，便于總剛度矩陣的局部分解。

在最終簡化模型的分析中，效率的提高主要在剛度矩陣的分解階段。下面采用Cholesky分解法對圖5～圖7的初始簡化模型和最終簡化模型的剛度矩陣分解進行比較。參與比較的剛度矩陣是從ANSYS的full文件中導出的Harwell－Boeing文件。分解的節點數以及時間比如表4所示，其中，可重用節點數為公共域中的節點數，局部分解的節點數為各新增特征的節點數。局部分解從抑制特征接觸面上的最小節點編號開始，參與局部分解的剛度矩陣元素數量大大減小；最終簡化模型的局部分解與完全分解相比，所需時間較短。

表4 圖5～圖7中最終簡化模型剛度矩陣的局部分解與完全分解

5 結語

提出了一種面向CAE的模型簡化誤差評估與信息重用方法。首先對初始簡化模型的網格節點編號進行預處理，確保再分析時的計算代價最小；然后對初始簡化模型進行有限元分析，并將分析結果作為被抑制特征的邊界條件，評估各抑制特征對模型分析結果的影響；最后通過重用初始簡化模型剛度矩陣的分解結果，對最終簡化模型剛度矩陣進行局部分解，實現了對最終簡化模型的快速分析。該方法既保證了最終簡化模型的分析可靠性，又提高了分析效率。同時，這種分析信息重用的方法也可用于模型修改設計的分析或同系列產品的分析。

作為將來的工作，一方面，繼續研究簡化模型局部求解的高效方法，研究迭代法求解中的信息重用方法；另一方面，對特征之間的依賴關系，特別是模型簡化時抑制特征對其他特征的局部影響進行研究，以滿足設計者對細節特征設計和評價的深層次需求。

［1］Foucalt G，Marin P，Léon J C.Mechanical Criteria for the Preparation of Finite Element Model［C］//Proceedings of the 13th International Meshing Roundtable.Williamsburg，2004：413－426.

［2］Lee J Y，Lee J H，Kim H，et al.A Cellular Topology－based Approach to Generating Progressive Solid Models from Feature－centric Models［J］.Computer－Aided Design，2004，36（3）：217－229.

［3］Lee S H.A CAD－CAE Integration Approach Using Feature－based Multi－resolution and Multi－abstraction Modelling Techniques［J］.Computer－Aided Design，2005，37（9）：941－955.

［4］Thakur A，Banerjee A G，Gupta S K.A Survey of Model Simplification Techniques for Physicsbased Simulation Applications［J］.Computer－Aided Design，2009，41（2）：65－80.

［5］Gopalakrishnan S H，Suresh K.A Formal Theory for Estimating Defeaturing－induced Engineering Analysis Errors［J］.Computer－ Aided Design，2007，39（1）：60－68.

［6］Gopalakrishnan S H，Suresh K.Estimating the Impact of Large Design Changes on Field Problems［C］//Proceedings of ACM Symposium on Solid and Physical Modeling.Beijing，2007：205－215.

［7］Ferrandes R，Marin P M，Léon J C，et al.A Posteriori Evaluation of Simplification Details for Finite Element Model Preparation［J］.Computers and Structures，2009，87（1/2）：73－80.

［8］龔志輝，李光耀，鐘志華，等.基于仿真誤差補償模型的回彈補償新方法［J］.中國機械工程，2008，19（11）：1358－1362.

［9］劉曉平，金燦，李書杰.有限元仿真分析軟件中建模的多態機理研究［J］.系統仿真學報，2007，19（3）：538－542.

［10］金燦，劉曉平.拋物問題中面向協調元的模型態誤差估計方法［J］.計算機輔助設計與圖形學學報，2011，23（4）：656－666.

［11］Householder A S.A Survey of Some Closed Form Methods for Inverting［J］.Matrices，1957，5（3）：155－169.

［12］Wu B S，Limb C W，Li Z G.A Finite Element Algorithm for Reanalysis of Structures with Added Degrees of Freedom［J］.Finite Elements in Analysis and Design，2004，40（13/14）：1791－1801.

［13］Wu B S，Li Z G.Reanalysis of Structural Modifications due to Removal of Degrees of Freedom［J］.Acta Mechanica，2005，180（1/4）：61－71.

［14］Kirsch U.A Unified Reanalysis Approach for Structural Analysis，Design，and Optimization［J］.Structural and Multidisciplinary Optimization，2003，25（2）：67－85.

［15］Kirsch U，Papalambros P.Exact and Accurate Reanalysis of Structures for Geometrical Changes［J］.Engineering with Computers，2001，17（4）：363－372.

［16］Erdogan M，Ibrahim G.The Finite Element Method and Applications in Engineering Using ANSYS［M］.Berlin：Springer，2006.

［17］周廣平，高曙明，崔秀芬，等.面向網絡基于分治策略的加工特征識別算法［J］.計算機輔助設計與圖形學學報，2005，17（6）：1152－1158.

［18］唐建國，高曙明，藺宏偉.面向CAE的模型簡化中的誤差評估與邊界補償［J］.計算機輔助設計與圖形學學報，2010，22（8）：1308－1315.