结构系统动力分析可以采用总体结构有限元法但该方法对于

《结构系统动力分析可以采用总体结构有限元法但该方法对于》由会员分享，可在线阅读，更多相关《结构系统动力分析可以采用总体结构有限元法但该方法对于（13页珍藏版）》请在装配图网上搜索。

1、结构系统动力分析可以采用总体结构有限元法，但该方法对于复杂大型结构进行分析存 在计算规模大，计算时间长，所用的磁盘空间、计算机系统太庞大，如飞机、车辆、船舶、 高层建筑等整体结构。特别是用有限元法进行较高频率振动分析时，要求结构被划分成非常 多的单元数以便获得详细的位移和应力特性。这时结构模型的节点自由度可能达到几十万甚 至上百万，直接求解如此庞大的模型是很困难。即使能够分析，也要耗费大量机时，效率极 低。模态综合法（Component Mode Synthesis ）就是在这样的背景下发展起来的一种缩减自 由度方法。它可以将大模型化小，先进行各个子结构的模态分析，然后进行模态综合。由于 仅采

2、用了各个子结构的低阶模态，因而使所建立整体结构动力模型的自由度数大大降低，而 且可以在不同的机器上对各子结构进行模态分析提高计算速度。 模态综合法的基本思想是根据复杂结构的特点将整体结构划分成若干子结构，对各个子 结构分别进行模态分析，得到其动力特性。再利用子结构间力平衡条件及位移协调条件将各 子结构部分低阶模态特性综合，由此得到整体结构的动力特性。 一.模态综合法基本原理 先将整个结构划分为s个子结构，每个子结构应该是容易分析的，子结构之间的连接尽 可能要“弱”些，使子结构之间耦合较小，每个子结构在力学上有较大的独立性。 1.建立各个子结构的模态矩阵［帖（厂=1,2, s），其中砂］

3、=砂F B］ ； © f是交界面完 r r r r r 全固定时所算出的一部分（参加综合的）固有模态，即主模态矩阵；eb是依次释放每一边 r 界自由度，使其得到单位位移而产生的静位移分布所组成的约束模态矩阵。它可以由对子结 k — k jj - 构直接求解一系列静平衡方程（边界自由度发生单位强迫位移）而得到： k —ii k -ji 下标i与j分别代表该子结构的内部及界面自由度。结构作自由振动时，内部自由度上 作用力为零，而界面自由度依次产生单位位移，则相当于^是一个单位矩阵。于是由上式 展开第一行可以得到： Ik ]* L ]+ tc ]* ]= t)] \u L

4、 —Ik 11 * L ]* b] ii i ij i ii ij 从而得到约束模态矩阵，显然e b的列数等于界面自由度数。 r B ]_ 2•用模态矩阵［e i］对子结构的刚度阵冈］和质［mi］量阵做第一次坐标变换，得到各个子 结构的缩减刚度矩阵［k ］和质量矩阵［m ］，再形式地拼装成不耦合的形式: ii [© ]= [e ] 1 ) [e] i ) ) [e ] s [k ] 1 ) ) [k] i ) ) [k ] s [m ] 1 ；[m] = 0 0 [m ] i 0 0 0 [m ] s

5、 其中[k ]=[© ]t[k』e ]，[m ]=[© ]t[mz][© ] i i i i i i i i 一般来说，选定的前N模态组成的模态矩阵，做坐标变化后得到的刚度矩阵，质量矩阵 为NxN阶矩阵，远比变换前的刚度矩阵，质量矩阵的阶数小得多。 3.第二次坐标变换，实现子结构的连接。 假如第r号子结构与第q号子结构有公共界面，则他们的位移向量分别写成： rF rB qF qB 下标 B 表示公共界面上的约束模态广义位移。由于在界面处接口点的位移必须相等，这 就导致与该位移相对应大的约束模态广义位移亦必须相等， 也就有： rB qB prF prB

