桥梁结构振动与控制分析研究结题报告

1、.第五期SRTP结题报告 桥梁结构振动与控制分析研究第五期SRTP结题报告桥梁结构振动与控制分析研究指导老师：沈火明 教授项目组成员：张 徐 20074968 李 恒 20074977 张冰清 20085379;桥梁结构振动与控制分析研究一、课题的研究意义及研究方法1.1 课题的研究意义桥梁结构的振动是引起桥梁损坏(破坏)的一个重要因素，引起桥梁振动的因素主要有：地震引起的振动、荷载引起的振动及车 - 桥耦合作用引起的振动。传统的结构强度设计方法通过增强结构物自身抗力来抵御地震作用，即由结构本身储存和消耗地震能量。但由于人类测震技术的不成熟，尚不能准确估计振动的强度和特性。因此，可能会出现结构

2、不满足安全性的要求而产生安全事故。近年来发展起来的结构控制技术是建筑结构抗震领域内的一个新的研究热点, 它是通过采用结构振动控制的理论与方法改变结构系统的动力学性能或阻尼耗散性能来增加和改善结构的抗震能力，是一种积极主动的对策。因此，近年来桥梁结构的振动控制倍受学术界、工程界的广泛关注, 并获得了长足的进步。结构振动与控制的研究与应用有着广泛的前景，它的研究和发展将给结构工程抗震设计带来一张革命，其巨大的经济效益和社会效益已得到证明。1.2 本文的研究思路和方法本文以竖向弯曲振动时桥梁跨中挠度的振动幅度为控制目标，通过分析安装TMD前后桥梁跨中挠度的振动幅度变化量来讨论TMD对桥梁振动的控制效

3、果，并探索TMD的参数优化。 对于TMD控制下的车桥耦合系统的振动，已被采用的数值研究方法有两类：一类在建立系统耦合方程组的基础上，借助编程语言或数值计算软件（如MATLAB和VB等），利用数值积分方法编程求解耦合方程组；另一类，借助仿真分析软件（如有限元软件和Simulink等）实现对系统的仿真分析。本文将采取第一类方法，编程求解方程组。针对简支梁桥在列车匀速通过时的竖向弯曲振动，本文先建立车桥耦合振动理论模型，利用数值计算方法结合MATLAB软件，编程求解车桥时变系统振动微分方程组，获得列车过桥时桥梁竖向振动位移响应；再建立车桥-TMD耦合振动理论模型，求解获得单个以及多个TMD控制下的桥

4、梁竖向位移响应，分析TMD的控制效果，并讨论TMD参数优化对控制效果的影响。1、 理论模型及求解方法2.1 车辆-简支梁桥竖向振动模型2.1.1 模型简化模型中，车体、转向架、车轮均被认为是刚体，相互间通过弹簧阻尼系统连接；不考虑车轮与钢轨表面粗糙，认为钢轨固结于桥梁上作为桥梁的一部分，不考虑钢轨及轨下结构局部变形造成的影响，于是轮轴的竖向位移等于轮轨接触点（即车轮与桥梁接触点）的竖向位移。由于仅考虑系统的竖向振动，本文中车辆简化为二系弹簧悬挂系统，仅考虑车体沉浮、点头，前后转向架构架沉浮运动，每节车辆四个自由度。桥梁采用简支欧拉梁模型。系统简化模型及坐标系建立如图1所示。图1 系统简化计算模

5、型图中各物理量意义如下：v -行车速度； yc-车体中心竖向沉浮位移；i-第i节车体点头位移； Jc-车体点头惯量；mc-车体质量； cs2-二系悬挂阻尼；ks2-二系悬挂刚度； mt-转向架构架与轮对质量之和；yt-构架中心竖向沉浮位移； cs1-一系悬挂阻尼；ks1-一系悬挂刚度； a-同一节车转向架中心距离；EI-桥梁抗弯刚度； d-前后两节车相邻轮对间距；L-桥梁全长； w-桥梁挠度，以水平位置为坐标起点；mk-TMD质量； kk-TMD刚度；Ck-TMD阻尼； yz-TMD竖向位移；li-第i个轮对与第一个轮对之间的距离；2.1.2 车辆系统振动微分方程第i节车振动方程车体沉浮运动：

