結(jié)構(gòu)系統(tǒng)動力分析可以采用總體結(jié)構(gòu)有限元法,但該方法對于VIP

1、結(jié)構(gòu)系統(tǒng)動力分析可以采用總體結(jié)構(gòu)有限元法，但該方法對于復(fù)雜大型結(jié)構(gòu)進(jìn)行分析存在計算規(guī)模大，計算時間長，所用的磁盤空間、計算機(jī)系統(tǒng)太龐大，如飛機(jī)、車輛、船舶、高層建筑等整體結(jié)構(gòu)。特別是用有限元法進(jìn)行較高頻率振動分析時，要求結(jié)構(gòu)被劃分成非常多的單元數(shù)以便獲得詳細(xì)的位移和應(yīng)力特性。這時結(jié)構(gòu)模型的節(jié)點(diǎn)自由度可能達(dá)到幾十萬甚至上百萬，直接求解如此龐大的模型是很困難。即使能夠分析，也要耗費(fèi)大量機(jī)時，效率極低。模態(tài)綜合法（Component Mode Synthesis）就是在這樣的背景下發(fā)展起來的一種縮減自由度方法。它可以將大模型化小，先進(jìn)行各個子結(jié)構(gòu)的模態(tài)分析，然后進(jìn)行模態(tài)綜合。由于僅采用了各個子結(jié)構(gòu)的

2、低階模態(tài)，因而使所建立整體結(jié)構(gòu)動力模型的自由度數(shù)大大降低，而且可以在不同的機(jī)器上對各子結(jié)構(gòu)進(jìn)行模態(tài)分析提高計算速度。模態(tài)綜合法的基本思想是根據(jù)復(fù)雜結(jié)構(gòu)的特點(diǎn)將整體結(jié)構(gòu)劃分成若干子結(jié)構(gòu)，對各個子結(jié)構(gòu)分別進(jìn)行模態(tài)分析，得到其動力特性。再利用子結(jié)構(gòu)間力平衡條件及位移協(xié)調(diào)條件將各子結(jié)構(gòu)部分低階模態(tài)特性綜合，由此得到整體結(jié)構(gòu)的動力特性。一.模態(tài)綜合法基本原理先將整個結(jié)構(gòu)劃分為s個子結(jié)構(gòu)，每個子結(jié)構(gòu)應(yīng)該是容易分析的，子結(jié)構(gòu)之間的連接盡可能要“弱”些，使子結(jié)構(gòu)之間耦合較小，每個子結(jié)構(gòu)在力學(xué)上有較大的獨(dú)立性。1.建立各個子結(jié)構(gòu)的模態(tài)矩陣，其中；是交界面完全固定時所算出的一部分（參加綜合的）固有模態(tài)，即主模態(tài)矩

3、陣；是依次釋放每一邊界自由度，使其得到單位位移而產(chǎn)生的靜位移分布所組成的約束模態(tài)矩陣。它可以由對子結(jié)構(gòu)直接求解一系列靜平衡方程（邊界自由度發(fā)生單位強(qiáng)迫位移）而得到：下標(biāo)與分別代表該子結(jié)構(gòu)的內(nèi)部及界面自由度。結(jié)構(gòu)作自由振動時，內(nèi)部自由度上作用力為零，而界面自由度依次產(chǎn)生單位位移，則相當(dāng)于是一個單位矩陣。于是由上式展開第一行可以得到： 從而得到約束模態(tài)矩陣，顯然的列數(shù)等于界面自由度數(shù)。2.用模態(tài)矩陣對子結(jié)構(gòu)的剛度陣和質(zhì)量陣做第一次坐標(biāo)變換，得到各個子結(jié)構(gòu)的縮減剛度矩陣和質(zhì)量矩陣，再形式地拼裝成不耦合的形式： ；其中，一般來說，選定的前模態(tài)組成的模態(tài)矩陣，做坐標(biāo)變化后得到的剛度矩陣，質(zhì)量矩陣為階矩陣

