平面四節(jié)點等參單元分析程序(共27頁)

1、變分原理與有限元大作業(yè)平面四節(jié)點等參單元分析程序姓 名：潘 清學 號：SQ10018014033完成時間：2011-4-26一、 概述通常情況下的有限元分析過程是運用可視化分析軟件(如ANSYS、ABAQUS、SAP等)進行前處理和后處理，而中間的計算部分一般采用自己編制的程序來運算。具有較強數(shù)值計算和處理能力的Fortran語言是傳統(tǒng)有限元計算的首選語言。隨著有限元技術(shù)的逐步成熟,它被應(yīng)用在越來越復(fù)雜的問題處理中，但在實際應(yīng)用中也暴露出一些問題。有時網(wǎng)格離散化的區(qū)域較大,而又限于研究精度的要求，使得劃分的網(wǎng)格數(shù)目極其龐大，結(jié)點數(shù)可多達數(shù)萬個，從而造成計算中要運算的數(shù)據(jù)量巨大，程序運行的時間較

2、長的弊端，這就延長了問題解決的時間，使得求解效率降低。因為運行周期長，不利于程序的調(diào)試，特別是對于要計算多種運行工況時的情況；同時大數(shù)據(jù)量處理對計算機的內(nèi)存和CPU 提出了更高的要求，而在實際應(yīng)用中，單靠計算機硬件水平的提高來解決問題的能力是有限的。因此，必須尋找新的編程語言。隨著有限元前后處理的不斷發(fā)展和完善，以及大型工程分析軟件對有限元接口的要求，有限元分析程序不應(yīng)只滿足解題功能，它還應(yīng)滿足軟件工程所要求的結(jié)構(gòu)化程序設(shè)計條件，能夠?qū)Υ鎯M行動態(tài)分配，以充分利用計算機資源，它還應(yīng)很容易地與其它軟件如CAD 的實體造型，優(yōu)化設(shè)計等接口。現(xiàn)在可編寫工程應(yīng)用軟件的計算機語言較多，其中C語言是一個較

3、為優(yōu)秀的語言，很容易滿足現(xiàn)在有限元分析程序編程的要求。C語言最初是為操作系統(tǒng)、編譯器以及文字處理等編程而發(fā)明的。隨著不斷完善，它已應(yīng)用到其它領(lǐng)域，包括工程應(yīng)用軟件的編程。近年來，C語言已經(jīng)成為計算機領(lǐng)域最普及的一個編程語言，幾乎世界上所有的計算機都裝有C的編譯器，從PC機到巨型機到超巨型的并行機，C與所有的硬件和操作系統(tǒng)聯(lián)系在一起。用C 編寫的程序，可移植性極好，幾乎不用作多少修改，就可在任何一臺裝有ANSI、C編譯器的計算機上運行。C既是高級語言，也是低級語言，也就是說，可用它作數(shù)值計算，也可用它對計算機存儲進行操作。二、 編程思想本程序采用C語言編程，編制平面四邊形四節(jié)點等參元程序，用以求

4、解平面結(jié)構(gòu)問題。程序采用二維等帶寬存儲整體剛度矩陣，乘大數(shù)法引入約束，等帶寬高斯消去法求解位移，然后求中間高斯點的應(yīng)力，最后用繞節(jié)點平均法講單元應(yīng)力等效到節(jié)點上，再將結(jié)果寫到tecplot文件中。在有限元程序中，變量數(shù)據(jù)需賦值的可分為節(jié)點信息，單元信息，載荷信息等。對于一個節(jié)點來說，需以下信息：節(jié)點編號（整型），節(jié)點坐標（實型），節(jié)點已知位移（實型），節(jié)點載荷（實型），邊界條件（實型）等。同樣，對于一個單元來說，需以下信息：單元的節(jié)點聯(lián)接信息（整型），材料信息（彈性模量，泊松比等）（實型）等。在FORTRAN 程序中，以上這些變量混合在一起，很難辨認，使程序的可讀性不好，如需要進行單元網(wǎng)絡(luò)的自

