L组Project课程设计(改).doc

辽宁工程技术大学力学与工程学院Project报告书题 目工程中板壳的强度、刚度计算和稳定 性分析（弹塑性力学） 班 级 工程力学091班 姓 名李慧飞 李世全 李天舒 梁霜 指导教师 李利萍 成 绩 辽宁工程技术大学力学与工程学院 制project任务书Project题目：工程中板壳的强度、刚度计算和稳定性分析。Project主要内容：通过FORTRAN语言对正方形薄板进行了连续体平面问题的有限元分析；通过板壳理论基础对负重轮盘进行分析。并对初步的板壳理论的认识有了进一步的了解。学生姓名： 李慧飞 李世全 李天舒 梁霜 指导教师签字： 年 月 日摘 要早期的薄板断裂问题是根据Kirehhoff假设进行分析。1961年M.L.Williams利用特征值展开方法详细推导了弯曲时直线裂纹尖端的应力表达式。1965年件G.C.Sih将Hilbert边值问题推广至板内，求得了直线裂纹尖端的应力强度因子。应用复变函数理论和保角映射方法，苏联学者解决了板内含有各种曲线型缺陷的弯扭问题，这些成果已列入附录中我们也曾用复变函数方法计算唇型和翼型裂纹的应力强度因子。1971年W.K.Wilson和Thompson首次用有限元法分析薄板弯曲的断裂问题。与薄板断裂分析相仿，最初薄壳断裂分析亦普遍应用Kirehhoff理论，基于扁壳理论，E.S.Folias研究含裂纹的球壳问题应用叠加原理将球壳含裂纹问题分解为无裂纹问题和裂纹面受分布力问题。前者满足基本方程和远离裂纹面处的边界条件，而后者满足基本方程并与前者解叠加后在裂纹面上满足内力自由的条件通过Fourier积分变换，首先化为对偶积分, 由于核函数求积困难，故进一步转换成奇异积分方程，求得含裂纹球壳的裂尖应力表达式的首项和强度因子。为克服上述缺点, 近年来已有较多的研究者采用Reissner精化理论进行研究。J.K.Knowles等、R.J.Hartranft和G.C.Sih分别导得Reissner的奇异性，在文中得出对称工型问题裂纹尖端应力、应变场的展开式，解决了有限尺寸复合型弯曲断裂问题。综观自1961年至今的五十余年时间内, 含裂纹板壳的断裂分析已取得很大进展, 但是还有一些问题, 如板壳弹塑性断裂分析, 复合材料的层合壳断裂以及实用工程计算方法等方面需要进一步研究。工程上大量采用薄的板壳型结构，它们在压力作用下，会在内部应力远小于材料的屈服极限应力（见材料的力学性能）时，突然产生垂直于压力方向的位移而降低承载能力，甚至发生破坏，这种现象称为失稳、皱损或屈曲。板壳失稳是由侧向位移引起的，因而失稳属于刚度问题。由于研究板壳失稳问题就形成了板壳稳定性理论。关键词：板壳；强度；刚度；稳定性AbstractEarly sheet fracture problem is analyzed according to the Kirehhoff hypothesis. In 1961 M.L.W illiams using eigenvalue method is deduced in detail on bending line of the crack tip stress expression. In 1965 a G.C.S ih will Hilbert boundary value problem is to promote in the plate, get the linear stress intensity factor of the crack tip. Application of complex function theory and conformal mapping method, the Soviet union scholars solved plate contains various curve type defect bending and twisting problem, these results have been listed in the appendix, we also used the methods of complex function calculation lip type and profile crack stress intensity factor. In 1971, for the first time W.K.W ilson and Thompson by finite element method analysis sheet bending fracture problems.And sheet fracture analysis are similar, initial shell fracture analysis is also widely used Kirehhoff theory, based on the shallow shell theory, E.S.F olias research including crack of spherical shell applications of the superposition principle will be spherical shell with crack problem is decomposed into no crack and the crack surface by distribution force problem. The former meet basic equation and away from the crack surface in the boundary conditions, which meet the basic equation and the former solution stack in the crack surface after meet internal force free condition through the Fourier integral transform, first into the dual integral, due to the nuclear function quadrature difficulties, so further converted into a singular integral equation, get including crack spherical shell of the crack tip stress expression first term and intensity factor.In order to overcome the above shortcomings, in recent years many researchers using Reissner refine theory research. J.K.K nowles, R.J.H artranft and G.C.S ih respectively to guide Reissner singularity, in the work that symmetrical type problems crack tip stress and strain field expansion, solve the hybrid finite size bending fracture problems.Synoptic since 1961 twenty years time, including crack plate and shell of fracture analysis has made great progress, but there are still some problems, such as plate and shell elastic-plastic fracture analysis, composite laminated shell fracture and practical engineering calculation method for further study.Works on a large number of thin lamella-type structure, and they are under pressure, in the internal stress is far less than the yield stress limit of the material (see the mechanical properties of the material), and suddenly generate a displacement perpendicular to the direction of the pressure and reduce the carrying capacity, or even destroyed, a phenomenon known as instability, wrinkled loss or buckling. Plates and Shells instability is caused by the lateral displacement, and thus instability belongs stiffness problem. As the study of plate and shell instability problems on the plate and shell stability theory.Key words: shell;strength;stiffness;stability目 录1绪论12板壳力学基础知识22.1几项假设22.2 板壳的应力、应变以及应力与应变的关系22.3 薄板的内力22.4 薄壳的内力32.5 薄板的边界条件32.6 薄壳的边界条件32.7 板壳的应力计算公式33工程中实例分析53.1 连续体平面问题的有限元分析53.1.1 理论依据与分析53.1.2程序原理及实现63.1.3 程序的原理如框图63.1.4 算例原始数据和程序分析103.1.5算例结果213.2 负重轮轮盘强度和刚度的板壳理论分析234.结论27参考文献 28工程中板壳的强度、刚度计算和稳定性分析1绪论有限单元法作为一门课程在现实工程中的应用已经十分的广泛，我组先对板壳力学的基本知识进行简单的学习，然后结合简单的正方形薄板为例用软件进行计算和分析，最后在现实生活中发现实例对实例进行计算并与理论进行比较得出结论。本文第一个例子利用计算机，结合FORTRAN语言和有限单元法的理论知识，对平面应力问题的薄板进行内力值的计算和相应的结构分析，并可以将由计算机计算的实际结果与由用弹性力学理论知识计算出来的解答相比较。第二个例子对负重轮轮盘进行了模型简化，运用板壳理论建立并求解了变形微分方程，得出轮盘内各点的应变和应力分量，并进一步对轮盘进行强度、刚度校核和设计，其结果与有限元分析结果接近。从而有效的对板壳的强度，刚度稳定性进行了一定的研究。2. 板壳力学基础知识2.1 几项假设1.板壳是均匀的、连续的，并且是各向同性的；2.板壳是线弹性的；3.板壳的变形是微小的；4.直法线假设，即认为板壳变形前垂直于中面的法线线段在变形后仍保持为直线，并垂直于变形后的中面，且其长度不变。5.法向应力很小，可以忽略；6.板的中面没有变形。2.2 板壳的应力、应变以及应力与应变的关系薄板壳内任一点沿z方向的位移wA与坐标z无关，仅是坐标x、y的函数，横向剪应变yz和zx应为零， 几何方程 （1-1）物理方程 （1-2） 或 （1-3）2.3 薄板的内力 （1-4）2.4 薄壳的内力 （1-5）2.5 薄板的边界条件1.简支边 ， （1-6）2.固定边 ， （1-7）3.自由边 ， （1-8）2.6 薄壳的边界条件1.简支边 ， ， （1-9）2.固定边 ， ， （1-10）3.自由边 ， （1-11）2.7 板壳的应力计算公式1.薄板 （1-12）其中惯性矩，静面矩 2.薄壳 （1-13）注意，在板壳弯曲问题中，数值上最大的是法向应力、和切向应力，因而是主要的应力，横向剪应力、数值较小，是次要的应力，一般说来，无须对它们进行计算。3工程中实例分析3.1 连续体平面问题的有限元分析正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，P=1kN/m，同时在沿对角线y轴上受一对集中压力，载荷为2kN，板厚t=1，泊松比=0，见下图：Y 2kN P=1kN/mP X o P P 2kN 图2.1正方形薄板四周受均匀载荷的作用Figure 2.1 square sheet around by the uniform load role3.1.1 理论依据与分析此问题，为弹性力学里的平面应力问题，在板的内部，到处都有z=0,yz=0,xz=0;x=f1(x,y),y=f2(x,y),xy=f3(x,y)，应力具有这种性质的问题，称为平面应力问题。弹性薄板在工程中应用很广泛，对于一些简单的情况，如等厚、单跨、无大孔口，外形规则（如矩形，圆形等）的薄板，已有一些解答和表格可资利用盘。由于连续平板的连续性，仅需要取其在第一象限的四分之一部分研究计算，然后做出一些辅助线将平板分成若干部分，在为每个部分选择分子单元，采用此模型化为4个全等的三角形单元，利用其对称性，四分之一的边界约束，荷载可等效如图所示。 1kN/m图2.2 三角形单元的内力Figure 2.2 triangle unit internal force3.1.2 程序原理及实现 用FORTRAN程序的实现。有节点信息文件NODE.IN和单元信息文件ELEMENT.IN,经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件BASIC.IN生成。该程序的特点如下：问题类型：可用于计算弹性力学平面应力问题和平面应变问题。单元类型：采用常应变三角形单元。位移模式：用线性位移模式。载荷类型：节点载荷，非节点载荷应先转换为等效节点载荷。材料性质：弹性体由单一的均匀材料组成。约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有三个自由度的独立约束。方程求解：针对半带宽刚度方程的GUASS消去法。输出文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT,IN。结果文件：输出一般的结果文件DATA.OUT。3.1.3 程序的原理如框图开始 输入数据（子程序READ_IN）BASIC.IN(基本信息文件)NODE.IN(节点信息文件)ELEMEENT.IN(单元信息文件)形成单元刚度矩阵（子程序FORM_KE）以半带存储方式形成整体刚度矩阵（BAND K）形成节点载荷向量（子程序FORM_P） 处理边界条件（子程序DO_BC） 求解方程获得节点位移（子程序SOLVE）计算单元及节点应力（子程序） 输出文件DATA.OUT 结束说明：（1）主要变量：ID：问题类型码：，ID=1时的平面应力问题，ID=2时的平面应变问题N_NODE： 节点个数N_LOAD： 节点荷载个数N_DOF： 自由度，N_DOF=N_NODE*2（平面问题）N_ELE： 单元个数N_BAND： 矩阵半带宽N_BC： 有约束的节点个数PE： 弹性模量PR： 泊松比PT： 厚度LJK_ELE(I，3)：单元节点编号数组，LJK_ELE(I,1),LJK_ELE(I,2),LJK_ELE(I,3)分别放单元I的三个节点的整体编号。X（N_NODE）,Y(N_NODE):节点坐标数组，X(I),Y(I)分别存放节点I的x,y坐标值。LJK_U(N_BC,3)：节点载荷数组，P_LJK(I,1)表示第I个作用有节点载荷的节点的编号，P_LJK(I,2),P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值。AK(N_DOF,N_BAND): 整体刚度矩阵AKE(6,6): 单元刚度矩阵BB(3,6): 位移应变转换矩阵（三节点的几何矩阵）DD(3,3): 弹性矩阵SS(3,6): 应力矩阵RESULT_N(N_NOF):节点荷载数组，存放节点荷载向量，解方程后该矩阵存放节点位移DISP_E(6): 单元的节点位移向量STS_ELE(N_ELE,3): 单元的应力分量STS_ND(N_NODE,3): 节点的应力分量（2） 子程序说明：READ_IN: 读入数据 BAND_K: 形成半带宽的整体刚度矩阵FORM_FE: 计算单元刚度矩阵 FORM_P: 计算节点载荷CAL_AREA: 计算单元面积 DO_BC: 处理边界条件CLA_DD: 计算单元弹性矩阵 SOLVE: 计算节点位移CLA_BB: 计算单元位移.应变关系矩阵CAL_STS: 计算单元和节点应力（3） 文件处理：源程序文件：chengxu,for 程序读入的数据文件：BASIC.IN,NODE.IN,ELEMENT.IN（需要手工生成）程序输出的数据文件：DATA.OUT(4) 数据文件格式：需读入的模型、基本信息文件BASIC.IN的格式如下表：栏目格式说明实际需输入的数据基本模型数据第1行，每两个数之间用“，”号隔开问题类型，单元个数，节点个数，有约束的节点数，有载荷的节点数材料性质第2行，每两个数之间用“，”号隔开弹性模量，泊松比，单元厚度节点约束信息在材料性质输入行之后另起行，每两个数之间用“，”号隔开LJK_U(N_BC,3)位移约束的节点编号，该节点x方向约束代码，该节点y方向代码节点载荷信息在节点约束信息输入行之后另起行，每两个数之间用“，”号隔开P_IJK(N_LOAD,3)载荷作用的节点编号，该节点x方向载荷该节点y方向载荷， 需读入的节点信息文件NODE.IN的格式如下表栏目格式说明实际需输入的数据节点信息每行为一个节点的信息（每行三个数，每两个数之间用空格或“，”分开LJK_U(N_BC,3)节点号，该节点的x坐标，该节点y方向坐标需读入的单元信息文件ELEMENT.IN的格式如下表栏目格式说明实际需输入的数据单元信息每行为一个单元的信息（每行有14个整型数4个为单元节点编号，对于3节点编号，第4个节点编号与第3个节点编号相同，后10个数无用，可输入“0”，每两个整型数之间用至少一个空格分开NE_ANSYS(N_ELE,14)单元的节点号1（空格）单元的节点号2（空格）单元的节点号3（空格）单元的节点号4（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0输出结果文件DATA.OUT格式如下表栏目实际需输入的数据节点位移I RESULT_N(2*I_1) RESULT_N(2*I)节点号 x方向位移 y方向位移单元应力的三个分量IE STE_ELE(IE,1) STE_ELE(IE,2) STE_ELE(IE,3)单元号 x方向应力 y方向应力 剪切应力节点应力的三个分量I STS_ND(I,1) STS_ND(I,2) STS_ND(I,3)节点号 x方向应力 y方向应力 剪切应力3.1.4 算例原始数据和程序分析（1） 模型基本信息文件BASIC.IN的数据为1,4,6,5,31,0,1.1，1,0,2,1,0,4,1,1,5,0,1,6,0,11，-0.5，-1.5,3，-1，-1,6，-0.5，-0.5（2） 手工准备的节点信息文件NODE.IN的数据为1 0.0 2.02 0.0 1.03 1.0 1.0 4 0. 0.5 1.0 0.6 2.0 0.（3） 手工准备的单元信息文件ELEMENT.IN的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 0 4（4） 源程序文件chengxu,for为：PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)DIMENSION STS_ELE(500,100),STS_ND(500,3)OPEN(4,FILE=BASIC.IN)OPEN(5,FILE=NODE.IN)OPEN(6,FILE=ELEMENT.IN)OPEN(8,FILE=DATA.OUT)OPEN(9,FILE=FOR_POST.DAT)READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,=PLANE STRESS PROBLEM=)25 FORMAT(/5X,=PLANE STRAINPROBLEM=)CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,& IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK)CALL FROM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N& STS_ELE,STS_ND)C to putout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),&STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMAT(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1),REAL(IJK_ELE(I,2).&REAL(IJK_ELE(I,3)REAL(IJK_ELE(I,3),&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1,N_ELE)72 FORMAT(7F9.4)C CLOSE(4)CLOSE(5)CLOSE(8)CLOSE(9)ENDCC to get the original data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,&PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),& P_IJK(N_LOAD,3),NE,ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,3)READ(4,*)PE,PR,PTREAD(4,*)(IJK_U(I,J),J=1,3),I=1,N_BC)READ(5,*)(ND_ANSYS(I,J),J=1,3),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUE DO 11 I=1,N_ELEDO 11 J=1,3 IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELE DO 20 I=1,3 DO 20 J=1,3IW=IABS(IJKK_ELE(IE,I)-IJK_ELE(IE,J)IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUE N=BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFRETURN ENDCC to form the stiffness matrix of elementSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6),SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3 DO 10 J=1,6 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J) DO 20 I=1,6DO 20 J=1,6 AKE(I,J)=0.0 DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT RETURNENDC to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE,& PR,PT,AK) DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100) N_DOF=2*N_NODEDO 40 I=1,N_DOFDO 40 J=1,N_BAND40 AK(I,J)=0 DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I=1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL.LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 CONTINUERETURNENDCC to calculate the area of element SUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YJJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNENDCC to calculate the elastic matrix of element SUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3 DO 10 J=1,310 DD(I,J)=0.0 DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/(1.0+PR)*2.0)RETURNENDCC to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3 DO 10 JJ=1,310 BB(II JJ)=0.0 BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3 DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE) RETURNENDCC to form the global load matrix SUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,IJK_ELE,X,Y,P_IJK,& RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3),& RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3) RETURN ENDCC to deal with BC(u) (here only for fixed displacement) using 1-0 methodSUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO 30 I=1,N_BC IR=IJK_U(I,1) DO 30 J=2,3IF(IJK+U(I,J).EQ.0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N=BAND10 AK(II,JJ)=0.0 DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0 END IF30 CONTINUERETURNENDCC to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)DIMENTION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF,GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF,LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IM L=I-K+1 C=AK(K,L)/AK(K,1) IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K) RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1 I=N_DOF-I1 IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1 IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQ K=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1) WRITE(8,50)50 FORMAT(/12X,*RESULT BY FEM2D *,/8X,&-DISPLACEMENT OF NODE-/5X,NODE NO,8X,X-DISP,8X,Y-DISP)DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6) RESULTENDCC calculate the stress component of elemen