6、 (2-1)车体点头运动： （2-2）构架沉浮运动： (2-3) (2-4)对车厢整体可得 （2-5） （2-6）其中p2i-1，p2i分别为第i节车厢前后两个车轮与桥梁之间的作用力。2.1.3 简支梁桥振动方程本为采用简支欧拉梁模型，不考虑桥梁阻尼时，振动方程为：（2-7) 其中，A-桥梁单位长度的质量； F(x,t)-t时刻x处作用在梁上的外力，包括桥梁自重和轮轨相互作用力，即 （2-8）其中，N表示车辆节数； 表示dirac函数； 将式（2-8）代入式（2-7）得， （2-9）为求解方程（2-9），利用分离变量法设 （2-10）其中为简支梁的振型函数，n为模态截断数，Tj(t)形态振幅函

7、数。将式（2-10）代入式（2-9），各项自0到L积分，利用振型函数的正交性与dirac函数的性质，并令,得j=1,2,.,n (2-11)将式（2-5）与式（2-6）代入式(2-11),消去变量Pi，得 j=1,2,.,n (2-12)为简化表达，令，整理得， j=1,2,.,n (2-13)这样，式（2-1）式（2-4）与式(2-13)一起组成车桥时变系统耦合振动微分方程组，方程组共有（4N+n)个方程，以yci,yt2i-1,yt2i,i,Tj(t)共（4N+n）个未知量为求解变量，结合初值条件，可利用动力学连续数值积分方法联合求解。2.2 车桥TMD耦合模型（以跨中悬挂单个TMD为例）

8、设TMD质量为mk，弹簧刚度为kk，阻尼为ck，悬挂位置为跨中，以静平衡位置为坐标起点，振动位移为yz，则TMD的运动方程为 （2-14）简支梁跨中悬挂TMD时，车辆系统振动方程不受影响，简支梁所受外力应考虑梁与TMD之间的相互作用，即通过TMD弹簧与阻尼器传递的力，于是，式（2-8）变为 （2-15）其中，为dirac函数。此时，式（2-13）变为 j=1,2,.,n （2-16）这样，式（2-1）式（2-4）与式(2-14)、式（2-16）一起组成车桥TMD时变系统耦合振动微分方程组，方程组共有（4m+n+1）方程，以yci,yt2i-1,yt2i,i,Tj(t),yz共（4m+n+1）个

9、未知量为求解变量，结合初值条件，可利用动力学连续数值积分方法联合求解。2.3 求解步骤与方法（以单个TMD控制系统为例）车桥TMD系统耦合方程组写成矩阵形式为： （2-17）其中，x-以yci,yt2i-1,yt2i,i,Tj(t),yz 构成的未知向量，即分块表示为：； M,C,K-总体等效质量,阻尼,刚度矩阵; f -等效载荷向量；2.3.1 系数矩阵a)质量矩阵分块表示其中, ； ,I表示单位矩阵； mz=mk; ,; ;其中，限于篇幅，hij(t)简写为hij.b)阻尼矩阵分块表示为 ,; ; ;， c)刚度矩阵分块表示为 ,; ; ;，d)载荷向量分块表示为;,.2.3.2求解方法本

10、文采用Newmark-法，利用MATLAB软件编程求解矩阵方程组（2-17），获得系统的位移响应。三、实际算例 根据（二）中建立的计算模型与求解方法，本文给出以下实际算例。实例研究一辆10节编组的列车通过一座简支梁桥时，引起的桥梁振动，并利用TMD控制桥梁振动。列车采用德国ICE动车和拖车，前后2节动车中间8节拖车编组，桥梁为一全长32m简支梁桥。具体车辆与桥梁参数见附录。3.1 桥梁静挠度根据材料力学中简支梁在均布压力下跨中挠度的计算公式，本文算例中的桥梁在自重下跨中挠度为：3.2 实施控制前，列车过桥时引起的跨中振动响应 如图2，给出了列车以100km/h车速匀速通过时，简支梁跨中挠度的时

