L组Project课程设计

1、辽宁工程技术大学力学与工程学院Project报告书工程中板壳的强度、刚度计算和稳定题目性分析弹塑性力性班级工程力学091班姓名李慧飞李世全李天舒梁霜指导教师尝匡成绩辽宁工程技术大学力学与工程学院制project任务书Project题目：工程中板壳的强度、刚度计算和稳定性分析。Project主要内容：通过FORTRAN语言对正方形薄板进行了连续体平面问题的有限元分析；通过板壳理论基础对负重轮盘进行分析。并对初步的板壳理论的认识有了进一步的了解。学生：李慧飞李世全李天舒梁霜指导教师签字：年月日摘要早期的薄板断裂问题是根据Kirehhoff假设进行分析。1961年M.L.Williams利用特征值展

2、开方法详细推导了弯曲时直线裂纹尖端的应力表达式。1965年件G.C.Sih将Hilbert边值问题推广至板内，求得了直线裂纹尖端的应力强度因子。应用复变函数理论和保角映射方法，苏联学者解决了板内含有各种曲线型缺陷的弯扭问题，这些成果已列入附录中我们也曾用复变函数方法计算唇型和翼型裂纹的应力强度因子。1971年W.K.Wilson和Thompson!次用有限元法分析薄板弯曲的断裂问题。与薄板断裂分析相仿，最初薄壳断裂分析亦普遍应用Kirehhoff理论，基于扁壳理论，E.S.FoliasW究含裂纹的球壳问题应用叠加原理将球壳含裂纹问题分解为无裂纹问题和裂纹面受分布力问题。前者满足基本方程和远离裂

3、纹面处的边界条件，而后者满足基本方程并与前者解叠加后在裂纹面上满足内力自由的条件通过Fourier积分变换，首先化为对偶积分，由于核函数求积困难，故进一步转换成奇异积分方程，求得含裂纹球壳的裂尖应力表达式的首项和强度因子。为克服上述缺点，近年来已有较多的研究者采用Reissne精化理论进行研究。J.K.Knowles等、R.J.Hartranf的G.C.Sih分别导得ReissnerB勺奇异性，在文中得出对称工型问题裂纹尖端应力、应变场的展开式，解决了有限尺寸复合型弯曲断裂问题。综观自1961年至今的五十余年时间内，含裂纹板壳的断裂分析已取得很大进展，但是还有一些问题，如板壳弹塑性断裂分析，复

4、合材料的层合壳断裂以及实用工程计算方法等方面需要进一步研究。工程上大量采用薄的板壳型结构，它们在压力作用下，会在内部应力远小于材料的屈服极限应力见材料的力学性能时，突然产生垂直于压力方向的位移而降低承载能力，甚至发生破坏，这种现象称为失稳、皱损或屈曲。板壳失稳是由侧向位移引起的，因而失稳届于刚度问题。由于研究板壳失稳问题就形成了板壳稳定性理论。关键词：板壳；强度；刚度；稳定性AbstractEarlysheetfractureproblemisanalyzedaccordingtotheKirehhoffhypothesis.In1961M.L.Williamsusingeigenvaluem

5、ethodisdeducedindetailonbendinglineofthecracktipstressexpression.In1965aGC.SihwillHilbertboundaryvalueproblemistopromoteintheplate,getthelinearstressintensityfactorofthecracktip.Applicationofcomplexfunctiontheoryandconformalmappingmethod,theSovietunionscholarssolvedplatecontainsvariouscurvetypedefec

6、tbendingandtwistingproblem,theseresultshavebeenlistedintheappendix,wealsousedthemethodsofcomplexfunctioncalculationliptypeandprofilecrackstressintensityfactor.In1971,forthefirsttimeW.K.WilsonandThompsonbyfiniteelementmethodanalysissheetbendingfractureproblems.Andsheetfractureanalysisaresimilar,initi

7、alshellfractureanalysisisalsowidelyusedKirehhofftheory,basedontheshallowshelltheory,E.S.Foliasresearchincludingcrackofsphericalshellapplicationsofthesuperpositionprinciplewillbesphericalshellwithcrackproblemisdecomposedintonocrackandthecracksurfacebydistributionforceproblem.Theformermeetbasicequatio

