Fortran平面鋼架有限元解析總結(jié)計(jì)劃

Fortran平面鋼架有限元剖析總結(jié)計(jì)劃Fortran平面鋼架有限元剖析總結(jié)計(jì)劃PAGE/PAGE25Fortran平面鋼架有限元剖析總結(jié)計(jì)劃PAGE結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

有限元剖析軟件的開(kāi)發(fā)

1.1程序功能

該程序?yàn)槠矫鎰偧莒o力剖析程序，能針對(duì)平面剛架間問(wèn)題進(jìn)行有限元計(jì)算，

計(jì)算桿端位移及桿端力大小。程序從磁盤(pán)文件中讀取單元編號(hào)、節(jié)點(diǎn)編號(hào)及坐標(biāo)、

資料屬性、荷載、界限條件等信息；將桿端位移，桿端力等計(jì)算結(jié)果以磁盤(pán)文件

的形式輸出，采納等帶寬二維數(shù)組儲(chǔ)蓄整體剛度矩陣并使用高斯消去法進(jìn)行求

解。1結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

1.2程序結(jié)構(gòu)及流程

2結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

開(kāi)始

是

結(jié)束

標(biāo)題及數(shù)組說(shuō)明

（讀入題目序號(hào)NO）

子程序READNO能否為零?

否

讀入數(shù)據(jù)并打印

子程序MKE

子程序MR

子程序MAKE

子程序TARN形成整體剛度矩陣

子程序MULV

形成結(jié)點(diǎn)載荷

子程序CALM

子程序MK解方程并打印桿端位移

子程序PE

計(jì)算并打印各桿軸力子程序MULV6

子程序MF

子程序SOLV

子程序MADE

3結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

1.3程序的輸入與輸出

詳細(xì)介紹輸入輸出數(shù)據(jù)的格式。如：數(shù)據(jù)文件分幾個(gè)部分，各有幾行，分別

包括哪些內(nèi)容及其種類(lèi)、先后次序，等等。

輸入，共有九行。第一行：7,13,5,1,2,2。分別為，7個(gè)結(jié)點(diǎn)，13個(gè)自由度，

5個(gè)單元，1個(gè)種類(lèi)，2個(gè)結(jié)點(diǎn)荷載，2個(gè)非結(jié)點(diǎn)荷載。

第二行：1,2,3,0.0,0.0,0,0,,6.0,0.0。分別為：一號(hào)結(jié)點(diǎn)的位移序號(hào)，x方向?yàn)?/p>

1，y方向?yàn)?，轉(zhuǎn)角為3，坐標(biāo)為(0.0,0.0)，由于二號(hào)結(jié)點(diǎn)固結(jié)在地面，所以二號(hào)結(jié)點(diǎn)的位移序號(hào)，x方向?yàn)?，y方向?yàn)?，轉(zhuǎn)角為0，坐標(biāo)為(6.0,0.0)。

第三行：4,5,6,0.0,6.0,4,5,7,0.0,6.0。分別為：三號(hào)結(jié)點(diǎn)的位移序號(hào)，x方向?yàn)?，y方向?yàn)?，轉(zhuǎn)角為6，坐標(biāo)為(0.0,6.0),四號(hào)結(jié)點(diǎn)位移序號(hào)x方向和y同樣，轉(zhuǎn)角為7，坐標(biāo)為(,0.0,6.0)。

第四行：8,9,10,6.0,6.0,0,0,11,0.0,12.0五.號(hào)結(jié)點(diǎn)位移序號(hào)，x方向?yàn)?，y方

向?yàn)?，轉(zhuǎn)角為10，坐標(biāo)為（6.0,6.0）。由于六號(hào)結(jié)點(diǎn)鉸接在地面，所以六號(hào)結(jié)點(diǎn)的位移序號(hào)，x方向和y方向?yàn)?，轉(zhuǎn)角為11，坐標(biāo)為（0.0,12.0）。

第五行：12,0,13,6.0,12.0.由于七號(hào)結(jié)點(diǎn)與地面用滑動(dòng)支座固定，所以七號(hào)

結(jié)點(diǎn)的位移序號(hào)，x方向?yàn)?2，y方向?yàn)?，轉(zhuǎn)角為13，坐標(biāo)（6.0,12.0）.