11、程曲线。t=0s时，列车开始上桥，约9s时，列车刚好完全离开桥梁。图2 100km/h车速下跨中振动响应 从图中可以看出，当列车以100km/h的速度通过时，跨中最大挠度发生在列车刚上桥后，为60.4mm，相对于跨中静挠度27.9mm超出一倍以上，振幅为32.5mm，这对于桥梁结构本身和列车的安全都是有害的也是危险的，因此需要对其进行控制。图中，t=9s后列车离开桥梁，跨中位移并没有衰减，而是继续作以静挠度为均值，以列车刚离开桥瞬间振幅为幅值的简谐振动。这是因为算例中的简支梁没有考虑桥梁本身的阻尼，当列车离开桥梁后，桥梁自由振动，此后跨中响应仅决定于其初值条件即车离开桥时跨中的振动情况，由于不

12、存在阻尼这种振动将持续不会衰减，但这仅是算例中假设前提下的情况。为研究桥梁跨中位置在列车通过时相对于静平衡位置的振动，将图3.1中的时程曲线沿纵坐标向下平移27.9mm，得到图3如下：图3 跨中振幅时程曲线3.3 实施控制前，跨中振幅随车速的变化情况本为计算了列车以60 km/h-200 km/h不同车速通过桥梁时，跨中振幅的变化情况，如图4 所示。图4 桥梁的振幅随车速的变化曲线从图中可以看出，随着列车过桥车速的增加，简支梁跨中振幅基本呈增加趋势，个别车速下振幅出现局部峰值，车速为185km/h时振幅甚至接近42mm。3.4 利用单个TMD控制桥梁在列车通过时的振动本文以车速为100km/h

13、为例，讨论单个TMD对桥梁振动的控制作用。计算中的TMD质量、阻尼和刚度参数由Den Hartog 参数调整公式（式3-1）给出。， （3-1）其中，为TMD质量比，cz和cc分别为TMD的阻尼系数和临界阻尼系数。如图5所示，给出了车速100km/h，跨中悬挂质量比为0.08%的单个TMD时，桥梁跨中振幅的时程曲线。图5 单个TMD控制下的跨中振动曲线从图中可以看出，跨中悬挂质量比为0.08%的单个TMD时，跨中振幅为30.6mm。相对于控制前的32.5mm振幅，控制效果为5.85%。若改变TMD质量比，将获得不同的控制效果。于是，以控制效果最大为控制目标，可获得最佳TMD质量比参数。同时，图

14、中看出当列车离开桥梁后，跨中位移呈现向静平衡位置衰减趋势，这是由于模型中考虑了TMD阻尼的作用。3.5 单个TMD控制的最佳质量比如图6所示，本文给出了车速为100km/h时，单个TMD的控制效果随TMD质量比变化曲线。图6 单个TMD的控制效果随TMD质量比变化曲线从上图可以看出质量比=0.05%时，获得最佳控制效果为24.3%，此时跨中振幅为24.6mm，最大挠度52.5mm，该TMD减振效果显著。3.6 实现MTMD对桥梁振动的控制MTMD形式多样，涉及的参数也多如：TMD的悬挂位置，各TMD的质量，阻尼，刚度，TMD相互间的频率间隔等。本文算例仅讨论频率呈线性分布，等间距悬挂的5个TM

