八節(jié)點等參單元平面有限元(共37頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上八節(jié)點等參單元平面有限元分析程序土木工程學(xué)院2011.2專心-專注-專業(yè)目錄計算力學(xué)課程大作業(yè)1. 概述通常情況下的有限元分析過程是運用可視化分析軟件(如ANSYS、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è)計等接口?，F(xiàn)在可編寫工程應(yīng)用軟件的計算機語言較多，其中C語言是一個較為優(yōu)秀的語言，很容

3、易滿足現(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ù)值計算，也可用它對計算機存儲進行操作。2. 編程思想本程序采用C語言編程，編制平面四邊形八節(jié)點等參元程序，用以求解平面結(jié)構(gòu)問題。程

4、序采用二維等帶寬存儲整體剛度矩陣，乘大數(shù)法引入約束，等帶寬高斯消去法求解位移。在有限元程序中，變量數(shù)據(jù)需賦值的可分為節(jié)點信息，單元信息，載荷信息等。對于一個節(jié)點來說，需以下信息：節(jié)點編號（整型），節(jié)點坐標(biāo)（實型），節(jié)點已知位移（實型），節(jié)點載荷（實型），邊界條件（實型）和節(jié)點溫度（實型）等。同樣，對于一個單元來說，需以下信息：單元的節(jié)點聯(lián)接信息（整型），單元類型信息（桁架、梁、板、殼等）（整型） ，單元特性信息（厚度、內(nèi)力矩等）（實型），材料信息（彈性模量，泊松比等）（實型）等。在FORTRAN 程序中，以上這些變量混合在一起，很難辨認(rèn)，使程序的可讀性不好，如需要進行單元網(wǎng)絡(luò)的自適應(yīng)劃分，節(jié)點

5、及單元的修改將非常困難。在進行C語言編譯過程中，采用結(jié)構(gòu)struct 使每個節(jié)點信息存儲在一個結(jié)構(gòu)體數(shù)組中，提高程序的可讀性，使數(shù)據(jù)結(jié)構(gòu)更趨于合理。1.2.2.1. 八節(jié)點矩形單元介紹八節(jié)點矩形單元編號如Error! Reference source not found.所示圖 21八節(jié)點矩形單元的位移函數(shù)為： 其形函數(shù)為 式和式中，并且采用等參變換，則單元的坐標(biāo)變換式可取為 單元應(yīng)變矩陣為 式一般簡寫為 其中的子塊矩陣為 由于是、的函數(shù)，在中的、要按照復(fù)合函數(shù)來求導(dǎo)，即 從而有 因此，單元應(yīng)力矩陣為 單元剛度矩陣為 其中積分采用三點高斯積分， （高斯積分點的總數(shù)），和或是加權(quán)系數(shù)，和是單元內(nèi)

6、的坐標(biāo).。對于三點高斯積分，高斯積分點的位置： ，。單元等效節(jié)點荷載 結(jié)構(gòu)剛度矩陣 結(jié)構(gòu)結(jié)點荷載列陣 注意，對于式和式中的理解不是簡單的疊加而是集成?？倓偲胶夥匠?從式求出 將回代入式和式，得到和。1.2.2.1.2.2. 有限元分析的模塊組織一個典型的有限元分析過程主要包括以下幾個步驟：1) 讀輸入數(shù)據(jù)，定義節(jié)點及單元數(shù)組。2) 由邊界條件計算方程個數(shù)，賦值荷載列陣。3) 讀入在帶狀存儲的總剛度矩陣中單元和載荷信息。4) 定義總剛度陣數(shù)組。5) 組裝總剛度陣。6) 解方程組。輸入邊界條件(對稱條件)形成各荷載工況的節(jié)點荷載陣 總剛分解 回代求出位移及輸出 計算應(yīng)變、應(yīng)力形成單元剛陣單剛向總剛