5、適應(yīng)劃分，節(jié)點及單元的修改將非常困難。在進行C語言編譯過程中，采用結(jié)構(gòu)struct 使每個節(jié)點信息存儲在一個結(jié)構(gòu)體數(shù)組中，提高程序的可讀性，使數(shù)據(jù)結(jié)構(gòu)更趨于合理。三、平面四節(jié)點等參單元介紹四節(jié)點等參單元實際單元與基本單元的映射關(guān)系如圖 31所示圖 31坐標的映射關(guān)系為：其位移模式和坐標的映射有相同的插值函數(shù)，形函數(shù)為：單元應(yīng)變矩陣為：上式一般簡寫為：其中的子塊矩陣為由于是、的函數(shù)，在中的、要按照復(fù)合函數(shù)來求導，即從而有因此，單元應(yīng)力矩陣為單元剛度矩陣為其中積分采用三點高斯積分，（高斯積分點的總數(shù)），和或是加權(quán)系數(shù)，和是單元內(nèi)的坐標.。對于三點高斯積分，高斯積分點的位置： ，。結(jié)構(gòu)剛度矩陣結(jié)構(gòu)結(jié)

6、點荷載列陣 注意，對于上兩式中的理解不是簡單的疊加而是按照對應(yīng)的自由度集成?？倓偲胶夥匠虖氖缴鲜角蟪觯核摹⒂邢拊治龅哪K組織一個典型的有限元分析過程主要包括以下幾個步驟：1) 讀輸入數(shù)據(jù)，定義節(jié)點及單元數(shù)組。2) 由邊界條件計算方程個數(shù)，賦值荷載列陣。3) 讀入在帶狀存儲的總剛度矩陣中單元和載荷信息。4) 定義總剛度陣數(shù)組。5) 組裝總剛度陣。6) 解方程組。輸入邊界條件(力、位移)形成各荷載工況的節(jié)點荷載陣總剛分解 回代求出位移及輸出 計算應(yīng)變、應(yīng)力形成單元剛陣單剛向總剛投放坐標變換輸入原始參數(shù)計算總剛規(guī)模形成總剛方程向總節(jié)點荷載陣投放形成單元荷載陣調(diào)整幾何、彈性矩陣調(diào)整單元位移列陣其流程

7、圖可見下圖五、程序變量及函數(shù)說明1、控制信息np：結(jié)構(gòu)節(jié)點總數(shù)ne：結(jié)構(gòu)離散單元總數(shù)nr1，nr2：總的約束的節(jié)點數(shù)，nr1，x方向；nr2，y方向nd：每個單元的節(jié)點數(shù)nf：每個節(jié)點的自由度數(shù)ld：集中力載荷的個數(shù)nm：材料的種類nu1，nu2：非零位移邊界條件的節(jié)點數(shù)，nu1，x方向；nu2，y方向u1，u2：非零位移的大小，u1，x方向；u2，y方向n=nf*np ：結(jié)構(gòu)的節(jié)點位移總數(shù)ndf=nd*nf ：每個單元的節(jié)點自由度數(shù)2、輸入的原始數(shù)據(jù)x(np):節(jié)點的x方向坐標y(np):節(jié)點的y方向坐標me（nd，ne）：單元節(jié)點的總體編號nrr(nr1+nr2) ：約束為零的位移所對應(yīng)的

