有限元三角形單元程序設計PPT學習教案

1、會計學1 有限元三角形單元程序設計有限元三角形單元程序設計 6:19 (0,0)A (0, 2)B (4,0)D (4,1)C 1F q 考慮一個平面應力問題如圖所示，假設厚度h=1，材料為各 項同性，楊氏模量為E=1，泊松比為=0，相關力和位移邊界 條件如圖中所示，問題左端為固定約束。試用兩個三角形單 元分析此問題，三角形單元的網格劃分如圖所示。試求問題 各節(jié)點位移u、v和應力x，y和xy。 問題描述 第1頁/共21頁 6:19 三角形單元程序 (0,0)A (0, 2)B (4,0)D (4,1)C 1F 第2頁/共21頁 6:19 三角形單元程序 第3頁/共21頁 6:19 三角形單元程

2、序 第4頁/共21頁 6:19 三角形單元程序 (0,0)A (0, 2)B (4,0)D (4,1)C 第5頁/共21頁 6:19 三角形單元程序 第6頁/共21頁 6:19 三角形單元程序 第7頁/共21頁 6:19 三角形單元程序 1 2 iiii Nab xc y A , i x N , i y N TTe A hdxdyA KB DBB DB 第8頁/共21頁 6:19 三角形單元程序 fdK fKd 1 第9頁/共21頁 6:19 三角形單元程序 e e e e UdKdT 2 1 第10頁/共21頁 6:19 三角形單元程序 ( , )( , ) e x yx yBd ( , )

3、( , )x yx yD 1 2 ijij Ud 第11頁/共21頁 6:19 三角形單元程序 第12頁/共21頁 6:19 三角形單元程序 第13頁/共21頁 6:19 材料矩陣計算程序 %- % Purpose: % determine the constitutive equation for isotropic material % Variable Description: % elastic - elastic modulus % poisson - Poissons ratio % iopt=1 - plane stress analysis % iopt=2 - plane s

4、train analysis % iopt=3 - axisymmetric analysis % iopt=4 - three dimensional analysis %- 2 1 00 01 01 1 2 E D 第14頁/共21頁 6:19 材料矩陣計算程序 )1 (2 21 00 01 1 0 1 1 )21)(1 ( )1 ( E D matmtrx= elastic/(1+poisson)*(1-2*poisson)* . (1-poisson) poisson 0; poisson (1-poisson) 0; 0 0 (1-2*poisson)/2; matmtrx= ela

5、stic/(1+poisson)*(1-2*poisson)* . (1-poisson) poisson poisson 0; poisson (1-poisson) poisson 0; poisson poisson (1-poisson) 0; 0 0 0 (1-2*poisson)/2; 第15頁/共21頁 6:19 材料矩陣計算程序 G G G G G G 0 00 0002 0002 0002 稱 對 D )21 (1 , )2(1 E G E matmtrx= elastic/(1+poisson)*(1-2*poisson)* . (1-poisson) poisson po

6、isson 0 0 0; poisson (1-poisson) poisson 0 0 0; poisson poisson (1-poisson) 0 0 0; 0 0 0 (1-2*poisson)/2 0 0; 0 0 0 0 (1-2*poisson)/2 0; 0 0 0 0 0 (1-2*poisson)/2; 第16頁/共21頁 6:19 局部單元系統(tǒng)自由度計算程序 index=feeldof(nd,nnel,ndof) %- % Purpose: % Compute system dofs associated with each element % Variable Des

7、cription: % index - system dof vector associated with element iel“ % nnel - number of nodes per element % ndof - number of dofs per node %- k=0; i=1:nnel start = (nd(i)-1)*ndof; j=1:ndof k=k+1; index(k)=start+j; 第17頁/共21頁 6:19 B矩陣計算程序 kinmtx2=fekine2d(nnel,dhdx,dhdy) %- % Purpose: % determine the ki

8、nematic equation between strains and displacements % for two-dimensional solids % Variable Description: % nnel - number of nodes per element % dhdx - derivatives of shape functions with respect to x % dhdy - derivatives of shape functions with respect to y %- i=1:nnel i1=(i-1)*2+1; i2=i1+1; kinmtx2(

9、1,i1)=dhdx(i); kinmtx2(2,i2)=dhdy(i); kinmtx2(3,i1)=dhdy(i); kinmtx2(3,i2)=dhdx(i); 第18頁/共21頁 6:19 剛度矩陣組裝程序 kk=feasmb_2(kk,k,index) %- % Purpose: % Assembly of element matrices into the system matrix % Variable Description: % kk - system matrix % k - element matri % index - d.o.f. vector associated

10、with an element %- edof = length(index); i=1:edof ii=index(i); j=1:edof jj=index(j); kk(ii,jj)=kk(ii,jj)+k(i,j); 第19頁/共21頁 6:19 邊界條件施加程序（置一法） kk,ff=feaplyc2(kk,ff,bcdof,bcval) %- % Purpose: % Apply constraints to matrix equation kkx=ff % Variable Description: % kk - system matrix before applying constraints % ff - system vector before applying constraints % bcdof - a vector containging constrained d.o.f % bcval - a vector containing contained value %- For example, there are constraints at d.o.f=2 and 10 and their constrained values are 0.0 and 2.5, respectively. Then, bc