matlab.rar_matlab风_rayszmf_框架结构模态计算_质量矩阵 MATLAB_风荷载

function frame % %各参数意义很明确 %友好的界面 %格式易读模块化 global node element material K k t f F distant %参数说明 %---------node:结点，由各结点(x,y)坐标构成的矩阵 %---------element:单元，由(单元1号结点,2号结点,单元对应的材料种类号)构成的矩阵 %---------material:各种材料种类号所对应的的属性值 %---------K:总刚，矩阵 %---------k:单刚，三维向量，k(:,:,i)表示第i各单元的单刚 %---------t:转换矩阵，t(:,:,i)表示第i各单元的转换矩阵 %---------f:先用作单元的杆端力(整体坐标系下)，后用作单元的内力(局部坐标系下)，用f(:,i)表示第i个单元 %---------F:等效的结点荷载 %---------distant:各结点的位移 [filename,pathname]=uigetfile('.txt','请选择的数据文件(默认为input.txt)'); %获得输入文件名 filepath=strcat(pathname,filename); %路经＋文件名 if ~isempty(filepath) %缺省为input.txt filenames = strcat('数据文件:',filename); disp(filenames); else filename='input.txt'; end fid = fopen(filename,'r');%////////////////////////////////////////////////////////////////读取输入文件 %先读入控制信息 node_number = fscanf(fid,'%g',1);%结点数 node = fscanf(fid,'%g',[node_number 2]);%结点坐标x1,x2...;y1,y2... element_number = fscanf(fid,'%g',1);%单元数 element = fscanf(fid,'%g',[element_number 3]);%单元的组成：单元1的结点1，单元2的结点1...；单元1的结点2，单元2的结点2...；单元1对应的E、I、A，单元2对应的E、I、A... material_number = fscanf(fid,'%g',1);%不同E、I、A材料的种类数 material = fscanf(fid,'%g',[material_number 3]);% E1.E2...;I1.I2...;A1.A2... weiyi_number = fscanf(fid,'%g',1);%已知位移数 weiyi = fscanf(fid,'%g',[weiyi_number 3]);%结点1，结点2...；结点1的位移方向(1为水平、2为竖向、3为转角)，结点2的位移方向...；位移1，位移2... nf_number = fscanf(fid,'%g',1);%结点力 nf = fscanf(fid,'%g',[nf_number 3]); %结点1，结点2，...；作用于结点1上力的方向(1为水平、2为竖向、3为转角)，作用于结点2上力的方向...；力1的大小，力2的大小... ef_number = fscanf(fid,'%g',1);%单元力 ef = fscanf(fid,'%g',[ef_number 3]); %单元1，单元2，...；作用于单元1上力的种类(1为均布荷载、2为集中荷载、3为集中力偶)，作用于单元2上力的种类...；力1的大小，力2的大小... hhhh = fscanf(fid,'%g',1); %是否考虑轴向变形的标志 %是否考虑轴向变形 if hhhh == 1 disp('此例考虑轴向变形。');%////////////////////////////////////////////////////////考虑轴线变形与否A方法无数倍 else material(:,3)=1e6*material(:,3); %放大轴向变形 disp('此例不考虑轴向变形。'); end %初始化 K = sparse(node_number*3,node_number*3); %系统命令 F = sparse(node_number*3,1); f = zeros(6,element_number); %形成t阵－－转换矩阵 for ie = 1:1:element_number if(ie==1) t = formT(ie); %自定义函数fomT，获得转换矩阵 else t(:,:,ie) = formT(ie); end end %形成单刚k和总刚K for ie = 1:1:element_number if(ie==1) k = StiffnessMatrix(ie); %自定义函数 StiffnessMatrix，获得单刚阵 AssembleStiffnessMatrix(ie,k); %自定义函数 AssembleStiffnessMatrix，获得总阵 else k(:,:,ie) = StiffnessMatrix(ie); AssembleStiffnessMatrix(ie,k(:,:,ie)); end end %结点荷载计算 for inf = 1:1:nf_number n = nf(inf,1); d = nf(inf,2); F((n-1)*3+d) = nf(inf,3); %组装一下 end %disp(F); %单元载荷计算 for ief = 1:1:ef_number f(:,ef(ief,1)) = etn(ef(ief,:)); end for ie = 1:1:element_number Assembleef(ie,f(:,ie)); end %///////////////////////////////////////////////////////////////////////////////////////////乘大数法 for i=1:1:weiyi_number n = weiyi(i,1); d = weiyi(i,2); m = (n-1)*3+d; F(m) = weiyi(i,3)*K(m,m)*1e15; K(m,m) = K(m,m)*1e15; end %求出各位移 distant = K \ F;%/////////////////////////////////////////////////////////////////矩阵相除 % distant=inv(k)*F;也可 %disp(distant); %求出各杆杆端力 for ie = 1:1:element_number prtf(ie); end %形成输出文件 fid=fopen('output.txt','w');%/////////////////////////////////////////////////////////////输出文件 fprintf(fid,'结点位移矩阵▽=\n'); for in = 1:1:node_number for inn = 1:1:3 if abs(distant(3*(in-1)+inn)) > 1e-4 fprintf(fid,'结点%g的位移%g:%.2e\n',in,inn,distant(3*(in-1)+inn)); else ; end end end fprintf(fid,'\n\n'); for ie = 1:1:element_number fprintf(fid,'单元%g的杆端内力:\nF= %.2e , %.2e , %.2e , %.2e , %.2e , %.2e\n\n',ie,f(:,ie)); end fclose(fid); type output.txt return; %求出各杆杆端力的函数 function prtf(ie) global element t k distant f distanttemp = zeros(6,1); for i = 1:1:2 for p = 1:1:3 m = (i-1)*3+p; M = (element(ie,i)-1)*3+p; distanttemp(m) = distant(M); end end f(:,ie) = k(:,:,ie)*distanttemp+f(:,ie); f(:,ie) = t(:,:,ie)*f(:,ie); %disp(f(:,ie)); return %形成整体坐标系下单元固端力 function ft=etn(eft) global node element t ie = eft(1); xi = node(element(ie,1),1); yi = node(element(ie,1),2); xj = node(element(ie,2),1); yj = node(element(ie,2),2); L = ((xj-xi)^2+(yj-yi)^2)^(1/2); kind = eft(2); value = eft(3); ft=zeros(6,1); switch kind case 1 %%均布力 ft(2) = -(value*L)/2; ft(5) = -(value*L)/2; ft(3) = -(value*L^2)/12; ft(6) = (value*L^2)/12; case 2 %%集中力 ft(2) = -value/2; ft(5) = -value/2; ft(3) = -value*L/8; ft(6) = value*L/8; case 3 %%集中力偶 ft(2) = (3*value)/(2*L); ft(5) = -ft(2); ft(3) = value/4; ft(6) = value/4; end ft = transpose(t(:,:,ie))*ft; return %形成整体坐标系下的单刚矩阵 function kt = StiffnessMatrix(ie) global node element material t kt = zeros(6,6); E = material(element(ie,3),1); I = material(element(ie,3),2); A = material(element(ie,3),3); xi = node(element(ie,1),1); yi = node(element(ie,1),2); xj = node(element(ie,2),1); yj = node(element(ie,2),2); L = ((xj-xi)^2+(yj-yi)^2)^(1/2); kt = [ E*A/L 0 0 -E*A/L 0 0 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L -E*A/L 0 0 E*A/L 0 0 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L ];%/////////矩阵的赋值 kt = transpose(t(:,:,ie))*kt*t(:,:,ie);%///////////////////////////////////////////////////////矩阵相乘，转置 return %形成t阵 function tk = formT(ie) global node element xi = node(element(ie,1),1); yi = node(element(ie,1),2); xj = node(element(ie,2),1); yj = node(element(ie,2),2); L = ((xj-xi)^2+(yj-yi)^2)^(1/2); l1 = (xj-xi)/L; m1 = (yj-yi)/L; tk = [ l1 -m1 0 0 0 0 m1 l1 0 0 0 0 0 0 1 0 0 0 0 0 0 l1 -m1 0 0 0 0 m1 l1 0 0 0 0 0 0 1 ]; return %集成总刚矩阵 function AssembleStiffnessMatrix(ie,kt) global element K for i = 1:1:2 for j = 1:1:2 for p = 1:1:3 for q = 1:1:3 m = (i-1)*3+p; n = (j-1)*3+q; M = (element(ie,i)-1)*3+p; N = (element(ie,j)-1)*3+q; K(M,N) = K(M,N)+kt(m,n); end end end end return %形成P阵 function Assembleef(ie,ft) global element F for i = 1:1:2 for p = 1:1:3 m = (i-1)*3+p; M = (element(ie,i)-1)*3+p; F(M) = F(M)-ft(m); end end return