有限元程序設(shè)計(jì)-平面四邊形4結(jié)點(diǎn)等參有限單元法程序設(shè)計(jì)

1、有限元程序設(shè)計(jì)平面四邊形4結(jié)點(diǎn)等參有限單元法程序設(shè)計(jì)1、程序功能及特點(diǎn)a.該程序采用四邊形4節(jié)點(diǎn)等參單元，能解決彈性力學(xué)的平面應(yīng)力應(yīng)變問(wèn)題。b.前處理采用網(wǎng)格自動(dòng)劃分技術(shù)，自動(dòng)生成單元及結(jié)點(diǎn)信息。b.能計(jì)算受集中力、自重體力、分布面力和靜水壓力的作用。c.計(jì)算結(jié)點(diǎn)的位移和單元中心點(diǎn)的應(yīng)力分量及其主應(yīng)力。d.后處理采取整體應(yīng)力磨平求得各個(gè)結(jié)點(diǎn)的應(yīng)力分量。e.算例計(jì)算結(jié)果與ansys計(jì)算結(jié)果比較，并給出誤差分析。f.程序采用visual fortran 5.0編制而成。2、程序流程及圖框圖2- 程序流程圖圖2-子程序框圖其中，各子程序的主要功能為：input輸入原始數(shù)據(jù)huafen自動(dòng)網(wǎng)格劃分，形

2、成coor(2,np),x,y的坐標(biāo)值與單元信息cband形成主元素序號(hào)指示矩陣ma（*）sko形成整體剛度矩陣kconcr計(jì)算集中力引起的等效結(jié)點(diǎn)荷載rbodyr計(jì)算自重體力引起的等效結(jié)點(diǎn)荷載rfacer計(jì)算分布面力引起的等效結(jié)點(diǎn)荷載rdecop支配方程lu三角分解fobalu分解直接解法中的回代過(guò)程outdisp輸出結(jié)點(diǎn)位移分量stress計(jì)算單元應(yīng)力分量outstre輸出單元應(yīng)力分量stif計(jì)算單元?jiǎng)偠染仃噁dnx計(jì)算形函數(shù)對(duì)整體坐標(biāo)的導(dǎo)數(shù)，1,2,3,4。fun8計(jì)算形函數(shù)及雅可比矩陣jsfun 應(yīng)力磨平-單元下的k=ncnscn應(yīng)力磨平-單元下的右端項(xiàng)系數(shù)cnsumskn應(yīng)力磨平-單

3、元下的右端項(xiàng)集成到總體的psumstrs應(yīng)力磨平-單元下的集成到總體的kgaustrss高斯消元求磨平后的應(yīng)力3、輸入數(shù)據(jù)及變量說(shuō)明當(dāng)程序開(kāi)始運(yùn)行時(shí)，按屏幕提示，鍵入數(shù)據(jù)文件的名字。在運(yùn)行程序之前，根據(jù)程序中input需要的數(shù)據(jù)輸入建立一個(gè)存放原始數(shù)據(jù)的文件，這個(gè)文件的名字為indat.dat。數(shù)據(jù)文件包括如下內(nèi)容：（1）總控信息。共一條，4個(gè)數(shù)據(jù)。l，b，nnl，nnb，nm，nrl矩形體長(zhǎng)度b-矩形體寬度nnl-l方向上劃分的結(jié)點(diǎn)數(shù)nnb-b方向上劃分的結(jié)點(diǎn)數(shù)nm單元材料類(lèi)型數(shù)nr約束結(jié)點(diǎn)總數(shù)（2） 結(jié)點(diǎn)約束信息。共nr條，每條依次輸入：ip，ix，iyip結(jié)點(diǎn)號(hào)ix、iy分別為ip結(jié)點(diǎn)在

4、x，y方向的約束情況，如果約束填0，如果自由填1。（3）材料信息。共nm條，每條依次輸入：jj,（ae（i，jj），i=1，4）jj材料類(lèi)型號(hào)，（ae（i，jj），i=1，4）該材料的材料參數(shù)，共4個(gè)參數(shù)，排列順序?yàn)椋簭椥阅Ａ?、泊松比、容重、單元厚度。?） 荷載信息a. 荷載控制信息。共一條，3個(gè)數(shù)據(jù)ncp，izncp受集中力作用的結(jié)點(diǎn)數(shù)iz面力批數(shù)b. 若ncp0，輸入ip，px，pyip結(jié)點(diǎn)號(hào)px、py分別為ip結(jié)點(diǎn)x，y方向的集中力分量。c. 若iz0，輸入面力荷載信息，共iz批，按批輸入：js，nse，（wg（i）i=1，4）js面力批號(hào)nse第js批面力受到面力作用的單元個(gè)數(shù)，（w