8、nandawayfromthecracksurfaceintheboundaryconditions,whichmeetthebasicequationandtheformersolutionstackinthecracksurfaceaftermeetinternalforcefreeconditionthroughtheFourierintegraltransform,firstintothedualintegral,duetothenuclearfunctionquadraturedifficulties,sofurtherconvertedintoasingularintegraleq

9、uation,getincludingcracksphericalshellofthecracktipstressexpressionfirsttermandintensityfactor.Inordertoovercometheaboveshortcomings,inrecentyearsmanyresearchersusingReissnerrefinetheoryresearch.J.K.Knowles,R.J.HartranftandG.C.SihrespectivelytoguideReissnersingularity,intheworkthatsymmetricaltypepro

10、blemscracktipstressandstrainfieldexpansion,solvethehybridfinitesizebendingfractureproblems.Synopticsince1961twentyyearstime,includingcrackplateandshelloffractureanalysishasmadegreatprogress,buttherearestillsomeproblems,suchasplateandshellelastic-plasticfractureanalysis,compositelaminatedshellfractur

11、eandpracticalengineeringcalculationmethodforfurtherstudy.Worksonalargenumberofthinlamella-typestructure,andtheyareunderpressure,intheinternalstressisfarlessthantheyieldstresslimitofthematerial(seethemechanicalpropertiesofthematerial),andsuddenlygenerateadisplacementperpendiculartothedirectionofthepr

12、essureandreducethecarryingcapacity,orevendestroyed,aphenomenonknownasinstability,wrinkledlossorbuckling.PlatesandShellsinstabilityiscausedbythelateraldisplacement,andthusinstabilitybelongsstiffnessproblem.Asthestudyofplateandshellinstabilityproblemsontheplateandshellstabilitytheory.Keywords:shell;st

13、rength;stiffness;stability绪论1板壳力学基础知识22.1几项假设22.2板壳的应力、应变以及应力与应变的关系22.3薄板的内力22.4薄壳的内力32.5薄板的边界条件32.6薄壳的边界条件32.7板壳的应力计算公式3工程中实例分析53.1连续体平面问题的有限元分析53.1.1理论依据与分析53.1.2程序原理及实现-63.1.3程序的原理如框图-63.1.4算例原始数据和程序分析103.1.5算例结果213.2负重轮轮盘强度和刚度的板壳理论分析234.结论27参考文献28工程中板壳的强度、刚度计算和稳定性分析绪论有限单元法作为一门课程在现实工程中的应用已经十分的广泛，

14、我组先对板壳力学的基本知识进行简单的学习，然后结合简单的正方形薄板为例用软件进行计算和分析，最后在现实生活中发现实例对实例进行计算并与理论进行比较得出结论。本文第一个例子利用电脑，结合FORTRAN语言和有限单元法的理论知识，对平面应力问题的薄板进行内力值的计算和相应的结构分析，并可以将由电脑计算的实际结果与由用弹性力学理论知识计算出来的解答相比较。第二个例子对负重轮轮盘进行了模型简化，运用板壳理论建立并求解了变形微分方程，得出轮盘内各点的应变和应力分量，并进一步对轮盘进行强度、刚度校核和设计，其结果与有限元分析结果接近。从而有效的对板壳的强度，刚度稳定性进行了一定的研究。1. 板壳力学基础知

