程序设计课件1

有限元程式數據輸入(前處理)

單元幾何:單元編號;節點座標材料性質邊界條件(約束)

荷載(體積力,面積力)求解計算

單元剛度矩陣(形函數矩陣，B矩陣，D矩陣，高斯積分)

單元荷載列陣組集整體剛度矩陣和整體荷載列陣方程組求解結果整理（後處理）

等值線向量圖變形網格圖編程基本原則

結構化，模組化可讀性強

變數命名簡單直觀，適當注釋程式條理清楚，邏輯性強

節約記憶體，計算效率優化FEAP程式輸入基本變數

D(10,NUMAT)材料性質數組

F(NDF,NUMNP)節點力數組

ID(NDF,NUMNP)約束資訊數組

IE(NUMMAT)各組材料對應的單元類型

IX(NEN1,NUMEL)單元編號及材料組號

T(NUMNP)節點溫度

X(NDM,NUMNP)節點座標

NDF節點自由度數目

NDM問題的維數

NEN單元的節點數

NEN1NEN+1NUMEL單元總數NUMMAT材料性質組數NUMNP節點總數巨集命令COORELEMBOUNFORCMATEENDTANG計算整體剛度矩陣

FORM形成線性方程組的右端項SOLV求解方程

DISP輸出位移

STRE輸出應力TANG計算整體剛度矩陣單元數據定位

XL（NDM,NEN)單元內節點座標數組

UL(NDM,NEN)單元內節點位移數組

LD(NDF*NEN)方程編號單元數組的計算

以平面四邊形等參單元為例形函數形函數導數B矩陣D矩陣SUBROUTINESHAPEF(S,T,XL,XSJ,SHP)C---shapefunctionroutinefor4-nodeisoparametricquadrilateralCSHP(1,I)=X-direivativefoI–nodeshapefunctionCSHP(2,I)=Y-direivativefoI–nodeshapefunctionCSHP(3,I)=shapefunctionforI-nodeCXJ=JacobianarrayCXSJ=JacobiandeterminantCDIMENSIONSHP(3,4),XL(2,4),XS(2,2),SI(4),TI(4)DATASI,TI/-0.5,0.5,0.5,-0.5,-0.5,-0.5,0.5,0.5/CC---COMPUTESHAPEFUNCTIONSANDDERIVATIVESC---INNATURALCOORDINATESC

DO100I=1,4SHP(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---COMPUTEJACOBIANTRANSROMATIONC---FROMX,YTOS,TCDO200I=1,2DO200J=1,2XS(I,J)=0.0D0200K=1,4200XS(I,J)=XS(I,J)+XL(I,K)*SHP(J,K)

CC---COMPUTEJACOBIANDETERMINANTCXSJ=XS(1,1)*XS(2,2)-XS(1,2)*XS(2,1)CC---TRANSFORMNATURALDERIVATIVESTOC---X,YDERIVATIVESCDO300I=1,4TEMP=(XS(2,2)*SHP(1,I)–XS(2,1)*SHP(2,I))/XSJSHP(2,I)=(XS(2,2)*SHP(1,I)–XS(2,1)*SHP(2,I))/XSJ300SHP(1,I)=TEMPRETURNENDSUBROUTINEPGAUSS(L,LINT,R,Z,W)CC----GAUSSPOINTSANDWEIGHTSFORTWODIMENSIONSCDATALR/-1,1,1,-1,0,1,0,-1,0/,LZ/-1,-1,1,1,-1,0,1,0,0/LW/4/LINT=L+LGOTO(1,2,3)LC----11INTERGRATION1R(1)=0.Z(1)=0.W(1)=4.RETURN

C----22INTERGRATION2G=1./SQRT(3.)DO21I=1,4R(I)=G*LR(I)Z(I)=G*LZ(I)21W(I)=1.RETURNC----33INTERGRATIONG=1./SQRT(0.6)H=1./81.DO31I=1,9R(I)=G*LR(I)Z(I)=G*LZ(I)31W(I)=H*LW(I)RETURNEND

單元剛度矩陣D矩陣C----ISOPARAMETRICELEMENTSTIFFNESSCOMPUTATIONC----FORISOTROPCILINEARELASTICITYC----PLANESTRESSANDPLANESTRAINDIFFERONLYINC---VALUESOFTHECONSTANTSD(1),D(2),D(3)SUPPLIEDCC----FORPLANESTRESSC---D(1)=E/(1.-NU*NU)C----D(2)=NU*D(1)C----D(3)=E/(2.*(1.+NU))C----DVISANAREAWEIGHTINGC----LINTISTHENUMBERFOINTEGRATIONPOINTSC----NELISTHENUMBEROFNODESONTHEELEMENTC----SISTHEELMENTC----SG,TGAREINTERGRATIONPOINTSINC----NATRUALCOORDINATESC----WGISTHEINTERGARIONWEIGHTCC----FOREACHINTERGRATIONCOMPUTEC----CONTRIBUTIONTOELEMENTSTIFFNESSCDO100L=1,LINTCALLSHAPEF(SG(L),TG(L),XL,XSJ,SHP)DV=XSJ*WG(L)D11=D(1)*DVD22=D(2)*DVD33=D(3)*DV

CC----FOREACHJNODECOMPUTEDB=D*BCDO100J=1,NELDB11=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---FOREACHINODECOMPUTES=BT*DBCDO100I=1,JS(I+I-1,J+J-1)=S(I+I-1,J+J-1)+SHP(1,I)*DB11+SHP(2,I)*DB31S(I+I-1,J+J-1)=S(I+I-1,J+J)+SHP(1,I)*DB12+SHP(2,I)*DB32S(I+I,J+J-1)=S(I+I,J+J-1)+SHP(1,I)*DB31+SHP(2,I)*DB21100S(I+I,J+J)=S(I+I,J+J)+SHP(1,I)*DB32SHP(2,I)*DB22CC---COMPUTELOWERTRIANGULARPARTBYSYMMETRYCNL=NEL+NELDO200I=2,NELDO200J=1,IS(I,J)=S(J,I)RETURNEND

線性方程組求解特點:

稀疏對稱正定解法:

直接解法迭代解法

直接解法一維變帶寬儲存帶寬優化對應關係:

每行(每列)的最大寬度(最大高度)，對角線元素位置

PROFIL副程