有限元程序設(shè)計大作業(yè)

1、1. 不同板寬的孔邊應(yīng)力集中問題姓名：胡宇學(xué)號：21201201282. 摘要本文采用MATLAB和FOTRAN四節(jié)點平面單元，利用有限元數(shù)值解法對不同板寬的孔邊應(yīng)力集中問題進(jìn)行了數(shù)值模擬研究。對于不同的板寬系數(shù)（半板寬b/孔半徑r），得到了不同的應(yīng)力集中系數(shù)，并且與解析解進(jìn)行了比較，驗證了有限元解的正確性，并且得出了解析解的適用范圍。3. 引言通常情況下的有限元分析過程是運用可視化分析軟件(如ANSYS、ABAQUS、SAP等)進(jìn)行前處理和后處理，而中間的計算部分一般采用自己編制的程序來運算。具有較強數(shù)值計算和處理能力的Fortran語言是傳統(tǒng)有限元計算的首選語言。隨著有限元技術(shù)的逐步成熟,

2、它被應(yīng)用在越來越復(fù)雜的問題處理中。MATLAB是由美國MATHWORKS公司發(fā)布的主要面對科學(xué)計算、可視化以及交互式程序設(shè)計的高科技計算環(huán)境。它將數(shù)值分析、矩陣計算、科學(xué)數(shù)據(jù)可視化以及非線性動態(tài)系統(tǒng)的建模和仿真等諸多強大功能集成在一個易于使用的視窗環(huán)境中，為科學(xué)研究、工程設(shè)計以及必須進(jìn)行有效數(shù)值計算的眾多科學(xué)領(lǐng)域提供了一種全面的解決方案，并在很大程度上擺脫了傳統(tǒng)非交互式程序設(shè)計語言（如C、Fortran）的編輯模式，代表了當(dāng)今國際科學(xué)計算軟件的先進(jìn)水平。4. MATLAB部分1，計算模型本程序采用MATLAB編程，編制平面四邊形四節(jié)點等參元程序，用以求解近似平面結(jié)構(gòu)問題。本程序的研究對象為中央

3、開有小孔的長方形板，選取的材料參數(shù)為：板厚h=1、材料強度E=1.0e11 Pa、泊松比mu=0.3。此外，為方便網(wǎng)格的劃分和計算，本文所取板的長度與寬度相等。其孔半徑為r=1，板寬為2b待定。由于本程序的目的在于驗證有限元解的正確性和確定解析解的適用范圍，因此要求網(wǎng)格足夠細(xì)密，以滿足程序的精度要求。同時為了減小計算量，我采取網(wǎng)格徑向長度遞增的網(wǎng)格劃分方法。此種方法特點是，靠近小孔部分的網(wǎng)格細(xì)密，在遠(yuǎn)離小孔的過程中，網(wǎng)格逐漸變得稀疏。以下為網(wǎng)格節(jié)點坐標(biāo)計算和單元節(jié)點排序的MATLAB源程序：dr = (b-r)/(m+m2/2+m3/6) ;dfi= (pi/2)/n ;gNode = zer