15、识2.1几项假设1. 板壳是均匀的、连续的，并且是各向同性的；2. 板壳是线弹性的；3. 板壳的变形是微小的；4. 直法线假设，即认为板壳变形前垂直丁中面的法线线段在变形后仍保持为直线，并垂直丁变形后的中面，且其长度不变。5. 法向应力很小，可以忽略；6. 板的中面没有变形。2.2板壳的应力、应变以及应力与应变的关系薄板壳内任一点沿z方向的位移wA与坐标z无关，仅是坐标x、y的函数,横向剪应变yz和zx应为零,几何方程u1-1xxVyxy物理方程xyi(i(y2(1Ey)x)xy1-2xy9(2(x1E(2yE2(1)xyy)x)1-3h22.3薄板的内力h2xzdzh2MyyzdzMxyh2

16、h2xyzdzh21-4xzdzQyh2yzdz2.4薄壳的内力1xxJyh/2Mxyh/2Myxh/2h/2h/2xyzdzh/2QxxzdzQyyh/2h/2h/2h/2NxxdzNyh/2h/2h/2NxyNyxxydzh/2h/2h/2yzdzdzydz2.5薄板的边界条件1.简支边Mx2.固定边VxQx3.自由边Mxyy2.6薄壳的边界条件1.简支边xa0wxa0Nxxa0Mxvxa02.固定边3.自由边Nxxao,vxa°,w0xaVxxaQxa0Mxxa0uxaMxyMxySyTxxaNxyry2.7板壳的应力计算公式1.薄板Mxbxz,xIMxybxyIz,xzMyb

17、丁乙QxSI,yzQySI1-51-61-71-81-91-101-111-12其中惯性矩1bh3/12，静面矩Sh242.薄壳Nx12Mxx''2z,hhz0,xy6Qyyz6QxxzNxyhh4h24Ny12MyTh2z12MxyWz1-13注意，在板壳弯曲问题中，数值上最大的是法向应力x、y和切向应力xy,因而是主要的应力，横向剪应力yz、遂数值较小，是次要的应力，一般说来，无须对它们进行计算。工程中实例分析3.1连续体平面问题的有限元分析正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力,P=1kN/m，同时在沿对角线y轴上受一对集中压力，载荷为2kN,板厚t

18、=1,泊松比V=0见下歹0图:图2.1正方形薄板四周受均匀载荷的作用Figure2.1squaresheetaroundbytheuniformloadrole3.1.1理论依据与分析此问题，为弹性力学里的平面应力问题，在板的内部，到处都有定=0,yz=0,xz=0;x=f1(x,y),y=f2(x,y),xy=f3(x,y),应力具有这种性质的问题，称为平面应力问题。弹性薄板在工程中应用很广泛，对丁一些简单的情况，如等厚、单跨、无大孔口，外形规则如矩形，圆形等的薄板，已有一些解答和表格可资利用盘。由丁连续平板的连续性，仅需要取其在第一象限的四分之一部分研究计算，然后做出一些辅助线将平板分成假

19、设干部分，在为每个部分选择分子单元,采用此模型化为4个全等的三角形单元，利用其对称性，四分之一的边界约束，荷载可等效如下图图2.2三角形单元的内力Figure2.2triangleunitinternalforce3.1.2程序原理及实现用FORTRAN程序的实现。有节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件BASIC.IN生成。该程序的特点如下：问题类型：可用丁计算弹性力学平面应力问题和平面应变问题。单元类型：采用常应变三角形单元。位移模式：用线性位移模式。载荷类型：节点载荷，非节点载荷应先转换为等效节

20、点载荷。材料性质：弹性体由单一的均匀材料组成。约束方式：为“0位移固定约束，为保证无刚体位移，弹性体至少应有三个自由度的独立约束。方程求解：针对半带宽刚度方程的GUASS消去法。输出文件：由手工生成节点信息文件NODE.IN,和单元信息文件ELEMENT,IN。结果文件：输出一般的结果文件DATA.OUT。3.1.3程序的原理如框图开始输入数据子程序READ_INBASIC.IN（基本信息文件）NODE.IN（节点信息文件）ELEMEENT.IN（单元信息文形成单元刚度矩阵子程序FORMKE以半带存储方式形成整体刚度矩阵BANDK形成节点载荷向量子程序FORMP处理边界条件子程序DOBC求解方

21、程获得节点位移子程序SOLVE计算单元及节点应力子程序输出文件DATA.OUT结束说明：1主要变量:ID:问题类型码：,ID=1时的平面应力问题，ID=2时的平面应变问题节点个数NNODE:NLOAD:NDOF:自由度，N_DOF=N_NODE*2平面问题N_ELE：单元个数NBAND:矩阵半带宽NBC:有约束的节点个数PE:弹性模量PR:泊松比PT:厚度LJK_ELE(I,3):单元节点编号数组,LJK_ELE(I,1),LJK_ELE(I,2),LJK_ELE(I,3)分别放单元I的三个节点的整体编号XN_NODE，Y(N_NODE):节点坐标数节点荷载个数组，X(I),Y(I)分别存放节

22、点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):节点的应力分量(

23、2) 子程序说明:READIN:读入数据BANDK:形成半带宽的整体刚度矩阵FORM_FE:计算单元刚度矩阵FORM_P:计算节点载荷CAL_AREA:计算单兀面枳DO_BC:处理边界条件CLA_DD:计算单元弹性矩阵SOLVE:计算节点位移CLA_BB:计畀单兀位移.应尖关系矩阵CALSTS:算单兀和|小、应力(3) 文件处理：源程序文件：chengxu,for程序读入的数据文件:BASIC.IN,NODE.IN,ELEMENT.IN需要手工生成程序输出的数据文件：(4) DATA.OUT数据文件格式：需读入的模型、基本信息文件BASIC.IN的格式如下表：栏目格式说明实际需输入的数据基本模

24、型数据第1行，每两个数之间用:”号隔J1问题类型，单元个数，节点个数，有约束的节点数，有载荷的节点数材料性质第2行，抑个数之间用：”号隔J1弹性模量，泊松比，单元厚度节点约束信息在材料性质输入行之后另起行，每两个数之间用；”号隔升LJK_U(N_BC,3)位移约束的节点编号，该节点x方向约束代码，该节点y方向代码节点载荷信息在节点约束信息输入行之后另起行，每两个数之间用：”号隔开P_IJK(N_LOAD,3)载荷作用的节点编号，该节点x方向载荷该节点y方向载荷，需读入的节点信息文件NODE.IN的格式如下表栏目格式说明实际需输入的数据节点信息每行为一个节点的信息每行三个数，每两个数之间用空格或

25、：”分开LJK_U(N_BC,3)节点号，该节点的x坐标，该节点y方向坐标需读入的单元信息文件ELEMENT.IN的格式如下表栏目格式说明实际需输入的数据单每行为一个单元的信息每行有14个整NE_ANSYS(N_ELE,14)元型数4个为单"点编号，对丁3节点单元的节点号1空格信编号,第4个节点编号与第3个节点编单元的节点号2空格息号相同，后10个数无用，可输入“0；'每两个整型数之间用至少一个空格分开单元的节点号3空格单元的节点号4空格0空格0空格0空格0空格0空格0空格0空格0空格0空格0输出结果文件DATA.OUT格式如下表栏目实际需输入的数据节点位移IRESULT_N

26、(2*I_1)节点号x方向位移RESULT_N(2*I)y方向位移单元应力的三个分量IESTE_STE_ELE(IE,3)单元号x方向应力ELE(IE,1)y方向应力STE_ELE(IE,2)剪切应力节点应力的三个分ISTS_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的数据

27、为1.010.02.00.031.01.040.0.51.00.62.00.(3)手工准备的单元信息文件ELEMENT.IN的数据为12330000111101245500001111025322000011110335660000111104(4)源程序文件chengxu,for为：PROGRAMFEM2DDIMENSIONIJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)DIMENSIONSTS_ELE(500,100),STS_ND(500,3)OPEN(4,FILE=&#

28、39;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)20FORMAT(/5X,'=PLANESTRESSPROBLEM=')25FORMAT(/5X,'=PLANESTRAIN

29、PROBLEM=')CALLREAD_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,&IJK_ELE,X,Y,IJK_U,P_IJK)CALLBAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,&IJK_ELE,X,Y,PE,PR,PT,AK)CALLFROM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,RESULT_N)CALLDO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALLSOLVE(N_NODE,N_DOF

30、,N_BAND,AK,RESULT_N)CALLCAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N&STS_ELE,STS_ND)CtoputoutadatafileWRITE(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)(RE

31、AL(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)CCLOSE(4)CLOSE(5)CLOSE(8)CLOSE(9)ENDCCtogettheoriginaldatainordertomodeltheproblemSUBROUTINEREAD_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,&PT,IJK_ELE

32、,X,Y,IJK_U,P_IJK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),&P_IJK(N_LOAD,3),NE,ANSYS(N_ELE,14)REALND_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)DO10I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)CONTINUEDO11I=1,N_E