8、總體位移編號p（n）：載荷向量nu（nu1+nu2）：位移載荷mat（6，nm）：材料參數(shù)3、程序中的其他標識符LD（n）：存放結(jié)構(gòu)剛度陣所以主對角線元素在A（nn）中的序號IS（ndf）：單元節(jié)點位移和節(jié)點力在總體位移陣列和載荷陣列中對應(yīng)的序號EK（ndf，ndf）：總體坐標系下的單元剛度矩陣A（nn）：架構(gòu)剛度陣下三角變帶寬一維壓縮存儲的數(shù)組nn：數(shù)組A的元素個數(shù)RSTG(3):高斯積分點的值H(3):高斯積分點的加權(quán)系數(shù)S（6，ne）：各單元的應(yīng)力分量XJAC（2,2）：雅閣比矩陣RJAC（2,2）：雅閣比矩陣的逆PN（2,4）：4個節(jié)點處形函數(shù)對局部坐標的導數(shù)DNX（2,4）：4個節(jié)點

9、處形函數(shù)對整體坐標的導數(shù)FUN（4）：形函數(shù)的值六、計算流程圖輸出節(jié)點位移計算局部坐標系下的單元剛度矩陣坐標變換矩陣計算單元的IS(ndf)數(shù)組取出單元在總體坐標系下的節(jié)點位移單元剛度矩陣向結(jié)構(gòu)剛度矩陣疊加形成單元的IS(ndf)數(shù)組計算局部坐標下的單元剛度矩陣坐標變換矩陣總體坐標下的單元剛度矩陣形成LD（n）數(shù)組輸入原始數(shù)據(jù)n=nf·npndf=nd·nf開始 輸入控制信息結(jié)束計算局部坐標系下單元節(jié)點力應(yīng)力計算局部坐標系下單元節(jié)點位移求解線性方程組求得結(jié)構(gòu)節(jié)點位移進行約束處理 輸出節(jié)點位移和應(yīng)力七、計算結(jié)果與abaqus分析結(jié)果的比較1、中間帶圓孔平面應(yīng)力板的分析。寬40

10、m，長50m，圓孔位于板中心，半徑為4m，承受水平方向位移載荷，E=200Pa，取m。用abaqus建模離散，并計算。再講abaqus離散的節(jié)點和單元信息拷貝到本程序的輸入文件中，用該程序計算，結(jié)果輸出成tecplot文件，用tecplot可以查看結(jié)果。與abaqus的計算結(jié)果進行比較，位移和應(yīng)力云圖如下（左邊是程序計算結(jié)果，右邊是abaqus結(jié)果，下同）：2、纖維增強復(fù)合材料平面應(yīng)力板的分析。寬40m，長50m，增強纖維位于板中心，纖維半徑為8m，承受水平方向位移載荷，基體材料E=20000Pa，纖維材料E=8000Pa，取m。3、半無限大含裂紋板的應(yīng)力分析寬40m，長60m，裂紋位于板左側(cè)

11、中間位置，裂紋長10m，承受豎直方向位移載荷，E=200Pa，取m。結(jié)果分析比較：通過以上三個算例發(fā)現(xiàn)該程序可以用于計算單材料、雙材料、帶孔、含裂紋等各種平面問題，加載條件可以是加集中力和加位移，因此，該程序的適用范圍還是比較廣的。以上三個算例的自編程序所得結(jié)果與abaqus分析結(jié)果進行比較發(fā)現(xiàn)，兩者的計算結(jié)果很接近，而且自編程序?qū)τ诳走厬?yīng)力集中和裂尖應(yīng)力集中都能很好的表達，說明該程序有很好的精度。第三個算例在裂尖處數(shù)值上有些區(qū)別，但總的分布規(guī)律還是很吻合的，主要是因為本程序是用四節(jié)點等參單元，對于應(yīng)力的奇異性表達效果還不是很好。八、附錄附錄一：程序代碼#include<stdio.h&

12、gt;#include<math.h>#include<stdlib.h>#include<iomanip.h>#include<iostream.h>float *float_two_array_malloc(int m,int n) /實型二維動態(tài)數(shù)組分配 float *a; int i,j; a=(float *)malloc(m*sizeof(float *); for (i=0;i<m;i+) ai=(float *)malloc(n*sizeof(float);for (j=0; j<n; j+)aij=0; return