6、rF pqF rB qF pqB I 同样的，对于整个结构可以形式地写出， p ■[ I ] 0 _ F p - 0 [I ] B p ' B' 0 [丫 ]_ * [件]二[卩]* [ p] p B 式中，[p]是相应于所有子结构主模态的广义位移；[p]是所有独立的约束模态广义 FB 位移的集合；[p」是可以用[p」表示的非独立约束模态广义位移，[丫]表示其间转换关系。 B' B 利用[0 ]可以进一步将未耦合的总刚度阵[k]及总质量阵[m]缩减（耦合）为： [K]二[0]T[k][0] ； [

7、M]二[0]T[m][0] 最后，求解下列广义特征值问题： [K] * [x] 2[M][x] 它的各低阶特征值与原结构低阶特征值的误差，跟子结构剖分是否得当，截取模态是否 合理有很大关系。如果想要得到物理坐标下的结构模态向量（振型），需要逆向做两次坐标 变换。 下面以悬臂板模态分析为例，分别用商用软件ANSYS和基于有限元算法的MATLAB程 序，演示固定界面模态综合法的应用。 二.ANSYS模态综合法求解 ANSYS模态综合法是子结构在动力分析中的应用，其基本过程包括三方面：生成过程、 使用过程、扩展过程。ANSYS提供了友好的模态综合法(Component Mode Synt

8、hesis)CMS 向导(Preprocessor>Modeling>CMS),可以方便的定义超单元和交界面，而且可以对模态综 合法分析生成的文件进行管理和组织。 1. 创建超单元 2. 选择CMS求解方法(CMSOPT)：固定界面法或自由界面法。同时确定提取模态数、 频率范围等。对于自由界面法还要确定刚体模态。 3. 命名超单元矩阵文件(SEOPT)。 4. 对于考虑阻尼时，输入集中质量矩阵公式。 5. 定义主自由度：在超单元的交界面定义主自由度。 6. 保持数据库：这是必须做的，因为在扩展模态时需要有相同的数据库文件。 7. 求解生成超单元矩阵文件。 2. CMS的使用和扩

9、展过程 CMS的使用和扩展部分与子结构基本相同，但是CMS只支持模态分析。在子结构分析 中，需要对整体结构中的每一个子结构进行生成和扩展，然后将各个子结构集合成整个模型。 在自由界面模态综合法中，使用MODOPT扩展模态数，应小于每个超单元求解的模态数。 若需要扩展更多的模态，需要CMS的超单元重新求解更多的阶数。 以长和宽2的正方形悬臂板为例： 材料参数：弹性模量E=2.1ell泊松比丫=0.3密度p=7.3e3厚度t=0.05 网格划分：在长度和宽度方向上均划分20个单元(为了与后面程序一致)

10、

11、

12、 X /PREP7 BLC4, , ,2,2 !前处理 画 2X2 矩形 MP,DENS,l,7.3e3 !密度 MP,EX,l,2.1ell !弹性模量 MP,PRXY,l,0.3 !泊松比 ET,1,SHELL63 !单元类型 R,l,0.05,0.05,0.05,0.05,

13、 , , !实常数，厚度 TYPE,1 MAT,1REAL,1 !材料及单元类型编号 ESIZE,O.1,O, !变长 0.1 的单元 AMESH,all !画网格 DL,4, ,ALL, !施加固定约束 FINISH 直接求解命令流！整体求解 /SOL ANTYPE,2 MODOPT,LANB,10 EQSLVSPAR MODOPT,LANB,20,0,99999999, ,OFF SOLVE 表一 ANSYS直接求解结果 阶数 1 2 3 4 5 6 7 8 9 10 频率 11.2 27.4 68.7 87.8 100.1 175.2

14、 197.8 207.3 229.5 281.3 按固定约束方向，在板长度方向的中点将方板一分为二，得到两个子结构。 固定界面模态综合法求解命令流: /clear,nostart !清空变量，启动一个新分析 finish /filnam,use !文件名use整体体结构 /solu /prep7 antype,modal !分析类型 et,1,matrix50 !定义超单元 modopt,lanb,10 !求解前10阶频率 type,1 solve se,part1 !选择子结构一 finish se,part2 !选择

15、子结构一 save 表二ANSYS模态综合法求解结果 子结构一 阶数 1 2 3 4 5 6 7 8 9 10 频率 287.5 302.3 355.3 458.9 624.8 792.3 812.7 857.9 879.8 995.8 子结构二 阶数 1 2 3 4 5 6 7 8 9 10 频率 45.1 69.2 131.8 247.1 282.1 318.7 407.1 439.6 557.5 682.2 整体结构 阶数 1 2 3 4 5 6 7 8

16、9 10 频率 11.2 27.5 68.7 87.8 100.1 175.2 197.8 207.3 229.5 281.5 /filnam,partl ！ 文件名parti子结构一 /solu antype,substr ！ 分析类型子结构 seopt,part1,2 !子结构一 cmsopt,fix,20 !固定界面法求前20阶频率 cmsel,s,part1 !选定parti单元集合 cmsel,s,interface ! interface节点集合 m,all,all !选择所有 nsle solve finish save

17、 !求解并保存 /filnam,part2 !文件名part2子结构二 /solu antype,substr !分析类型子结构 seopt,part2,2 !子结构二 cmsopt,fix,20 !固定界面法求前20阶频率 cmsel,s,part2 !选定part2单元集合 cmsel,s,interface ! interface节点集合 m,all,all !选择所有 nsle solve finish save !求解并保存 与表一中直接求解结构对比，可以看出本例中，ANSYS模态综合法具有较高的精度（保 留一位小数），前10阶频率

18、基本上是一致的。在计算大型复杂结构，模态综合法可以节省大 量的计算时间和计算机资源，提高效率；同时可以灵活修改大系统的子系统设计。因此，对 于复杂大型结构，如飞机、车辆、船舶、高层建筑等结构，采用ANSYS模态综合法来对结 构进行模态分析，可以在精度和计算速度上得到较好的解决方案。 三.基于有限元方法的模态综合法MATLAB实现 1.薄板弯曲理论 等厚度薄板的坐标系如图(a)所示，板厚1/2平面，即xy平面为板的中面。从板中任取 一微元体dxdydz，在每一个面上作用的正应力和剪应力见图(b)。 薄板弯曲有如下三个假设： (1)垂直于中面方向的正应变*z极微小，可以忽略。

19、取£二0，由S =-° = 0得 z Z -z y)，说明板的任何一点的挠度O只与坐标x和y有关。 (2)应力分量 zx T 一 , 和zy远小于 xy， 因此可以忽略不计它们产生的正 Y 、剪应变 zx zy， 则由物理方程有 -po 2(1 + p) xy xy (3)薄板弯曲时，中面内各点不产生平行于中面的应变，即 x z=0 -x y z=0 xy z=0 综合上述三个假设, 可用得到: -v -u =——+—— -x -2O -2O -2O - 2O -x 2 -x 2 xy -2O ) 2(1 - p)' -x