4、os( (n+1)*(m+1), 2 ) ;for i=1:1:m for j=1:1:n+1 gNode( (i-1)*(n+1)+j, : ) = cos(dfi*(j-1)*(r+. dr*(i-1)+(i-1)2/2+(i-1)3/6),sin(dfi*(j-1)*(r+dr*(i-1)+(i-1)2/2+(i-1)3/6) ; endend for i=1:1:(n/2+1) gNode(n+1)*m+i, : )=b,b*tan(dfi*(i-1);endfor i=1:1:n/2 gNode(n+1)*m+(n/2+1)+i, : )=b/tan(dfi*(i+n/2),b;en

5、dgElement =zeros( m*n, 5 ) ;for i=1:m for j=1:n gElement( (i-1)*n+j, 1:4) = (i-1)*(n+1)+j, . i*(n+1)+j,. i*(n+1)+j+1,. (i-1)*(n+1)+j+1 ; EndendgElement( :, 5 ) = 1 ;以上源程序中g(shù)Node為節(jié)點坐標(biāo)矩陣，gElement為單元節(jié)點排列矩陣。圖（1）為5x10的網(wǎng)格劃分圖，圖（2）為結(jié)構(gòu)的微變形圖。圖（1）圖（2）假定該板所受的外載為p=1.0e9 Pa。由于是取1/4板進(jìn)行計算，需要在分界面上加上合適的邊界條件。以小孔圓心為原點，板

6、長方向為x軸，板寬方向為y軸。我所取得邊界條件為：與x軸平行的分界面上的節(jié)點的y向固定，x向可以移動；與y軸平行的分界面上的節(jié)點的x向固定，y向可以移動。劃分好網(wǎng)格，約束和載荷加好后。就可以計算了。取b=10，m=5，n=10時的應(yīng)力圖。圖（3）為x向應(yīng)力分布圖，圖（4）為y向應(yīng)力分布圖，圖（5）為剪切力分布圖。圖（3）圖（4）圖（5）2，孔邊應(yīng)力研究無限大板寬的孔邊應(yīng)力集中問題的彈性力學(xué)的解析解為： 在孔邊的y軸上有分布:為了驗證有限元解的正確性，我取10x20的網(wǎng)格進(jìn)行了計算，得到了y軸上的x向應(yīng)力分布曲線，并且與解析解進(jìn)行了比較。如圖（6）。y軸上的x向應(yīng)力解析解為：圖（6）同時，小孔邊

7、緣的剪切力分布如圖（7）所示。剪切力的解析解為圖（7）圖中，紅色*點為有限元節(jié)點解，黑色x點為解析解點。由圖可知，有限元解和解析解基本相同，有限元解是有效的。同時，為了研究不同板寬對應(yīng)力集中系數(shù)的影響，我選取了多組數(shù)據(jù)進(jìn)行對比。如表（1）所示。表（1）b or（）b-r2116.028646.11143213.626464.03154323.228973.52485422.909483.31696522.712293.20987622.5751103.14708732.8527113.10699832.7762113.017010932.7102123.0604171642.6707163.0

8、071262552.6589202.9911373662.6545242.9851504972.6526282.9824表中，,為徑向單元的數(shù)目。,分別為與之對應(yīng)的應(yīng)力集中系數(shù)（）。為半板寬b與孔半徑r的比值，即。解析解中的取值為3.下圖（7）為解析解和有限元解的比較。圖（7）圖中黑色o點為理論應(yīng)力集中系數(shù)，紅色x點為較低網(wǎng)格密度下的應(yīng)力集中系數(shù)變化曲線，藍(lán)色*點為較高網(wǎng)格密度下的應(yīng)力集中系數(shù)變化曲線。由圖可知，同為有限元解，網(wǎng)格密度越高，精確度就越高。同時，我們還可以知道，解析解并不是適用于所有的情況。當(dāng)板寬與孔徑的比值小于5時，解析解與有限元解的差別變大，即解析解不再適用。5. FORTR

9、AN部分 為了進(jìn)一步驗證結(jié)果的正確性，我同時運用FORTRAN對相同的網(wǎng)格進(jìn)行了計算。計算結(jié)果表明，有限元法所得的結(jié)果是可信的。6. 結(jié)論有限元解和解析解各有優(yōu)劣，有限元解的優(yōu)點是基本上對任意情況都適用，但是一般情況下計算量都很大，結(jié)果獲取速度慢。解析解則求解迅速，但并不適用于所有情況。本文所研究的孔邊應(yīng)力集中問題就很好的體現(xiàn)了兩種方法的特點。解析解并不適用于半板寬和孔半徑之比小于5的情況，但是對于板寬遠(yuǎn)大于孔徑的情況具有很高的精確度。而有限元解的精確度依賴于網(wǎng)格的密度和算法，在算法相同的情況下，為增加精度而加密網(wǎng)格會使計算量增加，在精度要求不高的情況下適用。附錄一 MATLAB源程序func

10、tion huyu2120120128 global r b m n p E mu% 定義板的尺寸 r = 1 ; % 小孔半徑 b = 10 ; %半板寬 format short e% 定義網(wǎng)格密度 m = 1; % 徑向劃分的數(shù)目 n = 2; % 周向劃分的數(shù)目（必須為偶數(shù)）n=2*m為最優(yōu)% 定義材料性質(zhì) E = 1e11 ; % 彈性模量 mu =0.3; % 泊松比 % 定義均布荷載大小 p = 1.0e9; %單位Pa% FemModel ; % 定義有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 把計算結(jié)果顯示出來 return

11、function FemModel% 定義有限元模型% 全局變量及名稱% r-圓孔半徑% b-板寬% m - 半徑方向單元數(shù)目 % n - 圓周方向單元數(shù)目% E - 彈性模量% mu - 泊松比% p - 均布載荷大小% gNode - 節(jié)點坐標(biāo)% gElement - 單元定義% gMaterial - 材料性質(zhì)% gBC1 - 第一類約束條件% gBC3 - 斜約束% gDF - 分布力% 返回值：% 無 global gNode gElement gMaterial gBC1 gDF global r b m n p E mu % 計算結(jié)點坐標(biāo) dr = (b-r)/(m+m2/2+m

12、3/6) ; dfi= (pi/2)/n ; gNode = zeros( (n+1)*(m+1), 2 ) ; for i=1:1:m for j=1:1:n+1 gNode( (i-1)*(n+1)+j, : ) = cos(dfi*(j-1)*(r+. dr*(i-1)+(i-1)2/2+(i-1)3/6),sin(dfi*(j-1)*(r+dr*(i-1)+(i-1)2/2+(i-1)3/6) ; end end for i=1:1:(n/2+1) gNode(n+1)*m+i, : )=b,b*tan(dfi*(i-1); end for i=1:1:n/2 gNode(n+1)*m

13、+(n/2+1)+i, : )=b/tan(dfi*(i+n/2),b; end % 單元定義 gElement =zeros( m*n, 5 ) ; for i=1:m for j=1:n gElement( (i-1)*n+j, 1:4) = i*(n+1)+j,. i*(n+1)+j+1,. (i-1)*(n+1)+j+1,. (i-1)*(n+1)+j ; end end gElement( :, 5 ) = 1 ; % 定義材料性質(zhì) gMaterial = E, mu, 0 ; % 定義約束條件 gBC1 = zeros( (m+1)*2, 3 ) ; for i=1:1:m+1 g

14、BC1(i, 1: 2)=(n+1)*(i-1)+1,(n+1)*(i-1)+1)*2; end for i=1:1:m+1 gBC1(m+1+i,1 :2 )=(n+1)*i,(n+1)*i-1)*2+1; end gBC1( :, 3 )=0; % 定義分布壓力 gDF = zeros(n/2,4); for i=1:1:n/2 gDF(i,1:2)=(n+1)*m+i,(n+1)*m+i+1; end gDF( :, 3 )=p; gDF( :, 4 )=p; returnfunction SolveModel% 求解有限元模型% 說明：% 該函數(shù)求解有限元模型，過程如下% 1. 計算單

15、元剛度矩陣，集成整體剛度矩陣% 2. 計算單元的等效節(jié)點力，集成整體節(jié)點力向量% 3. 處理約束條件，修改整體剛度矩陣和節(jié)點力向量% 4. 求解方程組，得到整體節(jié)點位移向量gMaterial global gNode gElement gBC1 gDF gK gDelta gElementStress % step1. 定義整體剛度矩陣和節(jié)點力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number * 2, node_number * 2 ) ; f = sparse( node_number * 2, 1 ) ; % s

16、tep2. 計算單元剛度矩陣，并集成到整體剛度矩陣中 element_number,dummy = size( gElement ) ; for ie=1:1:element_number fprintf( sprintf( '計算單元剛度矩陣并集成，當(dāng)前單元: %dn', ie ) ) ; k = StiffnessMatrix( ie ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 計算分布壓力的等效節(jié)點力，并集成到整體節(jié)點力向量中 df_number,dummy = size( gDF ) ; for idf = 1

17、:1:df_number enf = EquivalentDistPressure( gDF(idf,1), gDF(idf,2),gDF(idf,3),gDF(idf,4) ) ; i = gDF(idf, 1) ; j = gDF(idf, 2) ; f( (i-1)*2+1:(i-1)*2+2, 1 ) = f( (i-1)*2+1:(i-1)*2+2, 1 ) + enf( 1:2 ) ; f( (j-1)*2+1:(j-1)*2+2, 1 ) = f( (j-1)*2+1:(j-1)*2+2, 1 ) + enf( 3:4 ) ; end % step4. 處理約束條件，修改剛度矩陣

18、和節(jié)點力向量。采用乘大數(shù)法 bc1_number,dummy = size( gBC1 ) ; for ibc=1:1:bc1_number g = gBC1(ibc, 2 ) ; f(g) = gBC1(ibc, 3)* gK(g,g) * 1e15 ; gK(g,g) = gK(g,g) * 1e15 ; end % step5. 求解方程組，得到節(jié)點位移向量 gDelta = gK f ; % step6. 計算每個單元的結(jié)點應(yīng)力 gElementStress = zeros( element_number,5, 3 ) ; delta = zeros( 8, 1 ) ; for ie

19、= 1:element_number xi = -1 1 1 -1 ; eta = -1 -1 1 1 ; for n=1:4 B = MatrixB( ie, xi(n), eta(n) ) ; D = MatrixD( ie, 1 ) ; delta(1:2:7) = gDelta( gElement(ie,1:4)*2-1 ) ; delta(2:2:8) = gDelta( gElement(ie,1:4)*2) ; sigma = D*B*delta; gElementStress( ie, n, :) = sigma ; end end returnfunction k = Sti

20、ffnessMatrix( ie ) % 計算平面應(yīng)力等參數(shù)單元的剛度矩陣% 輸入?yún)?shù)：% ie - 單元號% 返回值：% k - 單元剛度矩陣% 說明：% 用高斯積分法求解平面等參數(shù)單元的剛度矩陣 k = zeros( 8, 8 ) ; D = MatrixD( ie,1) ; % 3 x 3 高斯積分點和權(quán)系數(shù) %x = -0.774596669241483, 0.0,0.774596669241483 ; %w = 0.555555555555556,0.888888888888889,0.555555555555556 ; % 2 x 2 高斯積分點和權(quán)系數(shù) x = -0.577350