第六行：1,2,1,1,3,1,4,5,1,3,6,1,5,7,1,分別為，1號(hào)和2號(hào)結(jié)點(diǎn)構(gòu)成的單元為1

號(hào)種類(lèi)。1號(hào)和3號(hào)結(jié)點(diǎn)構(gòu)成的單元為1號(hào)種類(lèi)，4號(hào)和5號(hào)結(jié)點(diǎn)構(gòu)成的為1號(hào)種類(lèi)，3號(hào)和6號(hào)結(jié)點(diǎn)構(gòu)成的單元為1號(hào)種類(lèi)，5號(hào)和7號(hào)結(jié)點(diǎn)構(gòu)成的單元為1號(hào)種類(lèi)。

第七行：分別為，彈性模量為E＝2×108kN/m2，截面面積A=0.16m2,慣性

矩I=0.002m4。

第八行：1號(hào)結(jié)點(diǎn)轉(zhuǎn)角方向的會(huì)集力偶為-20.0kN，3號(hào)結(jié)點(diǎn)會(huì)集力為10.0KN。

第九行：1號(hào)單元，受會(huì)集力（會(huì)集力型號(hào)為3），大小為15.0kN，到始端

的距離為3.0。5號(hào)單元，受均布力（均布力型號(hào)為1），大小為5.0kN，到端點(diǎn)

的距離為5.0。

第十行：0為計(jì)算停止符。

輸出：第一部分為輸入的數(shù)據(jù)。RESULTSOFCALCULATION以下為輸出

結(jié)果，第二部分的第一段為4個(gè)結(jié)點(diǎn)的x,y方向的位移和轉(zhuǎn)角。第二段為1,2,3

4結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

號(hào)單元的軸力，剪力和彎矩。

5結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

1.4程序求解中遇到的問(wèn)題

對(duì)實(shí)例進(jìn)行計(jì)算時(shí)，坐標(biāo)原點(diǎn)采納不同樣的點(diǎn)，會(huì)以致整個(gè)題目的坐標(biāo)值發(fā)生改變，輸入的內(nèi)容會(huì)有所不同樣，最后的結(jié)果也不同樣

對(duì)結(jié)點(diǎn)荷載和非結(jié)點(diǎn)荷載的正負(fù)判斷不同樣，結(jié)點(diǎn)荷載的方向和整體坐標(biāo)有關(guān)，非結(jié)點(diǎn)荷載方向判斷和局部坐標(biāo)有關(guān)。

在非結(jié)點(diǎn)荷載中，均布荷載和會(huì)集力到始端的距離判斷不同樣。

6結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

有限元剖析算例

2.1算例說(shuō)明

已知圖示剛架，各桿的資料及截面均同樣，彈性模量E＝2×108kN/m2,A=0.16m2,慣性矩

I=0.002m4，q=5kN/m.，一號(hào)單元會(huì)集力為15KN,一號(hào)結(jié)點(diǎn)會(huì)集力偶為20KN*M,三號(hào)結(jié)點(diǎn)集

中力為10KN.試求剛架的內(nèi)力。節(jié)點(diǎn)編號(hào)如圖

7結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

2.2理論剖析

對(duì)所采納的力學(xué)問(wèn)題進(jìn)行理論剖析，要有詳細(xì)的推導(dǎo)過(guò)程和計(jì)算結(jié)果。

內(nèi)力計(jì)算

對(duì)結(jié)構(gòu)進(jìn)行剖析，可以看出1,2,4單元構(gòu)成的是二次超靜定結(jié)構(gòu)，3,5單元是靜定結(jié)構(gòu)。所以先對(duì)3,5單元構(gòu)成的結(jié)構(gòu)進(jìn)行剖析。

如上圖所示，可以依據(jù)x，y方向力均衡，對(duì)結(jié)點(diǎn)七力矩均衡算得支座反力。再畫(huà)出其彎矩，剪力軸力爭(zhēng)。此后對(duì)1,2,4單元構(gòu)成的結(jié)構(gòu)剖析。

8結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

用力法解超靜定，將結(jié)點(diǎn)六的拘束消除，加上支座反力x1=1，x2=1.畫(huà)出M1，M2,MP圖。