4、，遠(yuǎn)比變換前的剛度矩陣，質(zhì)量矩陣的階數(shù)小得多。3.第二次坐標(biāo)變換，實(shí)現(xiàn)子結(jié)構(gòu)的連接。假如第號子結(jié)構(gòu)與第號子結(jié)構(gòu)有公共界面，則他們的位移向量分別寫成：，下標(biāo)B表示公共界面上的約束模態(tài)廣義位移。由于在界面處接口點(diǎn)的位移必須相等，這就導(dǎo)致與該位移相對應(yīng)大的約束模態(tài)廣義位移亦必須相等，也就有： 同樣的，對于整個結(jié)構(gòu)可以形式地寫出，式中，是相應(yīng)于所有子結(jié)構(gòu)主模態(tài)的廣義位移；是所有獨(dú)立的約束模態(tài)廣義位移的集合；是可以用表示的非獨(dú)立約束模態(tài)廣義位移，表示其間轉(zhuǎn)換關(guān)系。利用可以進(jìn)一步將未耦合的總剛度陣及總質(zhì)量陣縮減（耦合）為：；最后，求解下列廣義特征值問題：它的各低階特征值與原結(jié)構(gòu)低階特征值的誤差，跟子結(jié)構(gòu)剖

5、分是否得當(dāng)，截取模態(tài)是否合理有很大關(guān)系。如果想要得到物理坐標(biāo)下的結(jié)構(gòu)模態(tài)向量（振型），需要逆向做兩次坐標(biāo)變換。下面以懸臂板模態(tài)分析為例，分別用商用軟件ANSYS和基于有限元算法的MATLAB程序，演示固定界面模態(tài)綜合法的應(yīng)用。二.ANSYS模態(tài)綜合法求解ANSYS模態(tài)綜合法是子結(jié)構(gòu)在動力分析中的應(yīng)用，其基本過程包括三方面：生成過程、使用過程、擴(kuò)展過程。ANSYS提供了友好的模態(tài)綜合法（Component Mode Synthesis）CMS向?qū)В≒reprocessor>Modeling>CMS），可以方便的定義超單元和交界面，而且可以對模態(tài)綜合法分析生成的文件進(jìn)行管理和組織。1.

6、創(chuàng)建超單元2.選擇CMS求解方法(CMSOPT)：固定界面法或自由界面法。同時確定提取模態(tài)數(shù)、頻率范圍等。對于自由界面法還要確定剛體模態(tài)。3.命名超單元矩陣文件(SEOPT)。4.對于考慮阻尼時，輸入集中質(zhì)量矩陣公式。5.定義主自由度：在超單元的交界面定義主自由度。6.保持?jǐn)?shù)據(jù)庫：這是必須做的，因?yàn)樵跀U(kuò)展模態(tài)時需要有相同的數(shù)據(jù)庫文件。7.求解生成超單元矩陣文件。2. CMS的使用和擴(kuò)展過程CMS的使用和擴(kuò)展部分與子結(jié)構(gòu)基本相同，但是CMS只支持模態(tài)分析。在子結(jié)構(gòu)分析中，需要對整體結(jié)構(gòu)中的每一個子結(jié)構(gòu)進(jìn)行生成和擴(kuò)展，然后將各個子結(jié)構(gòu)集合成整個模型。在自由界面模態(tài)綜合法中，使用MODOPT擴(kuò)展模態(tài)

7、數(shù)，應(yīng)小于每個超單元求解的模態(tài)數(shù)。若需要擴(kuò)展更多的模態(tài)，需要CMS的超單元重新求解更多的階數(shù)。以長和寬2的正方形懸臂板為例：材料參數(shù)：彈性模量E=2.1e11 泊松比g=0.3 密度r=7.3e3 厚度t=0.05網(wǎng)格劃分：在長度和寬度方向上均劃分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, , , ！實(shí)常數(shù)，厚度TYPE,1 MAT,1

