非线性有限元作业_老骆整理.docx

轴对称问题的弹塑性分析开始读入单元号，B矩阵位数，单刚位数，单元节点号，单刚等各项参数EN1存储单元节点号，局部坐标系转换为全局坐标JD0打印错误 调用子函数DEMATR求[D]调用子函数BMATR求[B][S]==切线刚度阵[EK]=[S][Q1]=·JD·RN·H(I1)H(J1)结束NYY流程图：R=1[C]=[De][B]NY[C]=[Dep][B]返回各值解析解。厚壁筒受内压，采用Mises屈服准则 (1)经计算知，当t=12.5（12.91）时，材料处于弹塑性交界面。弹性区为: (2) (3)塑性区： (4) (5) 交界处有：， (6)最后解得残余应力为： (7a) (7b) (8a) (8b)有限元网格信息图：图1有限元网格 输入数据文件内容(详细信息见附件):DATA(1)NNODE MELEM IFU IFW IPF IPR NPP NRM HAC MSF NULOAD EXPNM(1-MELEM)NNNN(1-NNODE) R ZNFU(1-IFU) FUNFW(1-IFW) FWMPQ(1-IPF) NPQ*PQNPRNRZ(1-IPR) PRNRZE EMU SSS HH UNLOAD对理想塑性材料厚壁筒，从初始状态开始，历经加载后完全卸载。这一过程中，厚壁筒内会产生残余应力。沿径向R的残余应力如图2-3所示。图2 径向残余应力-半径曲线图3 切向残余应力-半径曲线 图2-3中分别给出了径向残余应力和切向残余应力随半径的变化，并且和解析解进行了比较。从图中可以看出，程序解和解析解在数值上能够很好的吻合，只是在径向残余应力最大的地方有少许偏差，这验证了程序计算结果的正确性。最大误差发生在径向残余应力达到最大的地方，达3.5%。因此，程序计算结果能够得到精确度比较高的解。对于不同线性硬化材料的厚壁圆筒，硬度对应力分布的影响如图4-5所示。图4 硬度分别为0.5和0.8时的径向应力曲线图5 硬度分别为0.5和0.8时的径向应力曲线如图4所示，随着硬度的增加，径向应力随之减小，而最大径向应力不在厚壁筒的表面。图5中显示了不同硬度下，切向应力随半径的变化。两种不同硬度下的应力曲线有交点，这说明硬度越大，弹性应力范围也就越大。弹塑性的分界层位于应力突变点。变形首先是弹性变形，应力的增量基本呈线性而且斜率比较大，这是弹性区；塑性变形区，应力基本保持为水平线。弹塑性分界层的位置就是图5中曲线的突变点，不同的硬度的分界位置会有少许差别。 轴对称问题的几何非线性分析开始读入单元号，B矩阵位数，单刚位数，单元节点号，单刚，大变形解法选项等参数单元内部自由度转换为全局坐标JD0错误 EN程序终止Q1=JD*RN*H(I1)H(J1)调用子函数BMATR求[B]清零BLIMATRIX阵计算[C]=[DE][BB],[CC]=[M][G][S]=,SS= [EK]=*JD*RN*H(I1)H(J1)+*JD*RN*H(I1)H(J1)结束NY流程图：Tu=1调用子函数BLIMATR求大位移刚度阵BL1并存入BL1MATR1XYN调用子函数GMATR求[G]矩阵并存入GMATR1X调用子函数MMATR求[M]矩阵并存入MMATR1X将[BL]矩阵加上[B]累加至[BB]中 有限元网格信息图（输入文件在附件中给出）：图6 有限元网格输入数据文件(详细信息见附件)：DATA(3)NNODE MELEM IPU IFW IPF IPR NPP EXPNM(1-MELEM) NNNN(1-NNODE) R ZNFU(1-IFU) FUNFW(1-IFW) FWMPQ(1-IPF) NPQ*PQNPRNRZ(1-IPR) PRNRZE EMU周边固支受均布载荷作用的圆板，中心点的载荷-位移曲线如图5所示。图7 圆板中心点的载荷-位移曲线如图7中所示，线性解和非线性解在很大程度上具有一致性。在初始的小变形情况下，线性解和非线性解结果是一致的；随着变形的增大，小变形线性假设已不再适用，就会产生几何分线性问题，导致了线性解和非线性解的偏差。线性解要大于非线性解，随着变形的增大，非线性程度的增加，两种结果的差值也会增大。对均布载荷为10的情况，非线性解为-1.229，TL法和UL法的误差分别为1.06%和0.81%。附件第一题输入文件（data1）43 8 0 34 1 0 10 2 1 1.0 0.011e-61 1 4 6 7 8 5 3 22 6 9 11 12 13 10 8 73 11 14 16 17 18 15 13 124 16