21、269189626, 0.577350269189626 ; w = 1, 1 ; for i=1:1:length(x) for j=1:1:length(x) B = MatrixB( ie, x(i), x(j) ) ; J = Jacobi( ie, x(i), x(j) ) ; k = k + w(i)*w(j)*transpose(B)*D*B*det(J) ; end endreturnfunction D = MatrixD( ie, opt )% 計算單元的彈性矩陣D% 輸入?yún)?shù)：% ie - 單元號% opt - 問題選項% 1 => 平面應(yīng)力% 2 => 平面

22、應(yīng)變% 返回值：% D - 彈性矩陣D global gElement gMaterial E = gMaterial( gElement(ie, 5), 1 ) ; % 彈性模量 mu = gMaterial( gElement(ie, 5), 2 ) ; % 泊松比 if opt = 1 % 平面應(yīng)力的彈性常數(shù) A1 = mu ; A2 = (1-mu)/2 ; A3 = E/(1-mu2) ; else % 平面應(yīng)變的彈性常數(shù) A1 = mu/(1-mu) ; A2 = (1-2*mu)/2/(1-mu) ; A3 = E*(1-mu)/(1+mu)/(1-2*mu) ; end D =

23、 A3* 1 A1 0; A1 1 0; 0 0 A2 ;returnfunction B = MatrixB( ie, xi, eta )% 計算單元的應(yīng)變矩陣B% 輸入?yún)?shù)：% ie - 單元號% xi,eta - 局部坐標(biāo) % 返回值：% B - 在局部坐標(biāo)處的應(yīng)變矩陣B N_x,N_y = N_xy( ie, xi, eta ); B = zeros( 3, 8 ) ; for i=1:1:4 B(1:3,(2*i-1):2*i) = N_x(i) 0; 0 N_y(i); N_y(i) N_x(i); endreturnfunction N_x, N_y = N_xy( ie, xi

24、, eta )% 計算形函數(shù)對整體坐標(biāo)的導(dǎo)數(shù)% 輸入?yún)?shù)：% ie - 單元號% xi,eta - 局部坐標(biāo) % 返回值：% N_x - 在局部坐標(biāo)處的形函數(shù)對x坐標(biāo)的導(dǎo)數(shù)% N_y - 在局部坐標(biāo)處的形函數(shù)對y坐標(biāo)的導(dǎo)數(shù) J = Jacobi( ie, xi, eta ) ; N_xi,N_eta = N_xieta( ie, xi, eta ) ; A=JN_xi;N_eta ; N_x = A(1,:) ; N_y = A(2,:) ;returnfunction N_xi, N_eta = N_xieta( , xi, eta )% 計算形函數(shù)對局部坐標(biāo)的導(dǎo)數(shù)% 輸入?yún)?shù)：% ie -