15、D的控制作用。图7给出了等质量比，频率间线性分布的5个TMD等间距悬挂时，跨中振动的时程曲线，其中，MTMD参数为质量比=0.001%，频率间隔df=3Hz。图7 MTMD控制下的跨中振幅曲线从上图可以看出，5个TMD控制时，跨中振幅为24.2mm，控制效果为25.54%，相对于最佳质量比的单个TMD控制效果虽然只增加了1.24%，但总质量却只有单个TMD的1/10，这非常有利于降低悬挂TMD对桥梁静挠度的影响。同时，如果调节MTMD的参数，可以获得不同的控制效果。于是，以控制效果最大为目标，以质量比和频率间隔为优化变量，根据优化理论，利用二维优化搜索方法可以获得MTMD的最佳参数。 由于涉及

16、的程序较大，本文受条件限制未能利用该方法研究MTMD的最佳参数，只是给出了保持总质量比的前提下控制效果随频率间隔的变化情况（如图8）。图8 MTMD控制效果随频率间隔的变化曲线从上图可以看出，控制总质量为0.05%，当频率间隔较低时，随着df增加，控制效果呈增加趋势但变化不明显；当df=11.5Hz时，控制效果取最大值；当频率间隔高于11.5HZ时，控制效果明显迅速衰减。所以取频率间隔为11.5Hz，可以获得最佳控制效果为26.1%。3.7 算例结论3.13.6的算例中，本文运用（二）中建立的模型以及求解方法，实现了TMD和MTMD对简支梁桥在列车匀速通过时的振动响应的控制。四、总结 本文讨论

17、了列车过桥时引起的车桥耦合振动，建立了车辆简支梁桥与车辆简支梁桥TMD耦合振动微分方程，推导了方程的求解过程。讨论了TMD对桥梁振动的控制作用，利用数值计算研究了TMD的参数优化对控制效果的影响。在建立理论模型的基础上，运用该模型计算了实际车桥耦合振动响应，并实现了TMD对桥梁的控制作用。但本文中的模型还十分粗糙，与实际情况差距较大，若要指导实际工程应用尚需更加符合实际的细化模型。同时，本文没有验证数值计算结果的可靠性。参考文献1肖盐平.桥梁振动控制的初步研究.西南交通大学研究生学位论文，2005，12.2王均刚,马汝建.TMD 振动控制结构的发展及应用.济南大学学报(自然科学版),2006，

18、04.3杨宜谦,张煅,周宏业,孙宁.用调频质量阻尼器抑制铁路桥梁竖向共振的研究.中国铁道科学,1998 年3 月.4郭文华, 路萍.TMD对高速列车通过简支箱梁桥时的振动控制研究.振动与冲击,2008,12.5许磊平,刘伟庆等.TMD 抑制桥梁振动仿真分析.防灾减灾工程学报.2009,6.6黄维平,强士中.大跨度悬索桥的双向TMD振动控制.桥梁建设,2000年第4期7张红艳,白长青,许庆余.多自由度复杂结构的TMD调谐减震控制研究.应用力学学报,2008,12.8张晶,王志强.调谐质量阻尼器(TMD) 在大跨斜拉桥减震控制中的应用.工业建筑,20079肖艳平,沈火明,叶献辉.两种桥梁振动控制方

19、法的对比分析.噪声与振动控制,2006年8 月10李小珍,强士中.列车-桥梁耦合振动研究的现状与发展趋势.铁道学报，2002年10月。11郝超，强士中.移动荷载作用下桥梁的振动控制.国外桥梁，1999.12左一舟.关于列车-轨道( 桥梁) 时变系统空间振动方程的建立及其求解.铁道科学与工程学报,2005年2月13晋智斌,强士中,李小珍.高速列车-桥梁竖向随机振动的时域分析方法.地震工程与工程振动，2008,06.附录一：实际算例中车桥参数表桥梁参数全长L=32m 弯曲刚度EI=5.18e10 单位长度质量m=1.08e4车辆参数表动车参数拖车参数车体质量/kg5.88e44.55e4构架质量/

