ANSYS Workbench中高斯移动热源下温度与应力仿真

上篇详细介绍了AnsysWorkbench中的塑性，本文主要介绍超弹性。超弹性材料在工程应用上很广泛，可以说随处可见，例如橡胶圈等。3超弹性3.1超弹性的基本理论（1）超弹性的概念橡胶是最典型的超弹性物质，具有以下特点：能承受大弹性(可恢复)变形，任何地方都可达100-700%。体积变化很小，几乎不可压缩。它们的应力-应变关系是高度非线性的。通常，拉伸状态下，材料先软化再硬化，而压缩时材料急剧硬化。（2）对超弹性的基本假设为了便于计算，将真实橡胶行为的理想化，对超弹性材料基本假设如下：材料是各向同性的，等温的；变形是可完全恢复的；材料是完全或接近不可压缩的。3.2AnsysWorkbench中的超弹性本构3.2.1HyperelasticExperimentData超弹性材料试验数据对于超弹性材料，应该注意的是用于拟合曲线的测试的数据（除体积测试数据）需要工程应力-应变数据，体积测试数据需要真实应力-应变数据。这与基于vonMises,Tresca塑性曲线拟合中收集数据方式不同（塑性曲线需要真实应力-应变）。对于试验数据需要调整滞回和应力软化行为，采用稳定的曲线（应该偏移至零应力和零应变）用来进行曲线拟合。对于橡胶超弹体进行模拟仿真，首先要进行材料曲线的拟合工作，WB提供了曲线拟合工具，可以帮助把实验数据转化成各种超弹模型能使用的应变能量密度函数系数，曲线拟合操作见实例1。UniaxialTestData:单轴试验数据。BiaxialTestData:双轴试验数据。ShearTestData:剪力试验数据。VolumetricTestData:体积试验数据。SimpleShearTestData:简单剪力试验数据。UniaxialTensionTestData:单轴拉伸试验数据。UniaxialCompressionTestData:单轴压缩试验数据。3.2.2Hyperelastic超弹性材料本构超弹性材料一般不定义杨氏模量。超弹性材料的典型体积行为可分为两类。第一类为不可压缩或几乎不可压缩的模型，比如橡胶聚合物，这类本构有Neo-Hookean、Mooney-Rivlin、多项式形式、Ogden、Arruda-Boyce、Gent、Yeoh和ExtendedTube；第二类为可压缩材料，比如泡沫在变形过程中会经历很大的体积变化，这类本构有Blatz-Ko和OgdenFoam。Neo-Hookean：这是最简单的超弹性模型，采用常剪切模量。然而，它局限于单轴拉伸时应变为30-40%和纯剪时应变为80-90%的情况。需要定义材料的初始剪切模量μ，材料不可压缩性参数d（初始体积模量=2/d）。Arruda-Boyce：一种基于统计力学的模型，应变范围为300%。需要定义材料的初始剪切模量μ，极限网络伸展λL，材料不可压缩性参数d（初始体积模量=2/d）。当λL接近无穷大时，该模型等同于Neo-Hookean模型。Gent：应变范围为300%。需要定义材料的初始剪切模量μ，极限值L，材料不可压缩性参数d（初始体积模量=2/d）。当L接近无穷大时，该模型等同于Neo-Hookean模型。Blatz-Ko：描述可压缩泡沫橡胶材料的最简单形式。需要定义材料的初始剪切模量μ。可以参考泡沫橡胶数据——泊松比为0.463，密度为1150kg/m，μ=1040MPaMooney-Rivlin：分为N=2,3,5和9项MooneyRivlin模型。2项是最常用的，在拉伸应变大致90-100%时有效，但2项MR不能很好描述压缩行为特性，必须保证提供足够的高阶项的数据，捕捉任何拐点。5和9项可以用于应变高达100-200%的情况。需要定义材料参数C01、C10、C11、C20.....d等。可以参考橡胶的IRHD硬度指标经验公式计算各材料参数。例如氢化丁睛橡胶H-NBR75（硬度为75MPa）各参数计算如下：又比如橡胶数据：密度为1150kg/m，泊松比为0.499,C10=0.55MPa，C01=0.14MPa。Polynomial多项式：与Mooney-Rivlin本构类似。1阶等效于N=2的MR本构，2阶等效于N=5的MR本构，3阶等效于N=9的MR本构。更高的阶数可以更好地拟合精确解。然而，它可能会导致拟合材料常数的数值困难，并且它还需要足够的数据来涵盖您可能感兴趣的整个变形范围。由于这些原因，不建议使用非常高的阶数。需要定义材料参数C10、.....d1、.....等。Yeoh：也称为简化多项式形式，一阶等效于Neo-Hookean本构。需要定义材料参数C10、.....d1、.....等。Ogden：与多项式形式类似，可以用于应变达700%的情况，可能更精确但计算更耗时。对于N=1和α1=2，等价于Neo-Hookean。对于N=2、α1=2和α2=-2，等效为2参数Mooney-Rivlin模型。需要定义材料参数μ1、.....α1、.....d1、.....等。OgdenFoam：描述可压缩泡沫材料，较Blatz-Ko在大应变下精度更高，应变可达700%。对于N=1、α1=-2、μ1=–μ和β=0.5，Ogden选项等价于Blatz-Ko选项。需要定义材料参数μ1、.....α1、.....β1、.....等。ResoponseFunction：以试验数据描述超弹性材料本构的响应函数。ExtendedTube：描述填充加强筋的橡胶材料，支持材料曲线拟合。MullinsEffect：描述随载荷变化的超弹性材料的拉伸软化，例如聚合填充材料在承受循环载荷后，卸载时要比加载时更软，造成了应力应变中的滞留效应和材料的不可恢复变形。实例1超弹性材料曲线拟合，以Yeoh3rd本构为例。注意：建立至少两倍于需要计算材料本构参数数目的数据点，同时需要考虑滞后等效应；一般应有三种试验数据(如单轴拉伸、双轴拉伸及平面拉伸)，以便完整描述材料本构参数，如果仅有单轴拉伸数据，则不能生成承受大剪切的模型；试验数据应该包括全部关注的应变范围，如果只有0~100%应变的测试数据，则不能生成承受150%应变的模型，外推数据会产生极大的误差。以此数据建立一个应力应变试验数据文档，推荐使用Excel文件，选择Excel中的数据，复制到剪切板中。Step1准备数据。下载图惜准备的单轴拉伸试验数据（本数据无实际工程应用意义），读者也可使用自己试验数据。链接：https://pan.baidu.com/s/1ll1DcVMSICBrIopo6zhfmA提取码：ebebStep2数据导入。新建一个EngineeringData，双击左键打开，如图所示。分别点击HyperelasticExperimentalData→UniaxialTestData和Hyperelastic→Yeoh3rdOrder，按住左键将两个菜单栏都拖至新建的材料1框处。在Outline菜单栏中点击UniaxialTestData进行编辑，在Properties菜单栏先点击Temperature栏，输入试验测试时的温度，在Strain-Stress栏右键粘贴EXCEL中的数据，即可将试验数据导入。Step2曲线拟合。在Outline菜单栏中点击Yeoh3rdOrder进行编辑，必须注意在ErrorNormforFit(拟合误差)进行选择。默认为NormalizedError，表示对输入数据不加权，适合一般情况。AbsoluteError表示对输入数据中大的数据加权，适合大应变情况，本例选择AbsoluteError。然后右键点击CurveFitting(曲线拟合)，在出现的菜单中先点击SolveCurveFit，再点击CopyCalculatedValuesToProperty，即可把拟合好的数据赋值给材料参数。以上便完成了材料曲线的拟合。3.3收敛性控制由于超弹性体的近乎不可压缩特性，仿真计算中往往会出现收敛困难的情况，经常让工程师苦不堪言，属于难点问题，各类书籍文献也会相应给出推荐求解方法，以下做简单归纳总结。3.3.1尽量简化为平面模型在进行处理计算几何选择方面，如果能简化为2D几何计算模型的尽力避免3D模型几何。3D计算模型相对2D模型计算采用成数倍数量的单元节点，同时非线性计算的迭代求解和存储耗费，这都给计算资源和效率带来了挑战。2D模型几何边线更方便控制尺寸份数以细化几何拐角，有效防止挤压或者拉伸时候的单元扭曲，提高收敛可能性。更好处理橡胶大变形状态下，橡胶与金属件接触关系以及橡胶自接触；规避开3D模型实体六面体单元的可能“刚硬”。如图所示是两类形状密封圈挤压状态模拟计算案例，模型计算都不可避免设置模型自接触、拐角单元细化以及单元形状控制以防止过度畸形。3.3.2单元技术的应用单元技术在几何非线性一文中作了详细介绍，以下作简要回顾，并针对超弹性体的问题做针对性介绍。剪切自锁表现为在弯曲过程中，完全积分的低阶单元呈现过度刚化，计算出的剪切应力明显过大，变形结果明显过小。体积自锁表现为完全或接近不可压缩的材料在体积变形问题上，完全积分单元呈现出的过度刚化，造成结构难以变形，收敛困难。除平面应力模型以外，其余平面应变、轴对称和3D模型均可能发生体积自锁。为克服剪切自锁与体积自锁，引入相应单元技术。（1）缩减积分缩减积分可解决剪切自锁和接近不可压缩的材料的体积自锁问题，但是只有一个积分点的低阶单元在发生完全变形时，会出现零能模式（即沙漏）。为避免沙漏，可以采用高阶单元或加密的低阶单元（厚度至少3层）。当厚度方向上只有一层单元时，使用完全积分单元有利于提高精度。几何结构的单元控制ElementControl设为手动Manual后,用户可手动触发积分方案选项，它将影响单元内积分点的数量。（2）增强应变只适用于四边形低阶单元Plane182或六面体低阶单元Solid185，在单元模型中增加附加自由度，能很好地解决剪切自锁问题，但会花费较多的计算时间。单元质量越好，计算精度越高。如果要采用增强应变单元，需要在WB的Mechanical中，对相应的Body插入命令：Plane182增强应变命令为keyopt,matid,1,1，Solid185增强应变命令为keyopt,matid,2,2。（3）混合U-P也称杂交单元，在单元内插静水压力作为附加自由度。WB的大部分18X等于均可应用，可以很好地解决体积自锁问题。如果要采用增强应变单元，需要在WB的Mechanical中，对相应的Body插入命令：keyopt,matid,6,1，注意Solid285默认已采用混合U-P单元，无需插入命令。如果采用混合U-P技术后任难收敛，可以在求解栏插入Command：Solc,,,,vtolvtol为容差值，取值范围0~1，默认1e-5，推荐1e-5~1e-2。使用该Command后，原定义的时间步长将全部重置，再结合重启控制，可保证收敛计算的充分条件。总结：对于剪切自锁，推荐采用增强应变；对于体积自锁，推荐采用混合U-P。注意Plane182和Plane183单元中的平面应力模型在处理不可压缩模型时没有困难，所有不需要缩减积分、增强应变、混合U-P技术。在WB中，几种常见的网格单元可使用的单元技术如下，低阶单元和高阶单元都支持完全积分计算。——表示不可修改，用户的修改将无效。3.3.3自适应网格设置在超弹性问题上，有时会遇到应用了增强应变或混合U-P技术，计算过程仍旧遇到因为单元畸形发散的问题。此时可以应用自适应网格技术，这是WB16之后的功能。非线性自适应区域条件使您能够在求解阶段更改网格以提高精度，而不会产生大量计算损失。非线性自适应区域功能是完全自动的。它根据某些标准充当重新网格控制器。自适应网格技术支持Plane182、Plane183、Plane222、Solid187、Solid285、Solid227等单元。使用时一般要求用户打开大变形控制，以及存储所有求解结果计算点StoreResultsAt=AllTimePoints，读者也可以打开重启动设置以帮助发散后重新处理求解子步等。如果发现自适应网格不可用，读者可以检查一下分析设置——输出控制：StoreResultsAt是否已设置为AllTimePoints。在Mechanical中插入自适应网格方法如下：右击“静态结构”——插入——非线性自适应。自适应网格标准有3种准则。Energy：如果单元应变能大于等于组件中平均应变能乘以自定义参数，单元将会分割与重划。以获得高精度区域模拟结果。Box：在指定区域内，进行分割和重划。Mesh：如果是3D网格，有三种选项控制，其中倾斜度Skewness：默认值为0.9，取值范围0~1，1为最差，当网格倾斜度差于用户设定值时，将进行网格重新生成。雅可比率JacobianRatio：这个值定义了单元积分点处的最大雅可比矩阵比。取值范围为-1~1。默认值为0.1（理想值为1，大于0.7可接受，质量较好，小于0.5，准确性不能保证）。当网格雅可比率差于用户设定值时，将进行网格重新生成。倾斜度和雅可比率SkewnessandJacobianRatio：设置倾斜度和雅可比两项的基准值，差于任意一项都将重划网格。如果是2D网格，只能控制最大拐角，默认160°，读者按照需求更改，越小表示越严格也越耗费资源。检查CheckAt：设置等距点和重复率，默认为1，表示每一子步都检查。自适应网格的应用实例见的小鲸鱼的小鲸鱼的文章《非线性大变形解决之道》，传送门在此ttps://mp.weixin.qq.com/s/14kM-QfOMEMfFRuw-9zzeQ3.3.4非线性控制超弹性问题一般都是大应变问题，需要用到分步控制，有时往往需要我们反复细化分步。超弹性计算过程中，也经常会遇到算了很久结果报错不收敛，有经验的工程师往往不会从头开始计算，而是设置稳定性后，从不收敛往前倒推几步重启计算（比如算到第90步不收敛，开启稳定性控制后从第80步重启），前提是您设置了重启控制。关于分步、重启、收敛控制、稳定性控制，都是在超弹性问题中需要常用到的，由于图惜在《一线工程师总结AnsysWorkbench工程应用之——结构非线性（上）：几何非线性》3篇文章中对此已经做了详细介绍，此处就不再凑字数了，感兴趣的可以自行搜索。以下对非线性控制做一个极不权威的使用经验总结：遇到非线性不收敛时，打开稳定性控制，使用“常数”或“减少”的默认设置，然后回调数步重启，往往能解决大部分不收敛问题。遇到非线性不收敛时，可以增加子步数量或者减小子步时间，然后回调数步重启。遇到强烈非线性问题，可以将非线性控制中的牛顿选项改为非对称，手动关闭力矩、位移、旋转收敛、打开线搜索，然后回调数步重启。更多的收敛控制方法，等待您的探索。以下为文献1中的综合实例。实例2超弹性综合实例概述：本实例采用平面应变模型，分析将球体压入橡胶的过程，实例中综合应用接触、网格划分、单元技术等。Step1建模。新建静力学分析项目。右击Geometry——新DM结构，因为模型为二维平面模型，只能在XYPlane基准上建立草绘。在1个草图中绘制，圆形(直径为12mm)表示刚性面，方形(长度20mm、高度12mm)表示橡胶面，两者相切于方形上边中点。再点击概念——草图面，选择刚才草图，操作中选择添加冻结，点击生成，便一次创建了2个面体。Step2定义材料。双击EngineeringData,新增一个超弹性材料“xiangjiao”，采用Ogden3rdOrder本构方程,具体数值依次为43438Pa、1.3、82.74Pa、5、-689.5Pa、-2、2.9E-08Pa^-1、0Pa^-1、0Pa^-1Step3定义零件行为。在项目管理主界面，点击项目中的“几何结构”栏，设置分析类型为2D。注意，虽然不在此处设置2D分析也能在Mechanical中设置平面应力/应变/轴对称等2D模型，但是会影响一些功能的使用，比如自适应网格技术，所以还是建议读者设置此处。双击项目中的“模型”栏，进入Mechanical界面。设置模型为2D平面应变、长方形模型材料为xiangjiao。对长方形模型右键插入Command:keyopt,matid,6,1(为处理体积自锁，采用混合U-P)；圆形Surface所有设置默认。注意：本例如果采用平面应力问题，则不存在体积自锁问题，不需要插入Command即可求解。Step4接触设置。圆形刚性面与橡胶上面始终处于接触状态，设置0.1的摩擦接触（此处与文献1略有不同）。在Contact选择橡胶上边线，在Target选择圆形的圆周线(两者定义不可颠倒)，定义接触类型为绑定，设置接触剪切=关闭。Step5网格划分。设置全局网格为二阶平面网格，即Plane183网格。单元中节点保留可以极大地减小计算规模，用较小的单元数即可保证计算精度；但是如果同样多的单元数，单元中节点保留即意味节点数较多，计算所需内存也增多。例如，计算过程中如果出现内存不足时，可以通过取消单元中节点降低计算规模。另外在非线性计算过程中，单元中节点保留的收敛难度大于单元中节点取消。针对复杂接触、弯曲、不可压缩材料这类复杂非线性问题，可以采用取消单元中节点。本例采用单元中节点保留是为了避免接触穿透现象。为保证类似本例的收敛计算，网格最好采用正方形或正方体网格(三维模型)，因此，定义网格Size比定义Sweep网格更有意义。再右键插入Size定义,在Geometry选择橡胶平面，在Type处定义ElementSize，ElementSize定义为1mm，Behavior定义为Hard(强制规定胶平面的单元尺寸为1mm,即1mmX1mm的正方形)。Step6求解设置。打开大变形。总步数默认为1步不修改，总时间默认为1s也不修改。分步初始、最小、最大时间依次为1e-003s、1e-003s、1e-002s，这个设置可等效为分步初始、最小、最大数量依次为1000,100,1000(分步时间与分步数量互为倒数)。Step7边界条件设置。对橡胶平面下边线加载FixedSupport(完全约束)，对上圆形面加载RemoteDisplacement(远程位移),大小方向为Y方向偏移-9mm，加载点位于坐标(8mm,26mm)。注意：加载点并不在圆形区域内，而在圆形区域外，这是根据实际测试工况定义的。读者可以先定义远程点后，再使用远程点加载远程位移。Step8计算与后处理。点击Slove计算，经过选代计算、结果收敛、后处理插入TotalDeformation(总位移)和橡胶平面的VectorPrincipalElasticStrain(弹性应变主矢量)结果如下图。如果读者遇到收敛困难，可以在求解栏插入Command：Solc,,,,1e-5，如果还难于收敛，可以应用上文提到的自适应网格、重启+稳定性控制、修改为线性单元、增加子步数量等。写在最后至此，非线性材料已经完成了弹塑篇、蠕变粘塑篇、超弹篇的总结，下一篇将拾缀拾缀其余遗漏的本构模型，敬请关注。参考文献：[1]《AnsysWorkbench有限元分析实例详解》——周炬、苏金英[2]《计算方法|浅析橡胶超弹体与粘弹性仿真(超弹篇)》——付稣昇[3]《ANSYS高级结构非线性培训手册》——安世亚太[4]ANSYS2022帮助文件来源：CAE初行者