7、投放坐標(biāo)變換輸入原始參數(shù)計算總剛規(guī)模形成總剛方程向總節(jié)點荷載陣投放形成單元荷載陣調(diào)整幾何、彈性矩陣調(diào)整單元位移列陣圖 22其流程圖可見Error! Reference source not found.。3. 程序變量及函數(shù)說明3.3.3.1. 主要變量說明：主要程序變量說明 wide分析模型的寬 hight分析模型的高wsize寬度方向網(wǎng)格劃分尺寸hsize高度方向網(wǎng)格劃分尺寸npoin節(jié)點總數(shù)nelem單元總數(shù)struct node500節(jié)點結(jié)構(gòu)體變量struct element 500單元結(jié)構(gòu)體變量ifpre500節(jié)點約束信息posgp3高斯積分點weigp3權(quán)系數(shù)gpcod29高斯積分

8、點總體坐標(biāo)bmatx316單元變形矩陣dmatx33單元彈性矩陣xjacm22雅可比矩陣xjaci22雅可比矩陣的逆矩陣djacb雅可比矩陣行列式的值estif136單元剛度矩陣shape8單元形函數(shù)deriv28形函數(shù)對局部坐標(biāo)的導(dǎo)數(shù)cartd28形函數(shù)對整體坐標(biāo)的導(dǎo)數(shù)A30000總體剛度矩陣eload16等效節(jié)點荷載gpcod29高斯積分點的總體坐標(biāo)3.1.3.2. 主要函數(shù)說明主要函數(shù)說明void meshing（ ）網(wǎng)格劃分void coordinate( )節(jié)點整體坐標(biāo)計算void input( )信息輸入void stif( )形成單元剛度矩陣void sfr2（）計算形函數(shù)對當(dāng)前

9、積分點及形函數(shù)對局部坐標(biāo)的到數(shù)值void jacobi2( )形成雅可比矩陣void modps( )形成彈性矩陣void bmatps( )形成應(yīng)變矩陣void load( )計算等效節(jié)點荷載void asem( )形成整體剛度矩陣和整體荷載列陣void solve( )求解整體方程void stre( )計算單元應(yīng)力4. 程序流程圖任何一個有限元程序處理過程，都可以劃分為三個階段：1) 前處理：讀入數(shù)據(jù)和生成數(shù)據(jù)。2) 處理器：有限元的矩陣計算、組集和求解。3) 后處理：輸出節(jié)點位移、單元應(yīng)變等。為了更清楚地描述程序，我們給出了程序的流程圖。 5. 程序應(yīng)用與ANSYS分析的比較為了驗證程

10、序的正確性，我們?nèi)×艘粋€簡單框架結(jié)構(gòu)，分別用自編程序和ANSYS進行分析和比較。4.5.5.1. 問題說明Error! Reference source not found.所示一簡支梁，高3m，長18m，承受均布荷載，取m，作為平面應(yīng)力問題。圖 51圖 52由于結(jié)構(gòu)和荷載對稱，只對右邊一半進行有限單元計算，如Error! Reference source not found.所示，而在y軸各節(jié)點布置水平連桿支座。5.2. ANSYS分析結(jié)果ANSYS分析中，采用plane82單元，在板單元上邊施加均布荷載，在y軸上的各結(jié)點約束UX方向的自由度，在板單元右下角的結(jié)點約束UX自由度，結(jié)點布置、網(wǎng)

11、格劃分方案如Error! Reference source not found.所示。圖 53Error! Reference source not found.和Error! Reference source not found.分別給出了結(jié)構(gòu)的位移圖和應(yīng)力云圖。圖 54 位移圖圖 55 各單元圖從Error! Reference source not found.和Error! Reference source not found.可以看到，跨中的位移和正應(yīng)力最大，與簡支梁承受均布荷載，跨中撓度和正應(yīng)力最大的力學(xué)常識相符合，可以初步認(rèn)為模型是正確的。Error! Reference sou

