面向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.