matlab论文12010245327马文建.doc

MATLAB课程论文利用MATLAB语言解决力学中非线性问题姓名马文健学号12010245327专业电子信息工程班级2010级电子班指导老师汤全武学院物电学院完成日期201112141利用MATLAB语言解决力学中非线性问题(马文健120102453272010级电子班)摘要MATLAB因其强大的功能而在世界范围内得以流行，特别在工程领域的应用甚为广泛。本文从弹性力学中的一个问题出发，先简要叙述并分析了目前厚壁筒弹塑性应力分析中所遇到的问题及其MATLAB在其中应用的思考；后建立了厚壁筒弹塑性区的模型，进行了各应力分析，用MATLAB中介绍的相关方法求解了弹塑性分界面半径并绘制了各应力关于半径的变化曲线，最后根据模型求解的结果进行了简要的分析。关键词MATLAB弹塑性应力厚壁筒一、提出问题厚壁筒在承受内压载荷的作用下，随着压力的增加，筒壁应力不断增加。厚壁筒在承受逐渐增加压力的过程中，会经历弹性阶段、筒体部分屈服阶段、整体屈服阶段、材料硬化、筒体过度变形、直至爆破失效阶段。而在分析厚壁筒的弹塑性应力分布时，遇到两个问题【1】弹塑性区分界面的半径的精确确定较为困难，因此半径是一个非线性方程；【2】弹塑性应力的分布曲线绘制难以精确，不能真实反映不同半径处的应力状况。二、分析问题针对厚壁筒弹塑性应力分析时遇到的两个问题，可利用MATLAB软件中的相关知识解决。对于弹塑性区分界面半径的确定，可利用方程求根的相关方法，从二分法、开方法、NEWTON法、NEWTON下山法以及弦截法中选取精度和收敛速度均较佳的方程求根数值方法进行求解；同时还可以利用MATLAB中的SOLVE函数求其精确解。对于弹塑性应力分布曲线的绘制，可以根据求出的各应力分布的具体方程，用PLOT函数进行绘制，同时要注意曲线后期的处理工作，以便更好地从图中得到各应力随半径分布的情况。1、模型建立及求解2图1受内压厚壁圆筒图2弹性区图3塑性区为简化分析，假设厚壁圆筒为理想弹塑性体，不考虑材料在塑性变形过程中塑性强化，筒体仅受内压PI作用，筒体的内半径为RI，外半径为RO。初始假设厚壁圆筒的内半径为RI1525MM，外半径为RO254MM，内壁所受的压力为PI340MPA，取其屈服极限为750MPA先分析筒体仅受内压PI作用的情形。当内压PI大于弹性极限压力PE时，圆筒内壁的屈服区向外扩展，筒体沿壁可以分成塑性区和弹性区两个区域，其中内侧为塑性区，外侧为弹性区。弹性区的内外半径分别为RC和RO，承受内压力为PC，见图2；假想两区域的交界圆面的半径为RC，则塑性区的内外半径分别为RI和RC，承受的内外压力分别为PI和PC，如图3所示。11塑性区(RIRRC)材料处于塑性状态时，筒壁微元体的平衡微分方程成立，即。01RD假设材料塑性变形时符合TRESCA屈服条件，代入上式，则有2RS，积分上是既有。易知有边界条件3SRDLN4RSA。将第一个条件代入式4。可求出A，并代,5IRIRP,6CCRRP入TRESCA屈服条件，可得到塑性区各应力分量的表达式12ZRLN750LN34075LN4103711LL68205LN750N34075N379RSIISIIZSIIPRRRRPR将第二个边界条件代入式7中可得塑性区交界面压力为。L10CCSIIRPP12弹性区(CORRR)3121应力分析弹性区内壁面即为弹塑性区交界面，即弹性区内壁面呈塑性状态。由式7、8、9可知在RRC有弹性区内壁面处各应力表达式为2020113CCCRRCRCZRRP弹性区内壁面开始屈服时，假设应力符合TRESCA屈服条件，则有。204SCCP由连续性可得式10和式4应为同一值，即有22022 5470LN750LN34015254L96115N07835CCCCCSIICCRRRRPXXFX122计算流程EMGU1预设迭代次数初值KK1迭代次数赋值F1,D1FF,XKXK1XKUF1/D1F2,D2FF,XK1WHILEABSF2ABSF1UU/2预设迭代次数赋值XK1XKUF1/D1F2,D2FF,XK1结束10X01KXKN输出迭代失败标志输出IX输出奇异标志5ENDEND(2)编写文件名为FUNC1M的函数文件FUNCTIONF,DFUNC1XF2542X2/25422LOGX/152509067D1’2542X2/25422LOGX/152509067’DSUBSDIFFD1(3)在命令窗口输入FFUNC1X,KMENDNEWTONF,180,0001124求解结果迭代次数KX值119605592196059231960594可知当迭代次数为3次时，有19605941960592000020001成立，故X的取值已符合要求，迭代成功。125计算验证运用MATLAB系统中的求解各类代数方程精确公式解的函数SOLVE求解此方程的精确解。SOLVE函数的调用格式为SOLVE‘方程或方程组’，未知量。编写程序如下SYMSXSOLVE 2LOGX/1525090672542X2/2542 ,X执行后屏幕显示ANS196060371可以看出由NEWTON下山法求解所得的值与精确解的差值为1960594196060370000970001，可知NEWTON下山法所求的解的精度已足够高。取196CRMM。126应力求解利用与塑性区各应力分析相同的方法，可求出弹性区各应力分量表达式6222702222027519654140139975196314SCRSCSCZRRRRR至此弹性区和塑性区各应力的表达式已全部求出。13应力分布曲线的绘制从上述的各应力关于半径R的表达式中可以看出，此六个方程皆具有一般的显式方程YFX的形式，可以利用曲线绘制法中的PLOT函数绘制。PLOT函数的调用格式为,, PLOTXY，其中X为曲线上点的横坐标，Y为纵坐标，’’中可包含确定曲线颜色、线形、两坐标轴上的比例等参数。131程序编写X1152501196取区间X219601254Y1750LOGX1411037曲线公式Y2750LOGX1336037Y3750LOGX1373537Y422329144107/X22Y522329144107/X22Y622329PLOTX1,Y1, R ,X1,Y2, B ,X1,Y3, C ,X2,Y4, R ,X2,Y5, B ,X2,Y6, C 输出三条有色曲线GRIDONAXIS0,260,400,600XLABEL 半径 YLABEL 应力 TITLE 厚壁筒弹塑性区应力分布曲线 GTEXT 轴向应力 GTEXT 周向应力 GTEXT 径向应力 BOXON说明GRIDON表示在图形窗口中画上网格线，以便更清楚分析各量的具体数值。GTEXT‘字符串’表示将字符串放在鼠标所指定的位置。XLABEL‘字符串’和YLABEL‘字符串’分别表示在X轴下方中部位置和Y轴左边中部位置输出字符串作为标题。PLOT括号后的R,B,C分