20、y -2O -t L z zo dz Et 3 -2O 12(1 -p 2) 上述应力在板内分别合成为薄板内力 zo dz xy -2O -x-x 简写成矩阵的形式F = DK，其中K =[- -2O -2O -2O ]t = Lo -x-x zo dz 2. 薄板弯曲单元 基于薄板弯曲理论的板单元中每个节点有三个自由度：绕度°，法线绕X轴转角0 ,绕 x y轴转角0 y '即 [a ]二 i 0 二[① xi 0 y i — i ]t dy dx ae = [a a a ]t 1 i n 等厚度板单元中，有4个

21、节点 标系,即 单元宽度为a,高度为b,在单元中心(x0,y0)建立局部坐 O x = x + ag 0 其中, y = y0 + bn a =丄(x — x ) 2 2 1 b = 1( y —y ) 2 2丿1 1 x = (x + x ) 0 2 2 1 y = 1( y + y ) 0 2 2 1 则单元内任 ii w =》(N①+ N① i i ix ix 1 + N ①)= ^L[N ][a ] iy iy i i 1 其中, [N ]=[N ,N ,N ] i i ix iy n = ](i+g

22、g )(1+nn )(2+gg +nn —g 2 —n2) i 8 i i i i n =1 bn (i+盘)(i+nn )(i—n2) ix 8 i i i n =1 ag (i+盘)(i+nn )(i—g 2) iy 8 i i i i = i,2,3,4 写成矩阵的形式W = Nae，则单元刚度矩阵可写成Ke =J BTDBdV，其中B = LN 单元质量矩阵可写成Me = J p NTNdV ，最终得到的单元质量矩阵和刚度矩阵都是 i2 X12的矩阵。 3.有限元方法求解悬臂板模态MATLAB程序 %清空变量 %弹性模量 % 泊松比 %密度

