平面弹性力学有限元源程序FORTRAN.docx
- 文档编号:29636517
- 上传时间:2023-07-25
- 格式:DOCX
- 页数:25
- 大小:20.96KB
平面弹性力学有限元源程序FORTRAN.docx
《平面弹性力学有限元源程序FORTRAN.docx》由会员分享,可在线阅读,更多相关《平面弹性力学有限元源程序FORTRAN.docx(25页珍藏版)》请在冰豆网上搜索。
平面弹性力学有限元源程序FORTRAN
平面弹性力学有限元源程序(FORTRAN)
说明:
1基本控制参数信息:
NG,NE,MC,NX,NB,EO,VO,DENSITY,T(共计5个整形数,4个实型数)
NG:
结构的结点总数;
NE:
结构的单元总数;
MC:
平面问题的类型,MC=0,为平面应力,MC=1,为平面应变;
NX:
荷载工况数;
NB:
支承位移数;
EO:
材料弹性模量(Pa);
VO:
材料泊松比; DENSITY:
容重(N/m3) T:
材料厚度(m);
2打印输出控制参数:
NWA,NEW,NWK,NWP(4个整形数) 等于1时,输出,否则不输出。
3单元结点信息:
(K,(IJK(I,K),I=1,3),K=1,NE)(每行4个整形数,共计NE行) K:
单元号; IJK(1,K):
K单元I结点编号; IJK(2,K):
K单元J结点编号; IJK(3,K):
K单元K结点编号;
4 结点坐标信息:
((K,XY(1,K),XY(2,K)),K=1,NG)(每行3个整形数,共计NG行) K:
结点号 XY(1,K):
K结点X坐标; XY(2,K):
K结点Y坐标;
5支承信息:
((K,MB(1,K),MB(2,K),ZB(K)),K=1,NB)(每行3个整形数,1个实型数,共计NB行) K:
支承位移序号; MB(1,K):
第K个支承位移所在的结点号; MB(2,K):
第K个支承位移的坐标方向; ZB(K):
第K个支承位移的数值;
6按NX荷载工况数输入荷载信息:
每一荷载工况如下 :
(1)NF,NP,NM(3个整型数) NF:
集中荷载个数; NP:
分布荷载个数; NM:
计自重单元数;
(2)若NF≠0,则输入下面数据 K,MF(1,K),MF(2,K),ZF(K)(每行3个整形数,1个实型数,共计NF行) K:
集中荷载序号; MF(1,K):
第K个集中荷载作用的结点号; MF(2,K):
第K个集中荷载的坐标方向; ZF(K):
第K个集中荷载的数值;
(3)若NP≠0,则输入下面数据 K,MP(1,K),MP(2,K),ZP(K)(每行3个整形数,1个实型数,共计NP行) K:
分布荷载序号; MP(1,K):
第K个分布荷载作用的结点号; MP(2,K):
第K个分布荷载的坐标方向; ZP(K):
第K个分布荷载的数值;
(4)若NM≠0,则输入下面数据 若NM≥NE,则表示计所有单元的自重,不需输入计自重的单元号; 若NM $DEBUG PROGRAMPLANE IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N) ALLOCATABLE: : IJK(: : ),XY(: : ),BCA(: : ),SK(: : ),STR(: : ),MB(: : ),ZB(: ),B(: ) ALLOCATABLE: : DELD(: : : ),TOD(: : ),DELST(: : : ),TOST(: : ),DELSUP(: : ),TOTSUP(: ) DIMENSIONEK(6,6) CHARACTERPN*40,FN*12 WRITE(*,'(A)')'本程序为计算平面问题的有限元程序' WRITE(*,'(A)')'特点: (1)采用三结点三角形单元;' WRITE(*,'(A)')' (2)采用等带宽存贮技术;' WRITE(*,'(A)')' (3)采用高斯消元法解线性方程组。 ' WRITE(*,'(/A)')'输入计算问题名(PN): ' READ(*,'(A)')PN CALLFNAME(PN,'.DAT',FN) WRITE(*,'(2A)')' 输入数据文件名为: ',FN OPEN(5,FILE=FN,STATUS='OLD') CALLFNAME(PN,'.OUT',FN) WRITE(*,'(/2A)')'结果输出数据文件名为: ',FN OPEN(6,FILE=FN,STATUS='UNKNOWN') CALLFNAME(PN,'.OU1',FN) WRITE(*,'(/2A)')'参数输出数据文件名为: ',FN OPEN(7,FILE=FN,STATUS='UNKNOWN') READ(5,*)NG,NE,MC,NX,NB,EO,VO,DENSITY,T WRITE(6,120)NG,NE,MC,NX,NB WRITE(6,130)EO,VO,DENSITY,T READ(5,*)NWA,NWE,NWK,NWP NT=2*NG ALLOCATE(IJK(3,NE),XY(2,NG),BCA(7,NE),STR(3,NE),MB(2,NB),ZB(NB),B(NT)) ALLOCATE(DELD(2,NG,NX),TOD(2,NG),DELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB)) CALLCLEAR(2,NG,TOD) CALLCLEAR(3,NE,TOST) CALLCLEAR1(NB,TOTSUP) IF(NG.EQ.0)THEN STOP000 ENDIF CALLINPUT(NG,NE,NB,IJK,XY,MB,ZB) CALLCALND(NG,NE,IJK,ND) ALLOCATE(SK(NT,ND)) IF(MC.EQ.0)THEN E=EO V=VO ELSE E=EO/(1-VO*VO) V=VO/(1-VO) ENDIF IP=0 NX1=NX A1=E/(1-V*V)/4.0 A2=0.5*(1-V) CALLABC(NG,NE,IJK,XY,BCA,NWA) CALLCLEAR(NT,ND,SK) DO100K=1,NE CALLKE(K,NE,T,V,A1,A2,IJK,BCA,EK,NWE) CALLSUMK(K,NE,NT,ND,IJK,EK,SK) 100 CONTINUE CALLCHECK(NT,ND,SK) 110 CONTINUE IP=IP+1 WRITE(6,'(///1X,A,I4)')"荷载工况=",IP READ(5,*)NF,NP,NM DOI=1,NT B(I)=0.0 ENDDO IF(NF.GT.0)THEN CALLPF(NF,NT,B) ENDIF IF(NP.GT.0)THEN CALLPP(NP,NG,NT,T,XY,B) ENDIF IF(NM.GT.0)THEN CALLPG(NM,NE,NT,DENSITY,T,IJK,BCA,B) ENDIF DOI=1,NB I1=2*(MB(1,I)-1)+2-MB(2,I) DELSUP(I,IP)=-B(I1) ENDDO IF(NF.NE.0.OR.NP.NE.0.OR.NM.NE.0)THEN CALLDBC(NB,ND,NT,NX,NX1,MB,ZB,SK,B) CALLGAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP) CALLSTRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR) CALLTRES(IP,NG,NE,NX,NT,B,STR,DELD,TOD,DELST,TOST) CALLSUPT(IP,NG,NE,NB,NX,T,V,A1,A2,IJK,MB,BCA,DELD,DELSUP, TOTSUP) CALLOUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP, TOTSUP) ELSE WRITE(*,'(2x,A,I4,A)')"荷载工况=",IP,"没有荷载! " WRITE(6,'(2x,A,I4,A)')"荷载工况=",IP,"没有荷载! " ENDIF NX1=NX1-1 IF(NX1.GT.0)GOTO110 120 FORMAT(1X,'结点总数=',I3,1X,'单元总数=',I3,1X,'问题类型=',& &I1,1X,'荷载工况数=',I2,1X,'支承位移数=',I2) 130 FORMAT(/1X,'弹性模量=',E10.3,1X,'泊松比=',F5.3,1X,'密度=',E10.3,& &1X,'厚度=',F6.3) END SUBROUTINEGAUSS(NT,ND,NX1,NX,SK,B,NWK,NWP) IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N) DIMENSIONSK(NT,ND),A(NT,NT),B(NT) CALLCLEAR(NT,NT,A) DOI=1,NT DOJ=1,ND IF((I+J-1).LE.NT)THEN A(I,I+J-1)=SK(I,J) ENDIF ENDDO ENDDO IF(NWK.EQ.1.AND.NX1.EQ.NX)THEN WRITE(7,'(A)')"结构总刚(列出右上三角元素)" DOI=1,NT WRITE(7,100)(A(I,J),J=1,NT) ENDDO ENDIF IF(NWP.EQ.1)THEN WRITE(7,'(1X,A,I3,A)')"第",NX-NX1+1,"荷载工况的荷载列阵" DOI=1,NT WRITE(7,'(1x,E11.4)')B(I) ENDDO ENDIF DO50K=1,NT-1 DO60I=K+1,NT C1=A(K,I)/A(K,K) DO70J=I,NT A(I,J)=A(I,J)-C1*A(K,J) 70 CONTINUE B(I)=B(I)-C1*B(K) 60 CONTINUE 50 CONTINUE B(NT)=B(NT)/A(NT,NT) DO80I=NT-1,1,-1 DO90J=I+1,NT B(I)=B(I)-A(I,J)*B(J) 90 CONTINUE B(I)=B(I)/A(I,I) 80 CONTINUE 100 FORMAT(1x,4000E10.3) END SUBROUTINECHECK(NT,ND,SK) IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N) DIMENSIONSK(NT,ND) M=0 DOI=1,NT IF(SK(I,1).LE.0.0)THEN M=M+1 ENDIF ENDDO IF(M.GT.0)THEN WRITE(*,*)"总刚矩阵的对角元素为负值! ! " STOP2222 ENDIF END SUBROUTINESTRESS(NE,NT,V,A1,A2,IJK,BCA,B,STR) IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N) DIMENSIOND1(6),D2(3),B(NT),S(3,6),IJK(3,NE),BCA(7,NE),STR(3,NE) CALLCLEAR(3,NE,STR) DO10K=1,NE DOI=1,3 II=IJK(I,K) D1(2*(I-1)+1)=B(2*(II-1)+1) D1(2*(I-1)+2)=B(2*(II-1)+2) ENDDO CALLCLEAR(3,6,S) C1=2*A1/BCA(7,K) DO20I=1,3 S(1,2*(I-1)+1)=C1*BCA(I,K) S(1,2*(I-1)+2)=C1*V*BCA(I+3,K) S(2,2*(I-1)+1)=C1*V*BCA(I,K) S(2,2*(I-1)+2)=C1*BCA(I+3,K) S(3,2*(I-1)+1)=C1*A2*BCA(I+3,K) S(3,2*(I-1)+2)=C1*A2*BCA(I,K) 20 CONTINUE CALLMATMUL(3,6,1,S,D1,D2) DOI=1,3 STR(I,K)=D2(I) ENDDO 10 CONTINUE END SUBROUTINEOUTPUT(IP,NG,NE,NB,NX,MB,DELD,TOD,DELST,TOST,DELSUP,TOTSUP) IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N) DIMENSIONMB(2,NB),DELD(2,NG,NX),TOD(2,NG) DIMENSIONDELST(3,NE,NX),TOST(3,NE),DELSUP(NB,NX),TOTSUP(NB) CHARACTERVE*6 WRITE(6,20)IP WRITE(6,30) DOI=1,NG WRITE(6,40)I,DELD(1,I,IP),TOD(1,I),DELD(2,I,IP),TOD(2,I) ENDDO WRITE(6,50) DOJ=1,NE WRITE(6,60)J,(DELST(L,J,IP),TOST(L,J),L=1,3) ENDDO WRITE(6,70) DOJ=1,NB IF(MB(2,J).EQ.1)THEN VE='x方向' ELSE VE='Y方向' ENDIF WRITE(6,80)MB(1,J),VE,DELSUP(J,IP),TOTSUP(J) ENDDO 20 FORMAT(/1X,'荷载工况=',I3,'的计算结果') 30 FORMAT(/,1X,'结点号',5X,'X位移增量',5X,'X累计位移',5x,& &'Y位移增量',5X,'Y累计位移') 40 FORMAT(1X,I3,6X,F10.4,4x,F10.4,4X,F10.4,4x,F10.4) 50 FORMAT(/,1X,'单元号',6x,'X应力增量',6x,'X累计应力',6x,& &'Y应力增量',6x,'Y累计应力',6x,'剪应力增量',6x,& &'累计剪应力') 60 FORMAT(2X,I3,6X,E10.3,5X,E10.3,5X,E10.3,5X,E10.3,6X,E10.3,& &6X,E10.3) 70 FORMAT(/,1X,'支座结点号',6x,'支反力方向',6x,'支反力增量',& &6x,'支反力累计量') 80 FORMAT(5X,I3,12X,A,6x,E10.3,8X,E10.3) END SUBROUTINEMATMUL(M,N,L,A,B,C) IMPLICITREAL*8(A-H,O-Z),INTEGER(I-N) DIMENSIONA(M,N),B(N,L),C(M,L) DO100I=1,M DO100J=1,L C(I,J)=0.0 DO100K=1,N 100 C(I,J)=C(I,J)+A(I,K)*B(K,J) END SUBROUTINECLEAR(M,N,A) IMPLICITREAL*8(A-H,O-Z),INTEGER
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 平面 弹性 力学 有限元 源程序 FORTRAN
![提示](https://static.bdocx.com/images/bang_tan.gif)