8、REAL,1 ！材料及單元類型編號ESIZE,0.1,0, ！變長0.1的單元AMESH,all ！畫網(wǎng)格DL,4, ,ALL, ！施加固定約束FINISH 直接求解命令流 !整體求解/SOL ANTYPE,2 MODOPT,LANB,10 EQSLV,SPAR MODOPT,LANB,20,0,99999999, ,OFF SOLVE 表一 ANSYS直接求解結(jié)果階數(shù)12345678910頻率11.227.468.787.8100.1175.2197.8207.3229.5281.3按固定約束方向，在板長度方向的中點(diǎn)將方板一分為二，得到兩個子結(jié)構(gòu)。/filnam,part1 ！文件名 par

9、t1 子結(jié)構(gòu)一/soluantype,substr ！分析類型 子結(jié)構(gòu)seopt,part1,2 ！子結(jié)構(gòu)一cmsopt,fix,20 ！固定界面法求前20階頻率cmsel,s,part1 ！選定part1單元集合cmsel,s,interface ！interface節(jié)點(diǎn)集合m,all,all ！選擇所有nsle solve finish save ！求解并保存/filnam,part2 ！文件名 part2 子結(jié)構(gòu)二/soluantype,substr ！分析類型 子結(jié)構(gòu)seopt,part2,2 ！子結(jié)構(gòu)二cmsopt,fix,20 ！固定界面法求前20階頻率cmsel,s,part2

10、！選定part2單元集合cmsel,s,interface ！interface節(jié)點(diǎn)集合m,all,all ！選擇所有nsle solve finish save ！求解并保存固定界面模態(tài)綜合法求解命令流：/clear,nostart ！清空變量，啟動一個新分析/filnam,use ！文件名 use 整體體結(jié)構(gòu)/prep7et,1,matrix50 ！定義超單元type,1 se,part1 ！選擇子結(jié)構(gòu)一se,part2 ！選擇子結(jié)構(gòu)一finish /solu antype,modal ！分析類型modopt,lanb,10 ！求解前10階頻率solvefinishsave表二 ANSYS

11、模態(tài)綜合法求解結(jié)果子結(jié)構(gòu)一階數(shù)12345678910頻率287.5302.3355.3458.9624.8792.3812.7857.9879.8995.8子結(jié)構(gòu)二階數(shù)12345678910頻率45.169.2131.8247.1282.1318.7407.1439.6557.5682.2整體結(jié)構(gòu)階數(shù)12345678910頻率11.227.568.787.8100.1175.2197.8207.3229.5281.5與表一中直接求解結(jié)構(gòu)對比，可以看出本例中，ANSYS模態(tài)綜合法具有較高的精度（保留一位小數(shù)），前10階頻率基本上是一致的。在計算大型復(fù)雜結(jié)構(gòu)，模態(tài)綜合法可以節(jié)省大量的計算時間和計算

12、機(jī)資源，提高效率；同時可以靈活修改大系統(tǒng)的子系統(tǒng)設(shè)計。因此，對于復(fù)雜大型結(jié)構(gòu)，如飛機(jī)、車輛、船舶、高層建筑等結(jié)構(gòu)，采用ANSYS模態(tài)綜合法來對結(jié)構(gòu)進(jìn)行模態(tài)分析，可以在精度和計算速度上得到較好的解決方案。三.基于有限元方法的模態(tài)綜合法MATLAB實(shí)現(xiàn)1. 薄板彎曲理論等厚度薄板的坐標(biāo)系如圖(a)所示，板厚12平面，即xy平面為板的中面。從板中任取一微元體dxdydz，在每一個面上作用的正應(yīng)力和剪應(yīng)力見圖(b)。薄板彎曲有如下三個假設(shè)：（1）垂直于中面方向的正應(yīng)變極微小，可以忽略。取，由得，說明板的任何一點(diǎn)的撓度只與坐標(biāo)x和y有關(guān)。（2）應(yīng)力分量、和遠(yuǎn)小于、和，因此可以忽略不計它們產(chǎn)生的正應(yīng)變、剪

