《matlab网格划分程序与matlab有限元的结合.docx》由会员分享，可在线阅读，更多相关《matlab网格划分程序与matlab有限元的结合.docx（16页珍藏版）》请在三一办公上搜索。
1、matlab网格划分程序与matlab有限元的结合1.distmesh是一个较好的网格划分程序，具体可以参考：2. matlab有限元可以参考徐荣桥的书3. 这里本人打算画一个园内包含一个椭圆的模型：具体程序如下：a.网格划分：fh = (p) 0.05+0.3*dellipse(p,0.5,0.2); fd = (p) ddiff(dcircle(p,0,0,1),dellipse(p,0.5,0.2);p,t=distmesh2d(fd,fh,0.05,-1,-1;1,1,-1,-1;-1,1;1,-1;1,1);b.在 distmesh2d.m最后插入语句，导出据格式（节点坐标，4,单元
2、材料属性，边界约束条件）fid = fopen(exam4_2.dat3, w);ilength,jlength二size( p );fprintf(fid,%dn,ilength);for i = 1:1:ilengthfprintf(fid,%dt%ft%fn,i,p(i,1),p(i,2);endilength,jlength二size( t );fprintf(fid,%dn,ilength);for i = 1:1:ilengthfprintf(fid,%dt%dt%dt%dt%dn,i,t(i,1),t(i,2),t(i,3),1);endfclose(fid);文件exam4_2
3、.dat3内容如下：5105(节点个数)1 -0.707106 -0.707106 (每个节点坐标)2 -0.707106 0.707107.9869(三角单元个数)1 5006 4951 4934 1(号, i，j，k，材料属性的编号)2 197 147 100 13 413 2 366 14413 2 323 15（5种材料属性）1 2.00e90.25 100.0 22.0e3蕾号，杨氏模量，泊松比，厚度平面成力填1，密度)2 2.60e90.201.023.0e33 2.85e100.201.025.0e34 1.85e100.201.023.0e35 2.85e100.201.022
4、.0e32（约束个R）1 2 0.0 （节点，方向【2,为y方向】，00（0为固定）1 1 0.04 2 0.0 （节点，方向【2,为y方向】，0.0（0为固定）41 0.0c.调用有限元程序exam_2.m（见徐荣桥书）exam4_2 exam4_2.dat3由后处理 exam4_2_post最大主应力云图如下：-0.02MPa-0.02MPa-0.01 MPa-0.01 MPa-0.01 MPa-O.OOMPa-O.OOMPa-O.OOMPc-O.OOMPc最大剪应力云图如下:-0.01 MPa-0.01 MPa-0.01 MPa-0.01 MPa-0.01 MPa-0.01 MPa-O.
5、OOMPa-OBOMPa-O.OOMPa最小主应力云图如下:-O.OOMPa-O.OOMPc-O.OOMPc-0.01 MPc-0.01 MPc-0.01 MPc-0.02MPc-0.02MPc-0.02MPcexam_2。m源程序如下：function exam4_2( file_in )%本程序为第四章的第二个算例，采用三角形单元计算隧道在自 重作用下的变形和应力% exam4_2( filename )%输入参数：% file_in 有限元模型文件%定义全局变量%gNode节点坐标%gElement单元定义%gMaterial材料性质%gBC1第一类约束条件%gK整体刚度矩阵%gDelt
6、a整体节点坐标%gNodeStress节点应力% gElementStress单元应力global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gNFif nargin 1file_in = exam4_2.dat;end%检查文件是否存在if exist( file_in ) = 0disp( sprintf(错误：文件 %s 不存在,file_in )disp( sprintf(程序终止)return ;end%根据输入文件名称生成输出文件名称 path_str,name_str,ext_str =
7、fileparts( file_in );ext_str_out = .mat;file_out = fullfile( path_str, name_str, ext_str_out);%检查输出文件是否存在if exist( file_out ) 二 0answer = input( sprintf(文件 %s 已经存在，是否覆盖？ ( Yes / No ): , file_out ), s);if length( answer ) = 0answer = no;endanswer = lower( answer );if answer = y | answer = yesdisp( sp
8、rintf(请使用另外的文件名，或备份已有的文件);disp( sprintf(程序终止);return ;endend%建立有限元模型并求解，保存结果FemModel( file_in ) ;%定义有限元模型SolveModel ;%求解有限元模型SaveResults( file_out ) ;% 保存计算结果%计算结束disp( sprintf(计算正常结束，结果保存在文件%s中, file_out ) ) ;disp( sprintf(可以使用后处理程序exam4_2_post.m显示计算 结果 ) ) ;return ;function FemModel(filename)% 定义有
9、限元模型%输入参数：%filename -有限元模型文件%返回值：% 无%说明：%该函数定义平面问题的有限元模型数据：%gNode节点定义%gElement 单元定义%gMaterial -材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩% gBC1约束条件global gNode gElement gMaterial gBC1 gNF%打开文件fid = fopen( filename, r);%读取节点坐标node_number = fscanf( fid, %d, 1 );gNode = zeros( node_number, 2 );for i = 1:node_numberdumm
10、y = fscanf( fid, %d, 1 );gNode( i,: ) = fscanf( fid, %f, 1, 2);end%读取单元定义element_number = fscanf( fid, %d, 1 );gElement = zeros( element_number, 4 );for i = 1:element_numberdummy = fscanf( fid, %d, 1 );gElement( i, : ) = fscanf( fid, %d, 1, 4);end%读取材料信息material_number = fscanf( fid, %d, 1);gMateria
11、l = zeros( material_number, 4 );for i = 1:material_numberdummy = fscanf( fid, %d, 1 );gMaterial( i,: ) = fscanf( fid, %f, 1,4);end%读取边界条件bc1_number = fscanf( fid, %d, 1 );gBC1 = zeros( bc1_number, 3 );for i = 1:bc1_numbergBC1( i, 1 ) = fscanf( fid, %d, 1 );gBC1( i, 2 ) = fscanf( fid, %d, 1 );gBC1( i
12、, 3 ) = fscanf( fid, %f, 1 );end%关闭文件fclose( fid );%集中力%节点号自由度号集中力值gNF=1,1,-800e3;1,2,-800e3;4,1,-800e3;4,2,-800e3; ;returnfunction SolveModel%求解有限元模型%输入参数：% 无%返回值：%无%说明：%该函数求解有限元模型，过程如下%1.计算单元刚度矩阵，集成整体刚度矩阵%2.计算单元的等效节点力，集成整体节点力向量%3.处理约束条件，修改整体刚度矩阵和节点力向量%4.求解方程组，得到整体节点位移向量%5.计算单元应力和节点应力global gNode g
13、Element gMaterial gBC1 gK gDelta gNodeStress gElementStress gNF% step1.定义整体刚度矩阵和节点力向量node_number,dummy = size( gNode );gK = sparse( node_number * 2, node_number * 2 );f = sparse( node_number * 2, 1 );% step2.计算单元刚度矩阵，并集成到整体刚度矩阵中element_number,dummy = size( gElement );for ie=1:1:element_numberdisp( s
14、printf(计算刚度矩阵，当前单元:%d, ie );k = StiffnessMatrix( ie );AssembleStiffnessMatrix( ie, k );end% step3 .把集中力直接集成到整体节点力向量中nf_number, dummy = size( gNF );for inf=1:1:nf_numbern = gNF( inf, 1 );d = gNF( inf, 2 );f( (n-1)*2 + d ) = gNF( inf, 3 );end% step3.计算自重产生的等效节点力% for ie=1:1:element_number% disp( sprin
15、tf(计算自重的等效节点力，当前单元：d, ie ) ) ;% egf = EquivalentGravityForce( ie );%i = gElement( ie, 1 );%j = gElement( ie, 2 );% m = gElement( ie, 3 );% f( (i-1)*2 + 1 : (i-1)*2+2 ) = f( (i-1)*2 + 1 : (i-1)*2+2 ) + egf( 1:2 ) ;% f( (j-1)*2 + 1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + egf( 3:4 );% f( (m-1)*2 +
16、 1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + egf( 5:6 );% end% step4.处理约束条件，修改刚度矩阵和节点力向量。采用乘大 数法bc_number,dummy = size( gBC1 );for ibc=1:1:bc_numbern = gBC1(ibc, 1 );d = gBC1(ibc, 2 );m = (n-1)*2 + d ;f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;gK(m,m) = gK(m,m) * 1e15 ;end% step 5.求解方程组，得到节点位移向量 gDelta
17、 = gK f ;% step 6.计算单元应力gElementStress = zeros( element_number, 6 );for ie=1:element_numberdisp( sprintf(计算单元应力，当前单元:%d, ie );es = ElementStress( ie );gElementStress( ie, : ) = es ;end% step 7.计算节点应力(采用绕节点加权平均) gNodeStress = zeros( node_number, 6 );for i = 1:node_numberdisp( sprintf(计算节点应力，当前结点:%d,
18、i );S = zeros( 1, 3 );A = 0 ;for ie=1:1:element_numberfor k=1:1:3if i = gElement( ie, k )area= ElementArea( ie );S = S + gElementStress(ie,1:3 ) * area ;A = A + area ;break ;endendendgNodeStress(i,1:3) = S / A ;gNodeStress(i,6) =0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2)人2 + 4*gNodeStress(i,3)人2 )
19、;gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2) + gNodeStress(i,6);gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)- gNodeStress(i,6);endreturnfunction k = StiffnessMatrix( ie )%计算单元刚度矩阵%输入参数：% ie -单元号%返回值：% k -单元刚度矩阵global gNode gElement gMaterialk = zeros( 6, 6 );E = gMaterial(
20、 gElement(ie, 4), 1 );mu = gMaterial( gElement(ie, 4), 2 );h = gMaterial( gElement(ie, 4), 3 );xi = gNode( gElement( ie, 1 ), 1 );yi = gNode( gElement( ie, 1 ), 2 );Xj = gNode( gElement( ie, 2 ), 1 );yj = gNode( gElement( ie, 2 ), 2 );xm = gNode( gElement( ie, 3 ), 1 );ym = gNode( gElement( ie, 3 ),
21、 2 );ai = xj*ym - xm*yj ;aj = xm*yi - xi*ym ;am = xi*yj - xj*yi ;bi = yj - ym ;bj = ym - yi ;bm = yi - yj ;ci = -(xj-xm);cj = -(xm-xi);cm = -(xi-xj);area = (ai+aj+am)/2 ;B = bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm; B = B/2/area ;D = 1-mu mu0mu 1-mu 000 (1-2*mu)/2;D = D*E/(1-2*mu)/(1 + mu);k
22、= transpose(B)*D*B*h*abs(area); returnfunction B = MatrixB( ie )%计算单元的应变矩阵B%输入参数：% ie -单元号%返回值：% B -单元应变矩阵 global gNode gElementxi = gNode( gElement( ie, 1 ), 1 ); yi = gNode( gElement( ie, 1 ), 2 ); xj = gNode( gElement( ie, 2 ), 1 ); yj = gNode( gElement( ie, 2 ), 2 ); xm = gNode( gElement( ie, 3
23、), 1 ); ym = gNode( gElement( ie, 3 ), 2 ); ai = xj*ym - xm*yj ;aj = xm*yi - xi*ym ;am = xi*yj - xj*yi ;bi = yj - ym ;bj = ym - yi ;bm = yi - yj ;ci = -(xj-xm);cj = -(xm-xi);cm = -(xi-xj);area = (ai+aj+am)/2 ;B = bi 0 bj 0 bm 00 ci 0 cj 0 cm ci bi cj bj cm bm; B = B/2/area ;returnfunction area = Ele
24、mentArea( ie )%计算单元面积%输入参数：% ie -单元号%返回值：% area 单兀面积global gNode gElement gMaterial xi = gNode( gElement( ie, 1 ), 1 ); yi = gNode( gElement( ie, 1 ), 2 ); xj = gNode( gElement( ie, 2 ), 1 ); yj = gNode( gElement( ie, 2 ), 2 ); xm = gNode( gElement( ie, 3 ), 1 ); ym = gNode( gElement( ie, 3 ), 2 );
25、ai = xj*ym - xm*yj ;aj = xm*yi - xi*ym ;am = xi*yj - xj*yi ;area = abs(ai+aj+am)/2);returnfunction AssembleStiffnessMatrix( ie, k )%把单元刚度矩阵集成到整体刚度矩阵%输入参数：%ie -单元号%k -单元刚度矩阵%返回值：% 无global gElement gKfor i = 1:1:3for j = 1:1:3for p=1:1:2for q = 1:1:2m = (i-1)*2 + p ;n = (j-1)*2+q ；M = (gElement(ie,i)-
26、1)*2+p ;N = (gElement(ie,j)-1)*2+q ;gK(M,N) = gK(M,N) + k(m,n);endendendendreturnfunction egf = EquivalentGravityForce( ie )%计算单元自重的等效节点力%输入参数% ie -单元号% 返回值% egf -自重的等效节点力global gElement gMaterialarea = ElementArea( ie );h = gMaterial( gElement( ie, 4 ), 3 );ro = gMaterial( gElement( ie, 4 ), 4 );w
27、= area * h * ro ;egf = -w/3 * 0; 1; 0; 1; 0; 1;returnfunction es = ElementStress( ie )%计算单元的应力分量%输入参数% ie -单元号%返回值% es单元应力分量列阵(1x6): sx, sy, txy, s1, s2,tmaxglobal gElement gDelta gMateriales = zeros( 1, 6 ) ; %单元的应力分量de = zeros( 6, 1 ) ; %单元节点位移列阵E = gMaterial( gElement(ie, 4), 1 );mu = gMaterial(
28、gElement(ie, 4), 2 );D = 1-mu mu0mu 1-mu000 (1-2*mu)/2;D = D*E/(1-2*mu)/(1 + mu);B = MatrixB( ie );for j = 1:1:3de( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 );de( 2*j ) = gDelta( 2*gElement( ie, j );endes(1:3) = D * B * de ;es(6) = 0.5*sqrt(es(1)-es(2)A2 + 4*es(3)人2 );es(4) = 0.5*(es(1)+es(2) + es(6);es(5) = 0.5*(es(1)+es(2) - es(6);returnfunction SaveResults( file_out )%显示计算结果%输入参数：% file_out -保存结果文件%返回值：% 无global gNode gElement gMaterial gBC1 gDelta gNodeStress gElementStresssave( file_out, gNode, gElement, gMaterial, gBC1, gDelta, gNodeStress, gElementStress);return
