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

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

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

3、阵；是依次释放每一边界自由度，使其得到单位位移而产生的静位移分布所组成的约束模态矩阵。它可以由对子结构直接求解一系列静平衡方程（边界自由度发生单位强迫位移）而得到：下标与分别代表该子结构的内部及界面自由度。结构作自由振动时，内部自由度上作用力为零，而界面自由度依次产生单位位移，则相当于是一个单位矩阵。于是由上式展开第一行可以得到： 从而得到约束模态矩阵，显然的列数等于界面自由度数。2.用模态矩阵对子结构的刚度阵和质量阵做第一次坐标变换，得到各个子结构的缩减刚度矩阵和质量矩阵，再形式地拼装成不耦合的形式： ；其中，一般来说，选定的前模态组成的模态矩阵，做坐标变化后得到的刚度矩阵，质量矩阵为阶矩阵

4、，远比变换前的刚度矩阵，质量矩阵的阶数小得多。3.第二次坐标变换，实现子结构的连接。假如第号子结构与第号子结构有公共界面，则他们的位移向量分别写成：，下标B表示公共界面上的约束模态广义位移。由于在界面处接口点的位移必须相等，这就导致与该位移相对应大的约束模态广义位移亦必须相等，也就有： 同样的，对于整个结构可以形式地写出，式中，是相应于所有子结构主模态的广义位移；是所有独立的约束模态广义位移的集合；是可以用表示的非独立约束模态广义位移，表示其间转换关系。利用可以进一步将未耦合的总刚度阵及总质量阵缩减（耦合）为：；最后，求解下列广义特征值问题：它的各低阶特征值与原结构低阶特征值的误差，跟子结构剖

5、分是否得当，截取模态是否合理有很大关系。如果想要得到物理坐标下的结构模态向量（振型），需要逆向做两次坐标变换。下面以悬臂板模态分析为例，分别用商用软件ANSYS和基于有限元算法的MATLAB程序，演示固定界面模态综合法的应用。二.ANSYS模态综合法求解ANSYS模态综合法是子结构在动力分析中的应用，其基本过程包括三方面：生成过程、使用过程、扩展过程。ANSYS提供了友好的模态综合法（Component Mode Synthesis）CMS向导（Preprocessor>Modeling>CMS），可以方便的定义超单元和交界面，而且可以对模态综合法分析生成的文件进行管理和组织。1.

6、创建超单元2.选择CMS求解方法(CMSOPT)：固定界面法或自由界面法。同时确定提取模态数、频率范围等。对于自由界面法还要确定刚体模态。3.命名超单元矩阵文件(SEOPT)。4.对于考虑阻尼时，输入集中质量矩阵公式。5.定义主自由度：在超单元的交界面定义主自由度。6.保持数据库：这是必须做的，因为在扩展模态时需要有相同的数据库文件。7.求解生成超单元矩阵文件。2. CMS的使用和扩展过程CMS的使用和扩展部分与子结构基本相同，但是CMS只支持模态分析。在子结构分析中，需要对整体结构中的每一个子结构进行生成和扩展，然后将各个子结构集合成整个模型。在自由界面模态综合法中，使用MODOPT扩展模态

7、数，应小于每个超单元求解的模态数。若需要扩展更多的模态，需要CMS的超单元重新求解更多的阶数。以长和宽2的正方形悬臂板为例：材料参数：弹性模量E=2.1e11 泊松比g=0.3 密度r=7.3e3 厚度t=0.05网格划分：在长度和宽度方向上均划分20个单元（为了与后面程序一致） /PREP7 BLC4, , ,2,2 ！前处理 画2×2矩形MP,DENS,1,7.3e3 ！密度MP,EX,1,2.1e11 ！弹性模量MP,PRXY,1,0.3 ！泊松比ET,1,SHELL63 ！单元类型R,1,0.05,0.05,0.05,0.05, , , ！实常数，厚度TYPE,1 MAT,1

8、REAL,1 ！材料及单元类型编号ESIZE,0.1,0, ！变长0.1的单元AMESH,all ！画网格DL,4, ,ALL, ！施加固定约束FINISH 直接求解命令流 !整体求解/SOL ANTYPE,2 MODOPT,LANB,10 EQSLV,SPAR MODOPT,LANB,20,0,99999999, ,OFF SOLVE 表一 ANSYS直接求解结果阶数12345678910频率11.227.468.787.8100.1175.2197.8207.3229.5281.3按固定约束方向，在板长度方向的中点将方板一分为二，得到两个子结构。/filnam,part1 ！文件名 par

9、t1 子结构一/soluantype,substr ！分析类型 子结构seopt,part1,2 ！子结构一cmsopt,fix,20 ！固定界面法求前20阶频率cmsel,s,part1 ！选定part1单元集合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

10、！选定part2单元集合cmsel,s,interface ！interface节点集合m,all,all ！选择所有nsle solve finish save ！求解并保存固定界面模态综合法求解命令流：/clear,nostart ！清空变量，启动一个新分析/filnam,use ！文件名 use 整体体结构/prep7et,1,matrix50 ！定义超单元type,1 se,part1 ！选择子结构一se,part2 ！选择子结构一finish /solu antype,modal ！分析类型modopt,lanb,10 ！求解前10阶频率solvefinishsave表二 ANSYS

11、模态综合法求解结果子结构一阶数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与表一中直接求解结构对比，可以看出本例中，ANSYS模态综合法具有较高的精度（保留一位小数），前10阶频率基本上是一致的。在计算大型复杂结构，模态综合法可以节省大量的计算时间和计算

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

