1、 數(shù)據(jù)輸入(前處理)求解計(jì)算 結(jié)果整理與輸出(后處理)統(tǒng)計(jì)結(jié)果：第1階段所用時(shí)間最長，第3階段次之，第2階段相對(duì)較少。有限元程序三大模塊數(shù)據(jù)輸入(前處理) 單元幾何: 單元編號(hào); 節(jié)點(diǎn)坐標(biāo) 材料性質(zhì) 邊界條件(約束) 荷載(體積力, 面積力)求解計(jì)算 單元?jiǎng)偠染仃?(形函數(shù)矩陣，B矩陣，D矩陣，高斯積分) 單元荷載列陣 組集整體剛度矩陣和整體荷載列陣 方程組求解結(jié)果整理（后處理） 等值線 矢量圖 變形網(wǎng)格圖編程基本原則 結(jié)構(gòu)化，模塊化 可讀性強(qiáng) 變量命名簡單直觀，適當(dāng)注釋 程序條理清楚，邏輯性強(qiáng) 節(jié)約內(nèi)存，計(jì)算效率優(yōu)化有限元分析軟件發(fā)展5個(gè)階段1.1950-1960，航空工業(yè)的發(fā)展促進(jìn)了有限元

2、程序的發(fā)展。2.人們致力于多功能通用有限元程序開發(fā)。如NASTRAN，ANSYS，SAP等。多采用FORTRAN語言開發(fā)。3.進(jìn)一步優(yōu)化程序代碼，采用面向?qū)ο蟮能浖O(shè)計(jì)方法，加強(qiáng)前后處理技術(shù)。4.適應(yīng)計(jì)算機(jī)系統(tǒng)的并行化改進(jìn)。5.20世紀(jì)90年代中期開始，與CAD/CAM進(jìn)行集成，向智能化方向發(fā)展。 近15-20年是有限元分析軟件商品化的迅速發(fā)展階段。整體剛度矩陣的形成單元數(shù)據(jù)定位 XL（NDM, NEN) 單元內(nèi)節(jié)點(diǎn)坐標(biāo)數(shù)組 UL (NDM, NEN) 單元內(nèi)節(jié)點(diǎn)位移數(shù)組 LD(NDF*NEN) 方程編號(hào)單元數(shù)組的計(jì)算 以平面四邊形等參單元為例形函數(shù)坐標(biāo)變換和位移函數(shù)都采用同樣的形函數(shù)和同樣多

3、的參數(shù)，所以稱為等參單元。iii11N41ii iiyNyxNxii iivNvuNu yi, xi, yi, xi, , , i, i,NNNNyxyxNNJ形函數(shù)導(dǎo)數(shù) i, i,1 yi, xi,NNNNJ于是B矩陣D矩陣33 11 12 1211D D DD D0000 x i,y i,y i,x i,iN NN NB00SUBROUTINE SHAPEF(S, T, XL, XSJ, SHP)C- shape function routine for 4-node isoparametric quadrilateralC SHP(1,I) = X-direivative fo I no

4、de shape function C SHP(2,I) = Y-direivative fo I node shape function C SHP(3,I) = shape function for I-nodeC XJ = Jacobian array C XSJ = Jacobian determinant CDIMENSION SHP(3,4), XL(2,4), XS(2,2), SI(4), TI(4)DATA SI,TI/-0.5, 0.5, 0.5, -0.5, -0.5, -0.5, 0.5, 0.5 /CC-COMPUTE SHAPE FUNCTIONS AND DERI

5、VATIVESC- IN NATURAL COORDINATESC DO 100 I = 1,4 SHP(3,I) = (0.5+SI(I)*S)*(0.5+TI(I)*T) SHP(1,I) = SI(I)*(0.5+TI(I)*T) SHP(2,I) = (0.5+SI(I)*S)*TI(I)CC-COMPUTE JACOBIAN TRANSROMATION C-FROM X,Y TO S,TC DO 200 I = 1,2 DO 200 J = 1,2 XS(I,J) = 0.0 D0 200 K = 1,4200 XS(I,J) = XS(I,J) + XL(I,K)*SHP(J,K)