25、 單元號% xi,eta - 局部坐標(biāo) % 返回值：% N_xi - 在局部坐標(biāo)處的形函數(shù)對xi坐標(biāo)的導(dǎo)數(shù)% N_eta - 在局部坐標(biāo)處的形函數(shù)對eta坐標(biāo)的導(dǎo)數(shù) x = -1, 1, 1, -1 ; e = -1, -1, 1, 1 ; N_xi = zeros( 1, 4 ) ; N_eta = zeros( 1, 4 ) ; N_xi(1) = x(1)*(1+e(1)*eta)/4; N_eta(1) = e(1)*(1+x(1)*xi)/4; N_xi(2) = x(2)*(1+e(2)*eta)/4; N_eta(2) = e(2)*(1+x(2)*xi)/4; N_xi(3)

26、= x(3)*(1+e(3)*eta)/4; N_eta(3) = e(3)*(1+x(3)*xi)/4; N_xi(4) = x(4)*(1+e(4)*eta)/4; N_eta(4) = e(4)*(1+x(4)*xi)/4;returnfunction J = Jacobi( ie, xi, eta )% 計算雅克比矩陣% 輸入?yún)?shù)：% ie - 單元號% xi,eta - 局部坐標(biāo) % 返回值：% J - 在局部坐標(biāo)(xi,eta)處的雅克比矩陣 global gNode gElement x = gNode(gElement(ie,1:4),1) ; y = gNode(gEleme

27、nt(ie,1:4),2) ; N_xi,N_eta = N_xieta( ie, xi, eta ) ; x_xi = N_xi * x ; x_eta = N_eta * x ; y_xi = N_xi * y ; y_eta = N_eta * y ; J = x_xi, y_xi; x_eta, y_eta ;returnfunction AssembleStiffnessMatrix( ie, k )% 把單元剛度矩陣集成到整體剛度矩陣% 輸入?yún)?shù):% ie - 單元號% k - 單元剛度矩陣% 返回值:% 無 global gElement gK for i=1:1:4 for j

