结构分析的有限元法与MATLAB程序设计，徐荣桥编源程序

function exam6_1 % 本程序为第六章的第一个算例，采用4～8结点四边形等参元计算受内压的旋转厚壁圆筒 % 输入参数： % 无 % 输出文件名称 file_out = 'exam6_1.mat' ; % 检查输出文件是否存在 if exist( file_out ) ~= 0 answer = input( sprintf( '文件 %s 已经存在，是否覆盖? ( Yes / [No] ): ', file_out ), 's' ) ; if length( answer ) == 0 answer = 'no' ; end answer = lower( answer ) ; if answer ~= 'y' | answer ~= 'yes' disp( sprintf( '请使用另外的文件名，或备份已有的文件' ) ) ; disp( sprintf( '程序终止' ) ); return ; end end FemModel ; % 定义有限元模型 SolveModel ; % 求解有限元模型 SaveResults( file_out ) ; % 保存计算结果 DisplayResults ; % 把计算结果显示出来 return function FemModel % 定义有限元模型 % 输入参数： % 无 % 返回值： % 无 % 全局变量及名称 % gNode ------ 节点坐标 % gElement --- 单元定义 % gMaterial -- 材料性质 % gBC1 ------- 第一类约束条件 % gBC3 ------- 斜约束 % gDF -------- 分布力 global gNode gElement gMaterial gBC1 gBC3 gDF r1 = 0.05 ; % 内半径 r2 = 0.1 ; % 外半径 theta = 10.0*pi/180.0 ; % 圆心角（弧度） % 节点坐标 node_number = 28 ; gNode = zeros( node_number, 2 ) ; j = 0 ; for i = 1:5:26 % 结点1,6,11,16,21,26的坐标 r = r1 + (r2-r1)/5 * j ; gNode(i,1) = r*cos(theta) ; gNode(i,2) = r*sin(theta) ; j = j + 1 ; end j = 0 ; for i = 4:5:24 % 结点4,9,14,19,24的坐标 r = r1+(r2-r1)/10+(r2-r1)/5*j ; gNode(i,1) = r*cos(theta) ; gNode(i,2) = r*sin(theta) ; j = j+1 ; end j = 0 ; for i=2:5:27 % 结点2,7,12,17,22,27的坐标 r = r1+(r2-r1)/5*j ; gNode(i,1) = r*cos(theta/2) ; gNode(i,2) = r*sin(theta/2) ; j = j+1 ; end j = 0 ; for i=3:5:28 % 结点3,8,13,18,23,28的坐标 r = r1+(r2-r1)/5*j ; gNode(i,1) = r ; gNode(i,2) = 0 ; j = j+1 ; end j = 0 ; for i=5:5:25 % 结点5,10,15,20,25的坐标 r = r1+(r2-r1)/10+(r2-r1)/5*j ; gNode(i,1) = r ; gNode(i,2) = 0 ; j = j+1 ; end % 单元定义 element_number = 5 ; gElement = zeros( element_number, 9 ) ; gElement(1,:)=[3 8 6 1 5 7 4 2 1] ; for i=2:5 gElement(i,1:8) = gElement(i-1,1:8)+5 ; gElement(i,9) = 1 ; end % 材料性质 gMaterial = [2.0e11 0.3 7800] ; % 第一类约束条件 bc1_number = 11; gBC1 = zeros( bc1_number, 3 ) ; j = 0 ; for i=3:5:28 % 结点3,8,13,18,23,28的边界条件(y方向约束) j = j + 1 ; gBC1(j,1) = i ; % 结点号 gBC1(j,2) = 2 ; % 自由度号 gBC1(j,3) = 0 ; % 约束值 end j = 6 ; for i=5:5:25 % 结点5,10,15,20,25的边界条件(y方向约束) j = j + 1; gBC1(j,1) = i ; % 结点号 gBC1(j,2) = 2 ; % 自由度号 gBC1(j,3) = 0 ; % 约束值 end % 斜约束条件 bc3_number = 17 ; gBC3 = zeros( bc3_number, 4 ) ; j = 0 ; for i=1:5:26 % 结点1,6,11,16,21,26的边界条件(y'方向约束) j = j + 1 ; gBC3(j,1) = i ; % 结点号 gBC3(j,2) = 2 ; % 自由度号 gBC3(j,3) = 0 ; % 约束值 gBC3(j,4) = theta ; % 倾斜角度 end j = 6 ; for i=4:5:24 % 结点4,9,14,19,24的边界条件(y'方向约束) j = j + 1; gBC3(j,1) = i ; % 结点号 gBC3(j,2) = 2 ; % 自由度号 gBC3(j,3) = 0 ; % 约束值 gBC3(j,4) = theta ; % 倾斜角度 end j = 11 ; for i=2:5:27 j = j+1 ; gBC3(j,1) = i ; % 结点号 gBC3(j,2) = 2 ; % 自由度号 gBC3(j,3) = 0 ; % 约束值 gBC3(j,4) = theta/2 ; % 倾斜角度 end % 分布压力 df_number = 1 ; gDF = [3 1 2 120e6 120e6 120e6] ; return function SolveModel % 求解有限元模型 % 输入参数： % 无 % 返回值： % 无 % 说明： % 该函数求解有限元模型，过程如下 % 1. 计算单元刚度矩阵，集成整体刚度矩阵 % 2. 计算单元的等效节点力，集成整体节点力向量 % 3. 处理约束条件，修改整体刚度矩阵和节点力向量 % 4. 求解方程组，得到整体节点位移向量 global gNode gElement gMaterial gBC1 gBC3 gDF gK gDelta1 gDelta gElementStress % step1. 定义整体刚度矩阵和节点力向量 [node_number,dummy] = size( gNode ) ; gK = sparse( node_number * 2, node_number * 2 ) ; f = sparse( node_number * 2, 2 ) ; % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中 [element_number,dummy] = size( gElement ) ; for ie=1:1:element_number disp( sprintf( '计算单元刚度矩阵并集成，当前单元: %d', ie ) ) ; k = StiffnessMatrix( ie ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 计算单元离心力的等效结点力，并集成到整体结点力向量中 for ie=1:1:element_number evf = EquivalentVolumeForce( ie ) ; node = gElement( ie, 1:8 ) ; f( (node-1)*2+1,2 ) = f( (node-1)*2+1,2 ) + evf(1:2:15) ; f( (node-1)*2+2,2 ) = f( (node-1)*2+2,2 ) + evf(2:2:16) ; end % step4. 计算分布压力的等效节点力，并集成到整体节点力向量中 [df_number, dummy] = size( gDF ) ; for idf = 1:1:df_number enf = EquivalentDistPressure( gDF(idf,1:3), gDF(idf,4:6) ) ; i = gDF(idf, 1) ; j = gDF(idf, 2) ; m = gDF(idf, 3) ; f( (i-1)*2+1:(i-1)*2+2, 1 ) = f( (i-1)*2+1:(i-1)*2+2, 1 ) + enf( 1:2 ) ; f( (j-1)*2+1:(j-1)*2+2, 1 ) = f( (j-1)*2+1:(j-1)*2+2, 1 ) + enf( 3:4 ) ; f( (m-1)*2+1:(m-1)*2+2, 1 ) = f( (m-1)*2+1:(m-1)*2+2, 1 ) + enf( 5:6 ) ; end % step5. 处理斜约束边界条件 [bc3_number,dummy] = size(gBC3) ; for ibc=1:1:bc3_number n = gBC3( ibc, 1 ) ; theta = gBC3( ibc, 4 ) ; c = cos(theta) ; s = sin(theta) ; T = [ c s -s c ] ; gK((n-1)*2+1:(n-1)*2+2,:) = T*gK((n-1)*2+1:(n-1)*2+2,:) ; gK(:,(n-1)*2+1:(n-1)*2+2) = gK(:,(n-1)*2+1:(n-1)*2+2)*transpose(T) ; f((n-1)*2+1:(n-1)*2+2,:) = T*f((n-1)*2+1:(n-1)*2+2,:) ; gBC1 = [gBC1; gBC3(ibc,1:3)] ; end % step6. 处理第一类约束条件，修改刚度矩阵和节点力向量。采用乘大数法 [bc1_number,dummy] = size( gBC1 ) ; for ibc=1:1:bc1_number n = 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 % step7. 求解方程组，得到节点位移向量 gDelta1 = gK \ f ; % step8. 把有斜约束的结点位移转换到整体坐标 gDelta = gDelta1 ; for ibc=1:bc3_number n = gBC3( ibc, 1 ) ; theta = gBC3( ibc, 4 ) ; c = cos(theta) ; s = sin(theta) ; T = [ c s -s c ] ; gDelta( (n-1)*2+1:n*2,:) = transpose(T)*gDelta( (n-1)*2+1:n*2,: ) ; end % step9. 计算每个单元的结点应力 gElementStress = zeros( element_number, 8, 3, 2 ) ; delta = zeros( 16, 2 ) ; for ie = 1:element_number xi = [ -1 1 1 -1 0 1 0 -1 ] ; eta = [ -1 -1 1 1 -1 0 1 0 ] ; for n=1:8 B = MatrixB( ie, xi(n), eta(n) ) ; D = MatrixD( ie, 1 ) ; delta(1:2:15, :) = gDelta( gElement(ie,1:8)*2-1, : ) ; delta(2:2:16, :) = gDelta( gElement(ie,1:8)*2, : ) ; sigma = D*B*delta ; gElementStress( ie, n, :, : ) = sigma ; end end return function k = StiffnessMatrix( ie ) % 计算平面应变等参数单元的刚度矩阵 % 输入参数： % ie -- 单元号 % 返回值： % k -- 单元刚度矩阵 %