13、应变和，则由物理方程有 （3）薄板弯曲时，中面内各点不产生平行于中面的应变，即，综合上述三个假设，可用得到： 上述应力在板内分别合成为薄板内力简写成矩阵的形式，其中2. 薄板弯曲单元基于薄板弯曲理论的板单元中每个节点有三个自由度：绕度，法线绕x轴转角，绕y轴转角，即 等厚度板单元中，有4个节点，单元宽度为a，高度为b，在单元中心(x0,y0)建立局部坐标系,即hxOyxO1234 其中， 则单元内任一点的绕度可用形函数和单元节点位移表示为： 其中，写成矩阵的形式，则单元刚度矩阵可写成，其中单元质量矩阵可写成，最终得到的单元质量矩阵和刚度矩阵都是12×12的矩阵。3. 有限元方法求解悬

14、臂板模态MATLAB程序clear all clc %清空变量tic %记时起点% 定义材料E=2.1e11; %弹性模量poisson =0.3; % 泊松比density=7.3e3; %密度t=0.05; %板的厚度lx=2; %x方向长度ly=2; %y方向长度jdx=21; %x向节点数jdy=21; %y向节点数%初始化%刚度和质量阵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)*

15、(jdy-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; endend%节点位移，约束以及自由度坐标disp(1:jdx*jdy,1:3)=1; % 位移阵constrai

16、nts=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; end endend %单元长度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

17、; % 标记，单元矩阵维数 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)+ek(jx,jy); m(tg(jx),tg(jy)=m(tg(jx),tg(jy)+

18、dm(jx,jy; end endend %总装刚度阵和质量阵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 %记时终点表三 MATLAB直接求解结果阶数12345678910频率11.20827.46968.76487.79799.961

19、174.77197.92207.16229.18299.24相对误差（%）0.070.250.090.030.130.240.060.060.136第十阶频率有很大差别，取出ANSYS分析结果的振型分析，可知第十节频率为平行于固支边的弯曲振型，在薄板弯曲单元中不包括这个方向上的自由度，因而取ANSYS分析结果中的第十一阶频率，其大小为300.1Hz，与MATLAB结果对比，误差为0.29%4. 基于有限元的模态综合法求解悬臂板模态MATLAB程序子结构的划分与ANSYS中处理一致，主模态分析程序如下：子结构一（只列出修改部分）：lx=2; %x方向长度 ly=1; %y方向长度jdx=21;

20、%x向节点数 jdy=11; %y向节点数 %y=0和1沿x轴去掉42个点，得到总自由度Tol_dof=3*jdx*(jdy-2);cons_Ori=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(add_num,master_mode);子结构二（

21、只列出修改部分）：lx=2; %x方向长度 ly=1; %y方向长度jdx=21; %x向节点数 jdy=11; %y向节点数 %总自由度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

22、子结构求解结果子结构一阶数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; %x方向长度 ly=1; %y方向长度jdx=21; %x向节点数 jdy=11; %y向节点数Tol_dof=3*jdx*(jdy-1);cons_Ori=1:1:jdx; %原有约束节点constraints=cons_Ori; %总约束 % 总刚

23、度矩阵中按是否在界面上分块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; %x方向长度 ly=1; %y方向长度j

24、dx=21; %x向节点数 jdy=11; %y向节点数 %总自由度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)*kij);%合并总模态fi_2=cat(2,fi_res2,fi_unres2);%用模态对刚阵和质量阵做第一次坐标变换k2=fi_2'*k*fi_2; m

25、2=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

26、;%非耦合模态综合矩阵，直接写成矩阵形式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_num);beta(master_mode+1+add_num:master_mode+2*add_num,master_mode+1:master_mode+add_num)=eye(add_num);b

27、eta(end-master_mode+1:end,end-master_mode+1:end)=eye(master_mode);%第二次坐标变换矩阵beta，由位移列车很容易求的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);%得到最终模态综

28、合方法求出频率表五 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,density与主程序中定义一致%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*pois

29、son*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/a2-

30、15*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*(-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

31、/(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-12*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*poiss

32、on*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*a-12*poisson*a-30*b2/a), -30*E*t3/(360-360*poisson2)*poisson, E*t3/(360-360*pois

33、son2)/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

34、*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*(-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*

35、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*a-12*poisson*a+15*b2/

36、a), E*t3/(360-360*poisson2)/a/b*(21-6*poisson-15*b2/a2-15*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*b-12*poisson*b+15*a2/b), E*t3/(360-360*poisson2)/a/b*(-8*b2+8*poisson*b2+20*a2), 0, E*t3/(

37、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*b

38、-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*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)/

39、a/b*(8*a2-8*poisson*a2+40*b2), E*t3/(360-360*poisson2)/a/b*(-3*a-12*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*poisso

40、n-15*b2/a2-15*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*(-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+1

41、5*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-30*b2/a2+15*a2/b2), E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b-15*a2/b

42、), 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*

43、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+12*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*a-3*poisson*a-15*b2

44、/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*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*

45、(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*poi

46、sson*a-15*b2/a), E*t3/(360-360*poisson2)/a/b*(21-6*poisson-15*b2/a2-15*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*(-21+6*poisson-30*b2/a2+15*a2/b2), E*t3/(360-360*poisson2)/a/b*(3*b+12*poisson*b-15

47、*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*(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