八節(jié)點(diǎn)等參元

1、計(jì)算力學(xué)課程大作業(yè)八節(jié)點(diǎn)等參單元平面有限元分析程序土木工程學(xué)院2011.2目錄1.概述12.編程思想22.1.八節(jié)點(diǎn)矩形單元介紹22.2.有限元分析的模塊組織53.程序變量及函數(shù)說明63.1.主要變量說明：63.2.主要函數(shù)說明74.程序流程圖85.程序應(yīng)用與ANSYS分析的比較95.1.問題說明95.2.ANSYS分析結(jié)果105.3.自編程序分析結(jié)果125.4.結(jié)果比較分析12參考文獻(xiàn)14附錄 程序源代碼15計(jì)算力學(xué)課程大作業(yè)1. 概述通常情況下的有限元分析過程是運(yùn)用可視化分析軟件(如ANSYS、SAP等)進(jìn)行前處理和后處理，而中間的計(jì)算部分一般采用自己編制的程序來(lái)運(yùn)算。具有較強(qiáng)數(shù)值計(jì)算和處

2、理能力的Fortran語(yǔ)言是傳統(tǒng)有限元計(jì)算的首選語(yǔ)言。隨著有限元技術(shù)的逐步成熟,它被應(yīng)用在越來(lái)越復(fù)雜的問題處理中，但在實(shí)際應(yīng)用中也暴露出一些問題。有時(shí)網(wǎng)格離散化的區(qū)域較大,而又限于研究精度的要求，使得劃分的網(wǎng)格數(shù)目極其龐大，結(jié)點(diǎn)數(shù)可多達(dá)數(shù)萬(wàn)個(gè)，從而造成計(jì)算中要運(yùn)算的數(shù)據(jù)量巨大，程序運(yùn)行的時(shí)間較長(zhǎng)的弊端，這就延長(zhǎng)了問題解決的時(shí)間，使得求解效率降低。因?yàn)檫\(yùn)行周期長(zhǎng)，不利于程序的調(diào)試，特別是對(duì)于要計(jì)算多種運(yùn)行工況時(shí)的情況；同時(shí)大數(shù)據(jù)量處理對(duì)計(jì)算機(jī)的內(nèi)存和CPU 提出了更高的要求，而在實(shí)際應(yīng)用中，單靠計(jì)算機(jī)硬件水平的提高來(lái)解決問題的能力是有限的。因此，必須尋找新的編程語(yǔ)言。隨著有限元前后處理的不斷發(fā)展

3、和完善，以及大型工程分析軟件對(duì)有限元接口的要求，有限元分析程序不應(yīng)只滿足解題功能，它還應(yīng)滿足軟件工程所要求的結(jié)構(gòu)化程序設(shè)計(jì)條件，能夠?qū)Υ鎯?chǔ)進(jìn)行動(dòng)態(tài)分配，以充分利用計(jì)算機(jī)資源，它還應(yīng)很容易地與其它軟件如CAD 的實(shí)體造型，優(yōu)化設(shè)計(jì)等接口?，F(xiàn)在可編寫工程應(yīng)用軟件的計(jì)算機(jī)語(yǔ)言較多，其中C語(yǔ)言是一個(gè)較為優(yōu)秀的語(yǔ)言，很容易滿足現(xiàn)在有限元分析程序編程的要求。C語(yǔ)言最初是為操作系統(tǒng)、編譯器以及文字處理等編程而發(fā)明的。隨著不斷完善，它已應(yīng)用到其它領(lǐng)域，包括工程應(yīng)用軟件的編程。近年來(lái)，C語(yǔ)言已經(jīng)成為計(jì)算機(jī)領(lǐng)域最普及的一個(gè)編程語(yǔ)言，幾乎世界上所有的計(jì)算機(jī)都裝有C的編譯器，從PC機(jī)到巨型機(jī)到超巨型的并行機(jī)，C與所有

4、的硬件和操作系統(tǒng)聯(lián)系在一起。用C 編寫的程序，可移植性極好，幾乎不用作多少修改，就可在任何一臺(tái)裝有ANSI、C編譯器的計(jì)算機(jī)上運(yùn)行。C既是高級(jí)語(yǔ)言，也是低級(jí)語(yǔ)言，也就是說，可用它作數(shù)值計(jì)算，也可用它對(duì)計(jì)算機(jī)存儲(chǔ)進(jìn)行操作。2. 編程思想本程序采用C語(yǔ)言編程，編制平面四邊形八節(jié)點(diǎn)等參元程序，用以求解平面結(jié)構(gòu)問題。程序采用二維等帶寬存儲(chǔ)整體剛度矩陣，乘大數(shù)法引入約束，等帶寬高斯消去法求解位移。在有限元程序中，變量數(shù)據(jù)需賦值的可分為節(jié)點(diǎn)信息，單元信息，載荷信息等。對(duì)于一個(gè)節(jié)點(diǎn)來(lái)說，需以下信息：節(jié)點(diǎn)編號(hào)（整型），節(jié)點(diǎn)坐標(biāo)（實(shí)型），節(jié)點(diǎn)已知位移（實(shí)型），節(jié)點(diǎn)載荷（實(shí)型），邊界條件（實(shí)型）和節(jié)點(diǎn)溫度（實(shí)型）