5、g（i），i1，4）該面力的特征參數(shù)共4個(gè)數(shù)據(jù)，第1個(gè)數(shù)為面力類(lèi)型，填1表示受靜水壓力作用，填2表示受均布法向壓力作用；第2個(gè)數(shù)為水壓密度，如果是均布?jí)毫η闆r，就填均布?jí)毫Φ募?；?個(gè)數(shù)為最高水位的y坐標(biāo)，如果是均布?jí)毫η闆r，可以填任意數(shù)；第4個(gè)數(shù)為面力作用的單元面的面號(hào)，單元面號(hào)的規(guī)定見(jiàn)圖2-3。（iew（m），m1，nse），iew（*）為受面力作用的單元的單元號(hào)，共nse個(gè)。 圖3-1 單元結(jié)點(diǎn)編碼與面號(hào)4、理論基礎(chǔ)和求解過(guò)程4.1、構(gòu)造插值函數(shù)：本有限元計(jì)算采取的是四邊形八結(jié)點(diǎn)等參元進(jìn)行插值計(jì)算的。直接調(diào)用教材115頁(yè)3.3.21的結(jié)果，寫(xiě)出所有插值函數(shù)：; 4.2位移插值函數(shù)及應(yīng)變

6、應(yīng)力求解：在有限元方法中單元的位移模式一般采用多項(xiàng)式作為近似函數(shù)，多項(xiàng)式的選取應(yīng)由低次到高次的完備多項(xiàng)式。位移模式的選取一般為：u。，為位移模式，為廣義坐標(biāo)向量。根據(jù)方程求解得出廣義坐標(biāo)，可將位移函數(shù)表示成結(jié)點(diǎn)位移的函數(shù)，即 ,，寫(xiě)成矩陣的形式為： n稱(chēng)為插值函數(shù)矩陣或形函數(shù)矩陣，為單元結(jié)點(diǎn)位移列陣。確定了單元位移后，可以很方便地利用幾何方程和物理方程求得單元的應(yīng)變和應(yīng)力。單元應(yīng)變?yōu)椋篵稱(chēng)為應(yīng)變矩陣，l是平面問(wèn)題的微分算子，其中：，根據(jù)物理方程可以求得單元應(yīng)力，其中，s稱(chēng)為應(yīng)力矩陣，b是應(yīng)變矩陣。由于是平面應(yīng)力問(wèn)題，e0和v0取為e和v，所以彈性矩陣。這部分內(nèi)容參考了教材第2.2節(jié)。4.3、等

7、參元變換：為了將局部（自然）坐標(biāo)中幾何形狀規(guī)則的單元轉(zhuǎn)換成總體坐標(biāo)中幾何形狀扭曲的單元，以滿(mǎn)足對(duì)一般形狀求解域進(jìn)行離散化的需要，必須建立坐標(biāo)轉(zhuǎn)換：,最方便的方法是將上式中表示成插值函數(shù)的形式，, 在有限元分析中，為建立求解方程，需要進(jìn)行各個(gè)單元面積內(nèi)的積分，其一般形式為：而g中包含場(chǎng)函數(shù)對(duì)于總體坐標(biāo)x，y的導(dǎo)數(shù)。由于場(chǎng)函數(shù)是用自然坐標(biāo)表述的，又因?yàn)樵谧匀蛔鴺?biāo)內(nèi)的積分限是規(guī)格化的，因此希望能在自然坐標(biāo)內(nèi)按規(guī)格化的數(shù)值積分方法進(jìn)行上述積分。為此需要建立兩個(gè)坐標(biāo)系內(nèi)導(dǎo)數(shù)之間的變換關(guān)系。設(shè), ，xi，yi，zi是結(jié)點(diǎn)在總體坐標(biāo)內(nèi)的坐標(biāo)值，ni為形狀函數(shù)，實(shí)際上它也是用局部（自然）坐標(biāo)表示的插值函數(shù)。對(duì)

