平面三角形單元有限元程序設(shè)計(jì)

1、、題目如圖1所示，一個厚度均勻的三角形薄板，在頂點(diǎn)作用沿板厚方向均勻分布的豎向載荷。已知：P=150N/m,E=200GPa,口=0.25,t=0.1m,忽略自重。試計(jì)算薄板的位移及應(yīng)力分布。1 .編寫有限元計(jì)算機(jī)程序，計(jì)算節(jié)點(diǎn)位移及單元應(yīng)力。（劃分三角形單元，單元數(shù)不得少于30個）；2 .采用有限元軟件分析該問題（有限元軟件網(wǎng)格與程序設(shè)計(jì)網(wǎng)格必須一致），詳細(xì)給出有限元軟件每一步的操作過程，并將結(jié)果與程序計(jì)算結(jié)果進(jìn)行對比（任選取三個點(diǎn)，對比位移值）；3 .提交程序編寫過程的詳細(xì)報(bào)告及計(jì)算機(jī)程序；4 .所有同學(xué)參加答辯，并演示有限元計(jì)算程序。有限元法中三節(jié)點(diǎn)三角形分析結(jié)構(gòu)的步驟如下：1）整理原始

2、數(shù)據(jù)，如材料性質(zhì)、荷載條件、約束條件等，離散結(jié)構(gòu)并進(jìn)行單元編碼、結(jié)點(diǎn)編碼、結(jié)點(diǎn)位移編碼、選取坐標(biāo)系。2）單元分析，建立單元剛度矩陣。3）整體分析，建立總剛矩陣。4）建立整體結(jié)構(gòu)的等效節(jié)點(diǎn)荷載和總荷載矩陣5）邊界條件處理。6）解方程，求出節(jié)點(diǎn)位移。7）求出各單元的單元應(yīng)力。8）計(jì)算結(jié)果整理。、程序設(shè)計(jì)網(wǎng)格劃分如圖，將薄板如圖劃分為6行，并建立坐標(biāo)系，則單元編號節(jié)點(diǎn)編號剛度矩陣的集成建立與總剛度矩陣等維數(shù)的空矩陣，已變單元剛度矩陣的集成。由單元分析已知節(jié)點(diǎn)、單元的排布規(guī)律，繼而通過循環(huán)計(jì)算求得每個單元對應(yīng)的節(jié)點(diǎn)序號。通過循環(huán)逐個計(jì)算：(1)每個單元對應(yīng)2種單元剛度矩陣中的哪一種；(2)該單元對應(yīng)總

3、剛度矩陣的那幾行哪幾列(3)將該單元的單元剛度矩陣加入總剛度矩陣的對應(yīng)行列循環(huán)又分為3層循環(huán)：(1)最外層：逐行計(jì)算(2)中間層：該行逐個計(jì)算(3)最里層：區(qū)分為第奇/偶數(shù)個計(jì)算Ke66k15656ke66k2'5656ke66kZ'5656Kke'ke'kZ'單元剛度的集成：邊界約束的處理：劃0置1法適用：這種方法適用于邊界節(jié)點(diǎn)位移分量為已知(含為0)的各種約束。做法：(1) 將總剛矩陣K中相應(yīng)于已知位移行主對角線元素置1,其他元素改為零；同時將載荷列陣R中相應(yīng)元素用已知位移置換。這樣，由該方程求得的此位移值一定等于已知量。(2) 將K中已知位移相應(yīng)的

4、列的非主對角成元素也置0,以保持K的對稱性。當(dāng)然，在已知位移分量不為零的情況下，這樣做就改變了方程左端的數(shù)值，為保證方程成立，須在方程右端減去已知位移對該方程的貢獻(xiàn)已知位移和相應(yīng)總剛元素的乘積。若約束為零位移約束時，此步則可省去。特點(diǎn)：(1) 經(jīng)以上處理同樣可以消除剛性位移(約束足夠的前提下)，去掉未知約束反力。(2) 但這種方法不改變方程階數(shù)，利于存貯。(3) 不過 元素。若是要求出約束反力，仍要重新計(jì)算各個劃去的總剛程序如下： 變量說明 NNODE NPION NELEM NVFIX FIXED NFORCE FORCE COORD LNODS YOUNG POISS THICKB DS

5、AESTIF ASTIF ASLOD ASDISP ELEDISP STRESS單元節(jié)點(diǎn)數(shù)總結(jié)點(diǎn)數(shù)單元數(shù)受約束邊界點(diǎn)數(shù)約束信息數(shù)組節(jié)點(diǎn)力數(shù)節(jié)點(diǎn)力數(shù)組結(jié)構(gòu)節(jié)點(diǎn)坐標(biāo)數(shù)組單元定義數(shù)組彈性模量泊松比厚度單元應(yīng)變矩陣(3*6)單元彈性矩陣(3*3)單元應(yīng)力矩陣(3*6)單元面積單元剛度矩陣總體剛度矩陣總體荷載向量節(jié)點(diǎn)位移向量單元節(jié)點(diǎn)位移向量單元應(yīng)力%*%初 始 化 clearformat short e clearNELEM=36 NPION=28 NVFIX=2 NFORCE=1YOUNG=2e11%設(shè)定輸出類型%清除內(nèi)存變量%單元個數(shù)(單元編碼總數(shù))%結(jié)點(diǎn)個數(shù)(結(jié)點(diǎn)編碼總數(shù))%受約束邊界點(diǎn)數(shù)%結(jié)點(diǎn)荷