23、%板的厚度 %乂方向长度 %丫方向长度 %乂向节点数 %丫向节点数 clear all clc tic %记时起点 % 定义材料 E=2.1e11; poisson =0.3; density=7.3e3; t=0.05; lx=2; ly=2; jdx=21; jdy=21; %初始化 %刚度和质量阵 Tol_dof=3*jdx*(jdy-1); %在y=0边界上施加固定约束，去掉21个节点 k(1:Tol_dof,1:Tol_dof)=0; %总刚阵 m(1:Tol_dof,1:Tol_dof)=0; %总质量阵 Tol_element=(jdx-1)*(j

24、dy-1); % 总单元数 %节点矩阵 en(1:Tol_element,1:4)=0; %每个单元四个节点，记录节点编号 for ni=1:jdx-1 for nj=1:jdy-1 en(ni+(nj-1)*(jdx-1),1)=ni+(nj-1)*jdx; en(ni+(nj-1)*(jdx-1),2)=ni+1+(nj-1)*jdx; en(ni+(nj-1)*(jdx-1),4)=ni+nj*jdx; en(ni+(nj-1)*(jdx-1),3)=ni+1+nj*jdx; end end %节点位移，约束以及自由度坐标 disp(1:jdx*jdy,1:3)=1; %

25、位移阵 constraints=1:1:jdx; % y=0 上被约束节点 disp(constraints,:)=0; dof=0; %自由度标记 for ni=1:jdx*jdy for nj=1:3 if disp(ni,nj)~=0 %有约束不编号 dof=dof+1; disp(ni,nj)=dof; %单元长度 el=lx/(jdx-1); %单元长度 eh=ly/(jdy-1); %单元高度 %四节点12自由度 单元刚度和质量矩阵 % km函数，输入材料参数和节点集中长度 [ek,dm]=km(el/2,eh/2,poisson,t,E,density

26、); %整体刚度矩阵和质量矩阵 tg(1:12)=0; % 标记，单元矩阵维数 for loopi=1:(jdx-1)*(jdy-1) %外循环，单元 for zi=1:4 %每个节点三个自由度 tg((zi-1)*3+1)=disp(en(loopi,zi),1); tg((zi-1)*3+2)=disp(en(loopi,zi),2); tg((zi-1)*3+3)=disp(en(loopi,zi),3); end for jx=1:12 for jy=1:12 if(tg(jx)*tg(jy)~=0) k(tg(jx),tg(jy))=k(tg(jx),tg(jy))