8、于等參變換，結(jié)點(diǎn)數(shù)、結(jié)點(diǎn)號(hào)和插值函數(shù)都不變。函數(shù)ni對(duì)的偏導(dǎo)數(shù)可以表示成： 集合成矩陣形式就是： 式中的j稱(chēng)為雅克比矩陣，利用, ，j可以表示為自然坐標(biāo)的函數(shù)，即：,ni對(duì)于x,y的偏導(dǎo)數(shù)可以用自然坐標(biāo)顯示地表示為:其中是的逆矩陣,它可以按照下式計(jì)算得到：,是j的伴隨矩陣，它的元素是j的元素的代數(shù)余子式。4.4、總體應(yīng)力磨平根據(jù)有限單元法第5章中的（5.3.15），即 其另一種矩陣形式為：令并利用以下表達(dá)式，和以及就可以像單元元?jiǎng)偠染仃嚨娇傮w剛度矩陣一樣，求出，只不過(guò)，這里有ke為12*12的矩陣，而“總剛”k為3*np階的矩陣，np為結(jié)點(diǎn)數(shù)。5、算例應(yīng)用本程序計(jì)算一個(gè)矩形土體受到均布荷載時(shí)的

9、位移和應(yīng)力，如下圖所示，三面約束的土體尺寸為40m*10m，取一半為20m*10m，e=1.0104 kn/m2，在右端受均布荷載 kn/m2，不考慮自重體力。給定網(wǎng)格大小為，有限元網(wǎng)格自動(dòng)劃分如圖3-2所示，單元總數(shù)2000，節(jié)點(diǎn)總數(shù)2121。由于對(duì)稱(chēng)性，右端約束為滑移支座，只限制x方向位移。雖然土體一般不為彈性，但是方法是一樣的，只是剛度矩陣改動(dòng)即可使用彈塑性模型。圖5-1 算例模型圖5-2 ansys計(jì)算模型圖實(shí)際計(jì)算結(jié)果與ansys分析結(jié)果的比較（詳細(xì)計(jì)算結(jié)果見(jiàn)后面）：（1） 位移比較比較項(xiàng)目x方向最大位移y方向最大位移本程序計(jì)算結(jié)果0.0863mm5.890mmansys分析結(jié)果0.

10、1200mm6.000mm誤差28.0%1.83%（2）應(yīng)力比較比較項(xiàng)目最大主應(yīng)力本程序計(jì)算結(jié)果9.996ansys分析結(jié)果9.940誤差0.56%誤差分析：本程序計(jì)算得到的y方向最大位移和ansys計(jì)算結(jié)果很相近，由于x方向上的位移不占主要部分，因此，其誤差雖然有些大，但對(duì)總體位移影響不大。 由于網(wǎng)格較密，因此，計(jì)算得到的單元應(yīng)力值（未磨平）與ansys結(jié)果相近，誤差小于1%，不必要應(yīng)力磨平也可以達(dá)到較好的精度。 本程序可以進(jìn)行總體應(yīng)力磨平，但是由于網(wǎng)格數(shù)較多，因此，磨平矩陣階數(shù)較大，可能計(jì)算時(shí)間也很長(zhǎng)，甚至無(wú)法計(jì)算。（3）實(shí)際計(jì)算結(jié)果matlab圖與ansys圖比較：(a)(b)圖5-3

11、x方向應(yīng)力圖，（a)ansys （b) matlab 圖5-4 主應(yīng)力計(jì)算結(jié)果的matlab 圖圖5-5 主應(yīng)力計(jì)算結(jié)果的ansys 圖圖5-6 剪應(yīng)力計(jì)算結(jié)果的matlab 圖圖5-7 剪應(yīng)力計(jì)算結(jié)果的ansys 圖圖5-8 ansys總位移圖圖5-9 ansys x方向位移圖圖5-10 ansys y方向位移圖圖5-11 ansys y方向應(yīng)力圖附錄：（1）輸入文件數(shù)據(jù)：平面四邊形四結(jié)點(diǎn)單元輸入數(shù)據(jù)l b nnl nnb nm nr20.0 10.0 101 21 1 141*約束信息數(shù)據(jù)結(jié)點(diǎn)號(hào) x向約束信息 y向約束信息1002003004005006007008009001000110

12、012001300140015001600170018001900200021002200430064008500106001270014800169001900021100232002530027400295003160033700358003790040000421004420046300484005050052600547005680058900610006310065200673006940071500736007570077800799008200084100862008830090400925009460096700988001009001030001051001072001093