33、LEDO11J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)10 CONTINUEN_BAND=0DO20IE=1,N_ELEDO20I=1,3DO20J=1,3IW=IABS(IJKK_ELE(IE,I)-IJK_ELE(IE,J)IF(N_BAND.LT.IW)N_BAND=IW20CONTINUEN=BAND=(N_BAND+1)*2IF(ID.EQ.1)THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFRETURNENDCCtoformthestiffnessmatrixofelementSUBROUTINEFORM_KE(IE,N_

34、NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),&AKE(6,6),SS(6,6)CALLCAL_DD(PE,PR,DD)CALLCAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO10I=1,3DO10J=1,6SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO20I=1,6DO20J=1,6AKE(I,J)=0.0DO20K=1,320AKE(I,J)=AKE(I,J)+SS(K,I)*BB

35、(K,J)*AE*PTRETURNENDCtoformbandedglobalstiffnessmatrixSUBROUTINEBAND_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_NODEDO40I=1,N_DOFDO40J=1,N_BAND40AK(I,J)=0DO50IE=1,N_ELECALLFORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,P

36、E,PR,PT,AKE)DO50I=1,3DO50II=1,2IH=2*(I=1)+IIDO50J=1,3DO50JJ=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)ENDIF50CONTINUERETURNENDCCtocalculatetheareaofelementSUBROUTINECAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSIONIJK_ELE(500,3),X(N_NODE),