27、+ek(jx,jy); m(tg(jx),tg(jy))=m(tg(jx),tg(jy))+dm(jx,jy; end end end %总装刚度阵和质量阵 end %求解特征值问题 [v,d] = eig(k,m); %eig函数求解特征值问题 tempd=diag(d); %记录特征值 [nd,sortindex]=sort(tempd); %对特征值排序 v=v(:,sortindex); %特征向量排序 mode_number=1:15; %求解阶数 frequency(mode_number)=sqrt(nd(mode_number) )/(2*pi);%记录角频

28、率 toc %记时终点 end end end 表三MATLAB直接求解结果 阶数 1 2 3 4 5 6 7 8 9 10 频率 11.208 27.469 68.764 87.797 99.961 174.77 197.92 207.16 229.18 299.24 相对误差(％) 0.07 0.25 0.09 0.03 0.13 0.24 0.06 0.06 0.13 6 DtEFLArSIEKT 第十阶频率有很大差别，取出ANSYS分析结果的振型分析，可知第十节频率为平行于 固支边

29、的弯曲振型，在薄板弯曲单元中不包括这个方向上的自由度，因而取ANSYS分析结 果中的第^一阶频率，其大小为300.1Hz,与MATLAB结果对比，误差为0.29% 畑 34 U/ia 4.基于有限元的模态综合法求解悬臂板模态MATLAB程序 子结构的划分与ANSYS中处理一致，主模态分析程序如下： 子结构一(只列出修改部分)： 子结构二(只列出修改部分): lx=2; %乂方向长度ly=l; %丫方向长度 jdx=21; %*向节点数 jdy=ll; %丫向节点数 %y=0和1沿x轴去掉42个点，得到总自由度 Tol_dof=3*jdx*(jdy-2);

30、 cons_0ri=1:1:jdx; %原有约束节点 %固定界面法人为施加约束 cons_Arti=jdx*(jdy-1)+1:1:jdx*jdy; %总约束 constraints=[cons_Ori,cons_Arti]; master_mode=15;add_num=jdx*3; %取前master_mode阶振型作为模态综合 fi_unres1=cat(1,v(1:master_mode,:)',zeros(a dd_num,master_mode)); lx=2; %乂方向长度 ly=1; %丫方向长度 jdx=21; %*向节点数 jdy=11; %丫向节点数 %总

31、自由度 Tol_dof=3*jdx*(jdy-1); %固定界面法人为施加约束 cons_Arti=1:1:jdx; %单边约束节点编号 constraints=cons_Arti; %总约束 master_mode=15;add_num=jdx*3; fi_unres2=cat(1,zeros(add_num,master_mode), v(1:master_mode,:)'); %主模态中，认为施加固定约束的节点，其位移为 零，补充完整加入模态综合 表四MATLAB子结构求解结果 子结构一 阶数 1 2 3 4 5 6 7 8 9 10 频

32、率 287.8 302.1 353.9 455.9 620.3 794.8 813.7 852.3 877.7 988.8 子结构二 阶数 1 2 3 4 5 6 7 8 9 10 频率 45.1 69.1 131.5 246.4 282.3 318.5 405.4 438.8 553.6 681.2 约束模态分析程序如下： 子结构一（只列出修改部分） lx=2; %乂方向长度ly=l; %丫方向长度 jdx=21; %*向节点数 jdy=ll; %丫向节点数 Tol_dof=3*jdx*(jdy-1); cons_O

33、ri=1:1:jdx; %原有约束节点 constraints=cons_Ori; %总约束 % 总刚度矩阵中按是否在界面上分块 add_num=jdx*3;master_mode=15; kii=k(1:end-add_num,1:end-add_num); kij=k(1:end-add_num,end-add_num+1:end); %提取约束模态 fi_res1=cat(1,-inv(kii)*kij,eye(add_num )); %合并总模态 fi_1=cat(2,fi_unres1,fi_res1); %用模态对刚阵和质量阵做第一次坐标变换 k1=fi_1'*k*fi_1;

34、m1=fi_1'*m*fi_1; 子结构二(只列出修改部分)： lx=2; %乂方向长度 ly=1; %丫方向长度 jdx=21; %*向节点数 jdy=ll; %丫向节点数 %总自由度 Tol_dof=3*jdx*jdy; %子结构没有约束 % 总刚度矩阵中按是否在界面上分块 add_num=jdx*3;master_mode=15; kii=k(add_num+1:end,add_num+1:end); kij=k(add_num+1:end,1:add_num); %提取约束模态 fi_res2=cat(1,eye(add_num),-inv(kii)*ki j); %合

35、并总模态 fi_2=cat(2,fi_res2,fi_unres2); %用模态对刚阵和质量阵做第一次坐标变换 k2=fi_2'*k*fi_2; m2=fi_2'*m*fi_2; 模态综合程序如下： master_mode=15;add_num=jdx*3;Tol_num=length(k1); %master_mode为主模态保留阶数，add_num为界面上自由度数，Tol_num为变换后刚度矩阵维数 K=zeros(2*Tol_num);M=zeros(2*Tol_num); %初始化模态综合矩阵 K(1:Tol_num,1:Tol_num)=k1;K(Tol_num+1:en

36、d,Tol_num+1:end)=k2; M(1:Tol_num,1:Tol_num)=m1;M(Tol_num+1:end,Tol_num+1:end)=m2; %非耦合模态综合矩阵，直接写成矩阵形式 beta=zeros(2*Tol_num,2*master_mode+add_num); beta(1:master_mode,1:master_mode)=eye(master_mode); beta(master_mode+1:master_mode+add_num,master_mode+1:master_mode+add_num)=eye(add_n um); beta(maste

37、r_mode+1+add_num:master_mode+2*add_num,master_mode+1:master_mode+add_num) =eye(add_num); beta(end-master_mode+1:end,end-master_mode+1:end)=eye(master_mode); %第二次坐标变换矩阵be ta,由位移列车很容易求的 KK=beta'*K*beta;MM=beta'*M*beta; %做第二次坐标变换 [v,d] = eig(KK,MM);tempd=diag(d);[nd,sortindex]=sort(tempd);v=v(:,sor

38、tindex); %求解特征值问题 mode_number=1:5;frequency_sum(mode_number)=sqrt(nd(mode_number))/(2*pi); %得到最终模态综合方法求出频率 表五 MATLAB 模态综合法求解结果 悬臂板 阶数 1 2 3 4 5 6 7 8 9 10 频率 11.5 29.1 87.7 106.8 152.5 260.1 309.6 496.1 645.2 818.3 function [k,m]=km(a,b,poisson,t,E,density) %变量a,b,poisson,t

39、,E,densit y与主程序中定义一致 %kk 为刚阵列向量，(太大了，没办法，定义出来只能是向量) kk=[E*t"3/(360-360*poisson"2)/a/b*(21-6*poisson+30*b"2/a"2+30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*a-12*poisson*a-30*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson-30*b"2/

40、a"2+15*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*b-12*poisson*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a-30*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(21-6*poisson-15*b"2/a"2T5*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*

41、a-3*poisson*aT5*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson+15*b"2/a"2-30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poisson*aT5*b"2/a); E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a

42、/b*(8*b"2-8*poisson*b"2+40*a"2), -30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(-3*bT2*poisson*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-8*b"2+8*poisson*b"2+20*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*bT5*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(2*b"2-2*pois

43、son*b"2+10*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b-30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-2*b"2+2*poisson*b"2+20*a"2), 0; E*t"3/(360-360*poisson"2)/a/b*(-3*aT2*poisson*a-30*b"2/a), -30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(8*a"2-8*poisson*a"2+4

44、0*b"2), E*t"3/(360-360*poisson"2)/a/b*(3*a-3*poisson*a+30*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-2*a"2+2*poisson*a"2+20*b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a+15*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(2*a"2-2*poisson*a"2+10*b"2), E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poi

45、sson*aT5*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-8*a"2+8*poisson*a"2+20*b"2); E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson-30*b"2/a"2+15*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*bT2*poisson*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a-3*poisson*a+30*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(

46、21-6*poisson+30*b"2/a"2+30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poisson*a+30*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson+15*b"2/a"2-30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*b+30*a"2/b), E*t"3/(360-360

47、*poisson"2)/a/b*(-3*aT2*poisson*a+15*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(21-6*poissonT5*b"2/a"2T5*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a+15*b"2/a); E*t"3/(360-360*poisson"2)/a/b*(-3*bT2*poisson*b+15*a"2/b), E*t"3/(

48、360-360*poisson"2)/a/b*(-8*b"2+8*poisson*b"2+20*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(8*b"2-8*poisson*b"2+40*a"2), 30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b-30*a"2/b), E*t"3/(360-360*poisso

49、n"2)/a/b*(-2*b"2+2*poisson*b"2+20*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*bT5*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(2*b"2-2*poisson*b"2+10*a"2), 0; E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a-30*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-2*a"2+2*poisson*a"2+20*b"2), E*t"3/

50、(360-360*poisson"2)/a/b*(3*a+12*poisson*a+30*b"2/a), 30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(8*a"2-8*poisson*a"2+40*b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*aT2*poisson*a+15*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-8*a"2+8*poisson*a"2+20*b"2), E*t"3/(360-360*poiss

51、on"2)/a/b*(3*a-3*poisson*a-15*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(2*a"2-2*poisson*a"2+10*b"2); E*t"3/(360-360*poisson"2)/a/b*(21-6*poissonT5*b"2/a"2T5*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*bT5*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a+15*b"2/a), E*t"3/(360-

52、360*poisson"2)/a/b*(-21+6*poisson+15*b"2/a"2-30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b-30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*aT2*poisson*a+15*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(21-6*poisson+30*b"2/a"2+30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*bT2*poisson*b-30*a"

53、2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poisson*a+30*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson-30*b"2/a"2+15*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*bT5*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a-3*poisson*a+30*b"2/a); E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson

54、*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(2*b"2-2*poisson*b"2+10*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-2*b"2+2*poisson*b"2+20*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(-3*bT2*poisson*b-30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(8*b"

55、2-8*poisson*b"2+40*a"2), -30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*bT5*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-8*b"2+8*poisson*b"2+20*a"2), 0; E*t"3/(360-360*poisson"2)/a/b*(3*a-3*poisson*aT5*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(2*a"2-2*poisson*

56、a"2+10*b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*aT2*poisson*a+15*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-8*a"2+8*poisson*a"2+20*b"2), E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poisson*a+30*b"2/a), -30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(8*a"2-8*poisson*a"2+40*b"2),

57、 E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a-30*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-2*a"2+2*poisson*a"2+20*b"2); E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson+15*b"2/a"2-30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b-30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poiss

58、on*a-15*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(21-6*poissonT5*b"2/a"2T5*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(3*b-3*poisson*bT5*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(3*a-3*poisson*aT5*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(-21+6*poisson-30*b"2/a"2+15*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*

59、(3*b+12*poisson*b-15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a-30*b"2/a), E*t"3/(360-360*poisson"2)/a/b*(21-6*poisson+30*b"2/a"2+30*a"2/b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*bT2*poisson*b-30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-3*aT2*poisson*a-30*b"2/a); E*t"3/(360-360*poisso

60、n"2)/a/b*(3*b-3*poisson*b+30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(-2*b"2+2*poisson*b"2+20*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(-3*b+3*poisson*b+15*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(2*b"2-2*poisson*b"2+10*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(3*b+12*poisson*bT5*a"2/b), E*t"3/(360-3

61、60*poisson"2)/a/b*(-8*b"2+8*poisson*b"2+20*a"2), 0, E*t"3/(360-360*poisson"2)/a/b*(-3*b-12*poisson*b-30*a"2/b), E*t"3/(360-360*poisson"2)/a/b*(8*b"2-8*poisson*b"2+40*a"2), 30*E*t"3/(360-360*poisson"2)*poisson; E*t"3/(360-360*poisson"2)/a/b*(3*a+12*poisson*aT5*b"2/a), 0, E*t"3/(360-360*poisson"

62、2)/a/b*(-8*a"2+8*poisson*a"2+20*b"2), E*t"3/(360-360*poisson"2)/a/b*(-3*a+3*poisson*a+15*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(2*a"2-2*poisson*a"2+10*b"2), E*t"3/(360-360*poisson"2)/a/b*(3*a-3*poisson*a+30*b"2/a), 0, E*t"3/(360-360*poisson"2)/a/b*(-2*a"2+2*poisson*a"2+20*b"2), E*t"3/(360-360*p

63、oisson"2)/a/b*(-3*aT2*poisson*a-30*b"2/a), 30*E*t"3/(360-360*poisson"2)*poisson, E*t"3/(360-360*poisson"2)/a/b*(8*a"2-8*poisson*a"2+40*b"2)]; % k将kk刚度向量用12*12的矩阵形式重新存储 k=zeros(12,12); for i=1:12 for j=1:12 k(i,j)=kk(12*(i-1)+j,1); end end % w为集中质量(平均分布在四个节点上) w=a*b*t*density; % 用符号变量求解形函数积

64、分求得单元质量阵 syms kx ytkxi yti real; ni=1/8*(1+kx*kxi)*(1+yt*yti)*(2+kx*kxi+yt*yti-kx"2-yt"2); %绕度 nix=T/8*b*yti*(1+kx*kxi)*(1+yt*yti)*(1-yt"2); %x轴转角 niy=1/8*a*kxi*(1+kx*kxi)*(1+yt*yti)*(1-kx"2); %y轴转角 n(1)=subs(ni,{kxi,yti},{-1,-1});n(2)=subs(nix,{kxi,yti},{-1,-1}); n(3)=subs(niy,{kxi,yti},{-1,-1

65、});n(4)=subs(ni,{kxi,yti},{1,-1}); n(5)=subs(nix,{kxi,yti},{1,-1});n(6)=subs(niy,{kxi,yti},{1,-1}); n(7)=subs(ni,{kxi,yti},{1,1});n(8)=subs(nix,{kxi,yti},{1,1}); n(9)=subs(niy,{kxi,yti},{1,1});n(10)=subs(ni,{kxi,yti},{-1,1}); n(11)=subs(nix,{kxi,yti},{-1,1});n(12)=subs(niy,{kxi,yti},{-1,1}); temp=n'*n;m1=int(temp,kx,-1,1); %局部坐标kx方向积分 m=int(m1,yt,-1,1); %局部坐标yt方向积分 m=m*w;m=double(m);% 符号变量转换成双精度格式数值格式 % %guyan 缩距 % m(i,j)=0; % for i=1:12 % end % for j=1:12 % end % if(m(i,j))<1e-1 % end