20、kg6.45e33.87e3车体的点头刚度/kg*m23.089e62.391e6二系阻尼系数/kg/s9e40.292e5一系阻尼系数/kg/s3e4 0.219e5二系弹簧刚度/N/m1.52e60.324e6一系弹簧刚度/N/m3.418e62.82e7轮对间隔/m11.46 17附录二：数值计算程序代码%用纽马克法计算结构的动力学运动方程，返回值为简支梁跨中挠度%单个TMD控制function w_mid=bz2t1(vv,t,mu)global L v ax1 ax2 dax1 dax2;rn=5;n=rn; %所考虑桥梁振型的前n阶nn=n+41;g=9.81; %重力加速度-m/

21、s2EI=5.18e10;m=10.8e3;L=32;v=vv/3.6; %列车移动荷载的速度-km/h(->m/s)Ms1=5.88e4; %车体质量-kgMp1=6.45e3; %构架质量（含轮对质量）-kgIs1=3.089e6; %车体的点头刚度-kg*m2Cs1=0.9e5; %二系阻尼系数-kg/sCp1=0.3e5; %一系阻尼系数-kg/sKs1=1.52e6; %二系弹簧刚度-N/mKp1=2.418e6; %一系弹簧刚度-N/max1=11.46; %轮对间隔- mMs2=4.55e4; %车体质量-kgMp2=3.87e3; %构架质量（含轮对质量）-kgIs2=2

22、.391e6; %车体的点头刚度-kg*m2Cs2=0.292e5; %二系阻尼系数-kg/sCp2=0.219e5; %一系阻尼系数-kg/sKs2=0.324e6; %二系弹簧刚度-N/mKp2=28.2e6; %一系弹簧刚度-N/max2=17; %轮对间隔- mdax1=3.67+3.5;dax2=3.67*2;c1=EI/m/L4*pi4; %桥梁自由振动频率-1/sc2=pi/L; %c2*v 车辆行进谐振频率-1/sc3=2*Mp1/m/L;c4=Ms1/m/L;c5=2*Is1/m/L/ax1;c6=(2*Mp1+Ms1)*g/m/L;c7=2*Mp2/m/L;c8=Ms2/m

23、/L;c9=2*Is2/m/L/ax2;c10=(2*Mp2+Ms2)*g/m/L;l=zeros(20,1); %第i个轮对与以一个轮对之间的距离li-ml(2)=ax1;l(3)=ax1+dax1;for j=4:18 if rem(j,2)=0 %判断奇偶 l(j)=l(j-1)+ax2; else l(j)=l(j-1)+dax2; endendl(19)=l(18)+dax1;l(20)=l(19)+ax1;k0=zeros(1,n);nt=5;f0=sqrt(c1); %基频Mk=mu*m*L;Kk=m*L*mu*c1/(1+mu)2;Ck=2*Mk*sqrt(3*mu/8/(1+

24、mu)3)*f0;m10=ones(1,n);m11=diag(m10);m21=zeros(41,n);m20=Mp1,Mp1,Ms1,Is1,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp2,Mp2,Ms2,Is2,. Mp1,Mp1,Ms1,Is1,Mk;m22=diag(m20);c11=zeros(n,n);c12=zeros(n,41);c221=Cp1+Cs1,0,-C

25、s1,0.5*ax1*Cs1; 0,Cp1+Cs1,-Cs1,-0.5*ax1*Cs1; -Cs1,-Cs1,2*Cs1,0; 0.5*ax1*Cs1,-0.5*ax1*Cs1,0,0.5*ax1*ax1*Cs1;c222=Cp2+Cs2,0,-Cs2,0.5*ax2*Cs2; 0,Cp2+Cs2,-Cs2,-0.5*ax2*Cs2; -Cs2,-Cs2,2*Cs2,0; 0.5*ax2*Cs2,-0.5*ax2*Cs2,0,0.5*ax2*ax2*Cs2;c220=zeros(4,4);c40=zeros(4,1);c42=zeros(1,40);c22=c221,c220, c220, c

