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

结构系统动力分析可以采用总体结构有限元法，但该方法对于复杂大型结构进行分析存 在计算规模大，计算时间长，所用的磁盘空间、计算机系统太庞大，如飞机、车辆、船舶、 高层建筑等整体结构。特别是用有限元法进行较高频率振动分析时，要求结构被划分成非常 多的单元数以便获得详细的位移和应力特性。这时结构模型的节点自由度可能达到几十万甚 至上百万，直接求解如此庞大的模型是很困难。即使能够分析，也要耗费大量机时，效率极 低。模态综合法（Component Mode Synthesis ）就是在这样的背景下发展起来的一种缩减自 由度方法。它可以将大模型化小，先进行各个子结构的模态分析，然后进行模态综合。由于 仅采用了各个子结构的低阶模态，因而使所建立整体结构动力模型的自由度数大大降低，而 且可以在不同的机器上对各子结构进行模态分析提高计算速度。模态综合法的基本思想是根据复杂结构的特点将整体结构划分成若干子结构，对各个子 结构分别进行模态分析，得到其动力特性。再利用子结构间力平衡条件及位移协调条件将各 子结构部分低阶模态特性综合，由此得到整体结构的动力特性。一.模态综合法基本原理先将整个结构划分为s个子结构，每个子结构应该是容易分析的，子结构之间的连接尽 可能要“弱”些，使子结构之间耦合较小，每个子结构在力学上有较大的独立性。1.建立各个子结构的模态矩阵帖（厂=1,2, s），其中砂=砂F B ； f是交界面完rrr rr全固定时所算出的一部分（参加综合的）固有模态，即主模态矩阵；eb是依次释放每一边 r界自由度，使其得到单位位移而产生的静位移分布所组成的约束模态矩阵。它可以由对子结kkjj -构直接求解一系列静平衡方程（边界自由度发生单位强迫位移）而得到：kiik-ji下标i与j分别代表该子结构的内部及界面自由度。结构作自由振动时，内部自由度上 作用力为零，而界面自由度依次产生单位位移，则相当于是一个单位矩阵。于是由上式 展开第一行可以得到：Ik* L +tc*=t)uL Ik11 * L*bii iijiiiij从而得到约束模态矩阵，显然e b的列数等于界面自由度数。rB_2用模态矩阵e i对子结构的刚度阵冈和质mi量阵做第一次坐标变换，得到各个子结构的缩减刚度矩阵k 和质量矩阵m ，再形式地拼装成不耦合的形式: ii =e 1)ei)e sk 1)ki)k sm 1；m =00m i000m s其中k = tke ，m = tmz iii iiii i一般来说，选定的前N模态组成的模态矩阵，做坐标变化后得到的刚度矩阵，质量矩阵 为NxN阶矩阵，远比变换前的刚度矩阵，质量矩阵的阶数小得多。3.第二次坐标变换，实现子结构的连接。假如第r号子结构与第q号子结构有公共界面，则他们的位移向量分别写成：rFrBqFqB下标 B 表示公共界面上的约束模态广义位移。由于在界面处接口点的位移必须相等，这就导致与该位移相对应大的约束模态广义位移亦必须相等，也就有：rB qBprFprBrFpqFrBqFpqB I同样的，对于整个结构可以形式地写出，p I 0 _Fp-0I Bp B0丫 _* 件二卩* ppB式中，p是相应于所有子结构主模态的广义位移；p是所有独立的约束模态广义FB位移的集合；p是可以用p表示的非独立约束模态广义位移，丫表示其间转换关系。BB利用0 可以进一步将未耦合的总刚度阵k及总质量阵m缩减（耦合）为：K二0Tk0 ； M二0Tm0最后，求解下列广义特征值问题：K * x2Mx它的各低阶特征值与原结构低阶特征值的误差，跟子结构剖分是否得当，截取模态是否 合理有很大关系。如果想要得到物理坐标下的结构模态向量（振型），需要逆向做两次坐标 变换。下面以悬臂板模态分析为例，分别用商用软件ANSYS和基于有限元算法的MATLAB程 序，演示固定界面模态综合法的应用。二.ANSYS模态综合法求解ANSYS模态综合法是子结构在动力分析中的应用，其基本过程包括三方面：生成过程、 使用过程、扩展过程。ANSYS提供了友好的模态综合法(Component Mode Synthesis)CMS 向导(PreprocessorModelingCMS),可以方便的定义超单元和交界面，而且可以对模态综 合法分析生成的文件进行管理和组织。1. 创建超单元2. 选择CMS求解方法(CMSOPT)：固定界面法或自由界面法。同时确定提取模态数、 频率范围等。对于自由界面法还要确定刚体模态。3. 命名超单元矩阵文件(SEOPT)。4. 对于考虑阻尼时，输入集中质量矩阵公式。5. 定义主自由度：在超单元的交界面定义主自由度。6. 保持数据库：这是必须做的，因为在扩展模态时需要有相同的数据库文件。7. 求解生成超单元矩阵文件。2. CMS的使用和扩展过程CMS的使用和扩展部分与子结构基本相同，但是CMS只支持模态分析。在子结构分析 中，需要对整体结构中的每一个子结构进行生成和扩展，然后将各个子结构集合成整个模型。 在自由界面模态综合法中，使用MODOPT扩展模态数，应小于每个超单元求解的模态数。 若需要扩展更多的模态，需要CMS的超单元重新求解更多的阶数。以长和宽2的正方形悬臂板为例：材料参数：弹性模量E=2.1ell泊松比丫=0.3密度p=7.3e3厚度t=0.05网格划分：在长度和宽度方向上均划分20个单元(为了与后面程序一致)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, , ,!实常数，厚度TYPE,1 MAT,1REAL,1!材料及单元类型编号ESIZE,O.1,O,!变长 0.1 的单元AMESH,all!画网格DL,4, ,ALL,!施加固定约束FINISH 直接求解命令流！整体求解/SOL ANTYPE,2 MODOPT,LANB,10 EQSLVSPARMODOPT,LANB,20,0,99999999, ,OFF SOLVE表一 ANSYS直接求解结果阶数12345678910频率11.227.468.787.8100.1175.2197.8207.3229.5281.3按固定约束方向，在板长度方向的中点将方板一分为二，得到两个子结构。固定界面模态综合法求解命令流:/clear,nostart!清空变量，启动一个新分析finish/filnam,use!文件名use整体体结构/solu/prep7antype,modal!分析类型et,1,matrix50!定义超单元modopt,lanb,10!求解前10阶频率type,1solvese,part1!选择子结构一finishse,part2!选择子结构一save表二ANSYS模态综合法求解结果子结构一阶数12345678910频率287.5302.3355.3458.9624.8792.3812.7857.9879.8995.8子结构二阶数12345678910频率45.169.2131.8247.1282.1318.7407.1439.6557.5682.2整体结构阶数12345678910频率11.227.568.787.8100.1175.2197.8207.3229.5281.5/filnam,partl！文件名parti子结构一/soluantype,substr！分析类型子结构seopt,part1,2!子结构一cmsopt,fix,20!固定界面法求前20阶频率cmsel,s,part1!选定parti单元集合cmsel,s,interface! interface节点集合m,all,all!选择所有nsle solve finish save!求解并保存/filnam,part2!文件名part2子结构二/soluantype,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阶频率基本上是一致的。在计算大型复杂结构，模态综合法可以节省大 量的计算时间和计算机资源，提高效率；同时可以灵活修改大系统的子系统设计。因此，对 于复杂大型结构，如飞机、车辆、船舶、高层建筑等结构，采用ANSYS模态综合法来对结 构进行模态分析，可以在精度和计算速度上得到较好的解决方案。三.基于有限元方法的模态综合法MATLAB实现1.薄板弯曲理论等厚度薄板的坐标系如图(a)所示，板厚1/2平面，即xy平面为板的中面。从板中任取 一微元体dxdydz，在每一个面上作用的正应力和剪应力见图(b)。薄板弯曲有如下三个假设：(1)垂直于中面方向的正应变*z极微小，可以忽略。取二0，由S =- = 0得 zZ -zy)，说明板的任何一点的挠度O只与坐标x和y有关。(2)应力分量zxT 一 , 和zy远小于xy，因此可以忽略不计它们产生的正Y、剪应变zxzy，则由物理方程有-po2(1+ p)xyxy(3)薄板弯曲时，中面内各点不产生平行于中面的应变，即x z=0-xy z=0xy z=0综合上述三个假设,可用得到:-v-u=+-x-2O-2O-2O- 2O-x 2-x 2xy-2O ) 2(1 - p) -xy-2O-t L zzo dzEt 3-2O12(1 -p 2)上述应力在板内分别合成为薄板内力zo dzxy-2O-x-x简写成矩阵的形式F = DK，其中K =-2O-2O-2Ot = Lo-x-xzo dz2. 薄板弯曲单元基于薄板弯曲理论的板单元中每个节点有三个自由度：绕度，法线绕X轴转角0 ,绕xy轴转角0 y 即a 二i0二xi0yi i tdydxae = aa a t1 i n等厚度板单元中，有4个节点标系,即单元宽度为a,高度为b,在单元中心(x0,y0)建立局部坐Ox = x + ag0其中,y = y0 + bna =丄(x x )2 2 1b = 1( y y )2 2丿11x =(x + x )02 21y = 1( y + y )0221则单元内任iiw =(N+ Ni i ix ix1+ N )= LN a iy iyi i1其中, N =N ,N ,N i i ix iyn = (i+gg )(1+nn )(2+gg +nn g 2 n2)i 8iii in =1 bn (i+盘)(i+nn )(in2)ix 8 iiin =1 ag (i+盘)(i+nn )(ig 2)iy 8 iiii = i,2,3,4写成矩阵的形式W = Nae，则单元刚度矩阵可写成Ke =J BTDBdV，其中B = LN单元质量矩阵可写成Me = J p NTNdV，最终得到的单元质量矩阵和刚度矩阵都是 i2X12的矩阵。3.有限元方法求解悬臂板模态MATLAB程序%清空变量%弹性模量% 泊松比%密度%板的厚度%乂方向长度%丫方向长度%乂向节点数%丫向节点数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)*(jdy-1); % 总单元数%节点矩阵en(1:Tol_element,1:4)=0;%每个单元四个节点，记录节点编号 for ni=1:jdx-1for 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;% 位移阵constraints=1:1:jdx; % y=0 上被约束节点disp(constraints,:)=0;dof=0;%自由度标记for ni=1:jdx*jdyfor nj=1:3if 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);%整体刚度矩阵和质量矩阵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);endfor jx=1:12for jy=1:12 if(tg(jx)*tg(jy)=0) k(tg(jx),tg(jy)=k(tg(jx),tg(jy)+ek(jx,jy); m(tg(jx),tg(jy)=m(tg(jx),tg(jy)+dm(jx,jy;endendend%总装刚度阵和质量阵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);%记录角频率 toc %记时终点endendend表三MATLAB直接求解结果阶数12345678910频率11.20827.46968.76487.79799.961174.77197.92207.16229.18299.24相对误差()0.070.250.090.030.130.240.060.060.136DtEFLArSIEKT第十阶频率有很大差别，取出ANSYS分析结果的振型分析，可知第十节频率为平行于 固支边的弯曲振型，在薄板弯曲单元中不包括这个方向上的自由度，因而取ANSYS分析结 果中的第一阶频率，其大小为300.1Hz,与MATLAB结果对比，误差为0.29%畑 34 U/ia4.基于有限元的模态综合法求解悬臂板模态MATLAB程序子结构的划分与ANSYS中处理一致，主模态分析程序如下：子结构一(只列出修改部分)：子结构二(只列出修改部分):lx=2;%乂方向长度ly=l; %丫方向长度jdx=21; %*向节点数 jdy=ll; %丫向节点数%y=0和1沿x轴去掉42个点，得到总自由度Tol_dof=3*jdx*(jdy-2);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; %丫向节点数%总自由度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子结构求解结果子结构一阶数12345678910频率287.8302.1353.9455.9620.3794.8813.7852.3877.7988.8子结构二阶数12345678910频率45.169.1131.5246.4282.3318.5405.4438.8553.6681.2约束模态分析程序如下： 子结构一（只列出修改部分）lx=2; %乂方向长度ly=l; %丫方向长度 jdx=21; %*向节点数 jdy=ll; %丫向节点数 Tol_dof=3*jdx*(jdy-1);cons_Ori=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; 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);%合并总模态 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:end,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(master_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(:,sortindex); %求解特征值问题mode_number=1:5;frequency_sum(mode_number)=sqrt(nd(mode_number)/(2*pi); %得到最终模态综合方法求出频率 表五 MATLAB 模态综合法求解结果悬臂板阶数12345678910频率11.529.187.7106.8152.5260.1309.6496.1645.2818.3function k,m=km(a,b,poisson,t,E,density)%变量a,b,poisson,t,E,densit y与主程序中定义一致%kk 为刚阵列向量，(太大了，没办法，定义出来只能是向量)kk=E*t3/(360-360*poisson2)/a/b*(21-6*poisson+30*b2/a2+30*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*a-12*poisson*a-30*b2/a),E*t3/(360-360*poisson2)/a/b*(-21+6*poisson-30*b2/a2+15*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*b-12*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a-30*b2/a),E*t3/(360-360*poisson2)/a/b*(21-6*poisson-15*b2/a2T5*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*aT5*b2/a),E*t3/(360-360*poisson2)/a/b*(-21+6*poisson+15*b2/a2-30*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*aT5*b2/a);E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(8*b2-8*poisson*b2+40*a2), -30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(-3*bT2*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(-8*b2+8*poisson*b2+20*a2),0,E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(2*b2-2*poisson*b2+10*a2),0,E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(-2*b2+2*poisson*b2+20*a2),0;E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a-30*b2/a),-30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(8*a2-8*poisson*a2+40*b2),E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*a+30*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-2*a2+2*poisson*a2+20*b2),E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a+15*b2/a),0,E*t3/(360-360*poisson2)/a/b*(2*a2-2*poisson*a2+10*b2),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*aT5*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-8*a2+8*poisson*a2+20*b2); E*t3/(360-360*poisson2)/a/b*(-21+6*poisson-30*b2/a2+15*a2/b2), E*t3/(360-360*poisson2)/a/b*(-3*bT2*poisson*b+15*a2/b), E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*a+30*b2/a), E*t3/(360-360*poisson2)/a/b*(21-6*poisson+30*b2/a2+30*a2/b2), E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*a+30*b2/a),E*t3/(360-360*poisson2)/a/b*(-21+6*poisson+15*b2/a2-30*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a+15*b2/a),E*t3/(360-360*poisson2)/a/b*(21-6*poissonT5*b2/a2T5*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a+15*b2/a);E*t3/(360-360*poisson2)/a/b*(-3*bT2*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(-8*b2+8*poisson*b2+20*a2),0,E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(8*b2-8*poisson*b2+40*a2),30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(-2*b2+2*poisson*b2+20*a2),0,E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(2*b2-2*poisson*b2+10*a2),0;E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a-30*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-2*a2+2*poisson*a2+20*b2),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*a+30*b2/a),30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(8*a2-8*poisson*a2+40*b2),E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a+15*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-8*a2+8*poisson*a2+20*b2),E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*a-15*b2/a),0,E*t3/(360-360*poisson2)/a/b*(2*a2-2*poisson*a2+10*b2);E*t3/(360-360*poisson2)/a/b*(21-6*poissonT5*b2/a2T5*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a+15*b2/a),E*t3/(360-360*poisson2)/a/b*(-21+6*poisson+15*b2/a2-30*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a+15*b2/a), E*t3/(360-360*poisson2)/a/b*(21-6*poisson+30*b2/a2+30*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*bT2*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*a+30*b2/a),E*t3/(360-360*poisson2)/a/b*(-21+6*poisson-30*b2/a2+15*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*a+30*b2/a);E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(2*b2-2*poisson*b2+10*a2),0,E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(-2*b2+2*poisson*b2+20*a2),0,E*t3/(360-360*poisson2)/a/b*(-3*bT2*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(8*b2-8*poisson*b2+40*a2),-30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(-8*b2+8*poisson*b2+20*a2),0;E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*aT5*b2/a),0,E*t3/(360-360*poisson2)/a/b*(2*a2-2*poisson*a2+10*b2),E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a+15*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-8*a2+8*poisson*a2+20*b2),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*a+30*b2/a),-30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(8*a2-8*poisson*a2+40*b2),E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a-30*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-2*a2+2*poisson*a2+20*b2);E*t3/(360-360*poisson2)/a/b*(-21+6*poisson+15*b2/a2-30*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*a-15*b2/a),E*t3/(360-360*poisson2)/a/b*(21-6*poissonT5*b2/a2T5*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*aT5*b2/a),E*t3/(360-360*poisson2)/a/b*(-21+6*poisson-30*b2/a2+15*a2/b2),E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b-15*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a-30*b2/a),E*t3/(360-360*poisson2)/a/b*(21-6*poisson+30*b2/a2+30*a2/b2),E*t3/(360-360*poisson2)/a/b*(-3*bT2*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a-30*b2/a);E*t3/(360-360*poisson2)/a/b*(3*b-3*poisson*b+30*a2/b),E*t3/(360-360*poisson2)/a/b*(-2*b2+2*poisson*b2+20*a2),0,E*t3/(360-360*poisson2)/a/b*(-3*b+3*poisson*b+15*a2/b),E*t3/(360-360*poisson2)/a/b*(2*b2-2*poisson*b2+10*a2),0,E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*bT5*a2/b),E*t3/(360-360*poisson2)/a/b*(-8*b2+8*poisson*b2+20*a2),0,E*t3/(360-360*poisson2)/a/b*(-3*b-12*poisson*b-30*a2/b),E*t3/(360-360*poisson2)/a/b*(8*b2-8*poisson*b2+40*a2), 30*E*t3/(360-360*poisson2)*poisson;E*t3/(360-360*poisson2)/a/b*(3*a+12*poisson*aT5*b2/a),0,E*t3/(360-360*poisson2)/a/b*(-8*a2+8*poisson*a2+20*b2), E*t3/(360-360*poisson2)/a/b*(-3*a+3*poisson*a+15*b2/a), 0,E*t3/(360-360*poisson2)/a/b*(2*a2-2*poisson*a2+10*b2), E*t3/(360-360*poisson2)/a/b*(3*a-3*poisson*a+30*b2/a), 0,E*t3/(360-360*poisson2)/a/b*(-2*a2+2*poisson*a2+20*b2), E*t3/(360-360*poisson2)/a/b*(-3*aT2*poisson*a-30*b2/a), 30*E*t3/(360-360*poisson2)*poisson,E*t3/(360-360*poisson2)/a/b*(8*a2-8*poisson*a2+40*b2);% k将kk刚度向量用12*12的矩阵形式重新存储 k=zeros(12,12);for i=1:12for j=1:12k(i,j)=kk(12*(i-1)+j,1);endend% w为集中质量(平均分布在四个节点上)w=a*b*t*density;% 用符号变量求解形函数积分求得单元质量阵syms kx ytkxi yti real;ni=1/8*(1+kx*kxi)*(1+yt*yti)*(2+kx*kxi+yt*yti-kx2-yt2);%绕度nix=T/8*b*yti*(1+kx*kxi)*(1+yt*yti)*(1-yt2);%x轴转角niy=1/8*a*kxi*(1+kx*kxi)*(1+yt*yti)*(1-kx2);%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);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