MP圖M1圖M2圖

9結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

此后畫(huà)出其彎矩，剪力，軸力爭(zhēng)

彎矩圖剪力爭(zhēng)軸力爭(zhēng)

10結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

位移計(jì)算

計(jì)算結(jié)點(diǎn)1位移，x方向加單位力1

其剪力與彎矩圖都為零，軸力爭(zhēng)為

依據(jù)公式:

軸力爭(zhēng)圖乘1x=（12.92*6*1）/EA=（12.92*6*1）/3*10^7*0.16=1.615e-5

y方向加單位力1，忽視剪力的影響，

彎矩圖圖乘

1y=-1/EI(3*15.05*3/2+1/2*3*2.82*2/3*3+1/2*1.27*(3+1.27/3)*17.87)+

（略小于程序結(jié)果）

11結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

加單位力偶，剪力與軸力爭(zhēng)為零，

彎矩圖為

計(jì)算結(jié)點(diǎn)3的位移，x方向加單位力1.

彎矩軸力

*2.71/3)+（12.92*6*1）/EA=-0.675E-2.其他位移同理可得。

12結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

2.3輸入輸出數(shù)據(jù)

輸入：

輸出：

13結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

14結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

2.4剖析結(jié)果

理論剖析中，由于內(nèi)力計(jì)算應(yīng)用了力法，所以程序所得結(jié)果和理論結(jié)果一致。而對(duì)位移進(jìn)行理論剖析時(shí)忽視了剪力的影響，所以理論位移略小于程序所計(jì)算的結(jié)果?？梢钥闯鲕浖牡恼_性很高，但是此軟件只合用計(jì)算平面桿系結(jié)構(gòu)，不可以解決彈性力學(xué)問(wèn)題.

結(jié)點(diǎn)1位移uvθ理論值1.615e-51.611*E-50.6617e-3程序值1.615e-51.640*E-50.6617e-3

15結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

程序源代碼

附上圓滿(mǎn)的程序源代碼。

PROGRAMPFAP

CANALYSISPROGRAMFORPLANEFRAME

REALK(200,200),KE(6,6),AKE(6,6),X(100),Y(100),AL(100),#EAI(3,100),PJ(100),PF(2,100),R(6,6),P(100),FF(6),

#FE(6),D(100),ADE(6),DE(6),RT(6,6),AFE(6),F(3)

INTEGER

JE(2,100),JN(3,100),JPJ(100),JPF(2,100),M(6),JEAI(100),NO

OPEN(6,FILE='ht2.TXT')

OPEN(8,FILE='ht.txt',STATUS='NEW')

READ(6,*)NO

WRITE(8,'(/9X,A5,I3,A1)')'(NO=',NO,')'

CALLREAD(NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF)DO5I=1,N

P(I)=0.0

DO5J=1,N

K(I,J)=0.0

DO10IE=1,NEL

CALLMKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALLMR(R,IE,JE,X,Y)

CALLMAKE(KE,R,AKE)

CALLCALM(M,IE,JN,JE)

CALLMK(K,AKE,M)

CONTINUEDO20IP=1,NPF

16結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

CALLMR(R,JPF(1,IP),JE,X,Y)

CALLTRAN(R,RT)

CALLPE(FE,IP,JPF,PF,AL)

CALLMULV6(RT,FE,AFE)

CALLCALM(M,JPF(1,IP),JN,JE)

CALLMF(P,AFE,M)

20CONTINUE

DO30I=1,NPJ

P(JPJ(I))=P(JPJ(I))+PJ(I)CALLSOLV(K,P,D,N)WRITE(8,'(/2(26(1H*),A))')'RESULTSOFCALCULATION'WRITE(8,'(/28X,A)')'NODELDISPLACEMENT'

WRITE(8,40)

