《MATLAB 有限元计算程.docx》由会员分享，可在线阅读，更多相关《MATLAB 有限元计算程.docx（22页珍藏版）》请在三一办公上搜索。
1、MATLAB 有限元计算程function exam7_1% 本程序为第七章的第一个算例，采用薄板矩形单元计算矩形板的弯曲% 输入参数： % 无% 定义板的尺寸a = 1 ; % 板x方向宽度b = 1 ; % 板y方向长度h = 0.01 ; % 板z方向厚度% 定义网格密度 m = 16 ; % x方向单元数目 n = 16 ; % y方向单元数目% 定义材料性质E = 2.1e11 ; % 弹性模量mu = 0.3 ; % 泊松比% 荷载参数load = 均布力 ; % 也可以是 集中力%load = 集中力 ;p = 1e6 ; % 载荷大小% 边界条件%bc = 简支 ; % 也可以
2、是 固支 bc = 固支 ;% 输出文件名称file_out = exam7_1.mat ;% 检查输出文件是否存在file_prompt = 0 ;if exist( file_out ) = 0 & file_prompt = 1answer = input( sprintf( 文件 %s 已经存在，是否覆盖? ( Yes / No ): , file_out ), s ) ;if length( answer ) = 0answer = no ;endanswer = lower( answer ) ;if answer = y | answer = yes disp( sprintf(
3、 请使用另外的文件名，或备份已有的文件 ) ) ;disp( sprintf( 程序终止 ) );return ; endendFemModel(a,b,h,m,n,E,mu,load,p,bc) ; % 定义有限元模型SolveModel ; % 求解有限元模型SaveResults( file_out ) ; % 保存计算结果DisplayResults(a,b,h,m,n,E,mu,load,p,bc) ; % 把计算结果显示出来 returnfunction FemModel(a,b,h,m,n,E,mu,load,p,bc)% 定义有限元模型% 输入参数：% a - 板x方向宽度%
4、b - 板y方向长度% h - 板z方向厚度% m - x方向单元数目 % n - y方向单元数目% E - 弹性模量% mu - 泊松比% load - 载荷类型% p - 载荷大小% bc - 边界条件% 返回值：% 无global gNode gElement gMaterial gBC1 gDF gNF gK% 计算结点坐标dx = a / m / 2 ; % 取板的1/4进行分析dy = b / n / 2 ;gNode = zeros( (m+1)*(n+1), 2 ) ;for i=1:n+1for j=1:m+1gNode( (i-1)*(m+1)+j, : ) = dx*(j
5、-1),dy*(i-1) ;endend% 定义单元gElement = zeros( n*n, 5 ) ;for i=1:nfor j=1:mgElement( (i-1)*m+j, 1:4) = (i-1)*(m+1)+j, .(i-1)*(m+1)+j+1, .i*(m+1)+j+1,.i*(m+1)+j ;endendgElement( :, 5 ) = 1 ;% 定义材料gMaterial = E, mu, h ; % 确定边界条件if bc = 简支gBC1 = zeros( (n+1)*2+m*2+1, 3 ) ;for i=1:m+1gBC1( i, : ) = i, 2, 0
6、 ; % x=0的边界上绕x轴的转角等于零endfor i=1:n+1gBC1( (m+1)+i, : ) = (i-1)*(m+1)+1, 3, 0 ; % y=0的边界上绕y轴的转角等于零endfor i=1:n+1gBC1( m+n+2+i, : ) = i*(m+1), 1, 0 ; % x=a/2的边界上挠度为零endfor i=1:mgBC1( (n+1)*2+m+1+i, : ) = n*(m+1)+i, 1, 0 ; % y=b/2的边界上挠度为零endelseif bc = 固支gBC1 = zeros( m*4+(n+1)*3+n, 3 ) ;for i=1:mgBC1(
7、i, : ) = i, 2, 0 ; % x=0的边界上绕x轴的转角等于零endfor i=1:ngBC1( m+i, : ) = (i-1)*(m+1)+1, 3, 0 ; % y=0的边界上绕y轴的转角等于零endfor i=1:n+1gBC1( m+n+i, : ) = i*(m+1), 1, 0 ; % x=a/2的边界上挠度为零endfor i=1:n+1gBC1( m+n+n+1+i, : ) = i*(m+1), 2, 0 ; % x=a/2的边界上绕x轴的转角为零endfor i=1:n+1gBC1( m+n+(n+1)*2+i, : ) = i*(m+1), 3, 0 ; %
8、 x=a/2的边界上绕y轴的转角为零endfor i=1:mgBC1( m+n+(n+1)*3+i, : ) = n*(m+1)+i, 1, 0 ; % y=b/2的边界上挠度为零endfor i=1:mgBC1( m*2+n+(n+1)*3+i, : ) = n*(m+1)+i, 2, 0 ; % y=b/2的边界上绕x轴的转角为零endfor i=1:mgBC1( m*3+n+(n+1)*3+i, : ) = n*(m+1)+i, 3, 0 ; % y=b/2的边界上绕y轴的转角为零endend% 确定载荷if load = 均布力gNF = ;gDF = gElement ;gDF(:,
9、5) = p ;elseif load = 集中力gDF = ;gNF = 1, 1, p ;endreturnfunction SolveModel% 求解有限元模型% 输入参数：% 无% 返回值：% 无% 说明：% 该函数求解有限元模型，过程如下% 1. 计算单元刚度矩阵，集成整体刚度矩阵% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量global gNode gElement gMaterial gBC1 gDF gNF gK gDelta gElementMoment gNodeMoment
10、% step1. 定义整体刚度矩阵和节点力向量node_number,dummy = size( gNode ) ;gK = sparse( node_number * 3, node_number * 3 ) ;f = sparse( node_number * 3, 1) ; % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中element_number,dummy = size( gElement ) ;for ie=1:1:element_numberdisp( sprintf( 计算单元刚度矩阵并集成，当前单元: %d, ie ) ) ;k = StiffnessMatrix
11、( ie ) ;AssembleStiffnessMatrix( ie, k ) ;end% step3. 计算分布压力的等效节点力，并集成到整体节点力向量中df_number, dummy = size( gDF ) ;for idf = 1:1:df_numberedf = EquivalentDistPressure( gDF(idf,1:4), gDF(idf,5) ) ;node = gDF( idf, 1:4 ) ;f( node*3-2 ) = f( node*3-2 ) + edf( 1:3:10 ) ;f( node*3-1 ) = f( node*3-1 ) + edf(
12、2:3:11 ) ;f( node*3 ) = f( node*3 ) + edf( 3:3:12 ) ;end% step4. 计算分布压力的等效节点力，并集成到整体节点力向量中nf_number,dummy = size( gNF ) ;for inf = 1:nf_numbernode = gNF( inf, 1 ) ;dim = gNF( inf, 2 ) ;f( (node-1)*3+dim ) = f( (node-1)*3+dim ) + gNF( inf, 3 ) ;end% step5. 处理第一类约束条件，修改刚度矩阵和节点力向量。采用乘大数法bc1_number,dumm
13、y = size( gBC1 ) ;for ibc=1:1:bc1_numbern = gBC1(ibc, 1 ) ;d = gBC1(ibc, 2 ) ;m = (n-1)*3 + d ;f(m) = gBC1(ibc, 3)* gK(m,m) * 1e10 ;gK(m,m) = gK(m,m) * 1e10 ;end% step5. 求解方程组，得到节点位移向量gDelta = gK f ;% step6. 计算单元弯矩gElementMoment = zeros( element_number, 4, 3 ) ;for ie=1:element_number disp( sprintf(
14、 计算单元弯矩，当前单元: %d, ie ) ) ;em = ElementMoment( ie ) ;gElementMoment( ie, :, : ) = em ;end% step7. 计算节点弯矩(采用绕节点平均)gNodeMoment = zeros( node_number, 3 ) ; nodenum = zeros( node_number,1 ) ;for i=1:node_number disp( sprintf( 计算节点弯矩，当前结点: %d, i ) ) ;for ie=1:element_numberindex = find( gElement( ie, 1:4
15、) = i ) ;if isempty(index)nodenum(i) = nodenum(i) + 1 ;gNodeMoment(i,1) = gNodeMoment(i,1) + gElementMoment( ie, index, 1 ) ;gNodeMoment(i,2) = gNodeMoment(i,2) + gElementMoment( ie, index, 2 ) ;gNodeMoment(i,3) = gNodeMoment(i,3) + gElementMoment( ie, index, 3 ) ;endendendgNodeMoment(:,1) = gNodeMo
16、ment(:,1) ./ nodenum ;gNodeMoment(:,2) = gNodeMoment(:,2) ./ nodenum ;gNodeMoment(:,3) = gNodeMoment(:,3) ./ nodenum ;returnfunction k = StiffnessMatrix( ie ) % 计算平面应变等参数单元的刚度矩阵% 输入参数：% ie - 单元号% 返回值：% k - 单元刚度矩阵global gNode gElement gMaterialk = zeros( 12, 12 ) ;E = gMaterial( gElement( ie, 5 ), 1
17、) ;mu = gMaterial( gElement( ie, 5 ), 2 ) ;h = gMaterial( gElement( ie, 5 ), 3 ) ;D = E * h3 / 12 / ( 1 - mu2 ) ;x = gNode( gElement( ie, : ), 1 ) ;y = gNode( gElement( ie, : ), 2 ) ;a = abs( x(3) - x(1) ) / 2 ;b = abs( y(3) - y(1) ) / 2 ;ab2 = a2/b2 ;ba2 = b2/a2 ;H = D/60/a/b ;xi = -1, 1, 1, -1 ;et
18、a = -1, -1, 1, 1 ;for i=1:4for j=1:4xi0 = xi(i)*xi(j) ;eta0 = eta(i)*eta(j) ;a11 = 3*H*(15*(ba2*xi0+ab2*eta0) .+(14-4*mu+5*ba2+5*ab2)*xi0*eta0) ;a12 = -3*H*b*(2+3*mu+5*ab2)*xi0*eta(i) .+ 15*ab2*eta(i)+5*mu*xi0*eta(j) ;a13 = 3*H*a*(2+3*mu+5*ba2)*xi(i)*eta0 .+ 15*ba2*xi(i)+5*mu*xi(j)*eta0) ;a21 = -3*H
19、*b*(2+3*mu+5*ab2)*xi0*eta(j) .+ 15*ab2*eta(j)+5*mu*xi0*eta(i) ;a22 = H*b2*(2*(1-mu)*xi0*(3+5*eta0)+.+ 5*ab2*(3+xi0)*(3+eta0) ;a23 = -15*H*mu*a*b*(xi(i)+xi(j)*(eta(i)+eta(j) ;a31 = 3*H*a*(2+3*mu+5*ba2)*xi(j)*eta0 .+ 15*ba2*xi(j)+5*mu*xi(i)*eta0) ;a32 = -15*H*mu*a*b*(xi(i)+xi(j)*(eta(i)+eta(j) ;a33 =
20、H*a2*(2*(1-mu)*eta0*(3+5*xi0)+.+ 5*ba2*(3+xi0)*(3+eta0) ;k( i*3-2:i*3, j*3-2:j*3 ) = a11 a12 a13;a21 a22 a23;a31 a32 a33 ;endendreturnfunction AssembleStiffnessMatrix( ie, k )% 把单元刚度矩阵集成到整体刚度矩阵% 输入参数:% ie - 单元号% k - 单元刚度矩阵% 返回值:% 无global gElement gKfor i=1:1:4for j=1:1:4for p=1:1:3for q=1:1:3m = (i-
21、1)*3+p ;n = (j-1)*3+q ;M = (gElement(ie,i)-1)*3+p ;N = (gElement(ie,j)-1)*3+q ;gK(M,N) = gK(M,N) + k(m,n) ;endendendendreturnfunction edf = EquivalentDistPressure( node, press )% 计算分布压力的等效节点力% 输入参数:% node - 结点号% press - 跟结点号对应的压力值% 返回值:% edf - 等效节点力向量global gNodeedf = zeros( 12, 1 ) ;x = gNode( node
22、, 1 ) ;y = gNode( node, 2 ) ;a = abs( x(3)-x(1) ) / 2 ;b = abs( y(3)-y(1) ) / 2 ;xi = -1; 1; 1; -1 ;eta = -1;-1; 1; 1 ;for i=1:4edf(i*3-2:i*3) = press*a*b; -press*a*b2/3*eta(i); press*a2*b/3*xi(i) ;endreturnfunction em = ElementMoment( ie )% 计算单元的结点弯矩% 输入参数:% ie - 单元号% 返回值:% em - 单元结点弯矩global gEleme
23、nt gDelta gMaterialE = gMaterial( gElement( ie, 5 ), 1 ) ;mu = gMaterial( gElement( ie, 5 ), 2 ) ;h = gMaterial( gElement( ie, 5 ), 3 ) ;em = zeros( 4, 3 ) ;node = gElement( ie, 1:4 ) ;delta = zeros( 12, 1 ) ;delta(1:3:10) = gDelta( (node-1)*3+1 ) ;delta(2:3:11) = gDelta( (node-1)*3+2 ) ;delta(3:3:1
24、2) = gDelta( (node-1)*3+3 ) ;xi = -1, 1, 1, -1 ;eta = -1, -1, 1, 1 ;D = E/(1-mu2)*1 mu 0mu 1 00 0 (1-mu)/2 ;for i=1:4B = MatrixB( ie, xi(i), eta(i) ) ;em( i, : ) = transpose(h3/12*D*B*delta) ;endreturnfunction B = MatrixB( ie, x, e )% 计算单元的应变矩阵B% 输入参数:% ie - 单元号% 返回值:% B - 单元应变矩阵global gNode gElemen
25、tB = zeros( 3, 12 ) ;xy = gNode( gElement(ie, 1:4), : ) ;a = abs(xy(3,1)-xy(1,1)/2 ;b = abs(xy(3,2)-xy(1,2)/2 ;xi = -1, 1, 1, -1 ;eta = -1, -1, 1, 1 ;for i=1:4x0 = xi(i)*x ;e0 = eta(i)*e ;B(:,(i-1)*3+1:(i-1)*3+3) = 3*b/a*x0*(1+e0) 0 b*xi(i)*(1+3*x0)*(1+e0)3*a/b*e0*(1+x0) -a*eta(i)*(1+x0)*(1+3*e0) 0x
26、i(i)*eta(i)*(3*x2+3*e2-4) -b*xi(i)*(3*e2+2*e0-1) a*eta(i)*(3*x2+2*x0-1) ;endB = B / a / b / 4 ;returnfunction SaveResults( file_out )% 保存计算结果% 输入参数：% 无% 返回值：% 无global gNode gElement gMaterial gBC1 gDF gDelta gNF gKsave( file_out, gNode, gElement, gMaterial, gBC1, .gDF, gDelta, gNF, gK ) ;returnfunct
27、ion DisplayResults(a,b,h,m,n,E,mu,load,p,bc)% 显示计算结果，并与解析解结果比较% 输入参数：% a - 板x方向宽度% b - 板y方向长度% h - 板z方向厚度% m - x方向单元数目 % n - y方向单元数目% E - 弹性模量% mu - 泊松比% load - 载荷类型% p - 载荷大小% bc - 边界条件% 返回值：% 无global gDelta gNodeMomentwmax = gDelta( 1, 1 ) ; % 中心挠度D = E*h3/12/(1-mu2) ; % 抗弯刚度mx = gNodeMoment( 1, 1
28、 ) ; % 中心内力矩(Mx)my = gNodeMoment( 1, 2 ) ; % 中心内力矩(My)mxy = gNodeMoment( 1, 3 ) ; % 中心内力矩(Mxy)fprintf( 板的尺寸abh = %g%g%gn, a, b, h ) ;fprintf( 弹性模量 = %gn, E ) ;fprintf( 泊松比 = %gn, mu ) ;fprintf( 载荷类型 = %sn, load ) ;fprintf( 载荷大小 = %gn, p ) ;fprintf( 边界条件 = %sn, bc ) ;fprintf( 网格密度 mn = %d%dn, m, n )
29、;if load = 均布力alpha = wmax * D / p / a4 ;beta = mx/p/a2 ;beta1 = my/p/a2 ;elseif load = 集中力alpha = wmax * D / p / 4 / a2 ;beta = mx/p/4 ;beta1 = my/p/3 ;endfprintf( 挠度系数 alpha = %.7fn, alpha ) ;fprintf( 弯矩系数 beta = %.7fn, beta ) ;fprintf( 弯矩系数 beta1 = %.7fn, beta1 ) ;% fid = fopen( exam7_1.out, a ) ;% fprintf( fid, %.1f %15.7f %15.7fn, b/a, alpha, beta ) ;% fclose( fid ) ;return
