汽车覆盖件冲压模具结构优化与分析

《汽车覆盖件冲压模具结构优化与分析》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。

西安科技大学硕士学位论文汽车覆盖件冲压模具的结构优化与分析姓名：徐自立申请学位级别：硕士专业：固体力学指导教师：郭志勇;娄茂昆@ 论文题目：汽车覆盖件冲压模具的结构优化与分析专业：固体力学硕士生：徐自立（签名）指导老师：郭志勇（签名）娄茂昆（签名）摘要近年来我国汽车工业飞速发展，众多汽车企业都在为缩短新车型研发周期、降低生产成本而努力。冲压模具作为生产汽车覆盖件的重要工艺装备，其设计质量和成本对汽车的生产研发有着重要影响。随着能源和钢铁价格不断上涨，控制模具成本显得尤为重要。另外，高强度钢板等新材料、伺服压力机等新设备的应用，对传统的冲压模具提出了新的问题和挑战。要解决这些问题就必须加强对模具结构优化设计和模具结构分析的研究。本文在总结原有模具结构分析方法的基础上，从覆盖件模具入手，提出了基于板料成形模拟的冲压模具结构优化方法和分析方法。采用冲压成形模拟和结构优化分析解耦的思路，运用拓扑优化技术，以某行李箱盖外板和后备箱盖外板为例，设计出新的冲压模具结构，并对其在冲压过程中的受力和变形规律进行了分析。结构优化方面，首先运用专业的成形模拟软件PAMSTAMP2G对板料冲压成形过程进行模拟，得到冲压过程中板料对模具的作用力。这些荷载通过笔者编写的FORTRAN程序Pam2OptiStruct.exe转换后直接施加到模具结构优化模型上。然后利用OptiStruct的拓扑优化技术，对模具结构进行优化，最后借助CAD系统CATIA在优化结果的基础上设计出新的模具结构。经结构分析验证，新结构相对传统结构减重15%左右，而柔度降低5%左右，同时保证变形和等效应力在许可的范围内。这说明新设计的模具结构优于传统模具结构。结构分析方面，板料冲压成形过程看作一个准静态过程。本文将连续的冲压过程用n个时间点的静态结构分析序列来近似，通过分析一系列的模具受力和变形，得到整个冲压过程中模具的受力和变形规律。关键词：覆盖件；冲压模具；结构优化；模具结构分析；数值模拟研究类型：应用研究 Subject：StructuralOptimizationandAnalysisofStampingDieforAutomobilePanelSpecialty：SolidMechanicsNam：XuZili(Signature)Instructor：GuoZhiyong(Signature)LouMaurice(Signature)ABSTRACTInrecentlyyears,astherapiddevelopmentofautomobileindustry,manyautomobileenterprisesaremakinggreateffortstoshortendevelopmentcycleanddiminishproductioncost.Stampingdie,asanimportantprocessequipmentofautomobilepanelproduction,itscostandqualityhasasignificantimpactonR&D(ResearchandDevelopment).Withthepricesrisingofenergyandiron,savingcostofdiehasbecomemoreandmoreimportant.Ontheotherhand,traditionalautopaneldiearefacingnewchallengesbecauseoftheapplicationofnewmaterialsandnewequipmentssuchashigh-strengthsteelandservopress.Alloftheseneedfurtherresearchonitsstructuraloptimizationandanalysis.Thispapersummarizedtheexistingmethodsofdiestructureanalysis.Startingwithautomobilepaneldie,thestructuraloptimizationandanalysismethodsforstampingdiewereproposedbasedonthesheetmetalformingsimulation.Withthethoughofdecouplingthestampingsimulationandstructureanalysis,usingtopologyoptimizationtechniques,thenewstamping-diestructuresforahoodandadecklidouterpanelweredesignedasanexample.Andthenthelawofstressanddeformationduringstampingprocesswasanalyzed.Instructuraloptimization,PAMSTAMP2Gwasusedtogettheforcesthatblankactedondiesurfacebysimulatingthestampingprocess.TheseforcesweretransformedanddirectlyappliedonfiniteelementmodelbyaFORTRANprogram:Pam2OptiStruct.exe.ThenanoptimizationprocessofdiestructurewasmadeusingtopologyoptimizationtechniquesofOptiStruct.Basedontheoptimizedresult,anew-stylestructurewasdesignedusingCADsystemCATIA.Thenewstructurehasabatterperformancethantraditional,thisisconfirmedbycomparingstructureresponses.Ithasaweightreductioneffectabout15%andastiffnessincreasingeffectabout5%.Atthesametime,thedeformationandequivalentstressareinthepermittedrange. Astherespectofstructureanalysis,thesheetmetalformingprocessisquasistatic.Thiscontinuousstampingprocesswasapproximatedbyasequenceofstaticstructureanalysisatn(nisaninteger)timepoints.Bythisway,thelawofdeformationandstressduringthestampingprocesswasdescribed.KeyWords:AutomobilePanelStampingDieStructuralOptimizationDieStructureAnalysisNumericalSimulationThesis:ApplicationResearch 1绪论1绪论1.1研究背景[1]国家十一五规划提出了“发挥骨干企业作用，提高自主品牌乘用车市场占有率”的发展战略。进入21世纪，中国汽车工业蓬勃发展。有关统计数据显示：1994年至2009年的15年间，中国汽车产量增长将近10倍（见图1.1）。2010年1-3月，中国汽车产销分别达到455.45万辆和461.06万辆，同比增长76.99%和71.78%。在国家宏观政策调整和国民经济的发展下，中国汽车必将继续保持强劲的增长势头。图1.12000~2009年中国汽车产量与增长率汽车产业竞争激烈，突出表现在新车型的开发和制造能力方面，要求做到新车型能[2]及时推向市场、广泛适应客户需求和低生产成本，在很短的时间内研究、开发并制造出高质量的汽车。为了缩短新车型开发周期、降低研发成本，汽车厂商大多在同一个底盘上修改汽车外表件而产生新的车型，因此汽车的更新换代在很大程上取决于车身。这就要求加快汽车覆盖件模具的设计和制造。构成汽车车身或驾驶室、覆盖发动机和底盘[3]的薄金属板料制成的异形体表面和内部零件称之为汽车覆盖件。一般车身的覆盖件达[2]70～110件，在“车型设计→模具设计与调试→汽车批量生产”的整个周期中，覆盖[4]件模具的设计与制造工作约占2/3的时间。在制造费用上，汽车覆盖件占整个汽车制[5]造费用的70%以上。新车型的快速推出需要有大量的覆盖件模具重新开发。在美国、日本等汽车制造业发达的国家，汽车模具在整个模具产业中占有50%左右的份额。而在我国，仅有1/3左右的模具产品服务于汽车制造业。因此，汽车模具市场有相当大的发展空间。汽车覆盖件模具是整个汽车模具的重要组成部分，也是技术最密1 西安科技大学硕士学位论文[6]集、加工难度最大的部分。与一般冲压件相比：汽车覆盖件材料相对厚度小、结构尺[7]寸大、形状复杂、尺寸精度和表面质量要求较高，成形时坯料上各部分的变形复杂、差别较大，各处应力很不均匀，所以覆盖件模具设计和制造难度大、质量和精度要求高。国外模具企业大量采用CAD/CAE/CAM一体化技术，CAD应用中3D设计的比例高达70%~90%，普遍采用高速切削加工技术。CAD/CAE/CAM技术的采用，保守地估计，可使模具设计与制造周期缩短2/5，模具生产成本降低1/3，进而大大降低整车成本，[8]增强市场的竞争力。计算机仿真技术也已经得到汽车界的广泛重视和应用。如福特公司第一辆样车定型之前，已经完成了全部仿真分析的95%以上，可以在工程开发阶段节省四千多万美元，在制造阶段节省超过十亿美元。通用公司由于在轿车设计开发过程中[9]全面采用了计算机仿真技术，使得开发时间由原来的39个月减少到24个月。CAX技术的广泛、深入应用是模具设计制造的必然趋势。经过近几年的迅速发展，我国的汽车覆盖件模具技术有了长足的发展，生产能力大为提高。原来行业中公认的四大模具厂（一汽模具制造有限公司、东风汽车模具有限公司、天津汽车模具有限公司和四川成飞集成科技股份有限公司）都已有了生产大中型汽车覆盖件模具200万元左右工时的能力，模具年产值都超过1亿元，有的还超过了2亿[6]元。这与2000年相比，能力已增加一倍以上。部分模具企业也逐渐开始使用大型CAD/CAE/CAM软件UG、CATIA、Pro/E进行模具的三维结构设计，利用板料成形仿真软件Dynaform、AutoForm、PAM-STAMP等进行冲压仿真以指导覆盖件模具的型面设计，然而大部分企业的CAD/CAE/CAM技术还处在较低的水平，我国覆盖件模具在设计、制造、标准、管理等方面与国际先进水平有很大差距。目前，覆盖件模具普遍存在的问题是：模具设计的科学性差，对经验的依赖性强，试模周期长，重新设计的次数多，需要反复修改试制才能满足要求，增加了成本，降低了企业竞争力。在整个汽车工业中，冲压模具的发展至关重要。应用CAX、CAPP(ComputerAidedProcessPlanning)、PLM(ProductsLifecycleManagement)、VPD(VirtualProductDevelopment)、CIMS(ComputerIntegratedManufacturingSystems)技术，建立敏捷制造系[10]统，向智能化、信息化方向发展是汽车和模具工业的必然趋势。长期以来，国内外学者和工程技术人员在板料成形数值模拟和模具的分析优化方面作了大量的工作。但是对于复杂、质量要求高的覆盖件模具，设计出来的模具很大程度上取决于设计人员的经验和水平，往往不能满足产品的技术要求，特别是成形质量要求高的高级轿车外覆盖件更是如此。因此目前我国在大型汽车覆盖件模具的设计制造方面仍部分依赖外国的技术支持。其次，由于汽车工业对轻量化、安全、低排放、低成本及燃料经济性的要求越来越高，高强度钢板将更多的应用到汽车车身上来。根据美国新一代ULSAB－AVC（国际钢铁协会1994开始提出的超轻概念车项目）轿车车身的用材情况看，高级高强度钢板将占80%以上。高强度钢板的使用将给汽车冲压模具的设计、制造和使用带来诸多问题：2 1绪论[11]模具破裂、模具磨损、零件的精度、成形质量等等。再者，冲压设备方面重型压力机、[12,13]高速压力机和伺服压力机也正在或即将投入冲压生产线，模具在新型设备下和传统油压机的力学性态有何不同；各大汽车模具公司纷纷建立自己的CAD/CAM系统，大[14][15]力发展研究CBD(Case-BasedDesign)技术，特征技术和参数化设计，这些技术能否推动产生出新型的模具。综上所述，对覆盖件模具进行力学分析，研究模具结构的优化形式，探索新的模具结构模式有非常重要的实际意义和应用价值。1.2板料成形模拟的研究现状板料成形，从力学角度而言，是一个同时包含几何非线性（大位移、大应变）、材料非线性（弹塑性、弹粘塑性、各向异性）、边界条件非线性（接触和摩擦）的非常复[16]杂的力学过程。板料成形数值模拟研究始于上世纪60年代。早期研究采用有限差分法，所分析的问题都是象圆板液压胀形（图1.2a）、半球形冲头或平底冲头下的胀形（图[17]1.2b）和拉延成形（图1.2c）这样的较简单问题。(a)液压胀形(b)半球冲头下的胀形(c)半球冲头下的拉延图1.2板料成形数值模拟早期的研究对象1977年，在美国通用召开的一个关于板料冲压成形力学分析的研讨会上，Kobayashi[18]和N.M.Wang分别提交了一篇关于板料成形数值模拟的文章。Kobayashi用刚塑性有限元法分析了图1.2中的三个问题，将计算得出的应变值与采用有限差分法的计算值和[19]实验值做了比较。N.M.Wang用弹粘塑性有限元法仿真了平底圆凸模和半球凸模下的拉深过程，并研究了粘性、凸模速度对应变分布的影响。这两篇文章开创了板料成形有限元数值模拟的先河。从此，板料成形数值模拟沿着这两篇论文的思路发展起来了。这些早期研究仅针对二维和轴对称情况下的诸如胀形和拉深等一些简单的问题，还没有人能够完满地模拟一般的三维成形，尤其对于成形中的起皱现象仍然无能为力。[20]1978年，N.M.Wang和Budiansky基于非线性薄膜壳理论，用弹塑性T.L.格式分析了任意几何形状模具的冲压成形问题。首次考虑了坯料在模具表面滑动和粘着效应的接触摩擦现象，并提出了不能简单地把分析一般板料成形问题的方法直接用于分析汽车覆盖件这样复杂的冲压成形过程，而必须考虑复杂零件的成形特点。这有力地推动了有[21]限元方法在冲压成形分析中的应用。1985年，Toh和Kobayashi采用壳单元理论用刚塑性有限元法首次分析了方形盒的拉深过程，标志着冲压成形三维有限元模拟的开始。[22]1981年福特公司的S.C.Tang完成了汽车覆盖件的行李箱盖冲压成形的CAE仿真，3 西安科技大学硕士学位论文[23]接着又完成了轿车前翼子板冲压成形的CAE仿真，至此开始了覆盖件板料成形的仿真分析。80年代末，材料成形有限元数值模拟迎来了一个蓬勃发展的时期。在这一时期，各类单元（膜、壳、体)的性能研究，静力隐式积分(staticimplicitintegration)的加速收敛算法，静力显式积分(staticexplicitintegration)的步长确定方法以及动力显式积分(dynamicexplicitintegration)在板料成形模拟中的应用的理论基础，各类材料本构关系的比较，接触和摩擦的求解方法，拉延筋的模拟，回弹模拟，自适应有限元和网格重划分，[17]模具表面的几何表示等研究课题都得到了人们广泛的关注。到90年代中期，板料成形模拟技术逐渐走向成熟，已经能够模拟复杂的汽车覆盖件冲压成形，并利用仿真结果指导实际生产。许多汽车制造企业都建立了覆盖件冲压成形仿真分析系统，例如福特公司采用MTL-FRM软件，上海大众和东风公司采用AutoForm软件，丰田采用ITAS-3D，通用、本田、克莱斯勒采用LS-DYNA3D，宝马和菲亚特采用DYNAFORM。目前，冲压成形仿真技术已广泛应用于汽车工业，为企业降低成本、提高质量和市场竞争力发挥[24]了重大作用。然而，由于板料成形的数值模拟涉及力学、材料、数值方法、计算机科学、制造技术等众多学科，影响数值仿真技术的参数极其复杂，人们对覆盖件的成形问题的研究仍有许多问题尚待解决。目前的研究热点主要有：(1)摩擦与接触通常采用库仑摩擦定律确定冲压成形中的摩擦力，但冲压成形中的摩擦不仅涉及板[25][26]料与模具表面的相对运动，还涉及相对速度方向的变化。林忠钦等提出了修正的[27]库仑摩擦定律；B.H.LEE等以润滑油粘度和表面粗糙度为参数，建立起一个新的摩擦模型。接触问题可分为接触力的计算和接触搜索方法。接触力计算方法主要有罚函数法、拉格朗日乘子法和防御节点法；接触搜索方法主要有主从面法、级域法和一体休算法。(2)回弹分析和控制冲压件最终形状取决于回弹量，回弹量超过允许容差就是成形缺陷，严重的会影响整车装配。回弹预测方法主要有解析法、实验法和有限元数值模拟法，利用有限元模拟[28]回弹具有速度快，成本低，不受实验条件限制等优点。Bui等人利用改进的假定应变[29]单元(enhancedassumedstrainelements)对回弹进行了数值仿真；Papeleuxl等指出回量不仅与数值参数（厚向积分点数、单元类型等）有关，还与物理参数（材料属性、硬化准则等）有关。回弹预测后还要根据回弹量重新设计模具，这就是回弹补偿。较为有效的方法是改变模具的几何形状进行补偿，如位移调整法(DisplacementAdjustment,DA法)和控制表[30]面过度弯曲法(SurfaceControlledOverbending,SCO法)。尽管许多学者对回弹进行了4 1绪论大量的研究，但回弹的预测、控制和补偿精度仍有待提高。(3)材料屈服模型材料屈服准则是影响分析精度的主要因素，并且汽车工业中新材料不断出现，高强度钢、铝合金板、镁合金板已应用于车身冲压，如何建立合理的屈服模型是研究的重要[31]方向。常用的屈服准则有：Hill48，Hill79，Hill90，Barlat91，Corus-Vegter模型。Banabic在Barlat89屈服模型的基础上提出了新的正交各向异性屈服模型(BBC2000)，经过验证优于Hill48和Hill90模型。该模型的最新发展为Comsa-Banabic2007和BBC2008模型。[32]李玉强等介绍了常用的材料模型在CAE中的应用和发展情况。(4)坯料形状设计优化坯料形状设计是工艺设计的重要环节，合理的坯料形状不仅可以节省原材料、降低成本，还可以防止拉深件的破裂和起皱，提高成形性能。但精确计算坯料形状非常困难。[26,33][34]常用的算法主要有滑移线法、几何映射法、势场仿真法、比拟法等。Shimt和Son等提出了基于板料形状变化灵敏度分析的方法，可以确定任意形状拉深件成形时板料的[35]最优形状。随后他们又提出一种新的方法——INOV(initialnodalvelocity)，并计算了L形杯状拉延和油箱底的坯料形状，得到了满意的结果。除了这些研究课题之外，冲压成形模拟中模具系统（凸凹模、压边圈、拉延筋）和工艺参数（压边力、润滑条件、冲压工序）的研究，用数值实验方法即反求法(inversemethod)反求材料参数同样是冲压成形中研究的热点。1.3冲压模具结构分析及优化研究现状传统的冲压模具结构设计要依据各企业的压力设备、设计规范和准则，结合模具设计员的经验，并没有特别关注模具的受力和变形情况。实际中低碳钢冲压模具受力情况不恶劣，模具设计的比较保守（安全系数取得较大），再加上传统低强度钢板的成形性好，模具的微小变形对材料成形影响很小，所以关于冲压模具结构分析的研究国内外并不多见。现有的关于模具结构优化研究多集中在对工艺参数的优化上，如优化板料形状[36-38]和尺寸、压边圈形状、压边力大小、拉延筋、凸凹模具圆角，很少有对结构形式本身的优化设计。但在锻造、挤压领域，模具受力非常恶劣，模具的变形会直接影响制件的成形精度，所以一些学者对这一领域的模具结构进行了研究。他们研究的方法可以作为冲压模具研究的借鉴。[39]实验研究方面，Kojima等率先测量了棒拉拔过程中模具上的压力分布。[40]Matsubara等设计了一种传感器测量模具和工件接触面上的压力，并先后用这种传感[41][42]器研究了纯铝平面应变冷挤压过程中、铝罐轴对称反挤压过程中和轴对称冷锻过[43][44]程中，模具和变形工件接触界面间的压力分布和特性。Murata等采用电容式位移传5 西安科技大学硕士学位论文[45]感器测量锻造过程中模具的变形。Hillery等用嵌入式应变片研究了棒材拉拔中的应[46]力。MacCormack等用有限元法分析了冷锻中的工件破裂过程。英国Balendra教授领导的研究小组从上世纪90年代开始就采用分析、实验和数值模[47-50]拟的方法对锻模的变形进行了深入研究。Balendra教授分析了锻件尺寸偏差的来源并指出净成形工艺中制件尺寸误差的模拟、预测和补偿的关键是知道成形过程中模腔内的压力分布。由于工件和模具的变形相互耦合，而当时的计算机软硬件水平还不足以处理如此复杂的耦合问题。因此，直至本世纪初，Balendra研究小组主要采用基于实验实测接触力的研究模型。他们采用压电陶瓷材料开发了一种三轴力传感器，利用这种传感器采集到模腔内的压力作为载荷输入，用有限元法分析模具的弹性变形，为净成形制件尺寸补偿模具设计提供指导。[51]Balendra研究小组的Lu等1996年提出了三种有限元模型：①弹——塑性有限元模型，即将模具视为弹性体，工件视为弹塑性体，对模具和工件进行整体耦合分析；②刚——塑性模型，即将模具视为理想刚性体，工件视为弹塑性体进行分析；③基于实验实测接触力的分离模型，即建立模具的线弹性模型和工件的弹塑性模型，将实验测得的接触力分别施加到模具和工件上对他们分别进行分析。但第一种模型过于复杂，超出了[52,53]当时计算机的求解能力；第三种方法要先制造出模具，这个过程耗资、费时。Lu等采用刚——塑性模型成功地模拟了锻造成形中模具模腔压力分布。在当时的计算机水平下，他们提出一种解耦合的两步模拟法，使用分离的模型来预测制件的尺寸误差。结果表明，解耦的两步模拟法具有足够的精度。[54-56]随着有限元软硬件水平的发展，Lee等同时通过实验方法和数值模拟方法研究了冷锻模具的弹性变形。实验中利用应变仪测量了模具的弹性应变；数值模拟运用DEFORM-2D，在成形分析中分别把模具视为刚性体和弹性体，工件分别采用刚塑性和弹塑性模型。通过不同组合与实验对比发现：将模具视为变形体，采用耦合方法分析得到的模具上的应变与实验更符合，但是耦合分析方法耗时长的多。冲压模具，特别是汽车覆盖件冲压模具通常体积大而且结构复杂，冲压模具结构分[57]析面临如下难题：(1)在建模方面，冲压模具尺寸大、结构复杂，有限元模型的精度和求解规模的平衡是一个大问题；(2)冲压成形过程中变形板料和模具之间的接触状态复杂，冲压模具实际受到的载荷信息难以准确获得；(3)在模型确认方面，因为冲压模具价格昂贵，实验耗时长久、耗资巨大，而且到目前为止没有切实可行的测量方法；(4)实际生产中随机因素不可避免，工艺波动对模具结构受力的影响复杂。[11,58,59][52,53]张贵宝、陈军等借鉴Lu等提出的解耦合两步模拟法思想，先通过板料6 1绪论成形数值模拟得到板料变形对模具的作用力，然后通过边界力的映射算法将作用力施加到模具的有限元模型上，代替板料对模具的作用力对模具进行结构分析。应用这种方法模拟底面为曲面的带法兰盒形件拉深成形过程，材料为高强度钢板DP600。实验中在模具上开3mm~5mm的浅槽，在槽内贴三轴应变花测量模具受力，计算关键部位的等效应力。经过对比验证，数值模拟和实验结果等效应力最大误差为14.81%。说明数值模拟方法得到冲压模具受力是可信的，而且兼顾了计算效率。1.4本文研究的思路及主要内容1.4.1研究思路本文将按照解耦合两步模拟法的思想，研究覆盖件模具在冲压过程中的力学响应，并在此基础上利用拓扑优化技术优化传统的模具结构形式，探索一种新型的模具结构。图1.2基于板料成形数值模拟的模具结构优化分析流程[60]本文提出的基于板料成形数值模拟的模具结构分析和优化流程如图1.2所示。获得冲压过程中模具表面荷载是覆盖件模具结构分析和优化的前提。为了使获取模具荷载和结构分析优化无缝衔接起来，首先划分模面网格，在冲压成形模拟、结构分析和优化设计中，模具采用同样的表面网格。这样通过冲压成形模拟所获得的模具表面节点荷载7 西安科技大学硕士学位论文[59]可以直接施加到结构分析和优化的有限元模型上，避免了文献中复杂的节点力转换和映射。冲压成形模拟中要求模面网格尽可能精确描述模面的几何形状，而结构分析优化中又对模面网格的质量提出了一定的要求。用专业的前处理软件HyperMesh划分模面网格，适当减小网格尺寸以保证模面网格质量，从而一次划分的网格可以在冲压模拟和结构优化分析中两次使用。1.4.2研究的主要内容(1)汽车覆盖件冲压成形的有限元模拟。通过冲压模拟得到冲压过程中板料对模具表面的作用力。以材料成形模拟的专业CAE软件PAM-STAMP2G2008为数值模拟平台，对两个典型的外覆盖件（行李箱盖外板、引擎盖外板）的拉深成形过程进行模拟，获取冲压过程中板料对模具的作用力。讨论板料成形模拟中几个关键问题：算法选择、材料模型、单元类型和接触与摩擦处理。(2)冲压模具的结构优化方法。以冲压模拟得到的板料对模具的作用力为荷载，对模具结构进行优化。利用著名的优化软件OptiSturct在模具结构分析的基础上，采用拓扑优化技术对模具结构进行优化。并针对以上两种典型的覆盖件模具，在CAD系统CATIA中设计出新的模具结构。(3)冲压模具结构分析方法。分析模具在冲压过程中任一时刻的力学响应，并力图了解模具响应在整个连续冲压过程中的变化规律。将板料成形模拟得到的模具受力情况施加到模具有限元模型上，利用集成在Altair公司HyperWorks中的有限元求解器RADIOSS分析传统模具和优化后模具的受力情况。8 2板料成形的数值模拟2板料成形的数值模拟2.1引言要分析优化模具结构，首先要得到的是冲压过程中模具的受力情况，即板料对模具表面的作用力。冲压过程中板料对模具表面的作用力为摩擦力，可分解为法向接触力和切向接触力。因此本文首先要进行的工作是进行覆盖件冲压成形过程的模拟。冲压成形模拟的正确与否、精确与否直接决定了后续结构优化和分析的准确性和有效性。得益于计算机技术的迅速发展，融合了CAD（计算机辅助设计）、数值计算、CG（计算机图形学）、CAM（计算机辅助制造）的CAE技术广泛应用于汽车、钢铁和航空航天等工业。在汽车工业，美国的通用、福特，德国的大众、奔驰，日本的丰田、三菱、日产等著名厂商应用CAE软件已20年左右，用其指导和优化板料设计，模拟板料成形过程，明显减少了模具的试模次数，缩短了模具制造周期，降低了成本，加快了整车生产。板料成形又称为冲压，这种成形方法通常是在常温下对板料进行成形，所以也称冷冲压。对于某个金属薄板零件，一般要经过多个工艺过程，才能直接用于装配。常见的[61]工艺过程有如下几种：(1)冲裁(Blanking)为了生产和运输方便，板料多是卷料。首先要将板材分离为一定形状的坯料才能用于冲压生产。冲裁时板料的剖面放大如图2.1所示。这是一个复杂的塑性剪切和破坏过程，在板料边缘会有局部的硬化现象。冲孔和修边也是冲裁过程。图2.1剪切时板料的局部剖面(2)弯曲(Bending)最简单的情况是沿一条直线将板料折弯。仅在弯曲部分发生塑性变形。如材料延性不够，在弯曲部位外表面容易发生破裂。但由于回弹的原因，精确控制弯曲角度非常困9 西安科技大学硕士学位论文难。常见的几种直线弯曲工艺如图2.2所示。图2.2(a)是折弯，板料一边被夹持不动，另一边被夹紧绕边缘旋转一定角度；图2.2(b)是压弯，冲头向下运动将材料压入V形模具中；图2.2(c)为滚轧成形，在滚轧机械中有成套的轧棍，大型面板及复杂的槽型截面都可以用这种工艺制造；2.2(d)是翻边，通常用于对冲压件边缘进行处理。将板料边缘弯曲直至贴在板料背面，从而使零件形成圆润的边缘，叫做包边。(3)拉深成形(Drawing)拉深也叫拉延，是利用模具使冲裁后的坯料变成开口的空心零件的冲压工艺方法。通常情况下一副拉延模包含凸模、凹模和压边圈。拉延工艺广泛应用于汽车覆盖件，与其它工艺配合，可以制造出形状极为复杂的零件。按照变形力学特点可分为：直壁圆筒形零件、盒形件、曲面形状零件（指曲面旋转体）和非旋转体曲面形状零件等四种类型。图2.3是一个典型的杯形拉延。图2.2几种弯曲工艺图2.3杯形拉延此外还有胀形、液压成形、管胀成形、热成形等多种工艺方式。通常情况下一个工业零件的成形过程包含许多基本的成形工艺。汽车覆盖件是板材冲压成形加工中最复杂的，一般要经过拉深、修边、翻边、冲孔等多道工序才能得到一个可用于装配的零件。能否得到合格的覆盖件，关键在于冲压的可行性。冲压的质量直接决定了零件的最终品10 2板料成形的数值模拟质，本文就是基于覆盖件冲压过程模拟而对拉延模具进行结构优化和受力分析的。为了给后续模具结构优化和分析提供准确的荷载信息，覆盖件拉延成形过程模拟必须达到以下目标：(1)拉延成形模拟中需采用与模具结构优化和分析中相同的模具表面网格。这样由成形模拟得到的模具表面节点力就可以直接施加到结构优化和分析模型中；(2)拉延成形模拟需保证模拟出的零件成形性能良好，即零件可以生产；(3)拉延成形模拟中的材料参数、工艺条件要尽可能接近实际生产情况，以保证结果与实际的可比性。由于板料成形是包括几何非线性、材料非线性、边界条件非线性的大变形问题，有限元算法的选择、单元类型、材料模型以及接触和摩擦的处理都将影响模拟结果的精度，下面将对这些问题分别论述。2.2计算平台的选择目前应用较为广泛的成形CAE软件主要有：瑞士的AutoForm、美国ETA(EngineeringTechnologyAssociatesInc.)的Dynaform，法国ESI的PAM-STAMP。[62]PAM-STAMP对板料成形在破裂、起皱、轮廓形状等方面的模拟精度达到90%以上。本文将采用PAM-STAMP2G2008进行板料冲压成形的模拟。PAM-STAMP2G是法国领先的数值模拟软件供应商ESI公司虚拟制造解决方案中的一个软件，是板料成形全工业链的完整解决方案。具有模面设计、板件展开反求、快速验证分析、精确成形过程分析、回弹分析等模块，还具有模面回弹补偿、板料形状优化、液压成形、橡胶成形、热成形等高级功能。PAM-STAMP2G求解器兼有隐式和显式算法，可以为不同的模拟过程设置不同的算法。重力阶段和回弹模拟采用隐式算法，精确成形模拟、外部摩擦算法可以使用显式算法，既保证了求解的精度，又显著降低了模拟时间。PAM-STAMP2G具有1阶和2阶Belytschlko-Tsay(BT)、Belytschlko-Wong-Chiang(BWC)、BWCFULLYINTEGRATED(FBWC)、Hughes-Tezduyar(HT)等多种形式的壳单元，用户可以根据精度和计算机配置选择合适的单元类型，同时它还可以模拟六面体单元。材料模型方面，PAM-STAMP2G具有Hill48、Hill90、Barlart89、CORUS-VEGTER、MOONEY-RIVLIN等屈服准则，支持MATFEM、EWK破裂准则，可以模拟各向同性及正交各向异性材料等向强化及随动强化行为。在接触方面，PAM-STAMP2G采用主从面法判断接触，可以使用罚函数法、非线性函数法、拉格朗日法计算接触和摩擦。在PAM-STAMP2G中对零件冲压成形过程进行模拟，调整各个成形参数（如压边力，拉延筋，板料线等），检测模拟结果的变薄率（满足刚度要求）、应变状态（一般通过成形极限图判定）等，直至模拟结果符合工艺要求。利用其后处理功能即可得到模11 西安科技大学硕士学位论文面上任一节点处在某个时刻的法向和切向接触力。PAM中的计算流程如下：图2.4PAM-STAMP模拟基本流程2.3板料成形的塑性力学概念2.3.1真实应力和真实应变在单轴拉伸实验中，常用的是工程应力和工程应变（或称名义应力，名义应变），定义式为P(2.1)eng.A0lll0100%(2.2)eng.ll00式中P为拉伸荷载；A为试样的初始横截面积；l为试样拉伸后的长度；l为试样的伸0长量；l为试样的标距。0但在实验过程中，试样的横截面积是不断减小的，所以工程应力并不能真实反映试样的受力情况；工程应变的计算也是基于试样的标距，也不能反应试样真实的变形。为了克服这些缺点，定义真实应力12 2板料成形的数值模拟P(2.3)A式中A为试样变形后横截面积；根据体积不变假设（不可压缩性）AlAl(2.4)00故真实应力Pl(2.5)Al00假设试样长度有一个增量dl，则真实应变的增量为dld=(2.6)l在小应变情况下，ll，真实应变增量和工程应变增量相近；在大应变情况下，0两者有较大差别。如果应变在一个方向上连续单调变化，对应变增量积分即可得到真实应变ldlldln(2.7)l0ll0与工程应变相比，真实应变有以下优点：(1)真实应变为可比应变，工程应变为不可比应变。设某物长l拉长一倍后变为2l，缩短一倍后变为l2，这两种情况下的工程应变分000别为22llll0000100%100%100%50%(2.8)ee12ll00真实应变值为22ll00ln0.693ln0.693(2.9)12ll00拉伸一倍与缩短一倍，这两种情况下的变形程度应该是一样的，可是用工程应变表示的变形数值相差甚远，而用真实应变表示数值相等，符号相反，显然后者更合理。(2)真实应变有可加性。设某物原长l，经历l后变为l，用工程应变表示这一过程为012llllll102120=(2.10)ee12totallll010显然：totale12e用真实应变表示为lll122lnlnln(2.11)12totallll01013 西安科技大学硕士学位论文则有：12ln(ll10)ln(ll21)ln(ll20)total(3)三个正应变之和为体积应变。设一长方体长、宽、高各为：lwt,,，变形后尺寸为：lwt,,，则体积应变有000ln(vv)ln[(lwt)(lwt)]ln(ll)ln(ww)ln(tt)(2.12)v0000000lwt对于塑性成形问题，由于问题的非线性，一般将整个变形过程看作是多个小变形过程的累积，最终的应变是通过对各阶段的小应变相加得到的。显然，真实应变的可加性符合这种需求。各阶段的小应变直接相加所得到的应变为对数应变。因此在塑性成形理论中，普遍采用以对数应变表示的真实应力——应变曲线。将真实应力应变画在对数坐标系中，一般的退火低碳钢会表现出如图2.5所示的性质。在弹性阶段的小应变区，近似地有EE或logloglog(2.13)在塑性阶段，曲线可用以下模型表示nK或log=log+logKn(2.14)式中E为弹性模量；n为硬化指数；K为强化系数；式(2.14)称为幂指型硬化模型，常用来描述退火低碳钢板材的塑性性质。图2.5普通板料的真实应力-应变的对数坐标曲线2.3.2材料的厚向异性系数设试样在宽度方向的应变ln(ww)，在厚度方向的应变ln(tt)，定义板料w0t0的厚向异性系数r为ln(/ww)w0r(2.15)ln(/)ttt0r值一般在特定的工程应变时测得，例如：eng15%，对于各向同性材料，应有r1。r1，板料在宽度方向上变形比在厚度方向上变薄更容易；r1，板料在厚度方向上变薄比在宽度方向上变形更容易。材料r值的作用在拉延成形中最为显著。r越大材料抵14 2板料成形的数值模拟抗失稳变薄的能力越大，越能发挥拉伸失稳前的最大强度，拉动凸缘部分形成更深的压延件。实验表明：r值与极限拉延比的关系基本上是线性的。对于板料来说，通常在轧制方向、与轧制方向成45°以及垂直于轧制方向上取三个试样测定的材料厚向异性系数。用下标表示与板料扎制方向的角度（图2.6），r表示轧0制方向；r表示与轧制方向成45°；r表示与轧制方向成90°。4590图2.6板料取样示意图2.3.3板料的成形极限图单轴拉伸试验中，通常用伸长率和断面收缩率来表征材料的延性。然而在板料成形中，材料的应力应变状态比较复杂，如何才能简单准确地判断材料破裂是一个重要问题。板料成形中厚度方向的应力与板料平面内的应力相比很小，这样就可以假设板料法向应力为零，则板料成形中材料的应力状态简化为平面应力。许多板料成形工艺中，板料在板料平面承受双向拉伸，破坏发生在局部颈缩之后。局部颈缩与圆形试件拉伸试验中的颈缩不同。对于金属板料试件，颈缩时试件在宽度和厚度方向上都在减小，而局部颈缩时，试件仅在厚度方向上减小，而宽度方向上没有减小。试样破裂发生在局部颈缩之后。也就是说，局部颈缩可以作为试样破坏的判定条件。成形中板料某点的应力应变状态如图2.7所示，设应变比，应力比。2121图2.7板料成形中某点的应力应变状态[63]在应变比保持常数的情况下，可以导出板料发生颈缩和局部颈缩条件分别为2*2(1n)(2.16)12(1)(22)15 西安科技大学硕士学位论文*n(2.17)11将颈缩和局部颈缩的判定准则式(2.16)、(2.17)绘制成图如图2.8所示。板料成形中一个主应变方向与板面垂直，即厚向应变。另两个主应变均在板面内。网格分析法常用来测量板料成形中各点在板面内的主应变。试验时先在板料上印制圆形（或方形）网格。板料变形后，圆变形为椭圆，测量椭圆两个主轴的长度就可以很方便地得到该点的主应变。把板料破坏时各点的主应变值画在一张图上，就叫做成形极限图(FLD,FormingLimitDiagram)，也叫Keeler-Goodwin图。图2.9是典型的低碳钢成形极限图。图中曲线是由实验得到的材料破裂判断准则，叫做成形极限曲线(FLC,FormingLimitCurve)。成形极限图上横轴表示板料平面内的最小主应变，纵轴表示最大主应变，[64]只把零件成形后各点的主应变在FLD中绘出，就可以判断零件的应变状态，如图2.10所示。大致分为四个区域：安全区、破裂区、拉拉应变区、拉压应变区。成形极限图的出现可以使人们根据它评估不同板料在复杂应变状态下相对成形性能的优劣，对板料的冲压成形是一种有效的方法。材料的成形极限图通常通过钢模胀形试验测得。PAMSTAMP2G中根据成形极限图将各点的成形质量分为六类（图2.11)：①破裂(cracks)、②过度变薄(excessivethinning)、③安全(safe)、④拉伸不足(insufficientstretching)、⑤起皱趋势(wrinklingtendency)、⑥强烈起皱(strongwrinklingtendency)。[63][63]图2.8颈缩和局部颈缩的判定准则图2.9典型的低碳钢成形极限图16 2板料成形的数值模拟图2.10成形极限图FLD图2.11PAMSTAMP中的成形质量分类北美深冲协会研究发现各种低碳钢在力学性能没有异常变化的情况下，FLD形状基本都是一样的，只是整个曲线的高低有所不同。利用这个特点只要知道FLC上的一点，就可以根据已知低碳钢的FLC推出材料的FLC来。这个点一般取FLC的最低点，叫做FLC0。可以用下面的Keeler’s公式来表示：t2.54mmFLC(23.314.13)/0tn.21022.54mmt5.33mmFLC(2020.67t1.94)/0.21tn(2.18)0t5.33mmFLC75.125/0.2n0上述公式中，当n0.22时，取n0.22；t表示板料的厚度(mm)；n表示板料的硬化指数。这种方法仅适用于低碳钢板，对特殊钢、铝板及其他材料，要根据试验得出的曲线输入。2.4冲压成形模拟的关键问题2.4.1算法选择板料成形中为了得到加载过程中应力和变形的演变历史以及保证求解的精度和稳定，通常采用增量法求解。应用全Lagrange形式建立系统离散后的有限元方程为MuCuKuP(2.19)式中M、C、K、P分别为质量矩阵、阻尼矩阵、刚度矩阵和节点力向量，u、u、u分别为节点加速度、速度和位移向量。增量法求解中，利用上一时间步或几个时间步的计算结果，求解下一个时间步的结果，如此往复直至模拟结束。对于相邻几个时间步t内，u、u、u之间的关系做函数假设，以消除对时间的积分运算。由于这些假设关系的不同，衍生出两种不同的解法：隐式算法和显示算法。又根据求解过程中是否考虑系统的动力效应（忽略加速度和速度的影响)，细分为静力隐式算法(StaticExplicitalgorithm)、静力显式算法(StaticImplicitalgorithm)、动力隐式算法(DynamicImplicitalgorithm)、动力显式算法(DynamicExplicit17 西安科技大学硕士学位论文algorithm)。下面说明在板料成形模拟中常用的动力隐式算法和动力显式算法。[64](1)动力隐式算法以Newmark方法为例，在t时间域内，假设uu[(1)uu]ttttttt12(2.20)uuutuutttttttt2式中和取不同值代表不同的数值方案。将式(2.20)代入(2.19)式后有1111KMCuPMuu1utt22ttttttttttt2(2.21)Cu11ututttt2其中K是u的函数，即式(2.21)两端都有未知解，这就是说每一个时间步长内，都要进行一定量的迭代计算才能获得当前时间步长下的位移增量utt。动力隐式算法有如下特点：①对于每一载荷增量步都需进行迭代，直至收敛为止；②每一节点都需满足平衡条件及接触条件；③对于复杂非线性问题，难以保证迭代计算的收敛性，即不一定能获得问题的解；④计算所需存储空间和CPU时间随求解问题规模的增大急剧增大（问题规模与系统自由度大致呈平方关系）。(2)动力显式算法将系统平衡方程式(2.19)中的质量矩阵和阻尼矩阵对角化，对时间进行离散处理，并用中心差分方法进行时间积分，假设1ut2utt2ututtt(2.22)1uuttutt2t将(2.22)式代入系统的运动方程(2.19)，则得到中心差分法的递推公式11211MCuPKMuMCu(2.23)2tttt2t2tttttt2t上式左边的系数是个常数矩阵，右边在得出t和tt时的计算结果后，也是一常数向量，整个方程就变成了一个简单的线性方程组。动力显式算法有如下特点：①不需要形成刚度矩阵，不需要求解联立方程组；②对于每一载荷增量步无需迭代，不存在收敛性问题；18 2板料成形的数值模拟③计算所需存储空间和CPU时间随求解问题规模线性增加；④不检查增量步结束后的平衡条件，算法是条件稳定的，要求数值模拟所用时间步长小于临界时间步长。临界时间步长以可由下式估算ltK(2.24)cr1c其中Ec(2.25)21v式中K是为保证显式算法计算稳定而采用的比例系数，一般取0.9；l为单元的最小尺1寸；c为材料中的声速；为材料密度；E为弹性模量；v为泊松比。由于隐式算法每一步迭代都必须求解复杂的非线性方程组，使得它在每个时间步上的花费（时间和空间上的）明显高于显式方法，并且隐式算法存在迭代收敛性的问题，一旦某一步收敛条件不满足，将导致收敛错误停止运算。因此隐式算法对于变化剧烈的问题（爆炸、高速冲击、复杂接触、剧烈变形）往往无法稳定求解。而对于显式算法来说，无需进行迭代，从而不存在收敛性的问题。只要时间步长足够小，总能获得最终结果。鉴于显式算法的优点，动力显式有限元方法成为了模拟复杂的板料成形过程（汽车覆盖件拉延）的主流算法。PAM-STAMP2G将板料拉延过程人为分为三个阶段重力(Gravity)、压料(Holding)和拉延(Stamping)，各阶段定义将在后文阐述。对于重力过程，板料受力情况简单，变形较小，采用隐式算法。对于压料和拉延阶段，板料变形剧烈，接触条件复杂，采用动力显式算法求解。-6-7显式算法时间步长取得很小（按(2.24)式求得的临界时间步长约为10~10s），虽然这对处理冲压成形中复杂的接触问题有利，但会导致模拟时间大大加长，通常采用速度加速法(Velocityacceleration)和质量放大法(Massscalingacceleration)来增大时间步长，减少板料成形模拟时间。实际冲头速度一般介于0.1~1m/s，模拟中一般取为3~10m/s，速度放大了10~30倍。质量放大通过人为增加材料密度而增大时间步长来实现，通常质量放大系数取1~10倍。这两种方法可以显著降低模拟时间。2.4.2材料模型(1)屈服准则板料有限元模拟中常用的屈服准则有：Hill48，Hill90，Barlat模型。由于Hill48模型有良好的理论基础且计算简单，钢板的成形分析中，本构模型一般选用Hil148模型。19 西安科技大学硕士学位论文Hill48准则的屈服函数如下12222Hill48H1122F2233G33112N12(2.26)2其中：F、G、H、N是与材料各向异性有关的系数，HG2。2r2(2r1)(rr)045900F;G;Nr(1r)1rr(1r)9000900当FG1,N3时，就等同于Mises屈服准则。Barlar89准则对多晶体材料及铝合金、钛合金等的通用性很好，准确性明显高于Hill90准则，其屈服函数为mmmmFKKKK22K2(2.27)121222hhxxyyxxyy22其中：K;Kp，x，y为板料轧制方向和其垂直方向。12xy221pmmhr0(1r90)(1rr0)90；2rr090(1r0)(1r90)；p()[2(2a2c)]；bsl式中是等双拉状态的Cauchy主应力；是单项拉伸状态的Cauchy主应力；p1b1为沿轧制方向的等效应力；是纯剪切状态进的屈服剪应力；对于面心立方晶格材料，slm取8；对体心立方晶格材料，m取6；公式中只包含x、y平面内的三个应力分量，因此上述屈服函数只能应用于平面应力状态。(2)硬化法则由于大多数材料都存在包辛格效应，材料在进入塑性状态后卸载，然后再加载，屈服面会随着前面加载的历史有所改变，通常表现为应变强化现象。硬化法则是用来规定材料进入塑性变形后的后继屈服函数（又称加载函数或加载曲面）在应力空间中变化的规则。一般来说后继屈服面不仅与材料所受的应力状态有关，而且与塑性应变和强化参数有关，其后继屈服函数为Fk(,)0(2.28)ij常用的硬化法则有各向同性硬化模型、随动硬化模型和混合硬化模型。①各向同性硬化模型。此法则规定，当材料进入塑性变形以后，加载曲面在各方向均匀地向外扩张，但其形状、中心、及其在应力空间中的方位均保持不变。例如对于0的情形，初始屈服轨迹和后继屈服轨迹如图2.12(a)所示。如采用VonMises屈服3条件，则后继屈服函数可表示为F(,)kfk0(2.29)ij20 2板料成形的数值模拟112其中：fss,k()。是现时的弹塑性应力，它是等效塑性应变的函数；ijijspsp23s是偏斜应力张量分量。ij对于等向硬化模型，屈服面在所有方向上均匀扩张。加载面仅由其曾经达到过的最大应力点决定，与加载历史无关。②随动硬化模型。材料进入塑性后，加载曲面在应力空间作一刚体移动，但其形状、大小和方位均保持不变。后继屈服函数表示为Fk(,,)0(2.30)ijij0其中是加载曲面的中心在应力空间内的移动张量，k为初始屈服条件中的材料参数。ij0与材料硬化特性以及变形历史有关。Prager运动硬化法则规定加载曲面中心的移动是ij沿着表征现时应力状态的应力点的法线方向，见图2.11(b)。如采用VonMises屈服条件，则Prager运动硬化法则的后继屈服函数为F(,,)kfk0(2.31)ijij00112其中：fss,k；是材料的初始屈服应力。ijijijijs0s023图2.12各种硬化模型示意图③混合硬化模型。试验证明弹塑性材料的屈服强化过程通常同时具有等向强化和随动强化特性，为了同时考虑各向同性硬化和随动硬化，Hodge首先提出了混合硬化模型。混合硬化模型的后继屈服函数可以表示为F(,,)kfk0(2.32)ijij112其中：fss,k(,M)；是移动张量的偏斜分量。ijijijijspij232.4.3单元类型板料成形模拟中常用的单元类型有：膜单元(Membrane)、体单元(Volume)和壳单元(Shell)。膜单元构造格式简单，对计算机内存要求小，计算效率高，但是薄膜单元忽略了板料成形过程中的弯曲效应，只考虑了沿厚度方向均匀分布且平行于中性面的应力，21 西安科技大学硕士学位论文忽略弯矩、扭矩和横向剪切，认为应变沿厚度方向是均匀分布的，无法实现回弹、起皱等现象的模拟。AutoForm软件使用一种近似弯曲刚度的膜单元(Bendingenhancedmembrane,BEM)，可以近似地模拟薄板(2.5mm)一般地拉延过程，可以在很短时间内[65]得到有相当精度的结果。所以AutoForm在汽车和模具企业得到了广泛应用。基于连续介质的实体单元的表达格式比薄膜单元更简洁，但为了准确模拟板料冲压成形问题，在板料厚度方向至少需要布置4层单元，因板料在厚度方向的尺寸相对于平面内的尺寸要小得多，为避免刚度矩阵的奇异，要求单元划分甚密，这往往导致计算量过大，模拟汽车覆盖件等复杂零件的成形问题难以忍受。通常只用在板厚大于4mm的情形下。壳单元大致分为二类：一类是基于Kirchhoff薄板理论的壳单元；另一类是基于Mindlin厚板理论的壳单元。Kirchhoff单元假设变形前垂直于中面的法线在变形后仍然为直线且与中面保持垂直，忽略了板壳横向的剪切变形；还假设平行于板中面的面层之间没有压力。这种方法需要构造C1连续单元。C1连续单元是指在相邻单元在共有节点处不仅要保证位移一致，还要保证在该节点的转角一致。这样构造出的单元形式非常复杂，效率很低。基于Mindlin理论的壳单元认为变形前垂直于中面的法线变形后仍然保持直线，但不一定垂直于变形后的中面。Mindlin单元将板单元位移和转角分别插值，只需构造C0连续单元，使问题简化。因此在商用有限元软件中（如PAM-STAMP2G，Dynaform等）普遍采用Mindlin厚壳单元。Mindlin单元在分析料厚较薄（即t/l<<1）的冲压成形及回弹问题时，易发生剪切锁死和零能模式。采用缩减积分和沙漏能控制方法解决剪切锁死以及由缩减积分带来的零能模式。[66]PAM-STAMP2G中支持的壳单元有：(1)Belytschlko-Tsay(BT)单元。采用缩减积分的线性壳元，计算效率高。但BT单元对单元翘曲敏感，不适用于单元翘曲严重的分析中。在板成形计算中，只要单元不是太粗糙，翘曲现象一般较少发生。因此BT单元是板料成形分析中最常用的单元。(2)Hughes-Tezduyar(HT)壳单元。HT单元采用选择性缩减积分的方法，避免了沙漏现象。对于沙漏现象严重的振荡或冲击主导的过程，可以用来替代BT单元。HT单元直接退化于实体单元，计算精度高，尤其是在单元扭曲较大的状态下，明显高于BT单元。但其计算耗费约为BT单元的4倍左右。(3)Belytschlko-Wong-Chiang(BWC)单元。BWC单元也采用一点积分控制沙漏现象，是BT单元的改进形式。对于翘曲严重的分析仍然适用。计算耗时约为BT单元的1.2~2倍。(4)BWCFULLYINTEGRATED(FBWC)壳单元。完全积分的BWC单元，在单元平面内使用4点的完全高斯积分，避免了沙漏现象和剪切锁死现象。引入了人工剪切应变场和薄膜应变场，避免了完全积分导致的适应能力不强的缺陷。计算时间较BT单元22 2板料成形的数值模拟有明显增加，约为BT单元的3~4倍。2.4.4接触和摩擦在覆盖件冲压成形过程中，由于模具刚体型面的复杂性，使得板料与模具表面之间的接触问题变得非常复杂，接触边界条件处于连续变化之中，由此产生的接触面间的相互作用也是复杂的和动态变化的。应用有限元方法分析冲压成形问题时，考虑包含润滑现象的接触问题非常困难。处理接触问题包含相互关联的两个方面：接触搜索和接触力计算。(1)接触搜索算法接触搜索方面，最主要的算法有：级域法、一体化算法和主从面法。①级域法。通过对接触体分别归类，定义层次，接触搜索从上级子体系到下级子体系按约定方式逐级有系统地进行，若上级搜索排除了接触的可能，其下级接触也被除去，从而大大减少计算量。②一体化算法。采用“桶式分类法”的搜索策略。将网格节点所在的空间，用许多相同规格的“桶子”（立方体或方格）进行分隔，按约定的排列方式将这些桶子不留间隙地充满整个节点空间，使得所有节点都落在这个桶子内。当接触可能发生时，离某个节点最近的点必定在该桶或周围几个桶内，从而提高搜索速度。③主从面法。适用于两个接触面已知的接触搜寻。两个将要互相接触的面，一个被指定为主面(Master)，另一个被指定为从接触面(Slave)。从接触面上的节点不允许穿透主接触面，但是主接触面上的接触节点可以穿透从接触面。(2)接触力算法接触力计算方面常用的是罚函数法和拉格朗日乘子法。①罚函数法。接触力大小与穿透深度成正比。穿透越大，接触力越大，阻止其穿透的惩罚措施越强，因此该接触力又称为罚函数值，比例系数为罚因子，方法被称为罚函数法。②拉格朗日乘子法。通过引入拉格朗日乘子，将接触力作为未知量代入运动方程，联立求解方程组得出接触力，该算法得出的结果比较精确，但因为求解联立方程组的原因，计算耗时比罚函数长，对网格要求也比较高。PAM-STAMP2G中的接触处理方法主要有：罚函数接触算法(Penalty)、非线性函数接触算法(Accurate—Nonlinearpenalty)和拉格朗日接触算法(Lagrange)。Accurate接触算法本质上还是罚函数法，只是在接触力计算中罚因子不再是常数，而是非线性变化，从而小穿透也能迅速产生足够的接触反力，避免穿透的发生。Accurate算法计算精度高，是PAM-STAM2G的默认算法。但它对初始穿透和模具网格要求较高。初始穿透不能超过0.01，模具网格必须连续，且尽可能精细。23 西安科技大学硕士学位论文2.5实例分析对两个典型的汽车覆盖件进行冲压成形模拟，导出冲压过程中模具的节点力信息。一个是行李箱盖外板（如图2.13所示），另一个是引擎盖外板（如图2.14所示）。两个零件均采用单动拉延成形，成形装置的简图见图2.15。成形过程可分为三步：首先，将板料放置在凸模和压边圈上；然后，凹模向下运动至与压边圈贴合将板料压紧；最后，凹模、压边圈连同板料一起向下运动直至拉延到底。对应于这三个步骤，PAMSTAMP2G中分别称为：①重力阶段(Gravity)：板料落在凸模和压边圈上在重力作用下自然下垂的过程。板料发生弹性变形。采用隐式算法求解；②压料阶段(Holding)：凹模向下运动与压边圈闭合，将板料压紧的过程。板料边缘与压边圈接触部分被凹模和压边圈挤压，可能发生塑性变形，只要变形部位不在产品范围内，可以接受。中部可能被凸模顶起，发生弹性变形而不能有塑性变形，否则会在最终产品表面留下痕迹，影响成品质量。采用显示算法求解；③成形阶段(Stamping)：凹模、压边圈带动板料向下运动直至到底的过程。板料大部分发生塑性变形至设计形状的过程。采用显示算法求解。图2.13行李箱盖外板图2.14引擎盖外板图2.15典型的单动拉延示意图24 2板料成形的数值模拟2.5.1行李箱盖外板冲压成形模拟(1)模型建立为了保证节点的连续性和一致性，在冲压模拟中使用的模面网格应和结构优化分析中使用的表面网格相同。在冲压成形模拟中，模面网格要求准确描述模面的几何形状，以便精确判断和计算接触与摩擦。对于单元形状和质量要求并不是很高。由于CAD文件中不可避免地产生许多碎面和多余边界线（如图2.16所示）。如果用冲压模拟软件划分模面网格，将产生很多狭长的四边形网格和内角比较小的三角形网格，单元长宽比(aspectratio)甚至大于100，三角形单元最小内角小于20度（如图2.17所示），这在结构优化分析中是不能接受的。为此，使用专业的前处理软件HyperMesh划分模面网格。它可以自由控制网格大小和形状，选择合适的网格划分算法，直至得到质量满意的网格。图2.18是HyperMesh划分的模面网格。图2.16凹模模面几何Maximumaspectratio:181.7Maximumaspectratio:2.54Minimuminnertriangularangle:0.313°Minimuminnertriangularangle:24.23°图2.17PAMSTAMP2G划分的模面网格图2.18HyperMesh划分的模面网格将HyperMesh中划分的模面网格导入PAMSTAMP2G中建立完整的冲压模拟模型如图2.19所示。冲压成形模拟中模具被视为刚性体，只包含参与成形的模面部分。导25 西安科技大学硕士学位论文料销和定位块是辅助板料定位之用。板料为0.8mm的DC06钢板，材料性能参数根据实际生产情况在默认参数的基础上修改而来，材料参数见图2.20。其中E为弹性模量，r0、r45、r90为厚向异性系数，应力应变曲线为真实应力和真实应变曲线。屈服准则为Hill48准则；选择等向硬化模型；润滑条件为一般润滑，摩擦因数为0.12；拉延过程压边力为1400kN；凹模虚拟速度取5mm/ms；板料初始单元大小为32mm，单元细化级数为5，即成形过程最小单元为2mm，单元类型为Belytschlko-Tsay(BT)单元；接触算法选择非线性罚函数法(Accurate)。设置完毕即可提交求解器计算。图2.19行李箱盖外板冲压模拟模型材料参数'DC06'Hardeningcurve:E:210000MPaYieldstress:125.7MPa材料模型正交各向异性Hill48准则r01.86r451.72r902.28硬化准则:等向硬化图2.200.8mmDC06钢板的材料参数(2)成形质量分析重力阶段后板料形状及位移分布如图2.21所示，图中的位移是板料相对初始位置的位移，仅能反映相对移动情况，不能用来描述板料变形。压料阶段完成后板料形状如图2.22所示，板料的塑性应变分布如图2.23所示。对26 2板料成形的数值模拟[63]于没有明显屈服点的材料，通常取应变为0.2%时的应力为屈服应力。图中的塑性应变分布表示，大部分的塑性应变发生在工艺补充面上，不会对产品表面产生影响。图2.21重力阶段位移分布图2.22压料后板料形状图2.23板料塑性应变成形过程可以自定义输出几个冲压成形过程，如图2.24所示。最后得到板件的成形极限图（图2.25）及成形质量分布图（图2.26）。借助成形极限图可以判断任意单元的应变状态。由图2.26可知板料产品部分的应变状态均处于安全区域，图中破裂区域处于工艺补充面上，不影响产品性能，整个零件成形性良好。(a)合模(b)到底前44mm27 西安科技大学硕士学位论文(c)到底前28mm(d)到底图2.24行李箱盖外板成形过程图2.25零件的成形极限图图2.26零件的成形质量(3)模具荷载导出拉延成形过程共设置输出10步成形结果，可以得到每步结束时刻板料对模具的作用力。这个力可分解为沿模具表面的法向和切向的分量，即法向接触力和切向接触力。把模面每个节点每步的荷载连接起来，就是节点冲压过程的荷载曲线。以凹模运行到底，板料成形结束时刻为例，得到模具表面的法向和切向接触力如图2.27、2.28所示。图中箭头的方向代表了节点力的方向，箭头的长短和颜色代表节点力的大小。从模面荷载分布图可以看出，在成形完毕、模具到底时刻，凸模模面上节点的法向接触力最大为254.80N，最小为0；切向接触力最大为19.64N，最小为0。凸模表面受力主要集中在外圆角处。法向接触力和切向接触力在模面外圆角处都较大，而平缓区域（曲率较小部位）和内圆角处受力为零或基本不受力。利用PAMSTAMP菜单中Project->Export->Contour将法向接触力和切向接触力数据输出为ASCII码文件，分别为NCP.asc和TCP.asc，以供结构优化和分析使用。28 2板料成形的数值模拟图2.27到底时模具表面的法向接触力图2.28到底时模具表面的切向接触力2.5.2引擎盖外板冲压成形模拟(1)模型建立同样地，先在HyperMesh中划分好模面网格如图2.29所示，导入到PAMSTAMP2G中。建立冲压模拟模型如图2.30所示。板料为0.8mm的DC06钢板，材料参数见图2.20。屈服准则为Hill48准则；选择等向硬化模型；润滑条件为一般润滑，摩擦因数为0.12；拉延过程压边力为1800kN；凹模虚拟速度取5mm/ms；板料初始单元大小为32mm，单元细化级数为5，即成形过程最小单元为2mm，单元类型为Belytschlko-Tsay(BT)单元；接触算法选择非线性罚函数法(Accurate)。设置完毕即可提交求解器计算。图2.29引擎盖外板的模面网格图2.30引擎盖外板冲压模拟模型(2)成形质量分析重力阶段完成后板料各点的位移见图2.31。29 西安科技大学硕士学位论文图2.31重力阶段板料位移分布压料阶段完成后，板料形状及塑性变形分布如图2.32和2.33所示。由图2.33可知压料阶段，板料没有发生塑性应变。图2.32压料后板料形状图2.33板料塑性应变成形过程设置10个输出步，零件成形过程如图2.34所示。图2.35和2.36分别是成形极限图和成形质量分布图。从它们可以看出零件成形质量良好，工艺可行。(a)合模(b)到底前60mm30 2板料成形的数值模拟(c)到底前20mm(d)到底图2.34引擎盖外板成形过程图2.35引擎盖外板成形极限图图2.36引擎盖外板成形质量分布(3)模具荷载导出成形过程模拟完毕后就可以得到任意输出步模具表面节点的法向和切向接触力的数据。到底时刻，凸模模面上各节点的法向和切向接触力的矢量图如图2.37和图2.38所示。由图可知，模面节点所受到的法向接触力最大为405.91N，最小为0；所受切向接触力最大为8.54N，最小为0。模面主要受力部位在外圆角处，其余部分基本不受力或受力很小。同样地，分别将模面节点受到的法向接触力和切向接触力输出为ASCII码文件，分别为NCP.asc和TCP.asc，供模具优化分析和结构分析使用。31 西安科技大学硕士学位论文图2.37到底时模面的法向接触力图2.38到底时模面的切向接触力2.6小结本章说明了冲压成形模拟的方法，并介绍了如何借助正确的冲压成形模拟得到冲压过程板料对模具表面作用力的方法。为了保证获得荷载的准确性，必须了解板料冲压成形的塑性力学特性和冲压成形模拟中的相关理论。因此本章首先说明了板料成形中描述材料特性的真实应力、真实应变、厚向异性系数以及用来判断板料成形质量的成形极限图。然后阐述了冲压成形模拟中求解算法的选择、材料模型的定义、单元类型的选择以及对接触和摩擦的处理。在进行完整的冲压成形模拟后，就可以得到任一时刻板料对模面的作用力，可分解为切向接触力和法向接触力。用两个实例：行李箱盖外板的冲压成形模拟和引擎盖外板的冲压成形模拟说明了获取模具表面节点力的过程。32 3模具结构拓扑优化3模具结构拓扑优化3.1引言目前模具厂所制造的拉延模具，其内部的承载结构都是“栅格式”结构，如图3.1所示。正如第一章所提到的，传统的模具结构设计遵循既有的规范、准则和经验，缺少对模具受力和变形的关注，但对于高强度钢板冲压成形等新技术其适应性未知。再加上能源和钢铁价格飞涨，从企业降低生产成本的角度来说，也需要优化现有的模具结构。本章将在上文冲压成形模拟的基础上，对传统的模具结构进行优化，探索一种新的模具结构模式。图3.1传统的“栅格式”结构结构优化技术是CAE技术中一个非常重要的组成部分，如今已经被用于产品设计，它以天生的优势正在改变传统的产品设计流程，帮助设计工程师们设计出创新和可靠的产品，这种全新的产品设计过程就是优化驱动的设计过程ODDP(OptimizationDrivenDesignProcess)。结构优化可以根据设计变量的类型分为不同的层次：第一个层次是在给定结构类型、布局和外形几何的条件下，把截面尺寸及材料参数作为设计变量，通常是尺寸优化(SizeOptimization)或参数优化；第二层次是把结构几何作为设计变量，即节点的空间位置可以发生变化，即形状优化(ShapeOptimization)；第三层次是将结构的拓扑关系作为设计变量，通过节点之间是否有连接来优化材料布局，即为拓扑优化33 西安科技大学硕士学位论文(TopologyOptimization)；第四层次表示结构类型也可以作为设计变量。显然，层次越高，优化的难度越大，同时在目标（如重量、刚度等）上的优化效益也越大。拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的设计区域内寻求材料的最优分布问题。结构拓扑优化包括离散结构拓扑优化和连续结构拓扑优化，通常所采取的策略是对于前者一般以结构最小重量为目标，而对于后者往往以结构最小柔顺性为目标，两者均反映结构设计中“传力路径最短”的思想，体现了材料的最合理分配。如果以应力为约束，那么反映了材料的最佳强度分布；如果以位移为约束，则反映了材料的最佳刚度分布。为了说明拓扑优化的基本思想，以一座大桥为例，桥面受均布荷载，下面有两个支座，如图3.2(a)所示。拓扑优化的初步结果见图3.2(b)，图中蓝色区域表示材料密度为零，即可以没有材料的区域。这与实际大桥的材料分布非常吻合。设计区域非设计区域(a)(b)图3.2大桥的拓扑优化(a)荷载(b)拓扑优化结果3.2OptiStruct的结构优化技术[67]结构优化方面已出现了许多成功的商业软件，如德国FE-DESIGN公司的TOSCA，日本Quint公司的OptiShape，美国VR&D公司的Genesis，美国Altair公司的OptiStruct，等，在许多大型通用有限元分析软件如ANSYS、Nastran中都有优化模块。AltairOptiStruct以其成熟、全面的技术得到大家的认可，是当今世界应用最广泛的优化类软件。OptiStruct是以有限元法为基础的最佳优化工具，提供拓扑优化(TopologyOptimization)、形貌优化(TopographyOptimization)、尺寸优化(SizeOptimization)、形状优化(ShapeOptimization)以及自由尺寸优化(Free-sizeOptimization)和自由形状优化(Free-shapeOptimization)。在产品设计的各个阶段灵活应用各种优化技术，加上全球首先引入制造工艺约束，如最大/最小成员尺寸控制、拔模约束、挤压路径约束、对称及模式重复约束，可以驱动设计出优秀的创新产品。OptiStruct优化可以基于线性静态分析、[68]模态分析、线性屈曲分析、频响分析、惯性释放分析以及多体动力学分析。可以将定义多种结构响应作为优化的约束条件或目标函数，如结构柔度、频率、体积/体积分数、34 3模具结构拓扑优化质量/质量分数、位移、应力等，并且与前处理器HyperMesh无缝融合，从而可以快捷地进行问题设置、提交计算和后处理。[69]图3.3是利用OptiStruct优化SUV车架的过程，优化设计与原型相比，减重23%，减少50%焊点数量和焊接长度，大胆采用HSS（高强度钢）和AHSS（超高强度钢），降低了制造成本。OptiStruct中进行拓扑优化设计的主要步骤如下：①导入结构几何和划分结构网格；②施加荷载和约束，创建荷载步；③设置材料参数，单元属性；④模型检查计算；⑤设置优化设计变量，设置优化设计响应，设置优化设计约束，设置优化设计目标函数；⑥模型优化计算；⑦后处理，检查模型刚度、节点位移、应力；⑧导出优化设计得到的几何模型。图3.4说明了这一过程。图3.3某SUV车架的优化设计35 西安科技大学硕士学位论文图3.4OptiStruct拓扑优化流程3.3拓扑优化的数学模型结构拓扑优化问题就是确定材料在设计空间中哪一点是实体材料，哪一点没有材料。对于结构中任意一点，数学上定义如下的特征函数1ifxs()x(3.33)0ifxs式中是给定的设计区域；是实体材料所占的区域；是孔洞所占的区域。ss在有限元求解中这个特征函数就可以表示每个单元的有无。如果以单元的有无作为设计变量，结构拓扑优化问题就是基于0和1离散变量的组合优化问题。其数学模型为TFind=(,,,)12nMinf()f(,,,)12nst..()g0j1,,m(3.34)j()h0k1,,mkhi0,1in1,,式中为设计变量，取1和0，表示单元的有无；f()是目标函数；g()是不等式约i束函数；h()是等式约束函数，n为单元个数。组合优化方法对于小规模问题具有较强的全局寻优能力，同时也不存在其它连续变量拓扑优化中的“奇异最优解”问题。但当优化问题的规模较大时，优化求解效率低且易出现所谓的“组合爆炸”问题。因此，将拓扑变量取为单元尺寸或材料物理参数等连续变36 3模具结构拓扑优化量，使原问题优化模型转化为连续变量问题的优化模型，通过连续变量取下限值，实现结构拓扑变更。这就是基于几何或材料（物理）描述形式的拓扑优化方法，主要有：(1)均匀化方法(HomogenizationMethod)[70]1988年Bendose和Kikuchi提出基于均匀化方法的结构拓扑优化设计。在结构材料中引入微结构（单胞），以微结构单胞尺寸为设计变量，计算其正则密度，建立材料密度与材料特性之间的关系，用孔洞尺寸的消长实现微结构的增删，从而得到最优的拓扑结构。它适用于连续体基于应力、位移约束或频率约束的拓扑优化分析。OptiStruct中持均匀化方法。单元中的微结构如图3.5所示。对于壳单元(Shellelement)有1(1ab)(1)(3.35)对于体单元(Solidelement)有1(1abc)(1)(1)(3.36)式中为单元正则密度；a、b、c为单元内单胞的尺寸。111-b1-c111-a1-a1-b1(a)(b)图3.5单元中的单胞(a)壳单元(b)实体单元(2)变厚度法(Thickness-varyingMethod)变厚度法是基于几何描述方式的优化方法。以基结构中单元厚度为设计变量，将连续体拓扑优化问题转化为广义的尺寸优化问题，通过删除厚度为尺寸下限的单元实现结构拓扑的变更。方法简单，适用于平面结构（如膜、板、壳等），应用到三维问题有一定难度。(3)进化结构优化法(EvolutionaryStructuralOptimization)[71,72]ESO法是澳大利亚华裔学者谢忆民在1993年提出的一种拓扑优化方法。其基本思想是通过将无效的或低效的材料一步步去掉，结构也将逐步趋于优化。在优化迭代中，该方法采用固定的有限元网格，对存在的材料单元，其材料数编号为非零数，而对不存在的材料单元，其材料数编号为零。当计算结构刚度矩阵等特性时，不计材料数编号为零的单元特性（通过数据映射转换，建立固定有限元网格数据信息和计算结构刚度矩阵等特性所需的有效网格数据信息关系），通过这种零和非零模式实现结构拓扑优化。37 西安科技大学硕士学位论文目前ESO法已被成功应用于热应力优化、热传导优化、触应力优化、可变荷载结构优化等。(4)密度法(DensityMethod)密度法属材料（物理）描述方式，其基本思想是人为地引入一种假想的密度可变的材料，假定材料物理参数（如许用应力，弹性模量）与材料密度间存在某种函数关系。优化时以每个单元的密度为拓扑设计变量，设计变量的个数与单元个数相当。密度法不仅可以采用结构的柔顺度为优化的目标函数，也可以用于特征值优化、柔性机构的优化、多学科优化等领域。OptiStruct中支持密度法，这也是本文将要采用的方法。一般情况下，密度法拓扑优化的数学模型可以表示如下TFind=(,,,)12nTMinC()FUst..FKU(3.37)VfV或MfM0001(in1,...,)mini式中为设计变量，是材料相对密度；n为设计空间内单元的个数；C()目标函数，i通常取系统柔度，也可取应变能、特征值、频率等；F为外力向量；U为结构的位移向量；K为系统总刚矩阵；f为体积（质量）约束分数，即优化后模型体积（质量）与初始结构体积（质量）的比值；V优化后结构体积；V初始结构的体积；M优化后结构0质量；M初始结构的质量；为最小相对密度，为了避免总刚度矩阵奇异，取一个0min较小的值。密度法中所指的密度是反映材料密度和材料特性之间对应关系的一种伪密度。常见的材料插值模型：固体各向同性惩罚微结构模型(SolidIsotropicMicrostructureswithPenalization，SIMP)、材料属性的有理近似模型(RationalApproximationofMaterialProperties，RAMP)。这两种方法都是假定材料弹性模量与单元相对密度之间的函数关系，通过引入惩罚因子对中间密度值进行惩罚，使中间密度值向0-1两端聚集，使连续变量的拓扑优化模型能很好地逼近0-1离散变量的优化模型，这时中间密度单元对应一个很小的弹性模量，对结构刚度矩阵的影响将变得很小。SIMP模型材料弹性模量与相对密度之间的关系为pEE(EE-)(01)(3.38)imini0minmini式中E为第i个单元的弹性模量；为第i个单元的相对密度；p为罚因子；E为单元ii0充满材料时的弹性模量；为数值求解的稳定性，通常取EE/1000；罚因子p的取min0值大小决定了抑制中间密度单元的效果，不妨取E1，插值模型如图3.6所示。038 3模具结构拓扑优化图3.6SIMP材料插值模型(p=1,2,3,4,10)RAMP模型材料弹性模量与相对密度之间的关系为iEE(3.39)i01p(1)i式中E为第i个单元的弹性模量；为第i个单元的相对密度；p为罚因子；E单元充ii0满材料时的弹性模量；为数值求解稳定，通常取EE/1000；同样地，取E1，插min00值模型如图3.7所示。图3.7RAMP材料插值模型(p=1,3,5,10,20)从图3.6和图3.7的可以看出，SIMP法和RAMP法具有很大的相似性。在相等的39 西安科技大学硕士学位论文罚因子情况下，SIMP方法对中间密度的惩罚比RAMP方法要缓和。OptiStruct中采用SIMP方法对材料参数进行插值，一般情况下取p=2~4。应用SIMP材料插值方法，取结构柔度为目标函数，单元密度为设计变量，体积分析为约束条件，拓扑优化的数学模型为TFind=(,,,)12nnTTpTMinC()FUUKUiukuiiii1st..FKU(3.40)VfV001(in1,...,)mini3.4优化模型的数值算法结构优化设计是一个不断修正设计变量，而得到系统最优解或局部最优解的过程。[68]OptiStruct中的迭代流程如图3.8所示。这其中涉及到收敛判定、灵敏度分析、优化模型求解等问题。图3.8OptiStruct优化迭代流程3.4.1收敛准则OptiStruct中使用两种收敛准则，即规则收敛和软收敛。当相邻两次迭代结果的目标函数的改变小于容差，并且约束条件的违反率小于1%时达到规则收敛。至少要进行3次物理分析，因为规则收敛准则是基于结构在最新的设计变量下的目标函数值。默认目标函数容差为0.005，可以通过设置OBJTOL来更改。40 3模具结构拓扑优化当相邻两次迭代的设计变量没有或只有很小的改变时，认为达到软收敛。不必要对最后一次迭代的目标函数值或约束函数值进行估值，因为模型没有发生显著改变。因此，软收敛比规则收敛少一次物理分析过程。3.4.2灵敏度分析结构灵敏度反映的是设计变量或者设计参数的改变对目标函数的影响，一般结构响应g可以从位移矢量u求得Tgqu(3.41)式中q为系数向量。结构响应的灵敏度就是对设计变量的偏导，或者响应的梯度TgqTuuq(3.42)有两种方法进行灵敏度分析：直接法和伴随变量法。对于一般的有限元方程KuF(3.43)对设计变量求偏导KuFuK(3.44)从而可以求出位移矢量u的灵敏度uFKKu(3.45)再反代入(3.42)式，即可求得得到结构响应的灵敏度。这就是直接法。在伴随变量法中，引入伴随变量aKa=q(3.46)则响应的灵敏度为TgqTFKuau(3.47)直接法适合于约束个数较多而设计变量较少的情况，如尺寸优化和形状优化。伴随变量法用于设计变量很多而约束较少的情况，如拓扑优化问题。3.4.3优化准则算法对于(3.40)式建立的优化模型，比较常用的方法有优化准则法(OptimalityCriteria，OC)和移动渐近线法(MethodofMovingAsymptotes，MMA)。在OptiStruct中，对于典型的拓扑优化问题使用优化准则法。41 西安科技大学硕士学位论文利用拉格朗日乘子法将(3.40)式中基于SIMP的约束最优化问题转化为无约束最优化问题进行求解。拉格朗日函数为LC(VfV)(KUF-)()(1)(3.48)1023min4式中：1是标量，2、3、4是向量。根据Kuhn-Tucker条件，拉格朗日函数在****(,,,)处取极值的必要条件12nLCV()KU0**1*2*3ii4iiiiVfV0FKU*(3.49)()03min*(1)040,0in1,,34ii*1min其中，3是由3i组成的列向量，4是由4i组成的列向量。*当min1时，设计变量的上下限约束均不起作用，有34ii0,0；*当imin时，设计变量的下限约束起作用，有34ii0,0；*当i1时，设计变量的上限约束起作用，有34ii0,0。故Kuhn-Tucker条件等价于下式LCV()KU*0if1**1*2*miniiiiiLCV()KU*0if1**12**iiiiiLCV()KU*(3.50)0if**1*2*iminiiiiVfV0FKU0,034iinTTVV*由CFUUKU，ii，式(3.50)中的第一式变为i1TUTTKUKUKUUUUKVUK0(3.51)***1i2*2*iiiii有平衡方程KUF(3.52)T将(3.52)两边转置，注意到刚度矩阵K为对称阵：KK，则有42 3模具结构拓扑优化TTUKF(3.53)*将式(3.52)和式(3.53)两边对i求偏导可得TKUTKUUK或UK(3.54)****iiii代入(3.51)式可得TTKU(2U)U(U)KV0(3.55)2**21iiip*TSIMP插值模型，可知：KKii0，取2U，代入上式得*1pTp()ukuv0(3.56)ii01ii*1pTkp()iukui0i令Di，优化准则法的相对密度的迭代公式如下v1iDkkifDkk1iiminiik1kkiminifDiimin(3.57)kk1ifDii1为阻尼因子，引入阻尼因子是为了保证数值计算的稳定性和收敛性。3.5模具结构的拓扑优化3.5.1模具荷载的转换前一章已经通过冲压成形模拟输出了覆盖件模具在冲压过程中某时刻板料对模具的作用力，即模具表面节点所受的法向接触力(NormalContactPressure)文件NCP.asc和切向接触力(TangentContactPressure)文件TCP.asc。现在只要将这些节点力映射到模具表面节点上就可以在拓扑优化和结构分析中使用了。由PAM-STAMP2G冲压模拟得到的是模面上所有节点处的法向和切向接触力，荷载数量是节点数的两倍，并且每个力的方向大小都不相同，显然不能通过手动的方法加载到有限元模型上。必须利用计算机自动实现，这就要把PAM-STAMP2G的荷载格式转化为优化软件OptiStruct能够识别的格式。需要注意，在PAM-STAMP2G中使用的是mm-kg-ms单位系统，得到力的单位为kN；而在OptiStruct中使用的是mm-t-s的单位系统，得到力的单位为N，二者之间相差1000倍。转换过程通过作者编写的FORTRAN程序Pam2OptiStruct.exe实现。程序运行时需要输入由冲压模拟软件PAM-STAMP得到的法向和切向接触力文件（NCP.asc和TCP.asc）的文件名以及前处理软件HyperMesh43 西安科技大学硕士学位论文输出的前处理文件的文件名，程序就可以将法向和切向接触力文件转换为OptiStruct可识别的荷载文件：NCP.INC和TCP.INC，并且插入到结构优化和分析的前处理文件中，这就相当于把模面节点所受的法向接触力和切向接触力施加到了结构优化和分析的模型上。图3.9是Pam2OptiStruct.exe运行时的界面。程序源代码见附录。图3.9Pam2OptiStruct.exe运行界面3.5.2荷载工况选择冲压过程中，模具受力是动态荷载，大小和方向都随时间变化。理论上来说，可以将冲压输出的步数设置较大，得到模具在许多个时间点的荷载。然后将每个节点在这些时间点上的数据连接起来形成该节点在冲压过程中的荷载曲线。最后将所有节点的荷载曲线输入有限元模型进行基于动力学分析的拓扑优化。但这样做既很复杂也无必要。冲压模拟文件的体积与输出的步数成正比，将冲压模拟的输出步数设大，将大大提高冲压模拟文件的体积。使得一个冲压模拟的文件达到5~10Gb，对计算机的存储空间是一个挑战。模具表面节点往往在上万个，也就是有数万条荷载曲线要载入到有限元模型上，这在目前一般的计算机上是无法进行运算求解的。拓扑优化反映的是结构的最佳传力路径。对于两种仅在荷载大小上有差异的工况，得到的结构形式是相同的。以第二章中行李箱盖外板的冲压过程为例，选择模面上10个关键节点（分布见图3.10），将它们的荷载曲线画出。图3.11是1至5号节点荷载曲线，图3.12是6至10号节点的荷载曲线。图中横轴表示凹模行程，行程等于0时Stamping阶段开始，行程等于71mm时凹模运行到底，成形完毕。44 3模具结构拓扑优化图3.10模面上选择的关键节点图3.111~5号节点的荷载曲线图3.126~10号节点的荷载曲线从图中可以看出，凸模外圆角处受力逐渐增大，到最后阶段上升到最大。1~2号节点受力从冲压开始就较大，并且在初始值附近波动；3~10号节点在冲压开始基本不受力，直到成形进行到三分之二以后才有接触力。其中3~8号节点在成形最后一步，节点力突然上升，在成形结束时刻达到最大值。这说明在成形最后一步节点力最大，模具的受力状况最为恶劣。因此可以选择成形结束时刻凸模表面荷载施加到凸模的有限元模型上，以此进行结构拓扑优化。设计出新的模具结构以后，再用冲压过程中的模具荷载对模具结构受力和变形进行验证计算。3.5.3目标函数和约束条件OptiStruct中可以定义多达26的结构内部响应（可由程序自动计算得出），如结构质量和体积、质量分数和体积分数、惯性矩、结构柔度、频率、VonMises应力、位移、屈曲因子、频率响应速度等。用户还可以用Fortran或C语言编写外部函数提取响应自定义为外部响应。这些响应有大部分都可作为目标函数和约束条件。本文的拓扑优化不涉及频率、屈曲、多荷载步等因素，而是在保证结构刚度的情况下，尽可能减小模具质量。由于离散系统的结构刚度不易比较，故以结构柔度为目标函45 西安科技大学硕士学位论文数，而以质量分数为约束条件。结构柔度C通过下式计算11TTCuKudV(3.58)22式中K为系统刚度矩阵；u为节点位移矢量；应变；为应力；V为设计材料的总体积；结构柔度是结构的应变能，可以认为是结构刚度的倒数。在一定的荷载下，结构柔度越小表示刚度越大。质量分数是优化后结构总质量与初始设计总质量的比。选择质量分数作为约束条件可以限制整个结构的质量。普通低碳钢模具的受力状况并不恶劣，其应力水平也不高，远没有达到材料的强度（这一点从后面的计算也可看出），故没有选用结构的应力作为约束条件。3.6模具拓扑优化实例3.6.1行李箱盖外板模具的拓扑优化(1)优化模型设置冲压模拟和结构优化使用同一模面网格，首先需要在上一章模面网格的基础上划分体网格。在模具设计标准中，模具顶部厚度一般为60mm。现将模面向下偏移50mm，这一部分为不可设计区域（图3.13中红色部分），其余部分为设计区域（图3.13绿色部分）。全部采用四面体单元，共有1317874个。凸模材料使用灰铸铁HT300。材料参数：33弹性模量105GPa，泊松比0.26，密度7.25×10kg/m。图3.13模具优化设计模型应用Pam2OptiStruct.exe将冲压模拟的荷载插入到模型文件中，然后重新用OptiStruct读入，荷载就已经施加到了模型上，如图3.14所示，图中箭头的方向代表了节点力的方向，箭头的长短代表节点力的大小。除此之外，还受到重力作用，在OptiStruct添加GRAV卡片实现。最后将模具底面节点全部约束。至此模具荷载施加完毕。设置优化参数如下：46 3模具结构拓扑优化目标函数：结构柔度最小化；约束条件：质量分数（设计后结构与初始质量的比）小于给定值(25%、30%、35%)；模面上特征点的位移小于允许值(0.1mm)；设计变量：设计区域单元的密度；模面上选取530个节点为位移约束的特征点，分布如图3.15所示。图3.14模具表面荷载图3.15模面位移约束特征点(2)结构设计对于质量分数约束，分别使用三种不同的水平(25%、30%、35%)进行优化设计，计算后得到单元的密度分布。密度为0表示该处材料为空，密度为1表示该处材料为实体。通常取0.3为密度阀值，在此基础上进行进一步的设计和验证。三种质量分数水平下模具结构优化结果如图3.16所示。从图中可以看出各质量分数下的优化结果基本相同。利用OptiStruct的OSSmooth功能，将优化结果导入CAD软件CATIA中，根据优化结果设计出新的模具结构见图3.17，传统的模具结构如图3.18所示。优化结果与传统的模具结构（图3.18）有明显不同。在模面受力较大的部位下有材料起支撑作用。明显地，二次设计出来的模具结构云去除了传统结构中的栅格，只在模具内部必要的地方利用材料支撑，在保证结构变形和刚度的前提下，甚至可以在模具侧面上适当开孔以降低模具质量。(a)47 西安科技大学硕士学位论文(b)(c)图3.16拓扑优化后密度大于0.3的单元分布(a)质量分数为25%(b)质量分数为30%(c)质量分数为35%图3.17二次设计的模具结构48 3模具结构拓扑优化图3.18传统的模具结构3.6.2引擎盖外板模具的拓扑优化(1)优化模型设置在HyperMesh中建立凸模的实体模型如图3.19所示。模型采用四面体网格和六面体网格混合，共有四面体单元181929个，六面体单元173040个。模面向下偏移50mm以上包括模面部分定义为不可设计区域（图中红色部分单元），将其余部分定义为设计区域（图中绿色部分单元）。凸模材料使用灰铸铁HT300。材料参数：弹性模量105GPa，33泊松比0.26，密度7.25×10kg/m。图3.19优化设计模型应用Pam2OptiStruct.exe施加模具表面法向和切向荷载如图3.20所示，图中箭头的方向代表了节点力的方向，箭头的长短代表节点力的大小。在OptiStruct添加GRAV卡片，约束模具底面节点。至此，模具荷载施加完毕。设置优化参数如下：目标函数：结构柔度最小化；约束条件：质量分数（设计后结构与初始质量的比）小于给定值(25%、30%、35%、40%)；模面上特征点的位移小于允许值(0.1mm)；49 西安科技大学硕士学位论文设计变量：设计区域单元的密度；模面上选取694个节点作为位移约束的特征点，分布如图3.21所示。图3.20模面荷载图3.21模面位移约束的特征点(2)结构设计同样地，四种不同质量分数约束水平下的优化结果如图3.22所示，利用OptiStruct的OSSmooth功能，将优化结果导入CAD软件CATIA中，根据优化结果设计出新的模具结构见图3.23，传统的模具结如图3.24所示。对比图3.23和图3.24可以明显看出，二次设计的模具结构去掉了中间的支撑结构，只在前后边缘处留有材料来支撑整个结构。不再是传统模具中的栅格式结构。50 3模具结构拓扑优化(a)(b)(c)(d)图3.22拓扑优化后密度大于0.3的单元分布(a)质量分数为25%(b)质量分数为30%(c)质量分数为35%(d)质量分数为40%51 西安科技大学硕士学位论文图3.23二次设计的模具结构图3.24传统的模具结构3.7小结拓扑优化技术如今已经广泛地应用到众多产品的设计研发当中，并且逐渐形成了以优化驱动的设计过程ODDP。本章首先介绍了OptiStruct的优化技术以及连续体拓扑优化的密度法的原理、思想和数学模型。最后提出了覆盖件模具的结构拓扑优化方法，设计出新型的模具结构。以两个实例：行李箱盖外板模具和引擎盖外板模具详细说明了这一方法，根据优化结果设计出了新型的模具结构。从结果可以看出，优化后的模型在结构上与传统结构有明显不同，去除了传统结构中栅格式承载部分，只在重要部分利用材料加以支撑。52 4模具结构分析4模具结构分析4.1引言模具结构分析对模拟获得的应力、变形准确度要求高，这决定了隐式算法更适合用于模具结构分析。如果采用动力显式算法，模具各部件本身质量相对较大，动力显式算法的虚拟速度使系统虚拟动能增加明显，模拟结果偏差较大，而且动力显式算法无须进行迭代，每步计算只能近似满足，大量的计算时间步导致累积误差较大，引起数值模拟结果的飘移，很难得到精确的模具应力和变形。本文的模具结构分析采用Altair公司HyperWorks中的RADIOSS求解器。RADIOSS是下一代的隐式和显式有限元求解器，可以进行线性静力学和动力学分析，也可用于复杂的非线性瞬态动力学和多体动力学分析。这一强大的多学科解决方案可以最大化产品的耐久性、NVH性能、碰撞性能、安全性能及可制造性，从而使得产品制造商们将创新的产品更快地推向市场。它与前处理器HyperMesh和结构优化软件OptiStruct无缝接合，为分析过程提供了很大的便利。上一章根据结构拓扑优化，已经重新设计出了一个新型的模具结构。然而，这个结构是否安全？是否优于原来的结构呢？这就需要对优化前后的模具结构进行分析，以比较二者的优劣。下面仍然以行李箱盖外板和引擎盖外板为例，说明分析方法。冲压过程是一个动力学过程，冲压过程中模具的受力和变形情况怎样？如果直接进行动力学分析，仍然会遇到1.3节中所指出的问题。同样利用解耦的思想，在结构分析的基础上，本章将提出冲压过程中模具结构分析方法。4.2行李箱盖外板模具结构的分析4.2.1传统模具结构分析与拓扑优化相同，在结构分析时模具有限元模型也使用与冲压模拟分析中相同的模面网格，这样就可以直接将冲压模拟得到的模具表面节点力通过3.5.1节中介绍的荷载转换程序Pam2OptiStruct.exe直接插入到有限元模型文件中。(1)模型建立行李箱盖外板的传统模具结构如图4.1所示。在HyperMesh划分好的模面网格（图2.18）基础上划分模具体网格如图4.2所示。由于模面采用三角形单元，为了避免用金字塔单元在四面体单元和六面体单元间的过渡，体网格也采用四面体单元。共有1752597个四面体单元。53 西安科技大学硕士学位论文划分好网格后，将模型导出为model.fem文件，利用Pam2OptiStruct.exe将冲压模拟得到凹模到底时刻的法向接触力文件和切向接触力文件插入模型文件中，重新读入模型，节点力就施加到模面节点上了。如图4.3所示，图中箭头的方向代表了节点力的方向，箭头的长短代表节点力的大小，模具受力主要集中在模面外圆角处。除模面节点力外，模具还受重力作用，底面节点应该全部约束。材料使用HT300。材料参数：弹性模33量105GPa，泊松比0.26，密度7.25×10kg/m。图4.1传统模具的几何模型图4.2传统模具的体网格图4.3模具表面荷载54 4模具结构分析(2)结果分析求解后得到节点的位移分布和VonMises应力分布如图4.4和图4.5所示（选取模具的顶视图、左视图、正视图和中间截面来观察，下同）。从位移分布图可以看出，最大-3位移为5.35×10mm。从模面向底面逐渐减小，位移较大的部位主要在接近模面中部节点力较大的区域以及模面两侧凸起的工艺补充面部分。VonMises等效应力云图显示，结构的最大等效应力为13.58MPa。等效应力较大的部位分布在节点力较大的部分，模具大部分应力水平较低，这也证明了在生产普通冲压件时模具的受力状况并不恶劣。图4.4传统模具的位移分布图4.5传统模具的VonMises应力分布55 西安科技大学硕士学位论文4.2.2优化后模具结构分析第三章中已经设定不同的质量分数约束得到了不同的优化结果，并且在此基础上设计出了新的结构形式。同样地，对设计出的新型模具（图3.17）进行结构分析如下：(1)模型建立用四面体单元划分结构网格如图4.6所示。共有1165352个四面体单元。用节点力转换程序Pam2OptiStruct.exe在模面上施加与图4.3相同的节点力。此外，还要加上重力，约束模具底面节点的自由度。材料使用HT300。材料参数：弹性模量105GPa，泊33松比0.26，密度7.25×10kg/m。有限元分析模型建立完毕。(2)结果分析得到二次设计的模具节点位移分布和VonMises应力分布如图4.7和图4.8所示。从-3位移分布（图4.7）可以看出，最大位移为5.41×10mm。从模面向底面逐渐减小。与传-3统结构的位移分布（图4.4）相比，模面上部位移分布较均匀，都在3.25~5.41×10mm之间，直至底面逐渐减小为0。这说明结构大部分材料都参与了变形，但变形在允许的范围内，充分利用了材料的承载性能。从VonMises等效应力分布可以看出，最大等效应力为12.55MPa，仍然在模面节点力较大的位置，与传统结构（图4.5）相比在同一水平上，略有降低。图4.6模具结构体单元56 4模具结构分析图4.7二次设计结构位移分布图4.8二次设计结构VonMises应力分布4.2.3对比分析结构优化是一个不断修改验证的过程，对每次优化的结果，都要与原结构进行分析对比，以确定修改的方向。第三章中按照不同的质量分数对行李箱盖外板的模具进行了优化，分别是25%、30%、35%，加上传统模具结构和二次设计的结构，一共有5种结构。对它们进行结构分析，并将结果汇总于表4.1。57 西安科技大学硕士学位论文表4.1优化设计结构和传统模具结构分析结果质量分数25%30%34.6%35%40.2%设计结构传统结构结构总质量(t)1.561.872.152.172.50结构柔度(N∙mm)297.54227.50205.83186.56212.95模面最大位移(×10-3mm)6.89e-35.31e-35.41e-34.15e-35.34e-3最大VonMises应力(MPa)16.3016.3012.5016.2013.60由表4.1可以看出，质量分数为25%和35%两种情况下的结构柔度比传统结构大，并且节点的最大位移和最大VonMises应力也高于传统结构。而二次设计的结构质量分数为34.6%，总质量为2.15t，相比传统结构的2.50t减少14%，与此同时柔度降低了3.3%，说明结构刚度增大。模面的最大位移和最大等效应力几乎没有变化。因此，二次设计的模具结构优于传统结构。4.3引擎盖外板模具结构分析4.3.1传统模具结构分析引擎盖外板的传统模具结构见图4.9。在模面网格的基础上划分结构体网格，全部采用四面体单元，共有402579个。与行李箱盖外板相同，利用Pam2OptiStruct.exe将模面节点力（同图3.19）施加到有限元模型上，同时施加重力并且约束模具底面节点，结构分析的模型就建立完毕。节点位移和VonMises等效应力分布如图4.10和图4.11所示。-3节点的最大位移为3.10×10mm。在模面前端和后半部节点力较大的部分位移较大。结构最大等效应力为22.69MPa，同样在模面前端节点力较大部位，其余部分应力很小。图4.9传统模具结构58 4模具结构分析图4.10模具的位移分布图4.11模具的VonMises应力分布4.3.2优化后的模具结构对图3.22中设计出的新型模具结构，进行结构分析。使用与冲压模拟中相同的模面网格，采用四面体单元对结构进行离散，如图4.12所示。共有567673个四面体单元。用节点力转换程序Pam2OptiStruct.exe在模面上施加与图3.19相同的节点力。此外，还要加上重力，约束模具底面节点的自由度。材料使用HT300。材料参数：弹性模量33105GPa，泊松比0.26，密度7.25×10kg/m。有限元分析模型建立完毕。利用RADIOSS求解计算后得到模具节点的位移和VonMises应力分布如图4.13和59 西安科技大学硕士学位论文-34.14所示。节点最大位移为3.13×10mm，与传统结构（图4.10）相比没有明显增加，但整个结构中参加变形有材料增多，更充分地利用了材料。从结构的VonMises等效应力分布可以看出，最大等效应力为22.67MPa，与传统模具相比没有太大差别。图4.12模具体单元划分图4.13模具节点位移分布60 4模具结构分析图4.14模具VonMises应力分布4.3.3对比分析上一章中按照不同的质量分数对行李箱盖外板的模具进行了优化，分别是25%、30%、35%、40%，加上传统模具结构和二次设计的结构，一共有6种结构。对它们进行结构分析并将结果汇总于表4.2。表4.2优化设计结构和传统模具结构分析结果质量分数25%30%30.1%35%39.5%40%设计结构传统结构结构总质量(t)1.621.941.992.272.562.59结构柔度(N∙mm)129.3987.5972.5570.7078.2948.14-3模面最大位移(×10mm)7.684.053.133.023.103.07最大VonMises应力(MPa)11.2011.2011.7011.2013.1011.30由表4.2可以看出，质量分数为25%和30%两种情况下的结构柔度比传统结构大，可认为刚度减小了，并且节点的最大位移和最大VonMises应力也高于传统结构，不能用来设计新的结构。随着质量分数的增加，结构柔度越来越小。这表示增大质量在一定程度上可以加大模具结构的刚度。但并不意味着质量越大，刚度一定增大。二次设计的结构质量分数为34.6%，总质量为1.99t，相对于传统结构的2.56t减少22%，与此同时柔度降低了7.3%，结构刚度反而增大，模面的最大位移和最大等效应力几乎没有变化。因此，二次设计的模具结构在保证模具变形和受力在可接受的情况下，使用较少的材料，设计出了具有更大刚度的结构。显然二次设计的模具结构优于传统结构。61 西安科技大学硕士学位论文4.4冲压过程的模具结构分析冲压过程中，模具的受力和变形是使用新材料、新技术的条件下关注的重点之一。在实际板料成形中，冲头速度一般0.1~1m/s，属于准静态的连续介质力学问题。因此，可以将一个冲压过程分为n步，提取第n步最后时刻的模面荷载，分别当作一个工况施加到模具结构上进行结构分析，然后将分析结果连接起来就可以得到整个冲压过程中模具的受力和变形情况。这就是基于板料成形模拟的模具结构分析方法，用流程图表示如图4.15所示。图4.15基于板料成形模拟的模具结构分析方法主要分为五个部分：(1)模具几何信息的提取。包括模具型面、模具三维实体数据。主要通过CAD软件完成，如CATIA、UG提取出模面几何为.igs文件，三维模型为.stp文件，供前处理程序划分网格使用。(2)冲压模拟和有限元分析的前处理。这里最主要的是划分网格，利用著名的专业前处理软件HyperMesh进行。划分好模面网格传递给冲压模拟程序PAMSTAMP，同时也在面网格的基础上划分体网格，以供结构分析使用。62 4模具结构分析(3)冲压过程模拟。主要是通过模拟得到冲压过程中各个关键时刻板料对模具的作用力，导出设定时刻的节点力数据文件，利用PAMSTAMP2G完成。利用HyperMesh划分的模面网格，进一步地设置好冲压模拟的模型，并且要保证模拟结果真实、成形效果良好，这样才能得到准确的模面节点力数据。(4)节点力荷载的施加。将冲压模拟得到的节点力数据分别施加到第2部分划分好的模具三维模型上，通过FORTRAN程序Pam2OptiStruct.exe修改有限元模型的数据文件实现。(5)有限元求解和分析。对于每个时刻施加了节点力荷载的模型继续设置，加上重力，约束必要的节点位移，对每个时刻的模型分别求解。然后将所有的求解结果连成一个连续过程进行分析。这部分通过前处理器HyperMesh、求解器RADIOSS和后处理软件HyperView实现。按照这种思想，将行李箱盖外板的冲压过程分为10步，输出每步结束时板料对模面的作用力，分别将这些作用力施加到模具有限元模型上进行结构分析，即可得到每一步最后时刻模具的受力和变形情况。图4.16和图4.17分别是冲压过程中某几个时刻的模具节点位移分布和VonMises等效应力分布。(a)Progression=0(b)Progression=36mm63 西安科技大学硕士学位论文(c)Progression=49mm(d)Progression=64mm(e)Progression=71mm图4.16冲压过程中模具的节点位移分布图4.16中是以凹模的行程(Progression)来划分的，Progression等于0代表压边过程结束，凹模、压边圈夹紧板料向下运动的开始时刻；Progression等于71mm表示凹模和板料完全到底，成形过程结束时刻。从节点位移分布图可以看出，拉延开始时节点位移以模具最高点为中心向四周和模具底面扩散（图4.16a）；随着拉延深度的扩大，模具表面上凸起的工艺补充面和外圆角部分开始与板料接触，这些地方的位移也逐渐增大（图4.16b、c）；随后凸模模面大部分与板料接触（图4.16d），此时模面上外圆角交汇处的节点位移明显增大（见图4.16d中红色圆圈处），比最开始接触的模具顶点区域还要大，并且一直持续到拉延过程结束（图4.16e），最终模面上两侧凸起的工艺补充面部分和外圆角部分的节点位移处在一个较高的水平。这与拉延模具凸模外圆角处受力通常最恶劣的经验一致。64 4模具结构分析(a)Progression=0(b)Progression=49mm(c)Progression=71mm图4.17冲压过程中模具的VonMises应力分布从VonMises等效应力分布图4.17可以看出，拉延开始时模具与板料接触部位应力比其它部位大，随着拉延的进行，模面上外圆角部分等效应力逐渐增大超过了最先与板料接触的模面顶端，与冲压过程中模具的节点位移分布规律相似。将冲压过程中凸模上最大节点位移和最大VonMises应力随凹模行程的变化曲线画出来如图4.18所示。图中表明凸模上最大变形和最大等效应力出现在最后一个时间步。最大等效应力在凹模行程的前14mm内较小，此时板料与凸模接触面积较小；凹模行程14~28mm阶段，等效应力加速增加，此时板料与凸模接触面积不断扩大，作用力增大；65 西安科技大学硕士学位论文凹模行程28~56mm阶段，等效应力呈波浪状起伏，略有上升，此时板料与凸模基本全部接触，处于拉延变形时期；凹模行程56~71mm阶段，等效应力迅速上升，这主要是因为板料大部分参与塑性变形，正在形成零件的关键特征，再加上凹模到底有一定的冲击作用。(a)最大节点位移变化曲线(b)最大VonMises应力变化曲线图4.18冲压过程凸模上最大节点位移和最大VonMises应力变化4.5小结优化设计结果的验证和冲压过程模具变形受力分析都需要结构分析。本章首先对第三章中通过拓扑优化设计出的模具结构和传统模具结构在成形最后时刻的情况分别进行了结构分析。通过比较发现，根据优化结果设计出的新型模具结构在提高模具整体刚度的同时，能显著降低模具质量。对于行李箱盖外板和引擎盖外板模具分别降低了14%和22%，并且模具的最大节点位移和最大VonMises等效应力与传统模具保持在同一水平上。这说明本文提出的模具结构拓扑优化方法可行。对于冲压过程中的模具受力和变形情况，提出了冲压过程中的模具结构分析方法，可以得到冲压过程中模具的受力和变形规律，并给出的详细步骤。以一个实例说明了冲压过程中的模具结构分析方法，分析了冲压过程中模具的变形和应力变化规律。结果说明在拉延开始阶段，首先与板料接触部位的节点位移和VonMises应力较大，随着拉延的进行，模面外圆角部位及其交汇处的节点位移和等效应力都显著增大，一直持续到拉延结束。66 5结论5结论我国汽车和模具工业的发展方兴未艾，必然要求模具技术不断发展和进步。冲压模具是汽车制造成本的重要部分。一个车型的模具总费用动辄上亿元，随着能源和钢材价格的不断上升，减少模具成本对控制整车成本有重要意义；另一方面传统的模具结构是一成不变的“栅格式”结构，随着高强度钢板等新材料、伺服压力机等新的冲压设备应用于汽车覆盖件的生产，对传统的覆盖件模具提出了新的问题和挑战。基于这两点，本文提出了基于板料成形数值模拟的冲压模具结构优化方法和模具结构分析方法。主要结论如下：(1)利用冲压模拟软件得到冲压过程中板料对模具的作用力，然后再将这些荷载施加到模具有限元模型上进行结构优化和分析。在冲压模拟中把模具看作刚性体，在结构优化和分析中把模具看作弹性体，避免了这两个过程的耦合求解。该方法符合当前的计算机软硬件水平，切实可行。(2)结构拓扑优化是冲压模具结构优化设计的有效方法。本文提出的基于板料成形数值模拟的模具结构优化方法能够设计出新型模具结构，新结构相比传统模具结构能够减重10%~20%，降低柔度5%左右，同时结构的最大变形和最大VonMises应力保持在同一水平。优化后的结构优于原结构。(3)基于板料成形数值模拟的模具结构分析方法把成形过程与结构分析过程解耦合，将连续的冲压成形过程用n个时间点的静态结构分析序列来近似，也可以得到模具结构应力和变形的规律，降低了分析成本和难度。(4)本文的模具结构优化和分析方法中，在划分好模面网格以后，板料成形模拟和模具三维有限元模型的建立可以同步进行，大大节省了分析时间。(5)本文的冲压模具结构优化和分析方法可以推广到其它条件（如高强度钢板模具、合金材料）和锻造、挤压模具的结构优化和分析中。本文从覆盖件模具着手，虽然经过长时间的努力，对传统冲压模具结构优化和分析进行了部分研究，但由于时间和精力的限制，在一些方面还有待提高和完善。在后续研究中还应进一步完成如下工作：(1)文中的模具结构优化和分析方法步骤较多，需要人工干预的步骤多，可以将荷载提取和转换作为一个单独的过程嵌入到前处理软件HyperMesh中，变成程序自动操作的过程。这样既可以减少劳动量，也可降低出错的几率。(2)优化设计的新结构需要在实践中不断检验与磨合才能朝着工艺化、实用化的方向前进，这方面需要资金和设备的支持。67 西安科技大学硕士学位论文致谢论文完稿之际，首先向悉心指导我的导师郭志勇教授和娄茂昆老师表示衷心的感谢！郭志勇教授为人豁达，睿智博学，治学严谨。从第一次聆听郭老师授课到现在已将近七年，郭老师深厚的学术造诣，执着的敬业精神，一切为学生着想的责任心以及朴实的生活作风深深地印在我的心中。郭老师不仅传授我许多专业知识，更加教导我如何做人，我将永远铭记郭老师“勤奋，敬业，不计较”的人生教诲。在奇瑞学习期间，娄茂昆老师在生活和学习上都给我了极大的支持、鼓励和关怀。从论文的选题、资料收集、分析研究直至成稿，都倾注了娄老师的心血。娄老师知识渊博、行事果断，具有国际化的视野和前瞻性的思想，尤其是他坚定的人生追求和永不停息的奋斗精神给我留下了深刻的印象，也时刻鞭策着我在生活中不断前进，这将使我终生受益。特别感谢能源学院的来兴平教授在生活和学习上对我的关心、帮助、鼓励和指导！来老师学术成就斐然，能以严谨的治学态度和极具人格魅力的作风影响学生，是学生人生道路上的灯塔。在此衷心祝愿老师们身体健康！工作顺利！阖家幸福！感谢理学院力学系的各位老师七年来对我的鼓励与帮助！感谢奇瑞汽车工程研究院的吴沈荣博士，感谢奇瑞汽车规划设计院同步工程科的工程师：孟详新、武建敏、任承峰、魏风波、谢伟、张婧、袁剑湘、贡建秀、刘小专、沈左军在论文研究期间给予的帮助！感谢师兄曹小平、朱向会，师妹高娟，同学付团伟、张莹博、李立波、曹建涛，感谢2007级研究生2班的同学，衷心感谢他们在生活上的关心和学习过程中的探讨！特别感谢我的父母和弟弟妹妹，他们的关心和鼓励是我最大的动力，他们克服重重困难而给予我的无私奉献是我最强大的后盾！最后，谨向百忙之中抽出宝贵时间评审论文和参加论文答辩的各位专家致以最诚挚的谢意！徐自立二零一零年一月68 参考文献参考文献[1]中华人民共和国国民经济和社会发展第十一个五年规划纲要.2006;Availablefrom:http://www.gov.cn/ztzl/2006-03/16/content_228841.htm.[2]杨征宇,常思勤.汽车覆盖件成形过程及模具制造的数字模拟[J].中国制造业信息化,2006,35(23):66-70.[3]现代模具技术编委会编.现代模具技术[J].汽车覆盖件模具设计与制造,1998:1-5.[4]赵力平,陈炜,季国荣等.基于自顶向下方式的汽车覆盖件模具结构参数化设计[J].模具工业,2005,3(289):3-7.[5]唐镜,黄鹤辉.汽车覆盖件模具冲压有限元仿真分析[J].机电一体化,2008,(9):57-59.[6]王永锋,曹静,刘锐.我国汽车覆盖件模具现状及发展思路[J].现代零部件,2009,(10):58-59.[7]王家昆,李琦,胡东勇等.汽车覆盖件冲压成形的数值模拟[J].模具工业,2006,32(2):1-5.[8]徐丙坤,施法中,陈中奎.板料冲压成形数值模拟中的几个关键问题[J].塑性工程学报,2001,8(2):32-35.[9]熊光楞,王克明,陈斌元等.计算机仿真技术在轿车工业中的应用与发展[J].系统仿真学报,2004,16(1):73-78.[10]周雄辉,王辉.模具敏捷制造系统研究[J].模具技术,2002,(1):45-49.[11]陈亮,陈军,张贵宝等.高强度钢板拉深模具结构有限元分析与试验研究[J].塑性工程学报,2008,15(2):12-45.[12]黄才元.我国冲压设备的现状与发展[J].机械工人（热加工）,2007,(4):23-25.[13]徐刚,鲁洁,黄才元.金属板材成形技术与装备的现状与发展[J].锻压装备与制造技术,2004,(4):16-22.[14]陈炜,董洪智,林忠钦等.汽车覆盖件模具结构CBD研究[J].模具工业,2000,(1):65-67.[15]陈弘彬,李志刚,朱培玉.特征技术在覆盖件模具结构设计中的应用[J].模具技术,1998,(1):36-40.[16]霍同如,徐秉业.板材成形的计算机辅助工程系统[J].力学进展,1996,26(4):548-557.[17]霍同如,徐秉业,宋军.板材成形有限元数值模拟的研究进展[J].力学进展,1996,26(3):379-388.69 西安科技大学硕士学位论文[18]S.KOBAYASHI,J.H.KIM.Deformationanalysisofaxisymmetricsheetmetalformingprocessedbyrigid-plasticfiniteelementmethod[M].MechanicsofSheetMetalForming,ed.KoistinenDP,WangNM.1978,NewYork,N.Y.,USA:PlenumPress.[19]N.M.WANG,M.L.WENNER.Elastic-viscoplasticanalysisofsimplestretchingformingprocesses[M].MechanicsofSheetMetalForming.1978,NewYork,N.Y.,USA:PlenumPress.367-391.[20]N.M.WANG,B.BUDIANSKY.Analysisofsheetmetalstampingbyafiniteelementmethod[J].JournalofAppliedMechanics,ASME1978,(100):73-82.[21]C.H.Toh,S.KOBAYASHI.Deformationanalysisandblankdesigninsquarecupdrawing[J].InternationalJournalofMachineToolDesignEngineering,1985,(25):12-32.[22]S.C.Tang.Verificationandapplicationofthedeformedshapeofautomobiletrunkdeck-lidduringthebinder-wrapstage[M].ExperimentalVerificationfoProcessModels.1981,MetalsPark,Ohio,USA:ASM.[23]S.C.Tang,B.CHAPPUISL.Evaluationofsheetmetalformingprocessdesignbysimplemodels[M],ed.ProcessMaterialsinManufacturing.1988,Chicago,Ill.,USA:ASME.[24]兰凤崇,王玉超,陈吉清.汽车覆盖件冲压成形仿真研究进展[J].计算机集成制造系统,2007,13(8):1586-1596.[25]B.CHAPPUISL,S.C.Tang,X.M.CHEN.Anumericallystablecomputermodelforsheetmetalforminganalysisby2Dmembranetheory,inProceedingsofSheetMetalandStampingSymposium.1993,Washington,D.C.,USA,SAE,155-163[26]林忠钦,李淑慧,于忠奇.车身覆盖件冲压成形仿真[M].2005,北京:机械工业出版社.[27]B.LEE,Y.T.KEUM,R.H.WAGONER.Modelingofthefrictioncausedbylubricationandsurfaceroughnessinsheetmetalforming[J].JournalofMaterialsProcessingTechnology,2002,130-131(20):60-63.[28]Q.V.BUI,PAPELEUXL,P.PONTHOTJ.Numericalsimulationofspringbackusingenhancedassumedstrainelements[J].JournalofMaterialsProcessingTechnology,2004,153-154:314-318.[29]PAPELEUXL,P.PONTHOTJ.Finiteelementsimulationofspringbackinsheetmetalforming[J].JournalofMaterialsProcessingTechnology,2002,125-126(785-791).[30]G.WEI,R.H.WAGONER.Diedesignmethodforsheetspringback[J].InternationalJournalofMechanicalSciences,2004,46(7):1097-1113.70 参考文献[31]D.Banabic,D.S.Comsa,T.Balan.Anewyieldcriterionfororthotropicsheetmetalsunderplane-stresscondition,inProceedingsoftheColdMetalFormingConference.2000,Cluj-Napoca,217[32]李玉强,王勇,董剑安等.覆盖件冲压CAE分析材料模型的工程应用与发展[J].模具技术,2009,(4):49-54.[33]兰箭,董湘怀,李志刚.板料成形毛坯展开方法研究[J].锻压技术,2000,25(4):21-25.[34]H.SHIM,K.SON,K.KIM.Optimumblankshapedesignbysensitivityanalysis[J].JournalofMaterialsProcessingTechnology,2000,104(3):191-199.[35]K.Son,H.Shim.Optimalblankshapedesignusingtheinitialvelocityofboundarynodes[J].JournalofMaterialsProcessingTechnology,2003,134(1):92-98.[36]肖海峰,徐艳琴.基于PAM-STAMP2G的模具结构的优化设计[J].锻压装备与制造技术,2009,(1):67-68.[37]徐金波,董湘怀.基于有限元分析的汽车覆盖件模具设计及优化[J].锻压技术,2004,(1):63-67.[38]林新波,李荣先,柳百成等.利用有限元仿真技术优化模具结构[J].锻压技术,2005,(2):46-49.[39]Y.Kojima,H.morikawa,T.Misuno.Measurementofpressuredistributionindrawingofbars[J].JournalofJapanSocietyforPrecisionEngineering,1985,26(295):855-861.[40]S.Matsubara,H.Kudo.Devisingandcalibrationofasimplesensorfordeterminingpressuredistributionontool-workpieceinterfaces[J].JournalofJapanSocietyforTechnologyofPlasticity,1991,32(364):589-596.[41]S.Matsubara,H.Kudo.Contactpressurecharacteristicsontool-workpieceinterfaceduringplane-strainextrusionofpurealuminium[J].JournalofJapanSocietyforTechnologyofPlasticity,1991,32(366):862-867.[42]S.Matsubara,H.Kudo.Distributionofcontactpressureontool-workpieceinterfaceinaxisymmetricbackwardextrusionofaluminiumcans[J].JournalofJapanSocietyforTechnologyofPlasticity,1991,32(366):868-873.[43]S.Matsubara,H.Kudo.Nonuniformpressuredistributiononpunchnosesurfaceinaxisymmetriccoldforgingprocesses[J].JournalofJapanSocietyforTechnologyofPlasticity,1991,32(366):874-879.[44]Y.Murata,Y.Yabuki.Measurementofdiedeformationinforgingbycapacitivedisplacementtransducer.inProceedingsoftheFourghInternationalConferenceonTechnologyofPlasticity.1993.[45]M.T.Hillery,S.Griffin.Anembedded-strain-gaugetechniqueofstressanalysisinrod71 西安科技大学硕士学位论文drawing[J].JournalofMaterialsProcessingTechnology,1994,47(1-2):1-12.[46]C.MacCormack,J.Monaghan.FailureanalysisofcoldforgingdiesusingFEA[J].JournalofMaterialsProcessingTechnology,2001,117(1-2):209-215.[47]R.Balendra,Y.Qin,X.Lu.Analysis,evaluationandcompensationofcomponent-errorsinthenet-formingofengineeringcomponents[J].JournalofMaterialsProcessingTechnology,2000,106:204-211.[48]R.Balendra.Nett-shapeforming:state-of-the-art[J].JournalofMaterialsProcessingTechnology,2001,115(2):172-179.[49]R.Balendra.Overviewof"Die-cavityelasticityconsiderationsfornett-forming(DICAVITY)"[J].JournalofMaterialsProcessingTechnology,2001,115(2):180-186.[50]R.Balendra,Y.Qin.Researchdedicatedtothedevelopmentofadvancedmetal-formingtechnologies[J].JournalofMaterialsProcessingTechnology,2004,145(2):144-152.[51]X.Lu,R.Balendra.EvaluationofFEmodelsforthecalculationofdie-cavitycompensation[J].JournalofMaterialsProcessingTechnology,1996,58(2-3):212-216.[52]X.Lu,R.Balendra.Pressurecontoursonformingdies:PartII:Finiteelementsimulation[J].JournalofMaterialsProcessingTechnology,2001,115(2):220-226.[53]X.Lu,R.Balendra.Finiteelementsimulationfordie-cavitycompensation[J].JournalofMaterialsProcessingTechnology,2001,115(2):227-232.[54]Y.S.Lee,J.H.Lee,T.Ishikawa.Analysisoftheelasticcharacteristicsatforgingdieforthecoldforgeddimensionalaccuracy[J].JournalofMaterialsProcessingTechnology,2002,130-131:532-539.[55]Y.S.Lee,J.H.Lee,J.U.Choi,T.Ishikawa.Experimentalandanalyticalevaluationforelasticdeformationbehaviorsofcoldforgingtool[J].JournalofMaterialsProcessingTechnology,2002,127(1):73-82.[56]Y.S.Lee,J.H.Lee,Y.N.Kwon,T.Ishikawa.Analysisoftheelasticcharacteristicsatdieandworkpiecetoimprovethedimensionalaccuracyforcoldforgedpart[J].JournalofMaterialsProcessingTechnology,2004,153-154:1081-1088.[57]WeilongHu,AmberWang,Advancedtoolsandmethodstoimprovedurability,stabilityandtryoutrobustnessofstampingdiesystem,inFordtechnicalreports.2007.[58]张贵宝,高强度薄板拉深模具结构分析关键技术研究及应用[D].2008,上海交通大学.[59]张贵宝,陈军,王晓方.基于板料成形数值模拟的冲压模具结构分析方法[J].机械工程学报,2008,44(8):174-179.[60]徐自立,娄茂昆,郭志勇,朱向会.基于板料成形数值模拟的覆盖件模具结构优化[J].72 参考文献锻压技术,2010,35(1):99-102.[61]Z.Marciniak,J.L.Duncan,S.J.Hu.MechanicsofSheetMetalForming[M].2002,Oxford:Butterworth-Heinemann.[62]T.BURANATHITI,J.CAO.Benchmarksimulationresults:automotiveunderbodycrossmember(Benchmark2),inNumisheet2005.2005,Detroit-USA,[63]W.F.Hosford,R.M.CADDEL.MetalFormingMechanicsandMetallurgy[M].2007,Cambridge,UK:CambridgeUniversityPress.[64]李泷杲,王书恒,徐岩.金属板料成形有限元模拟基础——PAMSTAMP2G(Autostamp)[M].2008,北京:北京航空航天大学出版社.[65]AutoFormEngineeringGmbH,AutoFormSoftwareManual.2008.[66]ESIGroup.PAM-STAMP2G2008User'sGuide.2008.[67]周传月,腾万秀,张俊堂.工程有限元与优化分析应用实例教程[M].2005,北京:科学出版社.[68]AltairEngineeringInc.,RADIOSS,MotionSolve,andOptiStructUser'sGuide.2008.[69]AltairEngineeringInc.;Availablefrom:http://www.altairhyperworks.com.cn/ResLibSearchResult.aspx?keywords=OptiStruct&industry=All&product_service=All&category=All.[70]M.P.Bendsoe,N.Kikuchi.Generatingoptimaltopologiesinstructuraldesignusingahomogenizationmethod[J].ComputerMethodsinAppliedMechanicsandEngineering,1988,71(2):197-224.[71]Y.M.Xie,G.P.Steven.Asimpleevolutionaryprocedureforstructureoptimization[J].Computstruct,1993,(49):885-896.[72]Y.M.Xie,G.P.steven.Optimaldesignofmultipleloadcasestructuresusinganevolutionaryprocedure[J].EngComput,1994,(11):295-302.73 西安科技大学硕士学位论文附录攻读硕士学位期间发表的学术论文[1]徐自立,娄茂昆,郭志勇,朱向会.基于板料成形数值模拟的覆盖件模具结构优化[J].锻压技术,2010,35(1):99-102.[2]娄茂昆,徐自立,郭志勇.冲压模具结构分析和优化方法研究.inAlatir2009HyperWorks技术大会论文集.2009.[3]徐自立,李明,郭志勇,曹小平.某选矿厂富液池的强度分析与结构优化[J].矿山机械,2008,36(17):91-95.[4]郭志勇,曹小平,徐自立.黄土地区静压预制桩应用测试研究[J].建筑科学,2008,24(9):68-70.[5]高崇霖,刘宝许,王巍,来兴平,徐自立.小净距立交隧道开挖围岩稳定性声发射监测[J].西安科技大学学报,2009,29(3):274-277.74 附录荷载转换程序Pam2OptiStruct.exe的FORTRAN源码:PROGRAMPam2OptiStructIMPLICITNONEINTEGERFLAG,ICHARACTERBUFFER*100,SIGNAL*695INTEGER*4NodeNumberREAL*8FX,FY,FZ!!INITIALIZATIONFLAG=010SIGNAL='Number'!欢迎信息WRITE(*,*)'***********************************************************'WRITE(*,*)'ThisprogramchangethedataformatfromPAMtoOPTISTRUCT.'WRITE(*,*)'BUILTBY:XUZILI'15WRITE(*,*)'VERSION:2.0.0'WRITE(*,*)'DATE:09/25/2009'WRITE(*,*)'***********************************************************'!打开PAM输出的NormalContactPressure文件并转化为OptiStruct格式20WRITE(*,*)'InputfilenameofNormalContactPressure:'READ(*,'(A100)')BUFFERIF(BUFFER(1:1).EQ.'')BUFFER='NCP.asc'!若输入空格则为默认文件名：NCP.ASCOPEN(UNIT=11,FILE=BUFFER)OPEN(UNIT=12,FILE='NCP.INC')25!找到节点力开始的位置DOWHILE(FLAG==0)READ(11,'(1X,A20)')BUFFERDOI=1,6IF(BUFFER(I:I).EQ.SIGNAL(I:I))THEN30FLAG=1ELSEFLAG=0EXIT75 西安科技大学硕士学位论文ENDIF35ENDDOENDDO!读入节点力信息，忽略分量都小于1E-6的节点力FLAG=0100FORMAT(I9,3(E13.5))40200FORMAT('FORCE,10,',I9,',0,1000,',E13.5,',',E13.5,',',E13.5)DOWHILE(FLAG==0)READ(11,100)NodeNumber,FX,FY,FZIF(NodeNumber==0)THENFLAG=145ELSEIF(.NOT.(ABS(FX)<1E-6.AND.ABS(FY)<1E-6.AND.ABS(FZ)<1E-6))WRITE(12,200)NodeNumber,FX,FY,FZENDIFENDDO50WRITE(*,*)'NCP.INC创建完毕'CLOSE(11)CLOSE(12)!NCP.INC创建完毕，开始创建TCP.INCWRITE(*,*)'InputfilenameofTangentContactPressure:'55READ(*,'(A100)')BUFFERIF(BUFFER(1:1).EQ.'')BUFFER='TCP.asc'!若输入空格则为默认文件名：TCP.ASCOPEN(UNIT=13,FILE=BUFFER)OPEN(UNIT=14,FILE='TCP.INC')FLAG=060DOWHILE(FLAG==0)READ(13,'(1X,A20)')BUFFERDOI=1,6IF(BUFFER(I:I).EQ.SIGNAL(I:I))THENFLAG=165ELSEFLAG=0EXIT76 附录ENDIFENDDO70ENDDOFLAG=0300FORMAT('FORCE,15,',I9,',0,1000,',E13.5,',',E13.5,',',E13.5)DOWHILE(FLAG==0)READ(13,100)NodeNumber,FX,FY,FZ75IF(NodeNumber==0)THENFLAG=1ELSEIF(.NOT.(ABS(FX)<1E-6.AND.ABS(FY)<1E-6.AND.ABS(FZ)<1E-6))WRITE(14,300)NodeNumber,FX,FY,FZ80ENDIFENDDOWRITE(*,*)'TCP.INC创建完毕'CLOSE(13)CLOSE(14)85!打开OptiStruct导出的*.FEM文件，并在其中插入'INCLUDENCP.INC'和'INCLUDETCP.INC'SIGNAL='$$HyperMeshCommandsforloadcollectorsnameandcolorinformation$'WRITE(*,*)'InputOptiStructexportedfilename:'READ(*,'(A50)')BUFFER!若输入空格则为默认文件名：MODEL.FEM90IF(BUFFER(1:1).EQ.'')BUFFER='MODEL.FEM'OPEN(UNIT=15,FILE=BUFFER,PAD='YES')OPEN(UNIT=16,FILE='MODEL.fem')FLAG=0DOWHILE(FLAG==0)95READ(15,'(A)')BUFFERWRITE(16,'(A)')TRIM(BUFFER)DOI=1,69IF(BUFFER(I:I).EQ.SIGNAL(I:I))THENFLAG=1100ELSEFLAG=077 西安科技大学硕士学位论文EXITENDIFENDDO105IF(FLAG==1)THENREAD(15,'(A100)')BUFFERWRITE(16,'(A)')TRIM(BUFFER)WRITE(16,'(A)')'$HMNAMELOADCOL10"ncp"'!插入法向接触力的荷载集WRITE(16,'(A)')'$HWCOLORLOADCOL1015'110WRITE(16,'(A)')'$'WRITE(16,'(A)')'$HMNAMELOADCOL15"tcp"'!插入切向接触力的荷载集WRITE(16,'(A)')'$HWCOLORLOADCOL1520'WRITE(16,'(A)')'$'WRITE(16,'(A)')'INCLUDENCP.INC'115WRITE(16,'(A)')'INCLUDETCP.INC'WRITE(16,'(A)')'$'ENDIFENDDODOWHILE(FLAG==1)120READ(15,'(A100)',END=400)BUFFERWRITE(16,'(A)')TRIM(BUFFER)ENDDO400CLOSE(15)CLOSE(16)125WRITE(*,*)'NORMALTERMINATION'ENDPROGRAMPam2OptiStruct78