专利名称：有限元建模中基于特征抑制的模型态生成方法
技术领域：
本发明涉及计算机辅助建模方法，更具体地说是一种应用在有限元建模中的模型态生成方法。
背景技术：
计算机对于提高社会生产力发挥了越来越重要的作用，特别是计算机辅助设计CAD、辅助工程CAE、辅助制造CAM在工业界日益成熟和普及，极大地提高了工业设计和生产的效率。采用有限元分析及基于有限元分析的优化处理，能够改进产品的结构设计，使产品在满足强度和刚度的情况下具有最合理的结构。
有限元分析FEA早在50年代首先在连续体力学领域中，比如在飞机结构静态和动态特性分析中得到应用，随后很快应用于分析热传导、电磁场和流体力学等连续性问题。
有限元分析的物理实质是用由有限个在节点处相连接的单元组成的组合体近似替代一个连续体，从而把连续体的分析问题转化为单元分析和单元组合的分析问题。有限元分析和计算机的强大的数据处理能力的结合，使过去不可能进行的一些大型复杂结构的分析变成了常规的计算任务。
进行有限元分析首先要建立待分析对象的几何模型。实际模型通常十分复杂，主要体现在其包含的特征种类和数量巨大，除了常见的孔、槽特征，还含有大量的过渡特征。如果建模的精度太高，无疑会为计算带来极大的负担，尤其是过渡特征的存在会大大增加计算量，因此，需要对实际模型进行特征简化。在对实际模型中的特征进行简化时涉及以下两个核心问题一是如何对细小的孔、槽以及过渡特征简化后对原计算结果造成的影响进行估计；二是如何估计出简化后的模型所需消耗的计算时间，以判断简化后模型是否能够提升计算效率。
针对此问题，在期刊《Computer-Aided Design》、2001年33卷13期第925-934页上，A Sheffer.的文章《Model Simplification for Meshing Using Face Clustering》(以下简称文献1)中，分析了几种常见的过渡特征对划分网格造成的消极影响；但并未针对过渡特征进行处理，更没有估计特征简化后对结果造成的影响和简化后的模型所需消耗的计算时间。
《Computer-Aided Design》2002年34卷22期第109-123页上，H Zhu，C H Menq的另一篇文献《B-Rep Model Simplification by Automatic Fillet/Round Suppressing for EfficientAutomatic Feature Recognition》(以下简称文献2)和《Computer-Aided Design》2004年第5期第24-28页上崔秀芬的文章《An Efficient Algorithm for Recognizing and SuppressingBlend Features》(以下简称文献3)中虽然提出了较为完善的过渡特征识别和简化方法，但是也没有估计过渡特征被抑制后对结果造成的影响和简化后的模型所需消耗的计算时间。
浙江大学学报自然科学版2006年第7卷第9期第1535-1543页上LEE Sang Hun的文章《Feature-based multiresolution techniques for product design》(以下简称文献4)提出了产品制造中基于细节层次LOD的建模方法，在这种方法中，将几何模型中的特征简单地分为加特征和减特征，通过用户决定LOD的值来决定模型中各个特征是否存在，以达到减小计算量的目的。但这种方法过分依赖用户经验；同时由于模型层次间几何差别过大，选择模型层次后还需要反复实验以验证模型是否合适，过程太复杂，且没有给出层次模型对应的估计误差，也没有估计简化后的模型需消耗的时间。
合肥工业大学2005届硕士学位论文《科学计算中的几何多态模型初探》(以下简称文献5)中，以中子输运程序MCNP为研究对象，对科学计算领域中采用某个模型态替换原始模型产生的计算误差的影响要素做了初步研究，并以实验数据为基础，总结了计算误差与特征体积之间的规律。根据总结的规律，采用特征抑制的方法形成模型态。每个模型态对应一个估计误差，给用户提供了一个很好的辅助建模功能。其研究方法具有一定的可借鉴性，但是存在领域局限性，即误差要素的总结方法没有考虑有限元领域的自身特性，结论的可靠性和通用性不强，且并未给出模型态计算效率的估计方法。
合肥工业大学2007届硕士学位论文《有限元领域中模型态的选择研究》(以下简称文献6)中，以实验为基础，从实验数据中总结了热传导领域中采用某个模型态替换原始模型产生的计算误差规律，并对影响计算时间的因素进行了初步分析和总结；在给出计算精度和计算时间的形式化表达式的基础上构造了遗传算法的适应度函数，通过一定数量的迭代次数选择出适应度最高的模型态；最后根据模型态编码使用特征抑制算法对原始模型进行改造得到最终几何模型。论文中所记载的模型态生成的具体步骤是 a、针对有限元不同的应用领域，建立只包含一个待简化特征的简单模型，这里的待简化特征是加工零件中最常见的槽、孔或柱面过渡特征，建立包含它们的模型，并构造稳态热传导实验方案。
b、按照特征体积参数和位置参数构造实验方案，即对含有若干组不同特征体积比V和特征与载荷面的距离D的槽、孔特征的模型进行分析计算，根据计算结果分析误差比EP、特征体积比V和特征与载荷面的距离D的内在关系，对其关系进行拟合。
c、用户利用自身领域知识判断模型中的过渡特征存在与否对结果的影响程度，若某过渡特征对计算结果的影响程度可以被忽略，则以网格划分的质量为标准，对过渡特征确定其被删除的几何属性阈值，避免计算时出现不可预知的问题；不能删除的过渡特征将被保留。
d、根据拟合出的槽、孔特征的特征体积比V、特征与载荷面的距离D与误差比EP之间的函数关系来确定为进一步提高计算效率可以抑制的槽、孔特征。
e、构造反应精度要素的表达式P＝1-∑EPi和时间要素的表达式，并基于此构造遗传算法适应度函数，选定遗传算法其他参数。这里EPi是第i个特征被抑制后产生的误差比，
是平均特征体积比，EPmax是用户输入的可接受的最大误差。
f、根据编码对几何模型处理，得到最终选择的模型态。
在上述步骤中a-d中，是通过拟合实验数据来获得误差比EP、特征体积比V、特征与载荷面的距离D之间的函数关系。但是，这种方式存在有两方面的不足首先，误差比、特征体积比和特征与载荷面距离之间的关系仅通过有限次实验的数据得出，具有随机性；其次，不能保证影响误差要素总结的全面性。在步骤(e)中，根据模型态中含有的特征的平均体积比大小来估计模型态需消耗的时间，显然缺少客观性和精确性。
发明内容
本发明是为避免上述现有技术所存在的不足之处，提供一种有限元建模中基于特征抑制的模型态生成方法，全面考虑模型简化后对有限元计算结果影响的因素，从而实现模型中各特征的变化对有限元计算结果敏感性的客观评价，建立模型中各要素的取舍与变化对计算结果的敏感性评价体系，并给出一种客观精确的模型态消耗时间的估计方法，在此基础上生成模型态。
本发明解决技术问题采用如下技术方案 本发明有限元建模中基于特征抑制的模型态生成方法的特点是按如下步骤进行 a、由用户提供待简化模型，所述待简化模型具有在进行有限元分析时的载荷条件，以及网格剖分条件，所述网格剖分条件为最大网格直径hmax和最小网格直径hmin； b、识别所述待简化模型中的可简化特征，并存入特征链表； c、在所述特征链表中选择待简化特征，确定待简化特征数量n；按任意顺序对所选择的n个待简化特征进行编号，计算获得所选择的每个待简化特征的体积Vi和位置参数Pi； d、依据每个所选择的待简化特征的体积Vi、位置参数Pi，以及每个所选择的待简化特征在其被抑制后的网格直径改变量，获得对应的待简化特征在其被抑制后的敏感值Ei＝ViPiHi；以所有被抑制特征的敏感值之和表征模型态的误差水平； e、通过对模型态网格单元进行数量的估算，估计出模型态进行网格划分后的单元数量Ne，所述单元数量Ne与模型态刚度矩阵阶数线性相关，以所述单元数量Ne反映模型态所需的有限元计算时间； f、根据模型的误差水平和单元数量Ne，采用遗传算法进行模型态的选择，获得简化模型编码； g、应用特征抑制法生成对应于所述简化模型编码的简化模型。
本发明方法的特点也在于 所述步骤c中简化特征位置参数Pi通过如下方式获得 a、给定每个待简化特征简化后对有限元计算结果影响的权重，设为q1，q2，…qn，其中qi∈
，(i＝0，1，…n)，qi＝0表征第i个待简化特征存在与否对有限元计算结果没有影响；qi＝1表征第i个待简化特征存在与否对有限元计算结果影响最大，qi越靠近1表征特征对有限元计算结果的影响越大； b、令i＝1； c、设第i个待简化特征到载荷面的距离di； d、计算第i个待简化特征的位置参数Pi＝qi/di； e、令i＝i+1，重复步骤c～d，直至模型中所有待简化特征的位置参数Pi计算完毕。
所述步骤d中网格直径改变量按如下步骤获得 a、若待简化特征均由平面组成，获取此待简化特征的最小边长hi，计算网格直径改变量Htemp＝hmax-min(hmax，hi)，若Htemp＞0，则采用此种特征简化方案，对应的网格直径改变量否则转下一步骤； b、若待简化特征中存在曲面，获取当前待简化特征i包含的所有高曲率曲面的平均高斯曲率θavgi，计算网格直径改变量，若Htemp＞0，则采用此种特征简化方案，对应的网格直径改变量 所述步骤e中对模型态网格单元进行数量的估算是按如下步骤进行 a、根据用户提供的最大网格直径hmax，最小网格直径hmin，取网格单元边长h＝(hmax+hmin)/2； b、结合网格单元类型和网格单元边长，计算出一个网格单元所占用的实际体积ve，体积计算公式因单元而异； c、计算出模型态的实际体积v，得到模型网格划分的单元数量Ne＝v/ve； 所述步骤f中采用遗传算法进行模型态选择按如下步骤进行 a、使用二进制编码表示模型态，编码长度是特征数量，每个特征对应于个体中的一个二进制位；0代表某特征不存在模型态中，1代表某特征存在于模型态中； b、用户输入可容忍的误差水平级别m，m∈[1，n]，形成最大允许误差，其中n是特征数量； c、随机生成个体数量为的种群，其中NF是特征数量；设置交叉概率Pc∈
，变异概率Pm∈
，每个个体的适应度计算公式为，其中EM＝∑Ei，i是模型态编码中为0的特征序号；α是权重因子，这里α选取为0.3； d、迭代次数为
，每迭代一次，估计出这代中所有个体编码对应模型态的网格单元数量，以此为基础计算适应度函数；若连续三代没有出现适应度更大的个体或者到达最大迭代次数，则最后一代的种群中适应度最大的个体即所获得的简化模型编码。
与已有技术相比，本发明的有益效果体现在 1、误差因素分析的客观性 本发明方法针对使用某个模型态代替原始模型进行计算产生误差的原因进行分析，得出了与误差相关的要素体积、特征位置参数和特征简化后网格直径改变量；摒弃了从经验出发，通过实验数据来寻找误差相关要素方法的随机性。
2、计算消耗时间估计的精确性 已有的方法中，只有在文献7中初步总结了衡量计算耗时的要素，但其提出将模型中的平均特征体积作为计算耗时的衡量参数，精确性不足。本发明方法将可以精确反映计算耗时的要素，即模型划分网格后的单元数量Ne作为模型态评价要素中的时间要素，能准确地反映模型态对应的计算时间，有益于最优模型态编码的选取。
3、遗传算法适应度函数构造的合理性 适应度函数的构造充分地考虑了模型态计算时间和计算精度之间的平衡关系，为了平衡两个要素的影响因子，对适应度函数中误差水平EM进行了数量级修正。
图1为本发明方法中模型态的选择变迁数学模型示意图。
图2为网格顶点在曲面上对应点外法矢量示意图。
图3为模型表面曲率和网格之间关系示意图。
图4为模型态编码示意图。
图5为模型态选取示意图。
图6为槽、孔和圆角的特征组成面示意图。其中，图6(a1)为含有槽特征的模型；图6(a2)为图6(a1)所示模型对应的面边图；图6(b1)为含有孔特征的模型；图6(b2)为图6(b1)所示模型对应的面边图；图6(c1)为含有圆角特征的模型；图6(c2)图6(c1)所示模型对应的面边图。
图7(a)是含有一个盲孔、一个通孔和一个槽特征的模型，图7(b)为图7(a)对应的特征子图。
图8(a)、图8(b)、图8(c)、图8(d)和图8(e)分别为方通孔、方盲孔、圆通孔、圆盲孔和槽的特征子图分类图。
图9为模型特征树示意图。
图10为本发明方法中待简化的实际模型示意图。
图11为以本发明方法所获得的简化模型示意图。
图中标号1基座、2内圆柱A、3内圆柱B、4内圆柱C、5内六棱柱A、6内八棱柱A、7内八棱柱B、8内六棱柱B、9槽A、10槽B、11槽C、12槽D、13槽E、14槽F、15圆角A、16圆角B、17圆角C、18圆角D。
以下通过具体实施方式
，并结合附图对本发明作进一步描述
具体实施例方式 本实施例给出模型误差水平构建方法按如下步骤进行 1、由用户提供待简化的模型，所述待简化的模型具有在进行有限元分析时的载荷条件，以及网格剖分条件，所述网格剖分条件为最大网格直径hmax和最小网格直径hmin； 2、对其进行特征识别，将识别出的特征分类存入特征链表；遍历特征链表，计算获得每个特征的体积Vi，每个特征的位置参数Pi； 其中，链表是一种常规的计算机数据结构； 在这一步骤中，待简化特征的体积Vi按常规的体积公式进行计算，特征的位置参数Pi通过如下方式获得 a、给定每个待简化特征简化后对有限元计算结果影响的权重，设为q1，q2，…qn，其中qi∈
，(i＝0，1，…n)，qi＝0表征第i个待简化特征存在与否对有限元计算结果没有影响；qi＝1表征第i个待简化特征存在与否对有限元计算结果影响最大，qi越靠近1表征特征对有限元计算结果的影响越大； b令i＝1； c求取第i个待简化特征包围盒形心，若载荷为体载荷，则获得载荷体的形心；计算第i个待简化特征形心到载荷面的距离di，若为体载荷，则计算第i个待简化特征形心到载荷体形心的距离di；若有多个载荷面则计算第i个待简化特征形心到每个载荷的距离，再取其平均值； d计算第i个待简化特征的位置参数Pi＝qi/di； e令i＝i+1，重复步骤c～d，直至模型中所有待简化特征的位置参数Pi计算完毕； 3、计算获得所有待简化特征所有特征组成面的最大高斯曲率θmax； 4、根据网格直径改变量计算方法，计算特征被简化后网格直径的改变量Hi； 依据最大曲率θmax、最大网格直径hmax和最小网格直径hmin按如下步骤获得网格直径改变量； a、令i＝1，计算ε＝θmaxhmin；若特征组成面均为平面，则获取此待简化特征的最小边长hi，计算网格直径改变量Htemp＝hmax-min(hmax，hi)，若Htemp＞0，对应的网格直径改变量Hi＝Htemp，否则转下一步骤； b、若特征组成面中存在曲面，则获取当前待简化特征i包含的所有曲面的平均高斯曲率θavgi，计算网格直径改变量，若Htemp＞0，则对应的网格直径改变量Hi＝Htemp，否则转下一步骤； c、令i＝i+1，重复步骤a、b，直至模型中所有待简化特征与其简化方式相应的网格直径改变量Hi，i＝1，2，…，n，计算完毕； 5、由误差水平估计公式Ei＝ViPiHi估计出该特征被抑制后的敏感值，并以所有被抑制特征的敏感值之和表征模型态的误差水平； 6、构造遗传算法的各参数并进行模型态的选择； 遗传算法的构建步骤如下 a、使用常规的二进制编码表示模型态，编码长度是特征数量，每个特征对应于个体中的一个二进制位；0代表某特征不存在模型态中，1代表某特征存在于模型态中； b、用户输入可容忍的误差水平级别m，m∈[1，n]，形成最大允许误差，其中n是特征数量； c、随机生成个体数量为的种群，其中NF是特征数量；设置交叉概率Pc∈
，变异概率Pm∈
，每个个体的适应度计算公式为，其中EM＝∑Ei，i是模型态编码中为0的特征序号；α是权重因子，这里α选取为0.3； d、迭代次数为
，每迭代一次，估计出这代中所有个体编码对应模型态的网格单元数量，以此为基础计算适应度函数；若连续三代没有出现更优秀的个体或者到达最大迭代次数，则遗传算法结束，获得最优模型态编码； 在这一步骤中，种群中某一代所有个体编码对应模型态的网格单元数量通过如下方式获得 (1)、根据用户输入的最大网格直径hmax，最小网格直径hmin，取单元边长h＝(hmax+hmin)/2； (2)、结合单元类型element type和单元边长，根据几何知识计算出一个单元占用的实际体积ve，体积计算公式因单元而异； (3)、按常规方式计算出模型态的实际体积v，得到模型网格划分的单元数量Ne＝v/ve； 7、根据最优模型态编码，遍历特征链表，使用特征抑制技术取得编码对应的几何模型。
本发明有限元领域中基于特征抑制的模型态生成方法所依据的原理包括 本发明所依据的多态模型理论是从系统论的角度出发对科学计算中的模型进行分析研究，将模型作为一个系统，模型中所包含的特征作为构成系统的要素。要素本身的可变性即每个要素都具有多态特性，同时要素之间的构成方式的变化会导致系统状态发生变化，因此将模型视为一个系统时，其本身呈现出多态特性。通过对多态模型所包含的特征要素及其组成方式的研究，为实际模型建立相应的系统，以特征作为系统的组成要素，研究各种特征的敏感性，建立特征敏感性的量化评价标准，寻找模型的态及在其变化过程中所蕴含的规律性，建立多态理论体系来表述复杂模型由于模型精度、计算时间、计算精度的变化而产生的多种状态。多态理论体系的核心思想是态的变迁，即原模型态到达目标模型态的动态过程如图1所示，将系统模型态记为M，可描述为M{f1，f2，…，fn}，其中fi(i＝1，2，…n)为其所包含的特征，为方便起见记为。设为原模型态，为变迁模型态。记各态计算代价(例如计算时间)为Cc(Mj)，(j＝0，1，2，…，n)。若模型态Mi到模型态Mj(i，j＝0，1，…n，i≠j)存在转化路径(路径可能不唯一)，设转化路径集合为Pij，模型态Mi沿路径p(p∈Pij)转化为模型态Mj的转化代价(例如计算结果的误差)记为
，则模型态Mi到模型态Mj的最小转化代价记为
。根据给定的转化代价阈值cε，可以取定满足转化代价阈值的模型态集合Tminct＝{Mj|CTmin(M0→Mj)＜cε}，再根据Tminct取定满足最小计算代价的模型态Mf，即Mf满足条件 本发明方法中使用模型态代替原始模型进行计算而产生误差的误差水平估算依据的原理 当有限单元类型和插值方式确定时，计算模型中特征对简化的敏感值只需考虑三个因素ΔT、ΔH和ΔU 其中，ΔT为特征简化对网格剖分产生影响的区域测度反映的要素。
ΔH为特征简化后网格直径改变量反映的要素。
ΔU为特征简化对网格剖分产生影响的区域内场函数的半范数反映的要素。
在计算模型中每个特征对简化的敏感值时需要对三个要素ΔT、ΔH和ΔU进行计算；但是要素ΔT和ΔH的精确值需要在网格剖分后才可能知道，要素ΔU的精确值需要在有限元计算后才可能知道，若要在有限元计算之前获得它们的精确值是困难的，并且对它们的精确计算所需时间可能会大于模型简化后节省的计算时间，从代价上考虑不宜进行精确计算，因此需要寻找一种方法对它们进行估值。一个明显的事实是特征简化对网格剖分产生影响的区域测度反映的要素ΔT与模型特征体积成正比，因此可用特征体积V来代替要素ΔT。在本发明方法中用特征体积V来代替要素ΔT，用特征位置参数P来代替要素ΔU，因此误差估算公式计算方法为Ei＝ViPiHi，其中Ei，Vi，Pi，Hi分别表示特征在第i中简化方式下产生的简化误差、特征体积、位置参数和网格直径改变量。
本发明方法中网格直径改变量计算方法所依据的原理 已知的情况是，CAD系统建立的模型通常会包含一些复杂特征，如细小特征和高曲率特征等，而有限元分析进行网格剖分时会在这些复杂特征处加大网格剖分密度，从而这些复杂特征处的网格直径比较小。删除或用简单特征替换复杂特征的目的是为了减小这些特征处的网格密度，也即加大这些特征处的网格直径。在误差估算公式计算方法原理中以指出特征简化前后网格直径的改变量是影响特征对简化敏感值的因素之一，本发明方法使用替换和删除两种方式对特征进行简化，不同的简化方式对网格直径改变量计算方法分别是 1、不含高曲率组成面的特征删除后网格直径改变量的计算方法 如果特征中不含高曲率特征组成面，则有限元分析在进行网格剖分时，特征上的网格直径将取决于用户能够容忍的最大网格直径和符合此特征几何属性网格直径的最小值，而符合此特征几何属性的网格直径应为此特征组成面中的最小边长。当特征删除后，模型中此特征不存在，则有限元分析进行网格剖分时应按用户所指定的最大网格直径进行剖分，即用户输入的hmax，因此不含高曲率组成面的特征删除后网格直径改变量应为hmax-min(hmax，h)，其中h为此特征的特征组成面中的最小边长。
2、含高曲率组成面的特征删除后网格直径改变量的计算方法 如果特征中含高曲率特征组成面，则应首先估计在高曲率特征组成面处的网格直径。首先阐述如何估计高曲率特征处的网格直径。
在有限元网格生成过程中，网格的密度控制包括两种一种是根据分析对象的几何特征和物理特性进行网格疏密的先验控制，另一种是根据当前网格的计算结果，在计算数据变化较大的区域加大网格剖分密度并在计算数据变化平缓的地方减小网格剖分密度的后验控制。本方法针对在有限元分析前对模型的预处理，因此只考虑网格密度的先验控制。根据分析对象的几何特征进行网格控制，主要是指根据网格对模型表面的拟合程度决定网格直径的大小。网格对模型表面的拟合程度可用网格顶点在曲面上对应点外法矢量的夹角来近似。如图2所示，模型的某一块表面使用一个三角形单元来拟合，三角形的顶点对应的外法矢量指尖的夹角决定了单元对模型表面的拟合程度。对于给定对模型表面拟合程度的容许误差η，则要求每一对外法矢量都满足(1-NiNj)≤η，i≠j且i，j∈(1，2，3)。
当用网格剖分算法对模型表面某一局部区域ΔS进行剖分时，考虑的是其平均法曲率，设为ρ，曲率球半径为r，可在以r为半径的球面上考虑ΔS上的网格直径。设△ABC为一三角网格，因球面上处处法曲率相同，所以△ABC为等边三角形，不失一般性，以边AB讨论△ABC的边长。
如图3所示，设在坐标系XOY中，A，B坐标分别为(x1，y1，z1)，(x2，y2，z2)，当模型表面拟合程度容许误差ε一定时，有如下等式 (x1-x2)2+(y1-y2)2+(z1-z2)2＝h2 ρ＝1/r 从而有 有限元分析剖分网格时，当一旦检测到即停止细化网格，因此可以认为。设模型中所有待简化特征的特征组成面的最大曲率为ρmax，则在此特征上应有最小网格直径hmin，因此，当待简化特征的特征组成面的平均曲率为ρavg时，此待简化特征上的网格直径应估计为 同样，当此类特征删除后，模型中此特征不存在，则有限元软件进行网格剖分时应按用户所指定的最大网格直径进行剖分，即用户输入的hmax，因此不含高曲率组成面的特征删除后网格直径改变量应为 3、含高曲率组成面的特征替换后网格直径改变量的计算方法 对于此类型待简化特征，特征上的网格直径估计方法与上一类型中估计方法相同，不同的是此类型的特征将有多种简化方案。对于平均曲率为ρavg的曲面，其表面的网格直径为，用正k变形替换后网格直径应为，因此网格直径改变量为 本发明方法中特征的位置参数计算方法依据的原理 如误差估算公式计算方法所依据的原理中所述，在场函数变化剧烈的区域，简化误差的三个要素之一ΔU也会较大，若对这些区域内特征进行简化会导致较大的误差，因此模型中可被简化的特征必须局限于场函数变化平缓的区域，也即特征在模型中所属区域的场函数变化剧烈程度是特征简化后对有限元计算结果影响程度的重要因素。在实际的有限元计算之前，可通过下面两种方法获得模型中特征所属区域场函数变化的剧烈程度。
1、依据经验 由用户给定模型特征所属区域场函数变化的剧烈程度，这是一个非常重要的参考值。在有限元分析中，特征的物理属性，即外加载荷及边界条件赋予特征的属性——特征所处区域场函数变化的剧烈程度，与特征几何形状有相同重要的影响，而这些物理属性在有限元计算之前无法获得精确值，可以根据以往分析经验给出特征物理属性值的排序，对于分析特征简化后对有限元计算结果影响有重要的参考作用。
2、距离的作用 根据圣维南原理，在平衡力系中，离载荷越远场函数变化越平缓，因此本发明方法将特征到载荷的距离作为反映场函数变化的剧烈程度的因素之一，将距离和专家指定的权重乘积作为ΔU的估计值。
根据常规理论，特征对简化的敏感值与距离成反比、与用户指定的权值成正比，所以取用户指定的权值和特征到载荷的距离的比值作为特征的位置参数，以此来反映特征所属区域场函数u变化的剧烈程度。
选取模型划分网格后的单元数量作为时间要素所依据的原理 在有限元领域中，各问题的求解过程都是相同的，步骤如下 a、对问题的求解区域进行剖分，剖分的结果是单元，求解区域离散为单元的集合。
b、构造试探函数空间。
c、计算单元刚度矩阵和单元载荷向量。
d、计算总体刚度矩阵和总载荷向量。
e、处理约束条件并求解。
以一维有限元问题为例，将求解区域[a，b]中，首先把[a，b]剖分为n个部分，设节点x0，x1，…，xn满足a＝x0＜x1＜…＜xn＝b，这样得到的子区间ei＝[xi-1，xi](i＝1，2，…，n)，即单元。单元ei的长度hi＝xi-xi-1。考虑构造比较简单的基函数vh，它满足在[a，b]上连续，而且在每个ei上是线性函数，所有满足这些条件的函数构成试探函数空间。再计算单元刚度矩阵和单元载荷向量将单元刚度矩阵kei扩展为(n+1)阶矩阵即把kei的四个元素对应放在一个(n+1)阶矩阵
的第(i-1)行和i行的2阶对角块中，而其他元素均为0。同理，可以将单元载荷向量Fei扩展成为(n+1)阶的向量，那么总体刚度矩阵，总载荷向量，考虑约束条件的前提下解代数方程组vT(Ku-F)＝0，即可得到有限元解。若认为这个近似解不够精确，可以使得剖分更细，即节点取的更多。从整个求解过程可以看出，有限元求解就是解一个代数方程组，显然，这个代数方程组内方程的多少决定了求解所需要的时间，而方程的数量就是总体刚度矩阵K的阶数，即问题域剖分得到的节点数量。虽然节点数量和计算所需要的时间不一定成线性关系，但是它是最能够反映所需计算时间的参数。而模型划分网格后的单元数量与节点数量成线性关系，所以通过算法来估计出单元数量从而反映出模型态计算消耗时间。
本发明方法中使用遗传算法选择模型态依据的原理 本发明方法中，以特征作为区分模型态的要素，模型态的表示方式为M{f1，f2，…，fn}。由排列组合的原理可知，含有n个特征的模型有2n个可能的模型态。如果模型中含有的特征数量太多，选择过程将耗费大量时间。
选择产生模型态的过程实质上是以成本函数为向导，在可能的解空间搜索最佳或几乎最佳的解的问题。因此，多态模型理论的核心问题就是如何在这么大的解空间里快速找到最优解的问题。常用的搜索算法有贪婪算法、模拟退火算法及遗传算法等。前两种算法虽然收敛速度快，但存在易陷入局部最优化的问题；而遗传算法具有搜索空间大的优点，但收敛速度较慢。
遗传算法是一种模拟生物在自然环境中的遗传和进化过程而形成的一种自适应全局优化概率搜索算法，它具有简单通用、鲁棒性强、适于并行处理的特点。它将实际问题中的可能解模拟成个体的生存环境，将目标函数模拟成个体的生存能力，将可能解的编码模拟为染色体。这样，从任一个初始种群出发，经过选择、交叉、变异三种运算产生新一代种群，经过多次迭代后，使其收敛于全局最优解或次最优解。
遗传算法的流程是设定目标和适应度函数，经过一定次数的迭代之后，找到适应度最大的个体。其中，算法中染色体结构的二进制编码方式所具有的特点与特征抑制的模型态生成方式具有很好的一致性与协调性。编码方式如图4所示，编码的长度就是特征的数量，每个编码对应一个特征的存在状态，0表示某特征不存在于这个模型态中，1表示相反的意义。通过对这个编码构造适应度函数来评价模型态的质量。本发明针对原始模型中的独立特征，结合有限元分析的特点，以热分析领域为切入点，提出一个同时考虑计算时间和计算精度的模型态选择产生标准。采用遗传算法作为搜索算法，通过搜索最佳或几乎最佳的模型态来建立合适的计算模型。如图5所示，根据个体长度随机生成一定数量的初始种群，进行选择、交叉和变异操作进行迭代，最后得到最优模型态的编码。
模型态的评价要素是模型网格划分后的单元数量Ne和误差水平Ei。NE越大，计算耗时越多；Ei越大，精度越低。因此模型态对应的NE和Ei越小，它的质量越高。大多数情况下，这2个要素之间有相互制约的关系，即对同一个模型来说，NE越小说明Ei可能较高，因此这是一个双目标优化问题。
传统的多目标优化方法是将各个子目标聚合成一个带正系数的单目标函数，系数由决策者(Decision Maker，DM)决定，或者由优化方法自适应调整。其中，加权法是一种常见的古典方法，是通过对目标函数的线性组合将MOP问题转换成SOP问题。
y＝f(x)＝ω1f1(x)+ω2f2(x)+…+ωkfk(x) ωi称为权重，通常权重可以正则化后使得∑ωi＝1，求解上述不同权重的优化问题则能够输出一组解。
基于古典多目标问题的加权解法，基于NE和EM构造适应度函数。在定义中，α的值对选取的模型态质量有重要影响，采用基于实验统计的方式求得α的最优值。穷举法可以把所有可能的模型态都列举出来，从中选取最优解，即全局最优模型态，因此通过比较不同α值选出的模型态质量，根据统计规律选出α的最佳值。
本发明特征识别及抑制算法依据的原理 本发明方法中，采用属性邻接图的方法进行特征识别。特征是由面组成的，即特征是特定的面集。
如图6所示，在一个方体模型中，分别含有槽、孔和过渡特征。其中槽特征由三个平面组成，孔特征由一个圆柱面组成而过渡特征由1/4圆柱面组成。在计算机图形学中的面-边图中，以面作为节点，面与面之间如果邻接的话，则代表面的节点之间有线段相连。在分别含有槽、孔和过渡特征的方体模型中，分别有三个节点、一个节点和一个节点代表相应特征。在图7中的方体模型含有1个开槽特征、1个圆盲孔特征和1个圆通孔特征，假设对模型中的所有面进行编号，槽特征对应的面编号为2，3，4；圆盲孔特征对应的面编号是11，12；圆通孔特征对应的面编号是13，那么在其对应的面-边图中，编号为2，3，4的三个节点包括与其相连接的所有边构成了该特征对应的子图；同理，编号为11，12的节点连同与其相连的所有边构成了圆盲孔特征的子图；编号为13的节点连同与其相连的所有边构成了圆通孔特征的子图，如图8所示。不同的特征对应的子图具有不同的特点，在图9中，根据槽、孔、过渡特征子图的不同特点，可以将模型中的特征分为孔、槽、过渡特征加以识别并存入链表，以便于特征抑制。
在文献6中，对特征识别及抑制的技术做了详细说明。步骤如下 1、判断每条边的凹凸性，生成模型的AAG图。
2、遍历AAG，找出所有的凸面，构成面集fvex，对fvex进行分析找出单一型特征，具体分类 (1)圆通孔特征凸面，类型为柱面，并且只有两条边。
(2)圆角特征凸面，类型包括环面、柱面、球面、样条曲面。
3、从AAG图中删除fvex，得到子图，其中所有的结点构成面集fcav。找出这一子图所有的连通分量，每个连通分量即为一个特征子图。对每个特征子图按照上述特征子图的分类特征进行识别，判断具体的特征类型。
4、圆盲孔特征两个面，有一个是叶子节点。两条相邻的边，默认叶子节点对应是平面，另一个面是圆柱面。
5、方盲孔特征存在一个面，所有邻边都为凹边。图中含有一个回路，回路中的节点对应面均属于此特征。
6、方通孔特征图中含有一个回路，回路中的节点对应面均属于此特征。每个面都有四条相邻边，且满足两条为凸边，两条为凹边。
7、槽三个面，两条凹边。
8、根据识别出的特征类型，生成对应的特征对象，记录对应的特征组成面。
9、生成模型对应的特征树，用于后续处理。
实施例 1、由用户提供待简化模型，待简化的模型具有在进行有限元分析时的载荷条件，以及网格剖分条件，如图10所示，待简化模型有一个中心坐标为(0，0，0)、边长为10个单位的正方体基座1，减去一个底面圆心为(2，-2，0)、半径为0.8个单位、高为10个单位的内圆柱A2，再减去一个底面圆心为(-1，-1，3.5)、半径为0.8个单位、高为3个单位的内圆柱B3，再减去一个底面圆心为(3.5，1，0)、半径为0.5个单位、高为10个单位的内圆柱C4，再减去一个中心坐标为(2，0，0)、边长为0.8个单位、高为10个单位的内六棱柱A5，再减去一个中心坐标为(2.5，2，-2.5)、边长为0.5个单位、高为5个单位的内八棱柱A6，再减去一个中心坐标为(1，2，0)、边长为0.8个单位、高为10个单位的内八棱柱B7，再减去一个中心坐标为(-1，2，3.5)、高为3个单位的内六棱柱B8，再减去对角顶点为(2，4，-5)和(3，5，5)的槽A9，再减去对角顶点为(-1，4，-5)和(0.5，5，5)的槽B10，再减去对角顶点为(-3，3，-5)和(-2，5，5)的槽C 11，再减去对角顶点为(2，-4，-5)和(3，-5，5)的槽D12，再减去对角顶点为(-1，-4，-5)和(0.5，-5，5)的槽E13，再减去对角顶点为(-3，-3，-5)和(-2，-5，5)的槽F14，对立方体的右侧面的四条边进行半径为0.8个单位的圆角化操作，得到圆角A15，圆角B16，圆角C17，圆角D 18，载荷面为立方体的左侧面，网格剖分条件为最大网格直径hmax＝2和最小网格直径hmin＝1； 2、建立特征链表 模型的载荷情况，有限元单元属性与网格剖分条件如表1所示。用户使用商业建模软件建立格式为.SAT的3D几何模型，将模型文件导入并进行特征识别，将识别出的特征分类存入特征链表；再次遍历特征链表，针对每一个特征，根据Di计算出该特征的位置参数Pi；通过如下方式获得位置参数Pi a、给定每个待简化特征简化后对有限元计算结果影响的权重为q1＝q2＝…＝q8＝0.5； b令i=1； c设第i个待简化特征包围盒形心，形心坐标为(xi，yi，zi)，计算第i个待简化特征形心到载荷面的距离di，di＝xi； d计算第i个待简化特征的位置参数Pi＝qi/di； e令i＝i+1，重复步骤c～d，直至模型中所有待简化特征的位置参数Pi计算完毕； 获得所选择待简化特征所有特征组成面的最大高斯曲率θmax＝0.25； 2、依据每个所选择的待简化特征的体积Vi、位置参数Pi、每个所选择的待简化特征在其相应简化方案下的网格直径改变量按计算公式Ei＝ViPiHi获得此特征在其相应简化方案下的敏感值Ei，如表1所示； 本实施例中网格直径改变量按如下步骤获得 a、计算ε＝θmaxhmin＝0.25，令i＝1； b、令i＝1，计算ε＝θmaxhmin；若特征组成面均为平面，则获取此待简化特征的最小边长hi，计算网格直径改变量Htemp＝hmax-min(hmax，hi)，若Htemp＞0，对应的网格直径改变量Hi＝Htemp，否则转下一步骤； c、若特征组成面中存在曲面，则获取当前待简化特征i包含的所有曲面的平均高斯曲率θavgi，计算网格直径改变量，若Htemp＞0，则对应的网格直径改变量Hi＝Htemp，否则转下一步骤； d、令i＝i+1，重复步骤b、c，直至模型中所有待简化特征与其简化方式相应的网格直径改变量Hi，i＝1，2，…，n，计算完毕，得到表1； 本实施例中敏感值为对应特征的体积、位置参数和相应简化方式之下的网格直径改变量之积，计为Ei＝ViPiHi，i＝1，2，…，n，得到表1； 表1特征属性列表 4、本实施例中有限单元选取用于热分析的Solid Tet 10node 87单元，即10节点单元，网格剖分采用默认的划分密度。由于平均边长是1.5个单位，因此单元的平均体积是3.375，模型态的体积为v，则其对应的单元数量NE＝v/3.375； 5、使用二进制编码表示模型态，编码长度是特征数量，在热分析领域中，用户根据经验可知圆角特征对网格划分的负面影响很大而对计算结果的影响可以忽略，将该模型的4个圆角特征先行抑制，于是特征编码为13位，对应于可能存在的213个模型态，模型中特征的顺序与表中的特征序号一致。其中初始模型态编码为1111111111111。
6、用户输入误差级别为1，计算出可容忍的量化误差Emax＝11.55。
7、随机生成个体数量为N＝100的种群；交叉概率Pc＝0.8，变异概率Pm＝0.02，取α＝0.3，那么每个个体的适应度计算公式为 其中 8、迭代次数为
，每迭代一次，使用网格剖分算法估计出这代中所有个体编码对应模型态的刚度矩阵阶数，再计算适应度函数。最终模型态的编码为0011010011011，经过特征抑制，几何形状如图13所示，剩下正方体基座1、圆盲孔特征3、正六面体通孔5、正六面体盲孔8、槽10、11、14。
权利要求
1、有限元建模中基于特征抑制的模型态生成方法，其特征是按如下步骤进行
a、由用户提供待简化模型，所述待简化模型具有在进行有限元分析时的载荷条件，以及网格剖分条件，所述网格剖分条件为最大网格直径hmax和最小网格直径hmin；
b、识别所述待简化模型中的可简化特征，并存入特征链表；
c、在所述特征链表中选择待简化特征，确定待简化特征数量n；按任意顺序对所选择的n个待简化特征进行编号，计算获得所选择的每个待简化特征的体积Vi和位置参数Pi；
d、依据每个所选择的待简化特征的体积Vi、位置参数Pi，以及每个所选择的待简化特征在其被抑制后的网格直径改变量，获得对应的待简化特征在其被抑制后的敏感值Ei＝ViPiHi；以所有被抑制特征的敏感值之和表征模型态的误差水平；
e、通过对模型态网格单元进行数量的估算，估计出模型态进行网格划分后的单元数量Ne，所述单元数量Ne与模型态刚度矩阵阶数线性相关，以所述单元数量Ne反映模型态所需的有限元计算时间；
f、根据模型的误差水平和单元数量Ne，采用遗传算法进行模型态的选择，获得简化模型编码；
g、应用特征抑制法生成对应于所述简化模型编码的简化模型。
2、根据权利要求1所述的有限元建模中基于特征抑制的模型态生成方法，其特征是所述步骤c中简化特征位置参数Pi通过如下方式获得
a、给定每个待简化特征简化后对有限元计算结果影响的权重，设为q1，q2，…qn，其中qi∈
，(i＝0，1，…n)，qi＝0表征第i个待简化特征存在与否对有限元计算结果没有影响；qi＝1表征第i个待简化特征存在与否对有限元计算结果影响最大，qi越靠近1表征特征对有限元计算结果的影响越大；
b、令i＝1；
c、设第i个待简化特征到载荷面的距离di；
d、计算第i个待简化特征的位置参数Pi＝qi/di；
e、令i＝i+1，重复步骤c～d，直至模型中所有待简化特征的位置参数Pi计算完毕。
3、根据权利要求1所述的有限元建模中基于特征抑制的模型态生成方法，其特征是所述步骤d中网格直径改变量按如下步骤获得
a、若待简化特征均由平面组成，获取此待简化特征的最小边长hi，计算网格直径改变量Htemp＝hmax-min(hmax，hi)，若Htemp＞0，则采用此种特征简化方案，对应的网格直径改变量否则转下一步骤；
b、若待简化特征中存在曲面，获取当前待简化特征i包含的所有高曲率曲面的平均高斯曲率θavgi，计算网格直径改变量，若Htemp＞0，则采用此种特征简化方案，对应的网格直径改变量
4、根据权利要求1所述的有限元建模中基于特征抑制的模型态生成方法，其特征是所述步骤e中对模型态网格单元进行数量的估算是按如下步骤进行
a、根据用户提供的最大网格直径hmax，最小网格直径hmin，取网格单元边长h＝(hmax+hmin)/2；
b、结合网格单元类型和网格单元边长，计算出一个网格单元所占用的实际体积ve，体积计算公式因单元而异；
c、计算出模型态的实际体积v，得到模型网格划分的单元数量Ne＝v/ve；
5、根据权利要求1所述的有限元建模中基于特征抑制的模型态生成方法，其特征是所述步骤f中采用遗传算法进行模型态选择按如下步骤进行
a、使用二进制编码表示模型态，编码长度是特征数量，每个特征对应于个体中的一个二进制位；0代表某特征不存在模型态中，1代表某特征存在于模型态中；
b、用户输入可容忍的误差水平级别m，m∈[1，n]，形成最大允许误差，其中n是特征数量；
c、随机生成个体数量为的种群，其中NF是特征数量；设置交叉概率Pc∈
，变异概率Pm∈
，每个个体的适应度计算公式为，其中，EM＝∑Ei，i是模型态编码中为0的特征序号；α是权重因子，这里α选取为0.3；
d、迭代次数为
，每迭代一次，估计出这代中所有个体编码对应模型态的网格单元数量，以此为基础计算适应度函数；若连续三代没有出现适应度更大的个体或者到达最大迭代次数，则最后一代的种群中适应度最大的个体即所获得的简化模型编码。
全文摘要
有限元建模中基于特征抑制的模型态生成方法，其特征是由用户提供待简化模型；识别待简化模型中可简化特征；确定待简化特征数量；获得待简化特征在其被抑制后的敏感值；以所有被抑制特征敏感值之和表征模型态误差水平；估计模型态进行网格划分后的单元数量；根据模型的误差水平和单元数量Ne进行模型态的选择，获得简化模型编码；再应用特征抑制法生成对应于简化模型编码的简化模型。本发明方法实现模型中各特征的变化对有限元计算结果敏感性的客观评价，建立模型中各要素的取舍与变化对计算结果的敏感性评价体系，并给出了客观精确的模型态消耗时间的估计方法，在此基础上生成模型态。
文档编号G06F17/50GK101241520SQ20081002056
公开日2008年8月13日 申请日期2008年2月1日 优先权日2008年2月1日
发明者刘晓平, 灿 金, 慧 石, 强 路, 郑利平, 李书杰, 敏 吴, 罗月童, 徐本柱 申请人:合肥工业大学