6、載個數(shù)%彈性模量POISS=0.25%泊松比THICK=0.1%厚度LNODS=123;245;253;356;478;485;589;596;6910;71112;7128;81213;8139;91314;91410;101415;111617;111712;121718;121813;131819;131914;141920;142015;152021;162223;162317;172324;172418;182425;182519;251926;192620;202627;202721;212728%單元定義數(shù)組(單元結(jié)點(diǎn)號)%相應(yīng)為單元結(jié)點(diǎn)號(編碼)、按逆時針順序輸入COORD=0

7、0;-0.751.5;0.751.5;-1.53;03;1.53;-2.254.5;-0.754.5;0.754.5;2.254.5;-36;-1.56;06;1.56;36;-3.757.5;-2.257.5;-0.757.5;0.757.5;2.257.5;3.757.5;-4.59;-39;-1.59;09;1.59;39;4.59%結(jié)點(diǎn)坐標(biāo)數(shù)組%坐標(biāo)：x，y坐標(biāo)(共NPOIN組)FORCE=10-15%結(jié)點(diǎn)力數(shù)組(受力結(jié)點(diǎn)編號,x方向,y方向)FIXED=2211;2811%約束信息(約束點(diǎn)，x約束，y約束)%有約束為1，無約束為0%*%生成單元剛度矩陣并組成總體剛度矩陣ASTIF=z

8、eros(2*NPION,2*NPION);%生成特定大小總體剛度矩陣并置0%*fori=1:NELEM%生成彈性矩陣DD= 1POISSPOISS 1000;0;(1-POISS)/2*YOUNG/(1-POISSA2)%*%計(jì)算當(dāng)前單元的面積A=-det(1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)/2%*%生成應(yīng)變矩陣Bforj=0:2-COORD(LNODS(i,(reb(j+1)=COORD(

9、LNODS(i,(rem(j+1),3)+1),2)m(j+2),3)+1),2);c(j+1)=-COORD(LNODS(i,(rem(j+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1);end0;c(3);b(3)/(2*A);B=b(1)0b(2)0b(3)0c(1)0c(2)0c(1)b(1)c(2)b(2)c(3)B1(:,:,i)=B;%*%求應(yīng)力矩陣S=D*BS=D*B;%求 解 單 元 剛 度 矩 陣%臨 時 向 量 , 用 來 記 錄 當(dāng) 前ESTIF=B'*S*THICK*A;a=LNODS(i,:);單元的節(jié)點(diǎn)編號forj

10、=1:3fork=1:3ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF(a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);%根據(jù)節(jié)點(diǎn)編號對應(yīng)關(guān)系將單元剛度分塊疊加到總剛%度矩陣中endendend%*%將約束信息加入總體剛度矩陣(對角元素改一法)fori=1:NVFIXifFIXED(i,2)=1ASTIF(:,(FIXED(i,1)*2-1)=0;%一列為零ASTIF(FIXED(i,1)*2-1),:)=0;%一行為零ASTIF(FIXED(i,1)*2-1),(FI

11、XED(i,1)*2-1)=1;%對角元素為1end%*%生成單元剛度矩陣并組成總體剛度矩陣*ifFIXED(i,3)=1ASTIF(:,FIXED(i,1)*2)=0;%一列為零ASTIF(FIXED(i,1)*2,:)=0;%一行為零ASTIF(FIXED(i,1)*2,FIXED(i,1)*2)=1;%對角元素為1endend%*%生成荷載向量ASLOD(1:2*NPION)=0;%總體荷載向量置零fori=1:NFORCEASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%*%求解內(nèi)力ASDISP=ASTIFASLOD'%計(jì)

12、算節(jié)點(diǎn)位移向量ELEDISP(1:6)=0;%當(dāng)前單元節(jié)點(diǎn)位移向量fori=1:NELEMforj=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出當(dāng)前單元的節(jié)點(diǎn)位移向量endiSTRESS=D*B1(:,:,i)*ELEDISP'%求內(nèi)力end（程序計(jì)算結(jié)果和有限元軟件得出的結(jié)果稍有偏差，可能是程序某些地方數(shù)據(jù)輸入時出了問題，還在尋找具體原因）有限元軟件分析¥：ni)GraphicsPropertiesA-DefinitionSuppr蚪*dN。Stiffn,曰5BehaviorFlexibleCoo

13、rdinateSystemDefault'toordmateSystemReferenceTemperatureByEnvironrri'ertThickness100.mmThieknesl/odeManualOffsetTypeMiddle-MaterialAs5ignmntStructuralSteelMonlinearEffectsYesVof''SurfaceBody"建模HlProject0Model(A4/除Geonery76Sr三二圮Body,不CoordinateSystems-j,%GobaGoordiratesteri白螳Mesl

14、h,配AlInanglesktet*iod日白staticStHKtural(A5)*sE、0.Analysis烹ttig5B?每Sokjtkxi(A6)-JSdubcnIrifonratoriIrujniE©Ho惻網(wǎng)日、的Geonery二4埒知依已修叫j-yAMraxSarsj-j*GbbdWnWESystrn小博*di感打Tfiands米加!：L-A印落9利0臼StaticStructuralAS)l工、5al芯eSettngs三值Sohithn.J而UtwMrnatMDrtaiso;'AllTrianglesMethod,-MetnodScapeScopingMetho

15、dGecrrtr>SelectionGeometryBodyDefinntiofiSuppressedNoMethodTnanglesBementMidsideNodesU5fGotHl5dtnsT|)Maiknf'EdgeSizing*-Sizing事-5»peScopirgM近odfeem(t7SectionGecMry3EdgesDefi喻麗5叩P'MRCrjcTypeMurr陽dMiiarisNum:irjf-Mii3r$31BehaviorHardB«TypeNoBias網(wǎng)格劃分：添加載荷邊界約束A：StaticStructuralEquSontStressTypczEqjrvelcnt(xon-Miscs)$trc”-*op/DuttQinUnt:MPjiTimo:12017/1/711:140.03918910.0D31