词、曲和有限元-本构模型
诗词通过汉字组合赋予了时代感、意境美；音乐通过词曲组合释放了视听感、经典美；有限元通过点线面体、材料单元属性和边界条件展现了立体感、彩图Good。
从唐诗到宋词，再到毛泽东诗词，从孩提时代《射雕英雄传》主题曲，到小学时的《八仙过海》主题曲，中学时的《乙未豪客传奇》主题曲，大学时的《伤心太平洋》，再到《后来》，《小苹果》，直到如今的《沙漠骆驼》，诗词陪伴着我们成长，音乐激励着我们上进。
诗词、歌词好比有限元模型中的点、线、面、体、单元和节点；旋律宛如边界条件，诗词、音乐不同等同于分析项目也不同。因此，词、曲和有限元有异曲同工之妙。
言归正传，接下来总结有限元的两部分核心内容：本构模型和分析方法。并用经典命令流的方式对相关的内容举例论证，看此文章的朋友可以通过信息下载，在ANSYS软件经典界面中运行测试。
压力容器有限元分析类型也很多，结构分析、热分析、屈曲分析、疲劳分析…最基本的是分析方法和材料的本构模型。材料的本构模型包括：线弹性模型（各向同性材料、正交各向异性材料、各向异性材料）、弹性-理想塑性模型、塑性增量理论（材料的真实应力-应变曲线）。
分析方法和本构模型是相对应的，分析设计方法包括应力分类法和直接法。如图1所示。
应力分类法对应于线弹性模型即胡克定律；
极限载荷法对应于弹性-理想塑性模型；
弹塑性分析法对应于塑性本构关系----增量理论，输入材料的真实应力-应变曲线求解。
图1：压力容器的分析设计方法
实施一完整的压力容器分析设计的一般过程如下：
分析设计条件的提出----载荷分析----材料选择----结构设计----总体强度尺寸确定----局部结构强度尺寸初定----分析模型构建----应力分析----强度评定（通过）----技术说明----技术要求----图样绘制----设计文件整理。
压力容器分析设计时，首先进行载荷分析，确定结构的载荷形式、量级、作用区域和波动范围，诸如内压、外压、重力、惯性力及载荷循环次数。其次进行结构分析，确定应力分析部位和构造力学模型（包括力学和位移的边界条件），确定材料的本构关系，确定分析设计方法和结果处理方法。
一、 本构模型：
三种本构模型对应的曲线如下图2所示。线弹性是数学中的直线方程y=kx（k为直线的斜率，x为正值），有限元中为σ=Eε；弹性-理想塑性是双线性方程：一个是线弹性的直线方程；一个是y=c（常数），有限元中定义为结构的屈服强度，即结构的最大等效应力不会超过材料的屈服强度，而结构的应变在增加。塑性增量理论是指材料的真实应力应变曲线。
图2： 材料的本构模型曲线
A：应力分类法中的本构模型：
采用线弹性的胡克定律，软件中包含三种需要输入的材料参数：
1）Isotropic：各向同性弹性性能参数；应力分类法时对应于常见的塑性金属材料，不包括铸铁和铸钢等脆性材料。
2）Orthtropic：正交各向异性弹性性能参数；对应于非金属材料，如玻璃钢。
3）Anisotropic：各向异性弹性性能参数；对应于非金属材料。
软件中的材料树形结构如图3所示。那么这三种材料属性如何用命令流输入以及输入后软件的显示结果又是怎样？
图3：ANSYS软件中的线弹性材料模型
第一：如果是各向同性材料，这是材料力学和弹性力学的基本假设，此时只输入弹性模量和泊松比。
1）弹性模量E：材料在弹性变形阶段，其应力和应变成正比例关系（即符合胡克定律），其比例系数E称为弹性模量。关系式为：σ=Eε
2）切变模量G：指材料在弹性变形阶段内，剪切应力与对应剪切应变的比值。关系式为：τ=Gγ
3）泊松比μ：材料沿载荷方向产生伸长(或缩短)变形的同时，在垂直于载荷的方向会产生缩短(或伸长)变形。垂直方向上的应变ζ与载荷方向上的应变ε之比的负值称为材料的泊松比。以μ表示泊松比，则μ=-ζ/ε。在材料弹性变形阶段内，μ是一个常数。
针对各向同性材料，因G=E/[2(1+μ)]，因此，有限元分析时只要输入两个独立的弹性常数便可确定第三个，一般输入弹性模量EX和泊松比NUXY。
| ！各向同性材料的命令流如下： /PREP7 !材料1为各向同性材料 MP,EX,1,2e5 !定义材料1的弹性模量 MP,NUXY,1,0.3 !定义材料1的泊松比 |
图4：各向同性材料参数输入后软件的显示结果
第二：如果要做正交各向异性的材料，在相互垂直的两个方向上具有不同的性能指标。因此需要输入X、Y、Z三个方向的弹性模量以及XY、XZ、YZ三个面上的剪切模量以及泊松比。
例如深冷设备中内外罐体之间的非金属支撑材料----玻璃钢，需要输入9个参数。此时，如果模型中只有这一种材料，那么总体坐标系与你所定义的三个方向的参数要对应，X对应轴向；Y对应径向；Z对应环向。如果模型中既有各向同性又有各向异性材料，那么对于各向异性材料，建模时要单独定义局部坐标系，并且坐标系的方向要与所定义的材料属性对应。分析完成后进行评定时，需要评定轴向拉或压应力，剪切应力，弯曲应力等，分别小于其许用轴向拉压应力，许用剪切应力，许用弯曲应力。
正交各向异性材料的相关命令流如下：
| ！Orthtropic：正交各向异性弹性性能参数 /PREP7 !材料2为正交各向异性材料 MP,EX,2,2E5 MP,EY,2,4E5 MP,EZ,2,5E5 ! 设定弹性模量Ex,Ey,Ez MP,NUXY,2,.1 MP,NUYZ,2,.2 MP,NUXZ,2,.3 ! 设定泊松比υyx,υzy,υzx MP,GXY,2,2.5E5 MP,GYZ,2,2.5E5 MP,GXZ,2,2.5E5 ! 设定剪切模量Gxy,Gyz,Gxz MPLIST,1 ! 列显材料的数据 |
图5：正交各向异性材料参数输入后软件的显示结果
第三：Anisotropic：各向异性弹性性能参数。以下命令流中的数值仅为参考值。
| /PREP7 ! 采用<TB>/<TBTEMP>/<TBDATA>的ANEL方式设置材料3(不考虑温度变化) TB,ANEL,3,,,1 ! 设定材料3为各向异性,并用弹性系数矩阵 TBDATA,1,.5E-5,-.25E-6,-.6E-6 ! 从第1项开始设置数据 TBDATA,7,.25E-5,-.4E-6,,,,.2E-5 ! 从第7项开始设置数据 TBDATA,16,.4E-5,,,.4E-5,,.4E-5 ! 从第16项开始设置数据 TBLIST,ANEL,3 ! 列显各向异性材料的数据 |
图6：各向异性材料参数输入后软件的显示结果
另外，热分析中需要定义随温度变化的热导率参数，热应力间接耦合法分析时先定义材料的热导率求解温度场，然后进行由热到结构的单元类型的转换，定义与温度相关的材料的弹性模量和泊松比以及热膨胀系数参数。
定义不同温度下材料的热导率，命令流如下（以材料N06625为例）：
| /prep7 et,1,solid70 MPTEMP,,,,,,,, MPTEMP,1,20 MPTEMP,2,50 MPTEMP,3,100 MPTEMP,4,150 MPTEMP,5,200 MPTEMP,6,250 MPTEMP,7,300 MPTEMP,8,350 MPTEMP,9,400 MPTEMP,10,450 MPTEMP,11,500 MPTEMP,12,550 MPDATA,KXX,1,,9.8*1E-3 MPDATA,KXX,1,,10.2*1E-3 MPDATA,KXX,1,,10.9*1E-3 MPDATA,KXX,1,,11.7*1E-3 MPDATA,KXX,1,,12.4*1E-3 MPDATA,KXX,1,,13.2*1E-3 MPDATA,KXX,1,,13.9*1E-3 MPDATA,KXX,1,,14.7*1E-3 MPDATA,KXX,1,,15.4*1E-3 MPDATA,KXX,1,,16.1*1E-3 MPDATA,KXX,1,,16.9*1E-3 MPDATA,KXX,1,,17.6*1E-3 !N06625材料的热导率ASME II-D P769页 |
图7：热分析时热导率参数的输入
或者
ET,1,PLANE77 KEYOPT,1,3,1 MPTEMP,1,20,100,150,200,250,300 MPTEMP,7,350,400,450,500,550 MPDATA,KXX,1,1,41*1E-3,40.6*1E-3,40.4*1E-3,40.1*1E-3,39.5*1E-3,38.7*1E-3 MPDATA,KXX,1,7,37.8*1E-3,36.8*1E-3,35.8*1E-3,34.8*1E-3,33.9*1E-3 !属性一：SA387 Gr12 （15CrMoR及锻件）的导热系数 !!!!!!!!机械场分析!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!重新进入前处理!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! /PREP7 ETCHG,TTS KEYOPT,1,3,1 mptemp,1,21,38,93,149,204,260 mptemp,7,316,371,426,482,538 mpdata,ex,1,1,2.068e5,2.055e5,2.027e5,1.993e5,1.965e5,1.937e5 mpdata,ex,1,7,1.903e5,1.875e5,1.841e5,1.806e5,1.772e5 mpdata,prxy,1,0.3 !N06625材料的弹性模量与泊松比 MPTEMP,,,,,,,, MPTEMP,1,21 MPTEMP,2,38 MPTEMP,3,93 MPTEMP,4,149 MPTEMP,5,204 MPTEMP,6,260 MPTEMP,7,316 MPTEMP,8,399 MPDATA,ALPX,1,,12.06*1E-6 MPDATA,ALPX,1,,12.24*1E-6 MPDATA,ALPX,1,,12.78*1E-6 MPDATA,ALPX,1,,12.96*1E-6 MPDATA,ALPX,1,,13.14*1E-6 MPDATA,ALPX,1,,13.32*1E-6 MPDATA,ALPX,1,,13.32*1E-6 MPDATA,ALPX,1,,13.5*1E-6 !N06625材料的热膨胀系数 |
B：极限载荷法中的本构模型：
采用的弹性-理想塑性本构模型。
双线性等向强化（BISO）：命令流如下：
| TB,BISO,1,1,2, !双线性等向强化----弹性斜率与塑性斜率 TBTEMP,150 TBDATA,,183*1.5,,,,, !极限载荷法---弹性-理想塑性模型!强化模量为0 |
图8：双线性等向强化参数输入后软件的显示结果
双线性随动强化（BKIN）：命令流如下：
| TB,BKIN,1,3 !材料1----双线性随动强化 TBTEMP,20,1 TBDATA,1,315,2.09E4 TBTEMP,500,2 TBDATA,1,212,1.72E4 TBTEMP,800,3 TBDATA,1,153,1.33E4 |
图9：双线性随动强化参数输入后软件的显示结果
等向强化与随动强化区别：
如果材料在一个方向屈服强度提高（强化）在其它方向的屈服强度也同时提高，这样的材料叫等向强化材料。形象比喻：吹气球。
如果材料在一个方向的屈服点提高，其它方向的屈服应力相应下降，比如拉伸的屈服强度提高多少，反向的压缩屈服强度就减少多少，这样的材料叫随动强化材料。形象比喻：捏气球。
具体用那种强化模型要看你的材料是那种材料。如果你只是单向加载，（即没有加载到屈服，卸载，再反向加载到屈服）两种材料模型的效果是一样的。
C: 弹塑性分析法中的本构模型：
采用的塑性本构关系----增量理论，输入材料的真实应力-应变曲线求解。真实的应力-应变曲线的获得需要依据Asme VIII-2 附录3-D所提出的考虑材料塑性应变强化的本构模型，只要知道材料的屈服强度、抗拉强度、弹性模量及断后延伸率（标准保证值和实测值），就可以获得材料的真实应力-应变曲线，用于弹塑性分析的有限元数值求解中。这需要编程实现。 真实的应力-应变曲线程序如图9所示。
图10：材料的真实应力-应变曲线程序
多线性等向强化（MISO）：命令流如下：
| TB,MISO,1,1,13 tbtemp,20 TBPT,,0.000315,E*0.000315 !曲线的第一个点必须与材料弹性模量相对应; !不允许有大于弹性模量或小于零的斜率段 TBPT,, 0.001687, 170 TBPT,, 0.002902, 200 TBPT,, 0.005249,230 TBPT,, 0.010116,260 TBPT,,0.019603,290 TBPT,,0.03262,320 TBPT,,0.045473,350 TBPT,,0.058066,380 TBPT,,0.071604,410 TBPT,,0.086723,440 TBPT,,0.103685,470 TBPT,,0.122628,600 tbplot,miso,1 !绘制材料1的miso曲线 |
图11：多线性等向强化参数输入（应力-应变数据）后软件的显示结果
图12：多线性等向强化材料的真实应力-应变曲线
MISO可以使用多线性来表示使用Von Mises屈服准则的等向强化的应力－应变曲线，它适用于比例加载的情况和大应变分析。
但是，应用这个模型有两点是应当注意的：
1、曲线的第一个点必须与材料弹性模量相对应；
2、不允许有大于弹性模量或小于零的斜率段。
多线性随动强化（MKIN）：命令流如下：
| sts1=173 !应力 stn1=0.002111 !应变 sts2=187 !应力 stn2=0.002732 !应变 sts3=201 !应力 stn3=0.003533 !应变 sts4=215 !应力 stn4=0.004542 !应变 sts5=229 !应力 stn5=0.005758 !应变 sts6=243 !应力 stn6=0.007117 !应变 tb,mkin,1 !多线性随动强化 ---输入真实的应力应变曲线 tbtemp,,strain tbdata,,stn1,stn2,stn3,stn4,stn5,stn6 tbtemp,, tbdata,,sts1,sts2,sts3,sts4,sts5,sts6 tbplot,mkin,1 !绘制材料1的mkin曲线 |
图13：多线性随动强化参数输入（应力-应变数据）后软件的显示结果
图14：多线性随动强化材料的真实应力-应变曲线
第一部分完成，未完待续。。。。。。