13、應(yīng)變和，則由物理方程有 （3）薄板彎曲時，中面內(nèi)各點(diǎn)不產(chǎn)生平行于中面的應(yīng)變，即，綜合上述三個假設(shè)，可用得到： 上述應(yīng)力在板內(nèi)分別合成為薄板內(nèi)力簡寫成矩陣的形式，其中2. 薄板彎曲單元基于薄板彎曲理論的板單元中每個節(jié)點(diǎn)有三個自由度：繞度，法線繞x軸轉(zhuǎn)角，繞y軸轉(zhuǎn)角，即 等厚度板單元中，有4個節(jié)點(diǎn)，單元寬度為a，高度為b，在單元中心(x0,y0)建立局部坐標(biāo)系,即hxOyxO1234 其中， 則單元內(nèi)任一點(diǎn)的繞度可用形函數(shù)和單元節(jié)點(diǎn)位移表示為： 其中，寫成矩陣的形式，則單元剛度矩陣可寫成，其中單元質(zhì)量矩陣可寫成，最終得到的單元質(zhì)量矩陣和剛度矩陣都是12×12的矩陣。3. 有限元方法求解懸

14、臂板模態(tài)MATLAB程序clear all clc %清空變量tic %記時起點(diǎn)% 定義材料E=2.1e11; %彈性模量poisson =0.3; % 泊松比density=7.3e3; %密度t=0.05; %板的厚度lx=2; %x方向長度ly=2; %y方向長度jdx=21; %x向節(jié)點(diǎn)數(shù)jdy=21; %y向節(jié)點(diǎn)數(shù)%初始化%剛度和質(zhì)量陣Tol_dof=3*jdx*（jdy-1）; %在y=0邊界上施加固定約束，去掉21個節(jié)點(diǎn)k(1:Tol_dof,1:Tol_dof)=0; %總剛陣m(1:Tol_dof,1:Tol_dof)=0; %總質(zhì)量陣Tol_element=(jdx-1)*

15、(jdy-1); % 總單元數(shù)%節(jié)點(diǎn)矩陣en(1:Tol_element,1:4)=0; %每個單元四個節(jié)點(diǎn)，記錄節(jié)點(diǎn)編號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%節(jié)點(diǎn)位移，約束以及自由度坐標(biāo)disp(1:jdx*jdy,1:3)=1; % 位移陣constrai

16、nts=1:1:jdx; % y=0上被約束節(jié)點(diǎn)disp(constraints,:)=0;dof=0; %自由度標(biāo)記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); %單元高度%四節(jié)點(diǎn)12自由度 單元剛度和質(zhì)量矩陣% km函數(shù)，輸入材料參數(shù)和節(jié)點(diǎn)集中長度ek,dm=km(el/2,eh/2,poisson,t,E,density); %整體剛度矩陣和質(zhì)量矩陣 tg(1:12)=0