12、rce not found.給出了的截面上的正應(yīng)力和Error! Reference source not found.給出了處各橫截面方向位移，其中 的單位為，的單位為。表格 1考查點的y/m1.51.00.50-0.5-1.0-1.5正應(yīng)力-270.20-178.25-88.570.0088.57178.25270.21表格 2考查點的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3485-0.3380-0.3069-0.2571-0.1914-0.113305.3. 自編程序分析結(jié)果用自編程序進行分析時，采用了整個結(jié)構(gòu)模型進行計算，其他條件不變，因編程水平有限，在

13、后處理階段沒有給出節(jié)點處的位移與應(yīng)力，而只能得到高斯積分點處的位移與應(yīng)力，得到積分點處的應(yīng)力后，根據(jù)數(shù)值分析相關(guān)知識采用了插值進行處理，從而得到的截面上的正應(yīng)力和 處各橫截面方向位移，分別列出表格如下：表格 3考查點的y/m1.51.00.50-0.5-1.0-1.5正應(yīng)力-297.2-196.06-93.250.0093196.08297.23表格 4考查點的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3481-0.3376-0.3065-0.2568-0.1910-0.112505.4. 結(jié)果分析比較為了更直觀的比較自編程序和ANSYS的計算結(jié)果，將Error!

14、Reference source not found.和Error! Reference source not found.的數(shù)據(jù)繪入Error! Reference source not found.，將Error! Reference source not found.和Error! Reference source not found.的數(shù)據(jù)繪入Error! Reference source not found.。圖 56 應(yīng)力比較圖 57 位移比較 自編程序所得結(jié)果與ANSYS分析結(jié)果進行比較發(fā)現(xiàn)，應(yīng)力偏大而位移偏小。分析其中的原因，一方面是編程水平有限，自編程序有很多不完善之處，有很

15、多因素沒有考慮（溫度、變形、非線性等），所得結(jié)果可信度不是很高，只能得到應(yīng)力以及位移大概的分布情況；另一方面是在后處理階段，在對高斯積分點結(jié)果進行處理時，誤差難以避免，還有一方面的原因是在進行求解時保留數(shù)據(jù)精度不夠，造成誤差較大，并且進行數(shù)值運算時可能會造成大數(shù)吃小數(shù)，從而引起結(jié)果的差異。 參考文獻1 (美)史密斯(Smith,I.m.)著;王崧等譯.有限單元法編程(第三版)M.北京:電子工業(yè)出版社,20032 王勖成.有限單元法M.北京:清華大學(xué)出版社,20033 宋建輝,涂志剛.MATLAB語言及其在有限元編程中的應(yīng)用J.湛江師范學(xué)院學(xué)報.2003.6(24):101-1054 鄭大玉,

16、連宏玉, 周躍發(fā). 有限元編程與C語言J. 黑龍江商學(xué)院學(xué)報.1997.3(13):23-285 王偉,劉德富.有限元編程中應(yīng)用面向?qū)ο缶幊碳夹g(shù)的探討J.三峽大學(xué)學(xué)報.2001.2(23):124-1286 汪文立, 呂士俊.二級C語言程序設(shè)計教程M. 北京:中國水利水電出版社,20067 趙翔龍.FORTRAN 90學(xué)習(xí)教程M.北京:北京大學(xué)出版社,2002附錄 程序源代碼#include<stdio.h>#include<math.h>/*The gemotry of the model*/ float mesh2; float wide,hight; float

17、wsize,hsize,young,poiss,thick; int i,j,k,knelem,npoin,elem500,ielem; float coord21,props31,lnods8500,ifpre1; float posgp3,weigp3,estif136,elcod28; int npoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb; float gpcod29,bmatx316,dmatx33; float shape8,deriv28; float xjacm22,xjaci22,elcod28; float cartd28;