FORMAT(9X,'NO.N',4X,'X-DISPLACEMENT',2X,#'Y-DISPLACEMENT',3X,'ANG.ROT.(RAD)')DO60KK=1,NJ

DO50II=1,3

F(II)=0.0

I1=JN(II,KK)

50

WRITE(8,70)KK,F(1),F(2),F(3)

FORMAT(4X,I8,2X,3G16.5)WRITE(8,'(/30X,A)')'ELEMANTFORCES'WRITE(8,80)

FORMAT(2X,'NO.E',4X,'N(1)',9X,'Q(1)',9X,'M(1)',

#9X,'N(2)',9X,'Q(2)',9X,'M(2)')DO120IE=1,NEL

CALLMADE(IE,JN,JE,D,ADE)

CALLMKE(KE,IE,JE,JEAI,EAI,X,Y,AL)CALLMR(R,IE,JE,X,Y)17結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

CALLMULV6(R,ADE,DE)

CALLMULV6(KE,DE,FF)

DO100IP=1,NPF

IF(JPF(1,IP).EQ.IE)THEN

CALLPE(FE,IP,JPF,PF,AL)

DO90I=1,6

90FF(I)=FF(I)-FE(I)

ENDIF

CONTINUEWRITE(8,110)IE,(FF(I),I=1,6)

110FORMAT(I5,2X,6G13.7)

CONTINUEGOTO1

END

SUBROUTINEREAD(NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,#JPJ,PJ,JPF,PF)

REALX(100),Y(100),EAI(3,100),PJ(100),PF(2,100)

INTEGERJE(2,100),JN(3,100),JPJ(100),JPF(2,100),JEAI(100),

#TITLE(20)

READ(6,'(20A4)')(TITLE(I),I=1,20)

WRITE(8,'(/9X,20A4)')TITLE

READ(6,*)NJ,N,NEL,NM,NPJ,NPF

WRITE(8,'(/3(8X,A4,1H:,I2))')'NJ=',NJ,

#'N=',N,'NE=',NEL,'NM=',NM,'NPJ=',NPJ,'NPF=',NPF

WRITE(8,10)

FORMAT(/8X,'NO.N(1)(2)(3)',10X,'X',9X,'Y')READ(6,*)((JN(J,I),J=1,3),X(I),Y(I),I=1,NJ)

DO20I=1,NJ

WRITE(8,'(8X,1H(,I2,1H),3I6,4X,2F10.3)')I,JN(1,I),JN(2,I),

#JN(3,I),X(I),Y(I)

18結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

READ(6,*)(JE(1,I),JE(2,I),JEAI(I),I=1,NEL)

WRITE(8,25)

25FORMAT(/7X,'NO.E(1)(2)NO.MATNO.E(1)(2)NO.MAT')

N2=(NEL+1)/2

DO30I=1,N2-1

WRITE(8,40)I,(JE(J,I),J=1,2),JEAI(I),

#I+N2,(JE(J,I+N2),J=1,2),JEAI(I+N2)

#NEL,JE(1,NEL),JE(2,NEL),JEAI(NEL)

FORMAT(4X,4I6,2X,4I6)

READ(6,*)((EAI(I,J),I=1,3),J=1,NM)

WRITE(8,50)(J,(EAI(I,J),I=1,3),J=1,NM)

50FORMAT(/3X,'NO.MAT',6X,'ELASTICMODULUS',5X,#'AREA',7X,'MOMENTOFINERTIA'/(I6,5X,3G16.4))

70

WRITE(8,'(/20X,12HNODELLOADS)')

WRITE(8,'(16XA)')'NO.DISP.VALUE'

READ(6,*)(JPJ(I),PJ(I),I=1,NPJ)

DO60I=1,NPJ

WRITE(8,'(14X,I7,F16.3)')JPJ(I),PJ(I)

CONTINUE100WRITE(8,'(/20X,16HNON-NODELLOADS)')WRITE(8,'(7X,A,8X,A,9X,A)')'NO.EREAD(6,*)(JPF(1,I),JPF(2,I),PF(1,I),PF(2,I),I=1,NPF)

DO80I=1,NPF

WRITE(8,90)(JPF(J,I),J=1,2),PF(1,I),PF(2,I)

FORMAT(6X,I3,8X,I4,8X,2F10.3)

CONTINUE

19結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

RETURN

END

SUBROUTINEMKE(KE,IE,JE,JEAI,EAI,X,Y,AL)REALKE(6,6),X(100),Y(100),EAI(3,100),AL(100),LINTEGERJE(2,100),JEAI(100)II=JE(1,IE)

JJ=JE(2,IE)

MT=JEAI(IE)

L=SQRT((X(JJ)-X(II))**2+(Y(JJ)-Y(II))**2)

AL(IE)=L

A1=EAI(1,MT)*EAI(2,MT)/L

A2=EAI(1,MT)*EAI(3,MT)/L**3

A3=EAI(1,MT)*EAI(3,MT)/L**2

A4=EAI(1,MT)*EAI(3,MT)/L

KE(1,1)=A1

KE(1,4)=-A1

KE(2,2)=12*A2

KE(2,3)=6*A3

KE(2,5)=-12*A2

KE(2,6)=6*A3

KE(3,3)=4*A4

KE(3,5)=-6*A3

KE(3,6)=2*A4

KE(4,4)=A1

KE(5,5)=12*A2

KE(5,6)=-6*A3

KE(6,6)=4*A4

DO10I=1,6

DO10K=1,6

10KE(K,I)=KE(I,K)

20結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

RETURN

END

SUBROUTINEMR(R,IE,JE,X,Y)

REALR(6,6),X(100),Y(100),L,CX,CY

INTEGERJE(2,100)

I=JE(1,IE)

J=JE(2,IE)

L=SQRT((X(J)-X(I))**2+(Y(J)-Y(I))**2)

CX=(X(J)-X(I))/L

CY=(Y(J)-Y(I))/L

DO10J=1,6

DO10I=1,6

R(I,J)=0.0DO20I=1,4,3R(I,I)=CXR(I,I+1)=CY

R(I+1,I)=-CY

R(I+1,I+1)=CX

20R(I+2,I+2)=1

RETURN

END

SUBROUTINEMAKE(KE,R,AKE)

REALKE(6,6),R(6,6),RT(6,6),TMP(6,6),AKE(6,6)

CALLTRAN(R,RT)

CALLMULV(RT,KE,TMP)

CALLMULV(TMP,R,AKE)

RETURN

END

SUBROUTINECALM(M,IE,JN,JE)

INTEGERM(6),JN(3,100),JE(2,100),IE

21結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

DO10I=1,3

M(I)=JN(I,JE(1,IE))

M(I+3)=JN(I,JE(2,IE))

RETURNEND

SUBROUTINEMK(K,AKE,M)REALK(200,200),AKE(6,6)INTEGERM(6)

DO10I=1,6

DO10J=1,6

#K(M(I),M(J))=K(M(I),M(J))+AKE(I,J)

CONTINUERETURNEND

SUBROUTINEPE(FE,IP,JPF,PF,AL)REALFE(6),PF(2,100),AL(100),LINTEGERJPF(2,100)

A=PF(1,IP)

C=PF(2,IP)

L=AL(JPF(1,IP))

IND=JPF(2,IP)

DO5I=1,6

FE(I)=0.0

GOTO(10,20,30,40,50,60),IND

FE(2)=(7.*A/20.+3.*C/20.)*LFE(3)=(A/20.+C/30.)*L**2FE(5)=(3.*A/20.+7*C/20.)*LFE(6)=-(A/30.+C/20.)*L**2

RETURN

22結(jié)構(gòu)力學(xué)剛架問(wèn)題的有限元剖析

FE(5)=A*C**3*(2.*L-C)/2./L**3FE(2)=A*C-FE(5)FE(3)=A*C**2*(6.*L*L-8.*C*L+3.*C*C)/12./L/L

FE(6)=-A*C**3*(4.*L-3.*C)/12./L/L

RETURN

FE(2)=A*(L-C)**2*(L+2.*C)/L**3FE(3)=A*C*(L-C)**2/L**2FE(5)=A-FE(2)FE(6)=-A*C**2*(L-C)/L**2

RETURN

FE(2)=-6.*A*C*(L-C)/L**3FE(3)=A*(L-C)*(L-3.*C)/L**2FE(5)=-FE(2)FE(6)=A*C*(3.*C-2.*L)/L**2

RETURN

50FE(1)=A*(1.-C/L)

FE(4)=A*C/L

RETURN

FE(1)=C*L/2.

FE(4)=FE(1)

RETURNEND

SUBROUTINEMULV6(A,B,C)

REALC(6),A(6,6),B(6)

DO10I=1,6

C(I)=0.0

DO10J=1,6

C(I)=C(I)+A(I,J)