17、; % 標(biāo)記，單元矩陣維數(shù) for loopi=1:(jdx-1)*(jdy-1) %外循環(huán)，單元 for zi=1:4 %每個節(jié)點(diǎn)三個自由度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 %總裝剛度陣和質(zhì)量陣end%求解特征值問題v,d = eig(k,m); %eig函數(shù)求解特征值問題tempd=diag(d); %記錄特征值nd,sortindex=sort(tempd); %對特征值排序v=v(:,sortindex); %特征向量排序mode_number=1:15; %求解階數(shù)frequency(mode_number)=sqrt(nd(mode_number)/(2*pi);%記錄角頻率 toc %記時終點(diǎn)表三 MATLAB直接求解結(jié)果階數(shù)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分析結(jié)果的振型分析，可知第十節(jié)頻率為平行于固支邊的彎曲振型，在薄板彎曲單元中不包括這個方向上的自由度，因而取ANSYS分析結(jié)果中的第十一階頻率，其大小為300.1Hz，與MATLAB結(jié)果對比，誤差為0.29%4. 基于有限元的模態(tài)綜合法求解懸臂板模態(tài)MATLAB程序子結(jié)構(gòu)的劃分與ANSYS中處理一致，主模態(tài)分析程序如下：子結(jié)構(gòu)一（只列出修改部分）：lx=2; %x方向長度 ly=1; %y方向長度jdx=21;

20、%x向節(jié)點(diǎn)數(shù) jdy=11; %y向節(jié)點(diǎn)數(shù) %y=0和1沿x軸去掉42個點(diǎn)，得到總自由度Tol_dof=3*jdx*(jdy-2);cons_Ori=1:1:jdx; %原有約束節(jié)點(diǎn)% 固定界面法人為施加約束cons_Arti=jdx*(jdy-1)+1:1:jdx*jdy; %總約束constraints=cons_Ori,cons_Arti; master_mode=15;add_num=jdx*3;% 取前master_mode階振型作為模態(tài)綜合fi_unres1=cat(1,v(1:master_mode,:)',zeros(add_num,master_mode);子結(jié)構(gòu)二（

21、只列出修改部分）：lx=2; %x方向長度 ly=1; %y方向長度jdx=21; %x向節(jié)點(diǎn)數(shù) jdy=11; %y向節(jié)點(diǎn)數(shù) %總自由度Tol_dof=3*jdx*(jdy-1);%固定界面法人為施加約束cons_Arti=1:1:jdx; % 單邊約束節(jié)點(diǎn)編號constraints=cons_Arti; %總約束master_mode=15;add_num=jdx*3;fi_unres2=cat(1,zeros(add_num,master_mode),v(1:master_mode,:)');% 主模態(tài)中，認(rèn)為施加固定約束的節(jié)點(diǎn)，其位移為零，補(bǔ)充完整加入模態(tài)綜合表四 MATLAB

22、子結(jié)構(gòu)求解結(jié)果子結(jié)構(gòu)一階數(shù)12345678910頻率287.8302.1353.9455.9620.3794.8813.7852.3877.7988.8子結(jié)構(gòu)二階數(shù)12345678910頻率45.169.1131.5246.4282.3318.5405.4438.8553.6681.2約束模態(tài)分析程序如下：子結(jié)構(gòu)一（只列出修改部分）：lx=2; %x方向長度 ly=1; %y方向長度jdx=21; %x向節(jié)點(diǎn)數(shù) jdy=11; %y向節(jié)點(diǎn)數(shù)Tol_dof=3*jdx*(jdy-1);cons_Ori=1:1:jdx; %原有約束節(jié)點(diǎn)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);%提取約束模態(tài)fi_res1=cat(1,-inv(kii)*kij,eye(add_num);%合并總模態(tài)fi_1=cat(2,fi_unres1,fi_res1);%用模態(tài)對剛陣和質(zhì)量陣做第一次坐標(biāo)變換k1=fi_1'*k*fi_1; m1=fi_1'*m*fi_1;子結(jié)構(gòu)二（只列出修改部分）：lx=2; %x方向長度 ly=1; %y方向長度j

24、dx=21; %x向節(jié)點(diǎn)數(shù) jdy=11; %y向節(jié)點(diǎn)數(shù) %總自由度Tol_dof=3*jdx*jdy; %子結(jié)構(gòu)沒有約束% 總剛度矩陣中按是否在界面上分塊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);%提取約束模態(tài)fi_res2=cat(1,eye(add_num),-inv(kii)*kij);%合并總模態(tài)fi_2=cat(2,fi_res2,fi_unres2);%用模態(tài)對剛陣和質(zhì)量陣做第一次坐標(biāo)變換k2=fi_2'*k*fi_2; m

25、2=fi_2'*m*fi_2;模態(tài)綜合程序如下：master_mode=15;add_num=jdx*3;Tol_num=length(k1); %master_mode為主模態(tài)保留階數(shù)，add_num為界面上自由度數(shù)，Tol_num為變換后剛度矩陣維數(shù)K=zeros(2*Tol_num);M=zeros(2*Tol_num); %初始化模態(tài)綜合矩陣 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、;%非耦合模態(tài)綜合矩陣，直接寫成矩陣形式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);%第二次坐標(biāo)變換矩陣beta，由位移列車很容易求的KK=beta'*K*beta;MM=beta'*M*beta;%做第二次坐標(biāo)變換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);%得到最終模態(tài)綜

28、合方法求出頻率表五 MATLAB模態(tài)綜合法求解結(jié)果懸臂板階數(shù)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