有限元报告——温度场

有限元上机报告——温度场的有限元计算一．问题如图一平面结构在无热源情况下，给定热边界条件，用有限元分析温度分布。二．解决步骤1．对问题的分析采用简单的三角形单元，单元内温度假定为线性分布，即yaxaayxT321),(与平面结构一样，可用单元3个顶点nml、、的温度nmlTTT、、插值单元内部温度场，有eTTNyxT),(其中TnmleTTTT为e单元的节点温度列阵，而形状函数矩阵为nmlTNNNN简单三角形单元内假定的温度场是线性分布的，其形状函数应为2/)(ycxbaNllll对任一个单元e，如面积域为e，则单元泛函数为xy10000100ABDCdxdyyTxTyTxTdxdyyTxTUeee212122而eeTTFTNyxyTxTnmlnmlcccbbbF21所以，泛函数eTeeThTU21单元刚度矩阵nmlnmlnmlFFFcccbbbF2121所以nmlTnTmTlTFFFFFFFF241所以srsrssrrrsrseccbbcbcbhhh4141412.数据准备如图所示，划分单元格每节点有一个自由度，边界约束为1,2,3,4,5,6,7,12,13,18,19,24,25,30,31,33,34,35,36，温度相当于载荷分布，所以只有边界处有载荷。和之前分析步骤相同，可得数据文件INP.DAT。3.程序运行结果三．改变边界条件7123456xy13单元格不变，边界条件改变如下，则程序的运行结果为四．思考与讨论1.分析采用该种单元分析平面温度场时是否可以收敛于真实解。对同样的三角形单元，用节点温度插值单元内部温度的形状函数，与用节点位移插值单元内部位移的形状函数是完全一样的。与位移单元的分析一样，这种单元在单元交界处温度也是连续的，满足本问题的相容性要求。这种单元内部温度T(x,y)为完全的一次多项式，可以实现任意的常温度导数的温度状态，满足插值函数的完备性要求。因而，采用这种单元分析平面温度场时，有限元分析是可以收敛于真实解的。附录：1.input.TXT1,3,36,8,20,50,22.1,0.3,10.,0.,0.,20.,0.,40.,0.,60.,0.,80.,0.,100.,0.,0.,20.,20.,20.,40.,20.,60.,20.,80.,20.,100.,20.,0.,40.,20.,40.,40.,40.,60.,40.,80.,40.,100.,40.,0.,60.,20.,60.,40.,60.,60.,60.,80.,60.,100.,60.,0.,80.,20.,80.,40.,80.,60.,80.,80.,80.,100.,80.,0.,100.,20.,100.,40.,100.,60.,100.,80.,100.,100.,100.,0.,20.,40.,60.,80.,100.,20.,0.,0.,0.,0.,80.,40.,0.,0.,0.,0.,60.,60.,0.,0.,0.,0.,40.,80.,0.,0.,0.,0.,20.,100.,80.,60.,40.,20.,0.,1,2,8,2,3,9,3,4,10,4,5,11,5,6,12,1,8,7,2,9,8,3,10,9,4,11,10,5,12,11,7,8,14,8,9,15,9,10,16,10,11,17,11,12,18,7,14,13,8,15,14,9,16,15,10,17,16,11,18,17,13,14,20,14,15,21,15,16,22,16,17,23,17,18,24,13,20,19,14,21,20,15,22,21,16,23,22,17,24,23,19,20,26,20,21,27,21,22,28,22,23,29,23,24,30,19,26,25,20,27,26,21,28,27,22,29,28,23,30,29,25,26,32,26,27,33,27,28,34,28,29,35,29,30,36,25,32,31,26,33,32,27,34,33,28,35,34,29,36,35,1,2,3,4,5,6,7,12,13,18,19,24,25,30,31,32,33,34,35,36,2.PLANE.FORPROGRAMMAINDIMENSIONSK(300,30),EK(12,12),Q(300),MC(55),XY(2,100),XYE(2,4),QE*(12),NX(4,100)OPEN(7,FILE='input.TXT')REWIND7READ(7,*)NF,NE,NN,MB,ND,LE,LSREAD(7,*)E,UM,T10FORMAT(7I5)12FORMAT(3F15.2)WRITE(*,600)NF,NE,NN,MB,ND,LE,LS,E,UM,TME=NE*NFMS=NN*NFCALLINPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)WRITE(*,102)((XY(I,J),I=1,LS),J=1,NN)102FORMAT(10X,'XY'/,(2X,6F12.3))WRITE(*,101)(Q(I),I=1,MS)101FORMAT(10X,'Q'/,(2X,6F12.3))WRITE(*,500)((NX(I,J),I=1,NE),J=1,LE)WRITE(*,400)(MC(I),I=1,ND)500FORMAT(10X,'NX'/,(2X,12I6))600FORMAT(10X,'NFNENNMBNDLELSEUMT'/7(2X,I4),3(2X,F8.4))400FORMAT(10X,'MC'/,(2X,10I6))CALLSTIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,UM,T*)CALLSOLVE(SK,Q,MS,MB)OPEN(9,FILE='OUT.DAT')REWIND9WRITE(9,200)WRITE(9,250)(Q(I),I=1,MS,3)200FORMAT(5X,'DISPLACEMENT')250FORMAT(2X,5E14.5)WRITE(9,222)(Q(I),I=1,13,3)222FORMAT(2X,E14.5)cCALLSTRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)STOP1000ENDSUBROUTINEINPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)DIMENSIONXY(LS,NN),Q(MS),NX(NE,LE),MC(ND)READ(7,*)XYREAD(7,*)QREAD(7,*)NXREAD(7,*)MCCLOSE(7)10FORMAT(6F11.2)20FORMAT(12I5)RETURNENDSUBROUTINESTIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E*,UM,T)DIMENSIONSK(MS,MB),EK(ME,ME),Q(MS),NX(NE,LE),MC(ND),XY(LS,NN),XYE*(LS,NE)DO35I=1,MSDO35J=1,MB35SK(I,J)=0.DO200L=1,LEDO40J=1,NELJ=NX(J,L)DO40I=1,LS40XYE(I,J)=XY(I,LJ)DO50I=1,MEDO50J=1,ME50EK(I,J)=0.0CALLSTIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)IF(L.EQ.1)WRITE(*,70)EK70FORMAT(10X,'EK'/,(6E14.5))DO200I=1,NEDO200II=1,NFM=NF*(I-1)+IIM1=NF*(NX(I,L)-1)+IIDO200J=1,NEDO200JJ=1,NFN=NF*(J-1)+JJN1=NF*(NX(J,L)-1)+JJMN=N1-M1+1IF(MN)200,200,150150SK(M1,MN)=SK(M1,MN)+EK(M,N)200CONTINUEDO220I=1,NDM=MC(I)Q(M)=SK(M,1)*Q(M)*1E8220SK(M,1)=SK(M,1)*1E8RETURNENDSUBROUTINESOLVE(SK,Q,MS,MB)DIMENSIONSK(MS,MB),Q(MS)K1=MS-1DO125K=1,K1IF(K+MB-1-MS)105,106,106105N=K+MB-1GOTO110106N=MS110I1=K+1DO125I=I1,NL=I-K+1C=SK(K,L)/SK(K,1)J1=MB-L+1DO122J=1,J1M=J+I-K122SK(I,J)=SK(I,J)-C*SK(K,M)125Q(I)=Q(I)-C*Q(K)Q(MS)=Q(MS)/SK(MS,1)M=MS-1DO145I1=1,MI=MS-I1IF(MS-I+1-MB)135,136,136135N=MS-I+1GOTO140136N=MB140DO142J=2,NL=J+I-1142Q(I)=Q(I)-SK(I,J)*Q(L)145Q(I)=Q(I)/SK(I,1)WRITE(*,147)147FORMAT(5X,'DISPLACEMENT')WRITE(*,150)(Q(I),I=1,MS)150FORMAT(2X,15E14.5)RETURNENDSUBROUTINESTRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,*LS,E,UM,T)DIMENSIONQ(MS),QE(ME),NX(NE,LE),XY(LS,NN),XYE(LS,*NE)DO400L=1,LEDO160I=1,NEDO160J=1,NFN=NF*(I-1)+JN1=NF*(NX(I,L)-1)+J160QE(N)=Q(N1)WRITE(*,165)LWRITE(*,170)(QE(I),I=1,ME)165FORMAT(4X,'L=',I4)170FORMAT(6E14.5)DO200J=1,NELJ=NX(J,L)DO200I=1,LS200XYE(I,J)=XY(I,LJ)CALLSTE(XYE,QE,NE,LS,ME,E,UM,T)400CONTINUERETURNENDSUBROUTINESTIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)DIMENSIONEK(ME,ME),XYE(LS,NE)B(1)=XYE(2,2)-XYE(2,3)B(2)=XYE(2,3)-XYE(2,1)B(3)=XYE(2,1)-XYE(2,2)C(1)=XYE(1,3)-XYE(1,2)C(2)=XYE(1,1)-XYE(1,3)C(3)=XYE(1,2)-XYE(1,1)AE=(B(2)*C(3)-B(3)*C(2))/2DO30I=1,3DO30J=1,330EK(I,J)=B(I)*B(J)+C(I)*C(J)RETURNENDSUBROUTINESTE(XYE,QE,NE,LS,ME,E,UM,T)DIMENSIONQE(ME),XYE(LS,NE)A=(XYE(1,2)-XYE(1,1))/2B=(XYE(2,4)-XYE(2,1))/2CALLSTR(QE,ME,A,B,E,UM,-A,-B,T)CALLSTR(QE,ME,A,B,E,UM,A,-B,T)CALLSTR(QE,ME,A,B,E,UM,A,B