28、=1:1:4 for p=1:1:2 for q=1:1:2 s = (i-1)*2+p ; t = (j-1)*2+q ; S = (gElement(ie,i)-1)*2+p ; T = (gElement(ie,j)-1)*2+q ; gK(S,T) = gK(S,T) + k(s,t) ; end end end endreturnfunction enf = EquivalentDistPressure( i, j, pi, pj )% 計算線性分布壓力的等效節(jié)點力% 輸入?yún)?shù):% i,j - 結(jié)點號% pi,pj - 跟結(jié)點號對應(yīng)的壓力值% 返回值:% enf - 等效節(jié)點力向量

29、global gNode % 獲取結(jié)點坐標(biāo) xi = gNode( i, 1 ) ; xj = gNode( j, 1 ) ; yi = gNode( i, 2 ) ; yj = gNode( j, 2 ) ; sig(1:4,1)=0; % 2 x 2 高斯積分點和權(quán)系數(shù) x = -0.577350269189626, 0.577350269189626 ; w = 1, 1 ; for i=1:1:length(x) N1=(1-x(i)/2; N2=(1+x(i)/2; dN1=-0.5; dN2=0.5; dx=xi*dN1+xj*dN2; dy=yi*dN1+yj*dN2; p=pi

30、*N1+pj*N2; sig(1)=sig(1)+w(i)*N1*dy*p; sig(2)=sig(2)-w(i)*N1*dx*p; sig(3)=sig(3)+w(i)*N2*dy*p; sig(4)=sig(4)-w(i)*N2*dx*p; end enf = sig ;returnfunction DisplayResults% 顯示計算結(jié)果，并與解析解結(jié)果比較 global gNode gElement gMaterial gDelta gElementStress outfile gBC1 global gDF n m r p b Stress outfile=fopen('

31、huyu2120120128.txt','w'); Stress = zeros(m+1)*(n+1),3); %集成節(jié)點應(yīng)力矩陣 for k=1:1:3 for i=1:1:m Stress (n+1)*(i-1)+1,k)=gElementStress(n*(i-1)+1,4,k); end for j=1:1:n Stress (n+1)*m+j,k)=gElementStress(n*(m-1)+j,1,k); end Stress (n+1)*(m+1),k)=gElementStress(n*m,2,k); for i=1:1:m for j=1:1:n S

32、tress (i-1)*(n+1)+j+1,k)=gElementStress(i-1)*n+j,3,k); end end end xx_max=max(max(Stress(:,1); N_xx_max=find(Stress=xx_max); xx_min=min(min(Stress(:,1); N_xx_min=find(Stress=xx_min); yy_max=max(max(Stress(:,2); N_yy_max=find(Stress=yy_max)-(m+1)*(n+1); yy_min=min(min(Stress(:,2); N_yy_min=find(Stres

33、s=yy_min)-(m+1)*(n+1); xy_max=max(max(Stress(:,3); N_xy_max=find(Stress=xy_max)-2*(m+1)*(n+1); xy_min=min(min(Stress(:,3); N_xy_min=find(Stress=xy_min)-2*(m+1)*(n+1); fprintf(outfile, ' 表0 節(jié)點坐標(biāo)（mm) n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '結(jié)點號 X坐標(biāo) Y坐標(biāo) n' ) ; fprintf(

34、outfile, '-n' ) ; for i=1:size(gNode(:,1) fprintf(outfile, '%4d , %8.4f , %8.4f n',i,gNode(i,1),gNode(i,2) ; end fprintf( outfile,'=nnn' ) ; fprintf(outfile, ' 表1 節(jié)點位移、應(yīng)力 n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '結(jié)點號 x-位移 y-位移 Stress-xx Stress-yy

35、Stress-xy n' ) ; fprintf(outfile, '-n' ) ; for i=1:size(gNode(:,1) fprintf(outfile, '%4d %12.4e %12.4e %12.4e %12.4e %12.4en',. i, full(gDelta(2*i-1),full(gDelta(2*i),Stress(i,1),Stress(i,2),Stress(i,3) ; end fprintf( outfile,'=nnn' ) ; fprintf(outfile, ' 表2 單元節(jié)點序列 n&

36、#39; ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '單元號 材料組 節(jié)點號（逆時針排列） n' ) ; fprintf(outfile, ' 1 2 3 4 n' ) ; fprintf(outfile, '-n' ) ; for i=1:size(gElement(:,1) fprintf(outfile, '%4d , %4d , %4d , %4d , %4d , %4dn',. i, gElement(i,5), gElement(i,1:4) ; en

37、d fprintf( outfile,'=nnn' ) ; fprintf(outfile, ' 表3 材料性質(zhì) n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '材料號 彈性模量 泊松比 比重 n' ) ; fprintf(outfile, '-n' ) ; for i=1:size(gMaterial(:,1) fprintf(outfile, '%4d %8.3e %8.3e %8.3en',. i,gMaterial(i,1),gMater

38、ial(i,2),gMaterial(i,3) ; end fprintf( outfile,'=nnn' ) ; fprintf(outfile, ' 表4 邊界條件 n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, ' 序號 結(jié)點號 約束自由度 給定值 n' ) ; fprintf(outfile, '-n' ) ; for i=1:size(gBC1(:,1) fprintf(outfile, '%4d %4d %4d %4dn',. i,gB

39、C1(i,1),gBC1(i,2),gBC1(i,3) ; end fprintf(outfile, '=nnn' ) ; fprintf(outfile, ' 表5 分布力 n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '邊界 結(jié)點1 結(jié)點2 均布力1 均布力2 n' ) ; fprintf(outfile, '-n' ) ; for i=1:size(gDF(:,1) fprintf(outfile, '%4d %4d %4d %8.4e %8.4e

40、n', . i,gDF(i,1),gDF(i,2),gDF(i,3),gDF(i,4) ; end fprintf( outfile,'=nnn' ) ; fprintf( outfile,' 表6 xx向節(jié)點應(yīng)力極值 n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '極大點 Stress_xx_max 極小點 Stress_xx_min n' ) ; fprintf(outfile, '-n' ) ; fprintf(outfile,'%4d %12.4e %4d %12.4en', . N_xx_max,xx_max,N_xx_min,xx_min) ; fprintf(outfile, '=nnn') ; fprintf( outfile,' 表7 yy向節(jié)點應(yīng)力極值 n' ) ; fprintf(outfile, '=n' ) ; fprintf(outfile, '極大點 Stress_yy_max 極小點 Stress_yy_mi