26、220, c220, c220, c220, c220, c220, c220,c40; c220,c222, c220, c220, c220, c220, c220, c220, c220, c220,c40; c220, c220,c222, c220, c220, c220, c220, c220, c220, c220,c40; c220, c220, c220,c222, c220, c220, c220, c220, c220, c220,c40; c220, c220, c220, c220,c222, c220, c220, c220, c220, c220,c40; c22

27、0, c220, c220, c220, c220,c222, c220, c220, c220, c220,c40; c220, c220, c220, c220, c220, c220,c222, c220, c220, c220,c40; c220, c220, c220, c220, c220, c220, c220,c222, c220, c220,c40; c220, c220, c220, c220, c220, c220, c220, c220,c222, c220,c40; c220, c220, c220, c220, c220, c220, c220, c220, c22

28、0,c221,c40; c42,Ck;k11=zeros(n,n);ii=1:n;Pkk=sin(0.5*ii*pi);k11=(ii.4*c1);k11=diag(k11);k12=zeros(n,41);k221=Kp1+Ks1,0,-Ks1,0.5*ax1*Ks1; 0,Kp1+Ks1,-Ks1,-0.5*ax1*Ks1; -Ks1,-Ks1,2*Ks1,0; 0.5*ax1*Ks1,-0.5*ax1*Ks1,0,0.5*ax1*ax1*Ks1;k222=Kp2+Ks2,0,-Ks2,0.5*ax2*Ks2; 0,Kp2+Ks2,-Ks2,-0.5*ax2*Ks2; -Ks2,-Ks2,

29、2*Ks2,0; 0.5*ax2*Ks2,-0.5*ax2*Ks2,0,0.5*ax2*ax2*Ks2;k220=zeros(4,4);k40=zeros(4,1);k42=zeros(1,40);k22=k221,k220, k220, k220, k220, k220, k220, k220, k220, k220,k40; k220,k222, k220, k220, k220, k220, k220, k220, k220, k220,k40; k220, k220,k222, k220, k220, k220, k220, k220, k220, k220,k40; k220, k2

30、20, k220,k222, k220, k220, k220, k220, k220, k220,k40; k220, k220, k220, k220,k222, k220, k220, k220, k220, k220,k40; k220, k220, k220, k220, k220,k222, k220, k220, k220, k220,k40; k220, k220, k220, k220, k220,k220, k222, k220, k220, k220,k40; k220, k220, k220, k220, k220, k220, k220,k222, k220, k22

31、0,k40; k220, k220, k220, k220, k220,k220, k220, k220, k222, k220,k40; k220, k220, k220, k220, k220, k220, k220, k220, k220,k221,k40; k42,Kk;Y0=zeros(nn,1); %位移初始条件Y10=zeros(nn,1); %速度初始条件Y20=zeros(nn,1); %加速度初始条件X_mid=sin(pi/2),sin(pi),sin(3*pi/2),sin(2*pi),sin(5*pi/2); %前5阶振型在跨中取值if t=0 returnelse