13、001114001135001156001177001198001219001240001261001282001303001324001345001366001387001408001429001450001471001492001513001534001555001576001597001618001639001660001681001702001723001744001765001786001807001828001849001870001891001912001933001954001975001996002017002038002059002080002101002102012103

14、01210401210501210601210701210801210901211001211101211201211301211401211501211601211701211801211901212001212101*材料信息數(shù)據(jù)類(lèi)型號(hào) 彈性模量 泊松比 容重 單元厚度1 10000.0 0.3 0.0 1.0 *荷載信息數(shù)據(jù)受集中力作用的結(jié)點(diǎn)數(shù) 面力批數(shù)0 1集中力數(shù)據(jù)受集中力作用的結(jié)點(diǎn)號(hào) x向集中力分量 y向集中力分量*面力數(shù)據(jù)面力批號(hào) 受面力單元個(gè)數(shù) 面力類(lèi)型 面力集度 最大集度 面力作用面號(hào) 1 15 2 -10 -10 4*面力作用的單元的單元號(hào)20001980196019401

15、9201900188018601840182018001780176017401720（3）源程序： program fem dimension sk(200000),coor(2,10000),ae(4,100),mel(5,10000),wg(4),&jr(2,10000),ma(10000),r(10000),iew(10000),stre(3,10000),&toalfun(10000,10000),sums(10000),goodstr(10000) common /cmn1/ np,ne,nm,nr common /cmn2/ n,mx,nh common /cmn3/ rf(8)

16、,ske(8,8),nn(8) open(1,file=indat.dat,status=old) open(2,file=out.dat,status=new) call input(xl,b,ix,iy,ae,ncp,iz,nnl,nnb,jr) !輸入數(shù)據(jù) call huafen(xl,b,nnl,nnb,coor,mel) call cband (ma,jr,mel) !形成主元素序號(hào)指示矩陣ma do i=1,nh sk(i)=0.0 enddo call sk0(sk,mel,coor,jr,ma,ae) !形成整體剛度矩陣k do i=1,n r(i)=0.0 enddo if(

17、ncp.gt.0) call concr(ncp,r,jr) call bodyr(r,mel,coor,jr,ae) read(1,70) tl read(1,70) tl read(1,70) tl70 format(a) if(iz.gt.0)then do jj=1,iz read (1,*)js,nse,(wg(i),i=1,4) read(1,70) tl read(1,70) tl read(1,*)(iew(m),m=1,nse) call facer(iew,nse,r,mel,coor,jr,wg) enddo endif call decop (sk,ma) call f