5、等。同樣，對(duì)于一個(gè)單元來(lái)說，需以下信息：?jiǎn)卧墓?jié)點(diǎn)聯(lián)接信息（整型），單元類型信息（桁架、梁、板、殼等）（整型） ，單元特性信息（厚度、內(nèi)力矩等）（實(shí)型），材料信息（彈性模量，泊松比等）（實(shí)型）等。在FORTRAN 程序中，以上這些變量混合在一起，很難辨認(rèn)，使程序的可讀性不好，如需要進(jìn)行單元網(wǎng)絡(luò)的自適應(yīng)劃分，節(jié)點(diǎn)及單元的修改將非常困難。在進(jìn)行C語(yǔ)言編譯過程中，采用結(jié)構(gòu)struct 使每個(gè)節(jié)點(diǎn)信息存儲(chǔ)在一個(gè)結(jié)構(gòu)體數(shù)組中，提高程序的可讀性，使數(shù)據(jù)結(jié)構(gòu)更趨于合理。1.2.2.1. 八節(jié)點(diǎn)矩形單元介紹八節(jié)點(diǎn)矩形單元編號(hào)如圖 21所示圖 21八節(jié)點(diǎn)矩形單元的位移函數(shù)為： 其形函數(shù)為 式和式中，并且采用等參

6、變換，則單元的坐標(biāo)變換式可取為 單元應(yīng)變矩陣為 式一般簡(jiǎn)寫為 其中的子塊矩陣為 由于是、的函數(shù)，在中的、要按照復(fù)合函數(shù)來(lái)求導(dǎo)，即 從而有 因此，單元應(yīng)力矩陣為 單元?jiǎng)偠染仃嚍?其中積分采用三點(diǎn)高斯積分， （高斯積分點(diǎn)的總數(shù)），和或是加權(quán)系數(shù)，和是單元內(nèi)的坐標(biāo).。對(duì)于三點(diǎn)高斯積分，高斯積分點(diǎn)的位置： ，。單元等效節(jié)點(diǎn)荷載 結(jié)構(gòu)剛度矩陣 結(jié)構(gòu)結(jié)點(diǎn)荷載列陣 注意，對(duì)于式和式中的理解不是簡(jiǎn)單的疊加而是集成。總剛平衡方程 從式求出 將回代入式和式，得到和。1.2.2.1.2.2. 有限元分析的模塊組織一個(gè)典型的有限元分析過程主要包括以下幾個(gè)步驟：1) 讀輸入數(shù)據(jù)，定義節(jié)點(diǎn)及單元數(shù)組。2) 由邊界條件計(jì)算

7、方程個(gè)數(shù)，賦值荷載列陣。3) 讀入在帶狀存儲(chǔ)的總剛度矩陣中單元和載荷信息。4) 定義總剛度陣數(shù)組。5) 組裝總剛度陣。6) 解方程組。輸入邊界條件(對(duì)稱條件)形成各荷載工況的節(jié)點(diǎn)荷載陣 總剛分解 回代求出位移及輸出 計(jì)算應(yīng)變、應(yīng)力形成單元?jiǎng)傟噯蝿傁蚩倓偼斗抛鴺?biāo)變換輸入原始參數(shù)計(jì)算總剛規(guī)模形成總剛方程向總節(jié)點(diǎn)荷載陣投放形成單元荷載陣調(diào)整幾何、彈性矩陣調(diào)整單元位移列陣圖 22其流程圖可見圖 22。3. 程序變量及函數(shù)說明3.3.3.1. 主要變量說明：主要程序變量說明 wide分析模型的寬 hight分析模型的高wsize寬度方向網(wǎng)格劃分尺寸hsize高度方向網(wǎng)格劃分尺寸npoin節(jié)點(diǎn)總數(shù)nele