18、float bmatx316,dmatx33; float nload1; float press32,pgash2,dgash2;struct node int nodenu;/*The number of the node*/ float cor2;/*The coordinate of the node*/ float displacement2;/*The displacement of the node*/ float load2; /*The load of the node*/ float boundary2; node500; struct int elementnu;/*Th

19、e number of element*/ int localnu8;/*Local number*/ int localcorx8; int localcory8; /*Local coordinate of node*/ float qx,qy,t;/*The stress and strain*/ element500;void meshing(float wide,float hight,float mesh2) printf("Please input the meshing sizen"); scanf("%f%f",&wsize,&

20、amp;hsize); getchar(); mesh0=wide/wsize; mesh1=hight/hsize; /*The global coordinate of node*/ void coordinate() int i,j=1; float x,y; for(i=1;i<=2*mesh1+1;i+) x=0; while(x<=wide) if(i%2!=0) nodej.cor1+=wsize/2; nodej.cor2+=(i-1)*hsize; j+; else nodej.cor1+=wsize; nodej.cor2+=(i-1)*hsize; j+; m

21、ain() int top500,bottom500;/*The number of top and bottom element*/ void input(); void stif(); void jacobi2( ); void load(); void asem(); void solve( ); void stre(); npoin=8; for(i=1;i<=8;i+) lnodsi1=i;for(i=1;i<=3;i+) scanf("%f",&propsi1); printf("The EX is %enThe uo is %8f

22、nThe bt is %8f",props11,props21,props31); getch(); printf("Please input the gemotry of the modeln"); scanf("%f%f",&wide,&hight); getchar(); printf("The wide is %0.3f,the hight is %0.3f",wide,hight); getchar(); meshing(wide,hight,mesh); printf("The wide

23、 size is %f,the hight size is %fn",mesh0,mesh1); input(); getchar(); nelem=mesh0*mesh1; for(i=1;i<=nelem;i+) /*The element number*/ elementi.elementnu=i; for(i=1;i<mesh0;i+) npoin+=5; for(i=1;i<mesh1;i+) npoin+=3*mesh0+2; printf("%d",npoin); for(i=1;i<=npoin;i+) nodei.node

24、nu=i; for(i=1;i<=nelem;i+) for(j=1;j<=8;j+) elementi.localnuj=j; for (i=1;i<=mesh0;i+) bottomi=i; j=1; for(i=mesh0*(mesh1-1)+1;i<=nelem;i+) topj+=i; for(i=1;i<j;i+) printf("the top number is %dn",topi); printf("%dn",element1.localnu8); coordinate(); stif(); jacobi2

25、( ); load(); asem(); solve( ); stre();getchar();/*data input*/void input() int nnode=8; int notreadchar;/*input element node number*/printf("input element node numbern");for(ielem=1;ielem<=nelem;ielem+) for(i=1;i<=nnode;i+) scanf("%d",&lnodsiielem);/*Input restriction m

26、essage*/printf("double-digit is symbol of the restriction,1 representates stable,2 representates gived displacementn"); i=1; while(i) scanf("%d",&i); if(i!=0) scanf("%d",&j); scanf("%d",&ifprej); else break; /*Element stiffness matrix*/void stif()

27、int kevab,nstre,ievab,istre,lnode,jstre,jevab; float kgasp,exisp,etasp, dvolu; float btdbm,dbmat3; void sfr();/*Gauss point*/ posgp1=-sqrt(0.6); posgp2=0; posgp3=sqrt(0.6); /*weight coefficient*/ weigp1=5.0/9.0; weigp2=8.0/9.0; weigp3=5.0/9.0; /*number of stress*/ nstre=3; /*form elastic matrix*/ fo

28、r(ielem=1;ielem<=nelem;ielem+) printf("The number of element is %dn",ielem);for(i=1;i<=8;i+) lnode=lnodsiielem; for(j=1;j<=2;j+) coordjlnode=nodelnode.corj; elcodji=coordjlnode; young=props11; poiss=props21; thick=props31; modps(young,poiss); /*The initial value of k*/ kevab=0; fo

29、r(i=1;i<=16;i+) for(j=1;j<=i;j+) kevab+; estifkevab=0.0; /*The gauss point shape function and deriv*/ kgasp=0; for(igaus=1;igaus<=3;igaus+) for(jgaus=1;jgaus<=3;jgaus+) kgasp=kgasp+1; exisp=posgpigaus; etasp=posgpjgaus; printf("The number of gauss point is %dn",kgasp); sfr2(exi

30、sp,etasp); jacob2(ielem,djacb,kgasp,elcod); dvolu=djacb*weigpigaus*weigpjgaus*thick; bmatps()/*The initial value of elastic matrix D*/ kevab=0; for(ievab=1;ievab<=16;ievab+) for(istre=1;istre<=nstre;istre+) dbmatistre=0.0; for(jstre=1;jstre<=nstre;jstre+) dbmatistre=dbmatistre+bmatxjstreiev

31、ab*dmatxjstreistre; for(jevab=1;jevab<=ievab;ievab+) kevab=kevab+1; btdbm=0.0; for(istre=1;istre<=nstre;istre+) btdbm=btdbm+dbmatistre*bmatxistrejevab; estifkevab=estifkevab+btdbm*dvolu; /*float sfr2(float s,float t) float s2,t2,ss,tt,st,stt,sst,st2;/*Shape function*/*these variables are to si

32、mplify the formula */ s2=s*2.0; t2=t*2.0; ss=s*s; tt=t*t; st=s*t; stt=s*t*t; sst=s*s*t; st2=s*t*2.0; /*shape function*/ shape1=(-1.0+st+ss+tt-sst-stt)/4.0; shape2=(1.0-t-ss+sst)/2.0; shape3=(-1.0-st+ss+tt-sst+stt)/4.0; shape4=(1.0+s-tt-stt)/2.0; shape5=(-1.0+st+ss+tt+sst+stt)/4.0; shape6=(1.0+t-ss-s

33、st)/2.0; shape7=(-1.0-st+ss+tt+sst-stt)/4.0; shape8=(1.0-s-tt+stt)/2.0; /* differential coefficient to local coordinate*/ deriv11=(t+s2-st2-tt)/4.0; deriv12=-s+st; deriv13=(-t+s2-st2+tt)/4.0; deriv14=(1.0-tt)/2.0; deriv15=(t+s2+st2+tt)/4.0; deriv16=-s-st; deriv17=(-t+s2+st2-tt)/4.0; deriv18=(-1.0+tt

34、)/2.0; deriv21=(s+t2-ss-st2)/4.0; deriv22=(-1.0+ss)/2.0; deriv23=(-s+t2-ss+st2)/4.0; deriv24=-t-st; deriv25=(s+t2+ss+st2)/4.0; deriv26=(1.0-ss)/2.0; deriv27=(-s+t2+ss-st2)/4.0; deriv28=-t+st; */*Jacobi matrix*/void jacobi2( ) /* xjacm22:jacobi matrix;xjaci22:j-1;elcod28:local coordinate*/ float cart

35、d28,gpcod29; int i,j,k; for(i=1;i<=2;i+) gpcodikgasp=0.0; for(j=1;j<=8;j+) gpcodikgasp=gpcodikgasp+elcodij*shapej; for(i=1;i<=2;i+) for(j=1;j<=2;j+) xjacmij=0.0; for(k=1;k<=8;k+) xjacmij=xjacmij+derivik*elcodjk; /*Caculate the value of Jacobi matrix*/ djacb=xjacm11*xjacm22-xjacm12*xja

36、cm21; printf("The det of jacobi matrix is %fn",djacb); if(djacb<=1e-6) printf("The jacobi det of %d is less or equal 0n",ielem); /*The element of j-1*/ xjaci11=xjacm22/djacb; xjaci22=xjacm11/djacb; xjaci12=-xjacm12/djacb; xjaci21=-xjacm21/djacb; /*The deriv of shape function t

37、o global coordinate*/ for(i=1;i<=2;i+) for(k=1;k<=8;k+) cartdik=0.0; for(j=1;j<=2;j+) cartdik=cartdik+xjaciij*derivjk; /* Form elastic matris */*float modps() float bmatx316,dmatx33; float constant; int i,j,y,p; y=props11; p=props21; for(i=1;i<=3;i+) for(j=1;j<=3;j+) dmatxij=0.0; /*If

38、 Ntype=1,it is plane stress question*/constant=y/(1.0-p*p); dmatx11=constant; dmatx22=constant; dmatx12=constant*p; dmatx21=constant*p; dmatx33=constant*(1.0-p)/2.0; */*Form strain matrix*/*void bmatps() float cartd28,shape8,deriv28,bmatx316,dmatx33; int npoin,nelem,ngash,mgash; float gpcod29;/*The

39、initial value ofB*/ for(i=1;i<=16;i+) for(j=1;j<=3;j+) bmatxji=0.0; /*Form B*/ ngash=0; for(i=1;i<=8;i+) mgash=ngash+1; ngash=mgash+1; bmatx1mgash=cartd1i; bmatx1ngash=0.0; bmatx2mgash=0.0; bmatx2ngash=cartd2i; bmatx3mgash=cartd2i; bmatx3ngash=cartd1i; */ /*equal load of node*/void load()fl

40、oat nload500; float elcod23,eload16,posgp3,weigp3; float press32,pgash2,dgash2,noprs3; int iedge,nedge,kount; float sgmar,tau,s,dvolu,pxcom,pycom,n,m,t; int lnode,number,nloca; /*The number of load side*/ /*element load*/ /*The position of gauss point*/ printf("input the number of load side*/n&

41、quot;); scanf("%d",&number); posgp1=-sqrt(0.6); posgp2=0.0; posgp3=sqrt(0.6); /*The weight of gauss point*/ weigp1=5.0/9.0; weigp2=8.0/9.0; weigp3=5.0/9.0; /*The initial load of element*/ for(i=1;i<=nelem;i+) nloadi=0.0; for(i=1;i<=number;i+) /*The initial value of equal node loa

42、d*/ for(i=1;i<=16;i+) eloadi=0; scanf("%d",&nedge); for(iedge=1;iedge<=nedge;iedge+) for(i=1;i<=3;i+) scanf("%f%f",&sgmar,&tau); /*if sgmar=0,input node load q(1,2,3)and t(1,2,3)*/ if(fabs(sgmar)<1e-8)&&(fabs(tau)<1e-8) for(j=1;j<=3;j+) for(i

43、=1;i<=2;i+) scanf("%f",&pressji); else for(i=1;i<=3;i+) pressi1=sgmar; pressi2=tau; /*The coordinate of the load node*/ for(i=1;i<=3;i+) /*input load node*/ printf("input node load numbern"); scanf("%d",&noprsi); lnode=noprsi; for(j=1;j<=2;j+) coordj

44、lnode=nodelnode.corj; elcodji=coordjlnode; t=-1.0; for(igaus=1;igaus<=3;igaus+) s=posgpigaus; sfr(s,t); for(j=1;j<=2;j+) pgashj=0.0; dgashj=0.0; /*current point q,t and the value of deriv*/ pgashj=pgashj+pressij*shapei; dgashj=dgashj+elcodji*deriv1i; dvolu=weigpigaus; pxcom=dgash1*pgash2-dgash2*pgash1; pycom=dgash1*pgash1+dgash2*pgash2; for(i=1;i<=8;i+) nloca=lnodsiielem; if(nloca=noprs1) j=i+2; break; j=i+2; kount=0; for(k=i;k<=j;k+) kount=kount+1; n=(k-1)*2+1; m=n+1; /*if the local node numner >8,then it is 1*/ if(k>8)