13、 a;float *float_one_array_malloc(int n) /實型一維動態(tài)數(shù)組分配 float *a; int j; a=(float *)malloc(n*sizeof(float);for (j=0; j<n; j+)aj=0; return a;int *int_two_array_malloc(int m,int n) /整型二維動態(tài)數(shù)組分配 int *a; int i,j; a=(int *)malloc(m*sizeof(int *); for (i=0;i<m;i+) ai=(int *)malloc(n*sizeof(int);for (j=0;

14、 j<n; j+)aij=0; return a;int *int_one_array_malloc(int n) /整型一維動態(tài)數(shù)組分配 int *a; int j; a=(int *)malloc(n*sizeof(int);for (j=0; j<n; j+)aj=0; return a;/全局變量定義int np,ne,nr1,nr2,nd,nf,n,ndf,ld,nm,nu1,nu2;float *x,*y,*p,*pp,u1,u2; /x,y,me,nrr,p為輸入信息float coords24,det,fun4,pn24,xjac22,rjac22,dnx24; /

15、求單剛定義的變量int *LD,*me,*nrr,*nu,*IS,nn; /nn為變帶寬一維組中元素個數(shù) /IS表示單元節(jié)點位移在結(jié)構(gòu)位移陣列中對應(yīng)的序號float *EK,*s,*A,*ns,*mat; /EK為總體坐標系下單元的剛度矩陣 /A存儲結(jié)構(gòu)剛度矩陣的一維變帶寬數(shù)組 /s存放各單元應(yīng)力中心點的應(yīng)力/ns存放各節(jié)點的應(yīng)力void input() /輸入函數(shù)定義int i,j;FILE *fp1;fp1=fopen("input.dat","r"); /要輸入的內(nèi)容放在input.dat中if(fp1=NULL)printf("cann

16、't open the file !n");exit(1);/控制信息讀入 fscanf(fp1,"%d",&np);fscanf(fp1,"%d",&ne); fscanf(fp1,"%d",&nr1); fscanf(fp1,"%d",&nr2); fscanf(fp1,"%d",&nd); fscanf(fp1,"%d",&nf);fscanf(fp1,"%d",&ld); f

17、scanf(fp1,"%d",&nm); fscanf(fp1,"%d",&nu1); fscanf(fp1,"%f",&u1); fscanf(fp1,"%d",&nu2); fscanf(fp1,"%f",&u2); n=nf*np; ndf=nf*nd;/printf("%dt%ft%dt%fn",nu1,u1,nu2,u2);/動態(tài)數(shù)組生成x=float_one_array_malloc(np); y=float_one_arr

18、ay_malloc(np);nrr=int_one_array_malloc(nr1+nr2); /位移約束為零的節(jié)點，nr1-x方向約束，nr2-y方向約束nu=int_one_array_malloc(nu1+nu2); /位移約束非零的節(jié)點pp=float_two_array_malloc(ld,2); /非零載荷 （集中力）,ppld0對應(yīng)的自由度序號，ppld1集中力大小p=float_one_array_malloc(n); /載荷矩陣me=int_two_array_malloc(nd,ne); /單元節(jié)點的總體編號mat=float_two_array_malloc(4,nm)

19、;/單元材料數(shù)組，先定義為各向同性材料/原始數(shù)據(jù)讀入 for(i=0;i<np;i+)fscanf(fp1,"%f",&xi);fscanf(fp1,"%f",&yi); for(i=0;i<ne;i+) for(j=0;j<nd;j+)fscanf(fp1,"%d",&meji); meji=meji-1; /abaqus模型節(jié)點編號從1開始，C語言中從0開始for(i=0;i<nr1;i+) /讀入固定位移約束條件fscanf(fp1,"%d",&nrri

20、);/x方向約束的節(jié)點號nrri=(nrri-1)*2; /轉(zhuǎn)化為自由度序號，abaqus模型節(jié)點編號從1開始，C語言中從0開始for(i=0;i<nr2;i+) fscanf(fp1,"%d",&nrrnr1+i);/y方向約束的節(jié)點號，abaqus模型節(jié)點編號從1開始，C語言中從0開始nrrnr1+i=(nrrnr1+i-1)*2+1; for(i=0;i<ld;i+) /讀入集中力載荷fscanf(fp1,"%f%f",&ppi0,&ppi1); for(i=0;i<n;i+) pi=0.0;for(i=0

21、;i<ld;i+) j=int(ppi0); pj=ppi1; for(i=0;i<nm;i+) /讀入材料信息 for(j=0;j<4;j+)fscanf(fp1,"%f",&matji); mat0i=mat0i-1;mat1i=mat1i-1; /abaqus模型節(jié)點編號從1開始，C語言中從0開始for(i=0;i<nu1;i+) /讀入非零位移載荷fscanf(fp1,"%d",&nui);/x方向約束的節(jié)點號nui=(nui-1)*2; /轉(zhuǎn)化為自由度序號，abaqus模型節(jié)點編號從1開始，C語言中從0開

22、始for(i=0;i<nu2;i+) fscanf(fp1,"%d",&nunu1+i);/y方向約束的節(jié)點號，abaqus模型節(jié)點編號從1開始，C語言中從0開始nunu1+i=(nunu1+i-1)*2+1;fclose(fp1);printf("讀入信息完成.n");void output() /結(jié)果輸出函數(shù) int i,j;FILE *fp2;fp2=fopen("output.dat","w"); /運算結(jié)果放在output.dat中if(fp2=NULL)printf("cann&

23、#39;t open the file !n");exit(1); fprintf(fp2,"位移為：n"); /輸出位移for(i=0;i<n;i+=2) fprintf(fp2,"u%d=%f v%d=%fn",i/2+1,pi,i/2+1,pi+1);fprintf(fp2,"n各單元應(yīng)力：Sx,tSy,tSxy,tS1,tS2,tSmisesn"); /輸出單元應(yīng)力for(i=0;i<ne;i+) fprintf(fp2,"%dt",i+1);for(j=0;j<6;j+)fpr

24、intf(fp2,"%.4ft",sji);fprintf(fp2,"n");fclose(fp2);printf("寫入output.dat完成.n");void tecplot() /結(jié)果生成tecplot文件 int i,j;float bl=1.0; /位移顯示比例FILE *fp3;fp3=fopen("tecplot.plt","w"); /運算結(jié)果放在tecplot.plt中if(fp3=NULL)printf("cann't open the file !n&q

25、uot;);exit(1); fprintf(fp3," TITLE = ' test 'n"); fprintf(fp3," VARIABLES=x,y,Ux,Uy,Sx,Sy,Sxy,Smax,Smin,Smisesn"); fprintf(fp3," ZONE N=%d,E=%d,F=FEPOINT,ET=quadrilateraln",np,ne); for(j=0;j<np;j+)fprintf(fp3,"%.4ft%.4ft%.4ft%.4ft",xj+bl*pj*2,yj+bl*

26、pj*2+1,pj*2,pj*2+1);for(i=0;i<6;i+)fprintf(fp3,"%.4ft",nsij);fprintf(fp3,"n"); for(i=0;i<ne;i+)for(j=0;j<nd;j+)fprintf(fp3,"%dt",meji+1);fprintf(fp3,"n");fclose(fp3);printf("寫入tecplot.plt完成.n");void fld() /生成LD數(shù)組，存放組對角線元素在A中的位置 int k,i,j,j1,

27、ig,nw; LD=int_one_array_malloc(n); LD0=0; for(k=0;k<np;k+) ig=np; for(i=0;i<ne;i+) for(j=0;j<nd;j+) if(meji=k)for(j1=0;j1<nd;j1+)if(mej1i<ig) ig=mej1i; for(i=0;i<nf;i+)j=nf*k+i;nw=nf*(k-ig)+i+1;if(j!=0) LDj=LDj-1+nw; nn=LDn-1+1; /nn為變帶寬一維組中元素個數(shù),LD下標從零開始printf("生成LD數(shù)組完成.n"

28、);void fis(int ie) /形成單元ie的節(jié)點自由度號IS數(shù)組 int i,j,ii; /IS=int_one_array_malloc(ndf); for(i=0;i<nd;i+) for(j=0;j<nf;j+) ii=nf*i+j; ISii=nf*meiie+j; void fdnx(int ie,float gs,float gr)/-計算形函數(shù)對整體坐標的導數(shù)dnx24 int m,u,i,v; float g1,g2,xi4,eta4,rec; xi0=-1;xi1=1;xi2=1;xi3=-1;eta0=-1;eta1=-1;eta2=1;eta3=1;

29、/-讀入單元節(jié)點坐標for (i=0;i<4;i+)coords0i=xmeiie;coords1i=ymeiie;/-計算形函數(shù)和|J|for(m=0;m<4;m+)g1=(1+xim*gr);g2=(1+etam*gs);funm=0.25*g1*g2;pn0m=0.25*xim*g2;pn1m=0.25*etam*g1;for(u=0;u<2;u+)for(v=0;v<2;v+)det=0.0;for(m=0;m<4;m+)det=det+pnum*coordsvm;xjacuv=det;det=xjac00*xjac11-xjac01*xjac10;/-雅

30、格比矩陣求逆-rec=1.0/det;rjac00=rec*xjac11;rjac11=rec*xjac00;rjac01=-rec*xjac01;rjac10=-rec*xjac10;/-形函數(shù)對整體坐標的導數(shù)-for(u=0;u<4;u+)for(v=0;v<2;v+)dnxvu=0.0;for(m=0;m<2;m+)dnxvu=dnxvu+rjacvm*pnmu;void fek(int ie) /形成總體坐標系下四節(jié)點平面等參單元ie的剛度矩陣EK（8，8） int i,j,ii,jj,i1,i2,j1,j2,k,l; float D33,E,MU,d1,d2,d3,

31、dxx,dxy,dyx,dyy,gs,gsh,gr,grh,rstg3,h3,dnix,dniy,dnjx,dnjy;rstg0=-0.77459666924;rstg1=0;rstg2=-rstg0;h0=0.555555555556;h1=0.888888888889;h2=h0;/-形成D矩陣- for(i=0;i<nm;i+) /讀入材料if(ie<=mat1i)&&(ie>=mat0i)E=mat2i;MU=mat3i;D00=E/(1.0-MU*MU);D11=D00;D01=D00*MU;D10=D00*MU;D02=0;D12=0;D20=0;

32、D21=0;D22=E/(2*(1+MU);d1=D00;d2=D01;d3=D22;/-單元剛度矩陣清零- for(i=0;i<ndf;i+)for(j=0;j<ndf;j+)EKji=0.0;/-for(i=0;i<nd;i+) /-形成I1，I2，J1，J2-ii=2*i;i1=ii;i2=ii+1;for(j=0;j<nd;j+)jj=2*j;j1=jj;j2=jj+1;/-dxx,dxy,dyx,dyy- dxx=0;dxy=0;dyx=0;dyy=0; for(k=0;k<3;k+)gs=rstgk; /高斯點的值gsh=hk; /加權(quán)系數(shù)for(l=

33、0;l<3;l+)gr=rstgl; grh=hl; fdnx(ie,gs,gr);dnix=dnx0i;dniy=dnx1i;dnjx=dnx0j;dnjy=dnx1j;dxx=dxx+dnix*dnjx*det*gsh*grh;dxy=dxy+dnix*dnjy*det*gsh*grh;dyx=dyx+dniy*dnjx*det*gsh*grh;dyy=dyy+dniy*dnjy*det*gsh*grh;/-分塊形成單元剛度矩陣-EKi1j1=dxx*d1+dyy*d3;EKi2j2=dyy*d1+dxx*d3;EKi1j2=dxy*d2+dyx*d3;EKi2j1=dyx*d2+d

34、xy*d3;void fasum(float*A) /將單元剛度矩陣疊加到結(jié)構(gòu)剛度矩陣中并存到A（nn）中 int i,j,ni,nj,ij;for(i=0;i<ndf;i+)for(j=0;j<ndf;j+)ni=ISi;nj=ISj; if(nj<=ni) ij=LDni-(ni-nj); Aij=Aij+EKij; void fr(float*A)/約束處理子程序int i,j,t;for(i=0;i<nr1+nr2;i+) /固定位移約束處理j=nrri;t=LDj;At=(float)1.0E30;pj=0.0;for(i=0;i<nu1;i+) /非零

35、位移約束處理x方向j=nui;t=LDj;At=(float)1.0E30;pj=(float)1.0E30*u1;for(i=nu1;i<nu1+nu2;i+) /非零位移約束處理y方向j=nui;t=LDj;At=(float)1.0E30;pj=(float)1.0E30*u2;printf("約束處理完成.n");void cholesky(float*A) /用cholesky方法求解線性方程組 int i, j, k,ij,mi; float* z; /增廣矩陣 z = float_two_array_malloc(n, n+1); /* 將原始數(shù)據(jù)導入z

36、中 */ for(i=0; i<n; i+) for(j=0; j<=i; j+) if(i=0) mi=0; else mi=i-(LDi-LDi-1)+1;/第i行第一個非零元素的列號if(j<mi) zij=0;else ij=LDi-(i-j); zij=Aij; for(i=0;i<n;i+)for(j=i+1;j<n;j+)zij=zji; for(i=0;i<n;i+)zin=pi; float temp; /* 分解過程 */ for(k=0; k<n; k+) temp = 0; for(i=0; i<k; i+) temp =

37、 temp + zik * zik; zkk = (float)sqrt(zkk - temp); for(i=k+1; i<n+1; i+) temp = 0; for(j=0; j<k; j+) temp = temp + zjk * zji; zki = (zki - temp) / zkk; /* 回代過程 */ zn-1n = zn-1n / zn-1n-1; for(k=n-2; k>=0; k-) temp = 0; for(i=k+1; i<n; i+) temp = temp + zki * zin; zkn = (zkn - temp) / zkk;

38、 /*將位移陣列存入p中*/ for(i=0; i<n; i+) pi=zin; /*釋放z數(shù)組*/for (i=0;i<n;i+) delete zi; delete z;printf("求解線性方程組完成.n"); void stress(int ie)/ 求單元中心點的應(yīng)力 int i,j,ii,j1,j2;float d1,d2,d3,D33,B38,r8,sig3,E,MU,ss,rr,bi,ci,h1,h2,sx,sy,sxy;/-形成D矩陣- for(i=0;i<nm;i+) /讀入材料if(ie<=mat1i)&&(i

39、e>=mat0i)E=mat2i;MU=mat3i;D00=E/(1.0-MU*MU);D11=D00;D01=D00*MU;D10=D00*MU;D02=0;D12=0;D20=0;D21=0;D22=E/(2*(1+MU);d1=D00;d2=D01;d3=D22;ss=0.0;rr=0.0;fdnx(ie,ss,rr);/-形成B矩陣-for(i=0;i<4;i+)ii=2*i;j1=ii;j2=ii+1;bi=dnx0i;ci=dnx1i;B0j1=bi;B1j1=0.0;B2j1=ci;B0j2=0.0;B1j2=ci;B2j2=bi;/-找出該單元四個節(jié)點的八個位移，寫到r8中-for(i=0;i<4;i+