8、m單元總數(shù)struct node500節(jié)點(diǎn)結(jié)構(gòu)體變量struct element 500單元結(jié)構(gòu)體變量ifpre500節(jié)點(diǎn)約束信息posgp3高斯積分點(diǎn)weigp3權(quán)系數(shù)gpcod29高斯積分點(diǎn)總體坐標(biāo)bmatx316單元變形矩陣dmatx33單元彈性矩陣xjacm22雅可比矩陣xjaci22雅可比矩陣的逆矩陣djacb雅可比矩陣行列式的值estif136單元?jiǎng)偠染仃噑hape8單元形函數(shù)deriv28形函數(shù)對(duì)局部坐標(biāo)的導(dǎo)數(shù)cartd28形函數(shù)對(duì)整體坐標(biāo)的導(dǎo)數(shù)A30000總體剛度矩陣eload16等效節(jié)點(diǎn)荷載gpcod29高斯積分點(diǎn)的總體坐標(biāo)3.1.3.2. 主要函數(shù)說明主要函數(shù)說明void

9、meshing（ ）網(wǎng)格劃分void coordinate( )節(jié)點(diǎn)整體坐標(biāo)計(jì)算void input( )信息輸入void stif( )形成單元?jiǎng)偠染仃噕oid sfr2（）計(jì)算形函數(shù)對(duì)當(dāng)前積分點(diǎn)及形函數(shù)對(duì)局部坐標(biāo)的到數(shù)值void jacobi2( )形成雅可比矩陣void modps( )形成彈性矩陣void bmatps( )形成應(yīng)變矩陣void load( )計(jì)算等效節(jié)點(diǎn)荷載void asem( )形成整體剛度矩陣和整體荷載列陣void solve( )求解整體方程void stre( )計(jì)算單元應(yīng)力4. 程序流程圖任何一個(gè)有限元程序處理過程，都可以劃分為三個(gè)階段：1) 前處理：讀入數(shù)

10、據(jù)和生成數(shù)據(jù)。2) 處理器：有限元的矩陣計(jì)算、組集和求解。3) 后處理：輸出節(jié)點(diǎn)位移、單元應(yīng)變等。為了更清楚地描述程序，我們給出了程序的流程圖。 5. 程序應(yīng)用與ANSYS分析的比較為了驗(yàn)證程序的正確性，我們?nèi)×艘粋€(gè)簡(jiǎn)單框架結(jié)構(gòu)，分別用自編程序和ANSYS進(jìn)行分析和比較。4.5.5.1. 問題說明圖 51所示一簡(jiǎn)支梁，高3m，長(zhǎng)18m，承受均布荷載，取m，作為平面應(yīng)力問題。圖 51圖 52由于結(jié)構(gòu)和荷載對(duì)稱，只對(duì)右邊一半進(jìn)行有限單元計(jì)算，如圖 52所示，而在y軸各節(jié)點(diǎn)布置水平連桿支座。5.2. ANSYS分析結(jié)果ANSYS分析中，采用plane82單元，在板單元上邊施加均布荷載，在y軸上的各結(jié)

11、點(diǎn)約束UX方向的自由度，在板單元右下角的結(jié)點(diǎn)約束UX自由度，結(jié)點(diǎn)布置、網(wǎng)格劃分方案如圖 53所示。圖 53圖 54 位移圖和圖 55 各單元圖分別給出了結(jié)構(gòu)的位移圖和應(yīng)力云圖。圖 54 位移圖圖 55 各單元圖從圖 54 位移圖和圖 55 各單元圖可以看到，跨中的位移和正應(yīng)力最大，與簡(jiǎn)支梁承受均布荷載，跨中撓度和正應(yīng)力最大的力學(xué)常識(shí)相符合，可以初步認(rèn)為模型是正確的。表格 1給出了的截面上的正應(yīng)力和表格 2給出了處各橫截面方向位移，其中 的單位為，的單位為。表格 1考查點(diǎn)的y/m1.51.00.50-0.5-1.0-1.5正應(yīng)力-270.20-178.25-88.570.0088.57178.2

12、5270.21表格 2考查點(diǎn)的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3485-0.3380-0.3069-0.2571-0.1914-0.113305.3. 自編程序分析結(jié)果用自編程序進(jìn)行分析時(shí)，采用了整個(gè)結(jié)構(gòu)模型進(jìn)行計(jì)算，其他條件不變，因編程水平有限，在后處理階段沒有給出節(jié)點(diǎn)處的位移與應(yīng)力，而只能得到高斯積分點(diǎn)處的位移與應(yīng)力，得到積分點(diǎn)處的應(yīng)力后，根據(jù)數(shù)值分析相關(guān)知識(shí)采用了插值進(jìn)行處理，從而得到的截面上的正應(yīng)力和 處各橫截面方向位移，分別列出表格如下：表格 3考查點(diǎn)的y/m1.51.00.50-0.5-1.0-1.5正應(yīng)力-297.2-196.06-93.25