6、 CC-COMPUTE JACOBIAN DETERMINANTC XSJ = XS(1,1)*XS(2,2) - XS(1,2)*XS(2,1) CC- TRANSFORM NATURAL DERIVATIVES TO C- X,Y DERIVATIVESC DO 300 I = 1,4 TEMP = (XS(2,2)*SHP(1,I) XS(2,1)*SHP(2,I)/XSJ SHP(2,I) = (XS(2,2)*SHP(1,I) XS(2,1)*SHP(2,I)/XSJ 300 SHP(1,I) = TEMP RETURN ENDjTiijDBBk單元?jiǎng)偠染仃噅jDBQ xj,33 y

7、j,33 yj,11 xj,12 yj,12 xj,11jND NDND NDND NDQD矩陣33 11 12 1211D D DD D0000 11 -11 -ijeijdd| J | ktk K| J |WWkijijjieij t 32x i,22y i,31x i,21y i,32y i,12x i,31x i,11x i,ijQNQN QNQNQNQN QNQNk用高斯積分轉(zhuǎn)化為求和單元?jiǎng)偠染仃?SUBROUTINE PGAUSS(L,LINT,R,Z,W) 高斯積分權(quán)重CC- GAUSS POINTS AND WEIGHTS FOR TWO DIMENSIONSC DATA LR

8、/-1,1,1,-1,0,1,0,-1,0/,LZ/-1,-1,1,1,-1,0,1,0,0/ LW/4/ LINT = L + L GO TO (1,2,3) LC- 11 INTERGRATION 1 R(1) = 0. Z(1) = 0. W(1) = 4. RETURN C- 2 2 INTERGRATION2 G = 1./SQRT(3.) DO 21 I = 1,4 R(I) = G*LR(I) Z(I) = G*LZ(I)21 W(I) = 1. RETURNC- 3 3 INTERGRATIONG = 1./SQRT(0.6) H = 1./81. DO 31 I = 1,9

9、R(I) = G*LR(I) Z(I) = G*LZ(I)31 W(I) = H*LW(I) RETURN END C-ISOPARAMETRIC ELEMENT STIFFNESS COMPUTATIONC- FOR ISOTROPCI LINEAR ELASTICITY C-PLANE STRESS AND PLANE STRAIN DIFFER ONLY IN C- VALUES OF THE CONSTANTS D(1),D(2) ,D(3) SUPPLIEDCC-FOR PLANE STRESSC- D(1) = E/(1.-NU*NU)C- D(2) = NU*D(1)C- D(3

10、) = E/(2.*(1.+NU)C- DV IS AN AREA WEIGHTINGC- LINT IS THE NUMBER FO INTEGRATION POINTSC- NEL IS THE NUMBER OF NODES ON THE ELEMENT C- S IS THE ELMENT C- SG,TG ARE INTERGRATION POINTS INC- NATRUAL COORDINATESC- WG IS THE INTERGARION WEIGHTCC- FOR EACH INTERGRATION COMPUTE C- CONTRIBUTION TO ELEMENT S

11、TIFFNESSC DO 100 L = 1,LINT CALL SHAPEF(SG(L),TG(L),XL, XSJ, SHP) DV = XSJ*WG(L) D11 = D(1)*DV D22 = D(2)*DV D33 = D(3)*DV CC- FOR EACH J NODE COMPUTE DB = D*BC DO 100 J = 1, NEL DB11 = D11*SHP(1,J) DB12 = D12*SHP(2,J) DB21 = D12*SHP(1,J) DB22 = D11*SHP(2,J) DB31 = D33*SHP(2,J) DB32 = D33*SHP(1,J)CC

12、-FOR EACH I NODE COMPUTE S = BT*DBC DO 100 I = 1,J S(I+I-1,J+J-1) = S(I+I-1,J+J-1) + SHP(1,I)*DB11+SHP(2,I)*DB31 S(I+I-1,J+J-1) = S(I+I-1,J+J ) + SHP(1,I)*DB12+SHP(2,I)*DB32 S(I+I ,J+J-1) = S(I+I ,J+J-1) + SHP(1,I)*DB31+SHP(2,I)*DB21100 S(I+I , J+J) = S(I+I ,J+J ) + SHP(1,I)*DB32SHP(2,I)*DB22CC-COMP

13、UTE LOWER TRIANGULAR PART BY SYMMETRYC NL = NEL + NEL DO 200 I = 2,NEL DO 200 J = 1, IS(I,J) = S(J,I) RETURN END 線性方程組求解特點(diǎn): 大型 稀疏 對(duì)稱 正定 主元占優(yōu) 解法: 直接解法：高斯消元法，三角分解法，以及由此發(fā)展的分塊解法和波前法等。計(jì)算量確定，適用于中小規(guī)模問題。 迭代解法：雅可比迭代法，高斯-賽德爾迭代法和共軛梯度法等。計(jì)算時(shí)間與方程條件數(shù)有關(guān)，可避免誤差累積，適用于大規(guī)模問題。直接解法一維變帶寬儲(chǔ)存 帶寬優(yōu)化 對(duì)應(yīng)關(guān)系: 每行(每列)的最大寬度(最大高度)， 對(duì)角線元素位置 求解 三角分解Ka=rK=LU 1 . L L . . . . 1 L . 0 1 Lnn 1221.00nn2n22n1211U . 0 . . .