有限元方法与MATLAB程序设计 第5章 轴对称问题.pptxVIP

第5章 轴对称问题§5.1 位移模式和形函数§5.2 几何矩阵§5.3 应力及应力矩阵§5.4 单元刚度矩阵§5.5 等效结点力§5.6 程序设计§5.1 位移模式和形函数结点位移向量位移形函数矩阵三角形环单元§5.2 几何矩阵三角形面积三角形环单元形函数§5.2 几何矩阵几何矩阵记几何方程几何子矩阵§5.3 应力及应力矩阵§5.3 应力及应力矩阵弹性矩阵应力及应力矩阵近似取§5.4 单元刚度矩阵§5.5 等效结点力1 体积力自重需要坐标变换坐标变换§5.5 等效结点力1 体积力自重2 离心力等效结点力例5.1 厚壁圆筒受内压y221222324258m262830324m252729311.5161718192018202224171921234m1111213141510121416y9111315x0.5678910246813570122.5343.55x234§5.6 程序主程序function AxiaSymm轴对称问题[gxy,ndel,dofix,F,Em,mu,Th,nd,ne,U] = AxiSymExa; dofree = setdiff(1:2*nd,dofix);非位移约束自由度a1 = mu/(1-mu); a2 = (1-2*mu)/(1-mu)/2;D = Em/a2/2*[1,a1,a1,0; a1,1,a1,0; a1,a1,1,0; 0,0,0,a2];弹性矩阵K = zeros(2*nd,2*nd); for el = 1:neN(2:2:6) = 2*ndel(el,:); N(1:2:5) = N(2:2:6)-1;单元自由度xy = gxy(ndel(el,:),:); [B,A] = AxiaSymmStrain(xy);几何矩阵和单元面积%%15K(N,N) = K(N,N)+2*pi*B*D*B*Th*sum(xy(:,1))*A;%§5.6 程序主程序fmuintf(%4s%6s%8s%10s%12s\n,Node,r,z,ur,w)标题for j = 1:nd输出结点号，结点坐标，结点位移 fmuintf(%4i%8.2f%8.2f%12.3e%12.3e\n,j,gxy(j,:),U(2*j+(-1:0))) end%fmuintf(%4s%4s%4s%4s%8s%12s%12s%12s\n,Elem,i,j,k,Sr,Sc,Sz,Sr标题for el = 1:neN(2:2:6) = 2*ndel(el,:); N(1:2:5) = N(2:2:6)-1;单元自由度B = AxiaSymmStrain(gxy(ndel(el,:),:));几何矩阵和单元面积S = D*B*U(N);单元应力列向量%16Sm(el) = sqrt(S(1:2)*S(1:2)-S(1)*S(2)+3*S(3)^2);% 等效应力§5.6 程序算例数据function [gxy,ndel,dofix,F,Em,mu,Th,nd,ne,U]= AxiSymExa Th = 1; Em = 1; mu = 0.; ro = 1; w = 1;板厚，弹性模量，泊松比，密度，转速* gxy = [2,3;4,3;2,0;4,0];结点坐标*ndel = [1,4,2; 4,1,3];单元信息*nd = size(gxy,1); ne = size(ndel,1);结点数，单元数F = zeros(2*nd,1); U = zeros(2*nd,1); n = [1,2,3,1,2];for el = 1:neN = ndel(el,:); r = gxy(N,1); rc = sum(r)/3;A = abs(0.5*det([gxy(N,:),ones(3,1)]));单元面积for i = 1:3F(2*N-1) = F(2*N-1)+A*(9*rc^2+2*r(1)^2-r(n(i+1))*r(n(i+2)));%%%17§5.6 程序几何矩阵function [StrainM,A] = AxiaSymmStrain(xy)应变矩阵A = abs(0.5*det([xy,ones(3,1)]));面积rz = sum(xy)/3; n = [1,2,3,1,2];for i = 1:3% 单元% 单元a = xy(n(i+1),1)*xy(n(i+2),2)-xy(n(i+2),1)*xy(n(i+1),2);b = xy(n(i+1),2)-xy(n(i+2),2);c = xy(n(i+2),1)-xy(n(i+1),1); f = a/rz(1)+b+c*rz(2)/rz(1);StrainM(:,2*i-1:2*i) = [b,0;f,0;0,c;c,b]/A/2;矩阵% 几何18结