13、0.0093196.08297.23表格 4考查點(diǎn)的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的計(jì)算結(jié)果，將表格 1和表格 3的數(shù)據(jù)繪入圖 56，將表格 2和表格 4的數(shù)據(jù)繪入圖 57。圖 56 應(yīng)力比較圖 57 位移比較 自編程序所得結(jié)果與ANSYS分析結(jié)果進(jìn)行比較發(fā)現(xiàn)，應(yīng)力偏大而位移偏小。分析其中的原因，一方面是編程水平有限，自編程序有很多不完善之處，有很多因素沒有考慮（溫度、變形、非線性等），所得結(jié)果可信度不是很高，只能得到

14、應(yīng)力以及位移大概的分布情況；另一方面是在后處理階段，在對(duì)高斯積分點(diǎn)結(jié)果進(jìn)行處理時(shí)，誤差難以避免，還有一方面的原因是在進(jìn)行求解時(shí)保留數(shù)據(jù)精度不夠，造成誤差較大，并且進(jìn)行數(shù)值運(yùn)算時(shí)可能會(huì)造成大數(shù)吃小數(shù)，從而引起結(jié)果的差異。 參考文獻(xiàn)1 (美)史密斯(Smith,I.m.)著;王崧等譯.有限單元法編程(第三版)M.北京:電子工業(yè)出版社,20032 王勖成.有限單元法M.北京:清華大學(xué)出版社,20033 宋建輝,涂志剛.MATLAB語(yǔ)言及其在有限元編程中的應(yīng)用J.湛江師范學(xué)院學(xué)報(bào).2003.6(24):101-1054 鄭大玉, 連宏玉, 周躍發(fā). 有限元編程與C語(yǔ)言J. 黑龍江商學(xué)院學(xué)報(bào).1997.

15、3(13):23-285 王偉,劉德富.有限元編程中應(yīng)用面向?qū)ο缶幊碳夹g(shù)的探討J.三峽大學(xué)學(xué)報(bào).2001.2(23):124-1286 汪文立, 呂士俊.二級(jí)C語(yǔ)言程序設(shè)計(jì)教程M. 北京:中國(guó)水利水電出版社,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 wsize,hsize,young,poiss,thick; int i

16、,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; float bmatx316,dmatx33; float nload1

17、; 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;/*The number of element*/ int localnu8;/

18、*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,&hsize); getchar(); mesh0=wide/ws

19、ize; 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+; main() int top500,bottom500;/*The num

20、ber 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 %8fnThe bt is %8f",props11,props21

21、,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 size is %f,the hight size is %fn&qu

22、ot;,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.nodenu=i; for(i=1;i<=nelem;i+) for(j=

23、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( ); load(); asem(); solve( ); stre(

24、);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 message*/printf("double-digit is

25、 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() int kevab,nstre,ievab,istre,lnode,js

26、tre,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*/ for(ielem=1;ielem<=nelem;ielem+) pr

27、intf("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; for(i=1;i<=16;i+) for(j=1;j<=i;j

28、+) 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(exisp,etasp); jacob2(ielem,djacb,kgasp,

29、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+bmatxjstreievab*dmatxjstreistre; for(jevab=1;jeva

30、b<=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 simplify the formula */ s2=s*2.0; t2=t

31、*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-sst)/2.0; shape7=(-1.0-st+ss+tt+sst-s

32、tt)/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)/2.0; deriv21=(s+t2-ss-st2)/4.0; de

33、riv22=(-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 cartd28,gpcod29; int i,j,k; for(i=1;i<

34、;=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*xjacm21; printf("The det of jacobi

35、 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 to global coordinate*/ for(i=1;i<=

36、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 Ntype=1,it is plane stress question

37、*/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 initial value ofB*/ for(i=1;i<=16

38、;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()float nload500; float elcod23,eload16,

39、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"); scanf("%d",&nu

40、mber); 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 load*/ for(i=1;i<=16;i+) eloadi=0; s

41、canf("%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=1;i<=2;i+) scanf("%f",

42、&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+) coordjlnode=nodelnode.corj; elcodji=coordj

43、lnode; 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

44、*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) n=1; m=2; /*sum to get the equal node load*/ eloadn=eloadn+shapekount*pxcom*dvolu; eloadm=eloadm+shapekount*pycom*dvolu; print("The element load