18、oba (sk,ma,r) call outdisp(np,r,jr) call stress (coor,mel,jr,ae,r,stre) call gaustrss(mel,coor,ae,toalfun,sums,goodstr,stre,jr,r)!應(yīng)力磨平并輸出磨平后的應(yīng)力值 write(2,(a) program saff has been ended write(*,(a) program saff has been ended close(1) close(2) end!input為原始數(shù)據(jù)輸入子程序 subroutine input(xl,b,ix,iy,ae,ncp,iz

19、,nnl,nnb,jr) dimension ae(4,*),jr(2,*) character*400tl common /cmn1/ np,ne,nm,nr common /cmn2/ n,mx,nh read(1,70)tl read(1,70)tl read(1,*)xl,b,nnl,nnb,nm,nr write(2,10)xl,b,nnl,nnb,nm,nr10 format(5x,平面四邊形四結(jié)點(diǎn)單元有限元程序/5x,l=,f5.2,4x,&b=,f5.2,4x,nnl=,i3,4x,nnb=,i3,4x,nm=,i2,4x,nr=,i3) np=nnl*nnb do 15 i=

20、1,np do 15 j=1,2 15 jr(j,i)=1 read(1,70)tl read(1,70)tl read(1,70)tl write(2,110)110 format(/5x,約束信息/5x,結(jié)點(diǎn)號(hào),5x,x向約束,5x,y向約束) do 20 i=1,nr read(1,*)ip,ix,iy write(2,100)ip,ix,iy100 format(8x,i5,5x,i2,9x,i2) jr(1,ip)=ix jr(2,ip)=iy 20 continue n=0 do 30 i=1,np do 30 j=1,2 if (jr(j,i) 30,30,25 25 n=n+1

21、 jr(j,i)=n 30 continue read(1,70) tl read(1,70) tl read(1,70) tl do 55 j=1,nm read (1,*) jj,(ae(i,jj),i=1,4) write(2,910) jj,(ae(i,jj),i=1,4)55 continue910 format (/5x,材料信息/5x,類(lèi)型號(hào),5x,彈性模量,5x,泊松比,5x,容重& ,8x,單元厚度/5x,i5,4(5x,e8.3) read(1,70) tl read(1,70) tl read(1,70) tl read(1,*)ncp,iz write(2,920)nc

22、p,iz920 format (/5x,荷載信息/5x,受集中力作用的結(jié)點(diǎn)數(shù),8x,面力批數(shù)&/(2(12x,i5) read(1,70) tl read(1,70) tl70 format(a) return end !自動(dòng)劃分網(wǎng)格子程序，計(jì)算單元數(shù)，結(jié)點(diǎn)數(shù)，結(jié)點(diǎn)坐標(biāo)，每個(gè)單元包含的結(jié)點(diǎn)!單元的4個(gè)結(jié)點(diǎn)號(hào)與材料類(lèi)型號(hào)，從左下角逆時(shí)針編號(hào)順序； subroutine huafen(xl,b,nnl,nnb,coor,mel) dimension coor(2,*),mel(5,*) common /cmn1/ np,ne,nm,nr np=nnl*nnb ne=(nnl-1)*(nnb-1)

23、write(2,*)總結(jié)點(diǎn)數(shù) np=,np write(2,*)總單元數(shù) ne=,ne dl=xl/(nnl-1) db=b/(nnb-1) n=1 m=1 t=nnl*nnb do 10 k=nnb,t,nnb do 20 i=n,k ip=i coor(1,i)=dl*(k-nnb)/nnb20 continue n=k+110 continue do 30 j=nnb,t,nnb do 40 i=m,j coor(2,i)=db*(i-m)40 continue m=j+130 continue write(2,100) (j,(coor(i,j),i=1,2),j=1,np)100 f

24、ormat(/7x,結(jié)點(diǎn)號(hào),10x,x坐標(biāo),8x,y坐標(biāo)/(8x,i5,5x,2f12.6) n=1 do 50 k=1,(nnl-1) do 60 i=n,(nnb-1)*k nne=i mel(1,i)=i+k-1 mel(2,i)=mel(1,i)+nnb mel(3,i)=mel(1,i)+nnb+1 mel(4,i)=mel(1,i)+160 continue n=n+(nnb-1)50 continue write(2,111)(j,(mel(i,j),i=1,4),j=1,ne)111 format(/7x,單元號(hào),5x,結(jié)點(diǎn)1,5x,結(jié)點(diǎn)2,5x,結(jié)點(diǎn)3,5x,結(jié)點(diǎn)&4/(7x

25、,i5,5x,i5,5x,i5,5x,i5,5x,i5) do 70 ie=1,ne mel(5,ie)=170 continue end!*! lu三角分解子程序! * subroutine decop (sk,ma) !方程lu分解 dimension sk(*),ma(*) common /cmn2/ n,mx,nh do 50 i=2,n l=i-ma(i)+ma(i-1)+1 k=i-1 l1=l+1 if (l1.gt.k) goto 30 do 20 j=l1,k ij=ma(i)-i+j m=j-ma(j)+ma(j-1)+1 if (l.gt.m) m=l mp=j-1 if

26、 (m.gt.mp) go to 20 do 10 lp=m,mp ip=ma(i)-i+lp jp=ma(j)-j+lp sk(ij)=sk(ij)-sk(ip)*sk(jp) 10 continue 20 continue 30 if (l.gt.k) goto 50 do 40 lp=l,k ip=ma(i)-i+lp lpp=ma(lp) sk(ip)=sk(ip)/sk(lpp) ii=ma(i) sk(ii)=sk(ii)-sk(ip)*sk(ip)*sk(lpp) 40 continue 50 continue return end!*! lu三角分解后的回代求解! * subr

27、outine foba (sk,ma,r) !回代求解 dimension sk(*),ma(*),r(*) common /cmn2/ n,mx,nh do 10 i=2,n l=i-ma(i)+ma(i-1)+1 k=i-1 if (l.gt.k) goto 10 do 5 lp=l,k ip=ma(i)-i+lp r(i)=r(i)-sk(ip)*r(lp) 5 continue 10 continue do 20 i=1,n ii=ma(i)45 r(i)=r(i)/sk(ii)20 continue do 30 j1=2,n i=2+n-j1 l=i-ma(i)+ma(i-1)+1

28、k=i-1 if (l.gt.k) go to 30 do 25 j=l,k ij=ma(i)-i+j55 r(j)=r(j)-sk(ij)*r(i)25 continue30 continue return end!*！ 計(jì)算形函數(shù)對(duì)整體坐標(biāo)的導(dǎo)數(shù)! * subroutine fdnx (xyz,dnx,det,r,s,rjac,iven,nee) dimension xyz(2,*),dnx(2,4),rjac(2,2),iven(*) common /cmn5/ fun(4),pn(2,4),xjac(2,2) call fun8 (xyz,r,s,det)!求形函數(shù) if (det.l

29、t.1.0e-5)then !判斷j的大小 write(2,600) nee,r,s,det write(2,*) (iven(m),m=1,4) stop endif rec=1.00/det !求j的逆 rjac(1,1)=rec*xjac(2,2) rjac(2,2)=rec*xjac(1,1) rjac(2,1)=-rec*xjac(2,1) rjac(1,2)=-rec*xjac(1,2) do 30 k=1,4 !求形函數(shù)的導(dǎo)數(shù)，坐標(biāo)轉(zhuǎn)換后 do 20 i=1,2 dnx(i,k)=0. do 25 m=1,2 dnx(i,k)=dnx(i,k)+rjac(i,m)*pn(m,k)

30、 25 continue 20 continue 30 continue600 format(1x,err0r* negtive or zero jacobian determinant for &element/ele.=,i5,r=,f10.5,6x,s=,f10.5,det=,f12.5) return end!*！ !形函數(shù)子程序，計(jì)算雅可比行列式! * subroutine fun8 (xyz,r,s,det) dimension xyz(2,*),xi(4),eta(4) common /cmn5/ fun(4),pn(2,4),xjac(2,2) data xi/-1.0,1.0

31、,1.0,-1.0/!雙線(xiàn)性自然坐標(biāo)值 data eta/-1.0,-1.0,1.0,1.0/ do 10 i=1,4 g1=(1.0+xi(i)*r) g2=(1.0+eta(i)*s) fun(i)=0.25*g1*g2 pn(1,i)=0.25*xi(i)*g2 !自然坐標(biāo)下的形函數(shù)導(dǎo)數(shù) pn(2,i)=0.25*eta(i)*g1 10 continue do 80 i=1,2 do 75 j=1,2 det=0.00 do 70 k=1,4 det=det+pn(i,k)*xyz(j,k) !求j矩陣的各項(xiàng)值 70 continue xjac(i,j)=det 75 continue

32、 80 continue det=xjac(1,1)*xjac(2,2)-xjac(2,1)*xjac(1,2)!求j矩陣的行列式值 return end!*！ 計(jì)算單元?jiǎng)偠染仃囎映绦? *subroutine stif(xyz,ae,iven) dimension ae(4,*),dnx(2,4),xyz(2,*),iven(*),rjac(2,2) common /cmn1/ np,ne,nm,nr common /cmn2/ n,mx,nh common /cmn3/ rf(8),ske(8,8),nn(8) common /cmn4/ nee,nme common /cmn5/ fun

33、(4),pn(2,4),xjac(2,2) common /gauss/ rstg(3),h(3) do 40 i=1,8 rf(i)=0.00 do 30 j=1,8 ske(i,j)=0.00!賦初值 30 continue 40 continue e=ae(1,nme) !單元參數(shù)，e，v,重度 u=ae(2,nme) gamma=ae(3,nme) d1=e*(1.00-u)/(1.00+u)*(1.00-2.00*u) !剛度系數(shù)矩陣的一些系數(shù) d2=e*u/(1.00+u)*(1.00-2.00*u) d3=e*0.50/(1.00+u) do 120 i=1,4 ii=2*(i-

34、1) i1=ii+1 i2=ii+2 do 115 j=1,4 jj=2*(j-1) j1=jj+1 j2=jj+2 dxx=0 dxy=0 dyx=0 dyy=0 do 99 is=1,3 s=rstg(is) sh=h(is) do 98 ir=1,3 r=rstg(ir) rh=h(ir) call fdnx (xyz,dnx,det,r,s,rjac,iven,nee) dnix=dnx(1,i) dniy=dnx(2,i) dnjx=dnx(1,j) dnjy=dnx(2,j) dxx=dxx+dnix*dnjx*det*rh*sh !按三點(diǎn)高斯積分求解單元?jiǎng)偠戎?dxy=dxy+d

35、nix*dnjy*det*rh*sh dyx=dyx+dniy*dnjx*det*rh*sh dyy=dyy+dniy*dnjy*det*rh*sh 98 continue 99 continue ske(i1,j1)=dxx*d1+dyy*d3 !單元?jiǎng)偠?ske(i2,j2)=dyy*d1+dxx*d3 ske(i1,j2)=dxy*d2+dyx*d3 ske(i2,j1)=dyx*d2+dxy*d3115 continue120 continue return end!*！ 形成主元素序號(hào)指示矩陣ma（*）! * subroutine cband (ma,jr,mel) dimensio

36、n ma(*),jr(2,*),mel(5,*),nn(8) common /cmn1/ np,ne,nm,nr common /cmn2/ n,mx,nh do 65 i=1,n65 ma(i)=0 do 90 ie=1,ne !對(duì)桿件循環(huán) do 75 k=1,4 iek=mel(k,ie) !單元結(jié)點(diǎn)號(hào) do 95 m=1,2 jj=2*(k-1)+m nn(jj)=jr(m,iek) !nn八個(gè)元素，記錄單元每個(gè)結(jié)點(diǎn)自由度信息95 continue75 continue l=n do 80 i=1,8 nni=nn(i) if(nni.eq.0) goto 80 if(nni.lt.l)

37、 l=nni !找單元的最小自由度號(hào) 80 continue do 85 m=1,8 jp=nn(m) if(jp.eq.0) goto 85 !排除約束的 jpl=jp-l+1 !把單元自由度號(hào)從1開(kāi)始 if(jpl.gt.ma(jp) ma(jp)=jpl !公共結(jié)點(diǎn)找最大的自由度號(hào) 85 continue 90 continue mx=0 ma(1)=1 do 10 i=2,n if(ma(i).gt.mx) mx=ma(i) !找總體最大自由度號(hào)，也是半帶寬 ma(i)=ma(i)+ma(i-1) 10 continue nh=ma(n) write (*,500) n,mx,nh w

38、rite (2,500) n,mx,nh500 format (/5x,總自由度數(shù) n=,i5,3x,半帶寬 mx=,i5,3x,一維存儲(chǔ)位數(shù) &nh=,i7) return end!*！ !形成總體剛度矩陣子程序! * subroutine sk0(sk,mel,coor,jr,ma,ae) dimension sk(*),mel(5,*),coor(2,*),jr(2,*),ma(*),ae(4,*),&xyz(2,4),iven(4) common /cmn1/ np,ne,nm,nr common /cmn2/ n,mx,nh common /cmn3/ rf(8),ske(8,8),

39、nn(8) common /cmn4/ nee,nme common /gauss/ rstg(3),h(3) h(1)=0.5555555555555560 h(2)=0.8888888888888890 h(3)=h(1) rstg(1)=-0.7745966692414830 !三點(diǎn)高斯積分 rstg(2)=0.00 rstg(3)=-rstg(1) do 10 ie=1,ne nee=ie !單元號(hào) nme=mel(5,ie) !材料型號(hào) do 75 k=1,4 iek=mel(k,ie) !單元結(jié)點(diǎn)號(hào)碼 iven(k)=iek do 95 m=1,2 jj=2*(k-1)+m nn(

40、jj)=jr(m,iek) !nn八個(gè)元素，記錄單元每個(gè)結(jié)點(diǎn)自由度號(hào)95 xyz(m,k)=coor(m,iek) !結(jié)點(diǎn)坐標(biāo)75 continue call stif(xyz,ae,iven) !單剛 do 60 i=1,8 do 60 j=1,8 ii=nn(i) jj=nn(j) if (jj.eq.0).or.(ii.lt.jj) goto 60 jn=ma(ii)-(ii-jj) !自由度數(shù)相減，算出在一維數(shù)組中位置 sk(jn)=sk(jn)+ske(i,j) !形成總剛 60 continue 70 continue 10 continue return end!*！ 輸出位移計(jì)算結(jié)果子程序! * subroutine outdisp(np,r,jr) dimension r(*),jr(2,*),u(2) real xmax,ymax !x、y向最大位移值 xmax=0.0 ymax=0.0 jnxm