32、dt=0.005; %积分步长 nm=t/dt; Gama=0.50; Beta=0.25; a0=1/Beta/dt2; a1=Gama/Beta/dt; a2=1/Beta/dt; a3=0.5/Beta-1; a4=Gama/Beta-1; a5=dt/2*(Gama/Beta-2); a6=(1-Gama)*dt; a7=Gama*dt; w_mid=zeros(nm+1,1); %跨中挠度 diracl=zeros(1,20);for k=0:nm t=k*dt;for m=1:20 %diracli diracl(m)=diraci(l(m)/v,(L+l(m)/v,t);endi

33、i=1:n;mt=2*Mk/m/L*sin(ii*pi/2);ctt=-Ck*sin(ii*pi/2);ktt=-Kk*sin(ii*pi/2);P011=sin(ii*v*t*c2);P022=sin(ii*(v*t-ax1)*c2);P033=sin(ii*(v*t-ax1-dax1)*c2);P044=sin(ii*(v*t-ax1-dax1-ax2)*c2);P055=sin(ii*(v*t-ax1-dax1-ax2-dax2)*c2);P066=sin(ii*(v*t-ax1-dax1-ax2-dax2-ax2)*c2);P077=sin(ii*(v*t-ax1-dax1-ax2-d

34、ax2-ax2-dax2)*c2);P088=sin(ii*(v*t-ax1-dax1-ax2-dax2-ax2-dax2-ax2)*c2);P099=sin(ii*(v*t-ax1-dax1-ax2-dax2-ax2-dax2-ax2-dax2)*c2);P01010=sin(ii*(v*t-ax1-dax1-ax2-dax2-ax2-dax2-ax2-dax2-ax2)*c2);P01111=sin(ii*(v*t-ax1-dax1-4*ax2-4*dax2)*c2);P01212=sin(ii*(v*t-ax1-dax1-5*ax2-4*dax2)*c2);P01313=sin(ii*(

35、v*t-ax1-dax1-5*ax2-5*dax2)*c2);P01414=sin(ii*(v*t-ax1-dax1-6*ax2-5*ax2)*c2);P01515=sin(ii*(v*t-ax1-dax1-6*ax2-6*dax2)*c2);P01616=sin(ii*(v*t-ax1-dax1-7*ax2-6*dax2)*c2);P01717=sin(ii*(v*t-ax1-dax1-7*ax2-7*dax2)*c2);P01818=sin(ii*(v*t-ax1-dax1-8*ax2-7*dax2)*c2);P01919=sin(ii*(v*t-ax1-2*dax1-8*ax2-7*da

36、x2)*c2);P02020=sin(ii*(v*t-2*ax1-2*dax1-8*ax2-7*dax2)*c2);m12=-c3*diracl(1)*P011',-c3*diracl(2)*P022',-c4*(diracl(1)*P011'+diracl(2)*P022'),c5*(diracl(1)*P011'-diracl(2)*P022'),. -c7*diracl(3)*P033',-c7*diracl(4)*P044',-c8*(diracl(3)*P033'+diracl(4)*P044'),c9*

37、(diracl(3)*P033'-diracl(4)*P044'),. -c7*diracl(5)*P055',-c7*diracl(6)*P066',-c8*(diracl(5)*P055'+diracl(6)*P066'),c9*(diracl(5)*P055'-diracl(6)*P066'),. -c7*diracl(7)*P077',-c7*diracl(8)*P088',-c8*(diracl(7)*P077'+diracl(8)*P088'),c9*(diracl(7)*P077

38、9;-diracl(8)*P088'),. -c7*diracl(9)*P099',-c7*diracl(10)*P01010',-c8*(diracl(9)*P099'+diracl(10)*P01010'),c9*(diracl(9)*P099'-diracl(10)*P01010'),. -c7*diracl(11)*P01111',-c7*diracl(12)*P01212',-c8*(diracl(11)*P01111'+diracl(12)*P01212'),c9*(diracl(11)*P01

39、111'-diracl(12)*P01212'),. -c7*diracl(13)*P01313',-c7*diracl(14)*P01414',-c8*(diracl(13)*P01313'+diracl(14)*P01414'),c9*(diracl(13)*P01313'-diracl(14)*P01414'),. -c7*diracl(15)*P01515',-c7*diracl(16)*P01616',-c8*(diracl(15)*P01515'+diracl(16)*P01616'),

40、c9*(diracl(15)*P01515'-diracl(16)*P01616'),. -c7*diracl(17)*P01717',-c7*diracl(18)*P01818',-c8*(diracl(17)*P01717'+diracl(18)*P01818'),c9*(diracl(17)*P01717'-diracl(18)*P01818'),. -c3*diracl(19)*P01919',-c3*diracl(20)*P02020',-c4*(diracl(19)*P01919'+diracl(20)*P02020'),c5*(diracl(19)*P01919'-diracl(20)*P02020'),mt'c21=Cp1*P011;Cp1*P022;k0;k0;. Cp2*P033;Cp2*