37、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.0RETURNENDCCtocalculatetheelasticmatrixofelementSUBROUTINECAL_DD(PE,PR,DD)DIMENSIONDD(3,3)DO10I=1,3DO10J=1,310DD(I,J)=0.0DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD

38、(1,1)DD(3,3)=PE/(1.0+PR)*2.0)RETURNENDCCtocalculatethestrain-displacementmatrixofelementSUBROUTINECAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DIMENSIONIJK_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)DO10II=1,3DO10JJ=1,310BB(IIJJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y

39、(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)CALLCAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO20J1=1,620BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDCCtoformthegloballoadmatrixSUBROUTINEFORM_P(N_ELE,N_

40、NODE,N_LOAD,N_DOF,IJK_ELE,IJK_ELE,X,Y,P_IJK,&RESULT_N)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3),&RESULT_N(N_DOF)DO10I=1,N_DOF10RESULT_N(I)=0.0DO20I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20RESULT_N(2*II)=P_IJK(I,3)RETURNENDCCtodealwithBC(u)(hereonlyforfixeddisplace

41、ment)using"1-0"methodSUBROUTINEDO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)DIMENSIONRESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO30I=1,N_BCIR=IJK_U(I,1)DO30J=2,3IF(IJK+U(I,J).EQ.0)THENELSEII=2*IR+J-3RESULT_N(II)=0.0DO10JJ=2,N=BAND10AK(II,JJ)=0.0DO20JJ=2,II20AK(II-JJ+1,JJ)=0.0ENDIF30CONTINUERET

42、URNENDCCtosolvethebandedFEMequationbyGAUSSeliminationSUBROUTINESOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)DIMENTIONRESULT_N(N_DOF),AK(500,100)DO20K=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_DOFDO20I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO10J=1,IWM=J+I-K10AK(I,J)=AK(I,

43、J)-C*AK(K,M)20RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO40I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDK=J+I-130RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50FORMAT(/12X,

44、9;*RESULTBYFEM2D*',/8X,&'-DISPLACEMENTOFNODE-'/5X,'NODENO',8X,'X-DISP',8X,'Y-DISP)DO60I=1,N_NODE60WRITE(8,70)I,RESULT_N(2*I-1),RESULT_N(2*I)70FORMAT(8X,I5,7X,2E15.6)RESULTENDCCcalculatethestresscomponentofelementandnodeSUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJ

45、K_ELE,X,Y,RESULT_N,&STS_ELE,STS_ND)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6),&SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSIONSTS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10FORMAT(/8X,'-STRESSOFELEMENT-')CALLCAL_DD(PE,PR,DD)DO50IE=1,N_ELECALLCAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,A

