training.zip_一维梁_有限元_梁 有限元_梁振动_频响

%%用一维的平面梁单元建立梁模型，求解梁的谐响应 %%总体思路就是建立单元质量/刚度矩阵，再组件总质量/刚度矩阵，再求谐响应 % L 为梁的长度 % b 为梁的宽度 % h 为梁的厚度 % E 杨氏模量 % I 惯性矩 % ln 单元个数 % 定义全局变量 % gNode ------ 节点坐标 % gElement --- 单元定义 % gMaterial -- 材料性质 % gNF -------- 集中力的幅值 % gK --------- 整体刚度矩阵 % gM --------- 整体质量矩阵 % gDelta ----- 节点位移 global gNode gElement gMaterial gK gM gf gDelta freq L ln L=1.2; b=0.04; h=0.005; ln=30; dx = L/ln ; % 单元的宽度 nf = 1e1 ; % 激励的幅值（单位：N） % 节点坐标 gNode = zeros( ln+1, 2 ) ; for i=1:ln+1 k = i ; % 节点号 xk = (i-1)*dx ; yk = 0; % 节点的x坐标 gNode(k,:) = [xk, yk] ; end % 单元定义，一个梁单元有两个节点 gElement = zeros( ln, 2 ) ; for i=1:ln k = i ; % 单元号 n1 = i; % 第一个节点号 n2 = i+1 ; % 第二个节点号 gElement(k,:) = [n1, n2] ; end % 材料性质 % 每个原胞由三个单元组成，1,3单元为材料1，2单元为材料2 % 弹性模量 抗弯惯性矩I 截面积A 密度 gMaterial = [70e9, b*h^3/12, b*h, 2700 % 材料 1 70e9, b*h^3/12, b*h, 2700] ; % 材料 2 %激励力 节点 自由度号 集中力值 gNF = [1, 1, nf ] ; % 定义总体位移矩阵 gDelta = zeros( (ln+1)*2, 1 ) ; % 1. 定义整体刚度矩阵、质量矩阵和节点力向量 [node_number, dummy]=size(gNode); gK = sparse(node_number*2,node_number*2); gM = sparse(node_number*2,node_number*2); gf = sparse(node_number*2,1); % 2. 计算单元刚度矩阵、质量矩阵，并集成到整体刚度矩阵和质量矩阵中 [element_nember, dummy]=size(gElement); for ie=1:1:element_nember; ke = StiffnessMatrix(ie); me = MassMatrix( ie ) ; AssembleGlobalMatrix( ie, ke, me ) ; end % step2.1 把集中力直接集成到整体节点力向量中 [nf_number,dummy]=size(gNF); for inf=1:1:nf_number ielem = gNF( inf, 1 ) ; %节点 iedge = gNF( inf, 2 ) ; %自由度号 gf( (ielem-1)*2 + iedge ) = gNF( inf, 3 ) ; end % step3. 给定频率，计算稳态响应 % 假设激励f=p*sin(2*pi*freq*t),节点位移响应为gDisp=Disp*sin(2*pi*freq*t), % 则节点位移向量gDisp=(gK-omega^2*gM)\f； freq=2000; gDelta=zeros(node_number * 2, freq); displacement=zeros(1,freq); for i=10:1:freq omega=2*pi*i; gK=full(gK); gM=full(gM); gf =full(gf); gDelta(:,i)=abs((gK-omega^2*gM)\gf); displacement(1,i)=20*log10(gDelta((node_number-1)*2+1,i)/gDelta(1,i)); end fprintf( '节点位移\n' ) ; fprintf( ' 频率 z方向位移 ') ; [node_number,dummy] = size( gNode ) ; % 显示梁末端的位移 for i=1:freq fprintf( '%6d %16.8e\n',... i, gDelta((node_number-1)*2+1, i )); end plot(gDelta((node_number-1)*2+1, :) ); set(gca,'yscale','log') %传递率曲线 figure plot(displacement(1, :) ); function ke = StiffnessMatrix(ie) % 计算单元刚度矩阵 % 输入参数: % ie ------- 单元号 % icoord -- 坐标系参数，可以是下面两个之一 % 1 ---- 整体坐标系 % 2 ---- 局部坐标系 % 返回值: % k ---- 根据icoord的值，相应坐标系下的刚度矩阵 global gNode gElement gMaterial ln L % k = zeros( 4, 4 ) ; % 每个原胞由三个单元组成，1,3单元为材料1，2单元为材料2 % 先对ie的值进行判断，ie取3的余数若为2，即给单元赋予材料2的属性 if rem(ie,3)==2 E = gMaterial( 2, 1 ) ; I = gMaterial( 2, 2 ) ; A = gMaterial( 2, 3 ) ; else E = gMaterial( 1, 1 ) ; I = gMaterial( 1, 2 ) ; A = gMaterial( 1, 3 ) ; end % xi = gNode( gElement( ie, 1 ), 1 ) ; % yi = gNode( gElement( ie, 1 ), 2 ) ; % xj = gNode( gElement( ie, 2 ), 1 ) ; % yj = gNode( gElement( ie, 2 ), 2 ) ; % l = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ; l = L/ln; k = [ 12 6*l -12 6*l; 6*l 4*l^2 -6*l 2*l^2; -12 -6*l 12 -6*l; 6*l 2*l^2 -6*l 4*l^2]; ke = E*I/l^3 * k ; return end function me = MassMatrix( ie ) % 计算单元质量矩阵 % 输入参数: % ie ------- 单元号 % 返回值: % m ---- 整体坐标系下的质量矩阵 global gNode gElement gMaterial rho A ln L % me = zeros( 4, 4 ) ; % 每个原胞由三个单元组成，1,3单元为材料1，2单元为材料2 % 先对ie的值进行判断，ie取3的余数若为2，即给单元赋予材料2的属性 if rem(ie,3)==2 A = gMaterial( 2, 3 ) ; rho = gMaterial( 2, 4 ) ; else A = gMaterial( 1, 3 ) ; rho = gMaterial( 1, 4 ) ; end % xi = gNode( gElement( ie, 1 ), 1 ) ; % yi = gNode( gElement( ie, 1 ), 2 ) ; % xj = gNode( gElement( ie, 2 ), 1 ) ; % yj = gNode( gElement( ie, 2 ), 2 ) ; % l = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ; l = L/ln; me = rho*A*l/420*[ 156, 22*l, 54, -13*l 22*l, 4*l^2, 13*l, -3*l^2 54, 13*l, 156, -22*l -13*l, -3*l^2, -22*l, 4*l^2]; return end function AssembleGlobalMatrix( ie, ke ,me ) % 把单元刚度矩阵集成到整体刚度矩阵 % 输入参数: % ie --- 单元号 % k --- 单元刚度矩阵 % 返回值: % 无 global gElement gK gM for i=1:1:2 for j=1:1:2 for p=1:1:2 for q=1:1:2 m = (i-1)*2+p ; n = (j-1)*2+q ; M = (gElement(ie,i)-1)*2+p ; N = (gElement(ie,j)-1)*2+q ; gK(M,N) = gK(M,N) + ke(m,n) ; gM(M,N) = gM(M,N) + me(m,n) ; end end end end return end