46、E,BB)DO20I=1,3DO20J=1,6SS(I,J)=0.0DO20K=1,320SS(I,J)=SS(I,J)+DD(I,K)*BB(K,L)DO30I=1,3DO30J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO40J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE

47、,3)=TXY50WRITE(8,60)IE,STX,STY,TXY60FORMAT(1X,'ELEMENTNO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)CthefollowingpartistocalculatestresscomponentofnodeWRITE(8,55)55FORMAT(/8X,'-STRESSESOFNODE-')DO90I=1,N_NODEA=0.B=0.C=0.II=0DO70K=1,N_ELEDO70J=

48、1,3IF(IJK_ELE(K,J).EQ.I)THENA=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)ENDIF70CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75FORMAT(1X,'NODENO.='I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90

49、CONTONUERETURNENDCFEM2Dprogramend3.1.5算例结果chengxu,for所输出的数据文件DATA.OUT数据内容如下：=PLANESTRESSPROBLEM=RESULTSBYFEM2D-DISPLACEMENTOFNODE-NODENOX-DISPY-DISP.000000E+00-.52527E+01.000000E+00-.22527E+01-.10879E+01-.137363E+01000000E+00.000000E+00-.824176E+00.000000E+00-.182418E+01.000000E+00-STRESSOFELEMENT-E

50、LEMENTNO.=STX=-.108791E+01TXY=.439560E+00ELEMENTNO.=2STX=-.824176E+00TXY=.000000E+00ELEMENTNO.=3STX=-.108791E+01TXY=.307692E+00ELEMENTNO.4STX=-100000E+01TXY=-131868E+00-STRESSOFNODE-NODENO.1STX=-.108791E+01STY=-300000E+01TXY=.439560E+00NODENO.=TXY=.249084E+00NODENO.TXY=.205128E+00NODENO.TXY=.000000E

51、+00NODENO.=TXY=-586081E-01NODENO.=2STX=-.100000E+013STX=-.105861E+014STX=-.824176E+005STX=-970696E+006STX=-100000E+01STY=-220879E+01STY=-191575E+01STY=-225275E+01STY=-166667E+01STY=-137363E+01TXY=-131868E+003.2负重轮轮盘强度和刚度的板壳理论分析履带车辆负重轮是实心轮胎，目前负重轮的设计基本沿用经验方法，尤其是轮盘的设计，无理论可循。本文利用板壳理论对负重轮轮盘的强度和刚度进行分析计算，以

52、使负重轮的设计更加合理。负重轮的轮盘形状复杂多样，目的是增加轮盘的强度和刚度，同时也给设计计算带来很多麻烦，有的甚至成为工程上难以解决的问题。为此，对534式装甲运输车负重轮轮盘的计算模型进行以下简化：把中曲面是一个复杂形状的扁壳的轮盘简化为一圆板。由丁扁壳结构的应力分布比薄板更均匀，其数值都小丁薄板的相应应力，因而利用此模型进行强度和刚度分析使轮盘的设计更加可靠：1. 由丁两排单轮用中心部分的6个螺栓联接在一起，因此可以将螺栓以内的部分看作是绝对刚性的；以其中一排负重轮为研究对象，将另一排负重轮施加在中央刚性部分的外载简化为弯矩M;轮盘和轮圈通过橡胶轮缘同履带接触，其轴向位移很小。故可以将轮

53、盘部分看作是轴向固定的。简化后得到轮盘模型如图1所示。图2为负重轮轮圈受力图，图中bt为单排负重轮的宽度；p为橡胶轮缘与轮圈接触处作用力集度分布，峰值为P0O利用抛物线拟合轮圈与轮缘接触处作用力集度分布，可以求得作用力集度为：2-14P02P020xL2作用丁轮圈的总作用力为：P总L2LPo2%2dxL22PoL2bt32-2当负重轮所受的动载荷P总确定后，由式(2)可求出峰值poo作用在轮盘上的弯矩P总bt2在一定的载荷作用下，轮盘中面的挠度大小直接反映了轮盘的刚度。图3为承受集中力矩的圆板的变形图，0为在力矩M的作用下中央刚性部分的倾角根据板壳理论2.3，经过化简几何方程、物理方程、平衡方程，可以获得极坐标下的圆形板弯曲