耦合非线性系统（共8篇）
耦合非线性系统 篇1
1 小世界网络简介
在网络理论中, 小世界网络是一类特殊的复杂网络结构, 在这种网络中大部份的节点彼此并不相连, 但绝大部份节点之间经过少数几步就可到达。
用数学中图论的语言来说, 小世界网络就是一个由大量顶点构成的图, 其中任意两点之间的平均路径长度比顶点数量小得多。
本文利用蔡氏振荡器对小世界网络的运行情况在matlab中进行仿真。
2 小世界网络仿真研究
(1) 对蔡氏振荡器中各个元件赋值:
耦合强度c=0.01, 控制增益d=0.3, R1=0.3Ω。
(2) 为系统中各个变量赋值:
耦合强度为输入值:di1=d;延时为0:tao=0。
(3) 创建一个N*N的全0矩阵建立初始的环状的规则网络, 结点网络有N个节点, 每个结点向与它最近邻的m个结点连出边, 求出邻接矩阵。
(4) 逆时针的边重连, 从节点到N-m-1。
(5) 随机选取一个数r=rand (1) ;
取出邻接矩阵中的非0元素位置:
求出非0元素个数:M=length (unconect) ;
正向取整:r1=ceil (M*rand (1) ) ;
连接这一对点, 构成SW小世界网络:
(6) 恢复小世界网络的邻接矩阵:
(7) 去掉自身节点形成的环A (k, k) =0;
(8) 为连接矩阵对角线上元素赋值
(9) 耦合强度矩阵cm (i, j) 表示节点i和节点j之间的耦合强度。
(10) 定义系统中其余各个参数:
此时我们不难发现随着小世界网络中连接概率的大幅增加, 尽管状态量似乎还有着同步性, 但状态量之间的差异已经十分明显, 可认为此时无法实现牵制同步。
进一步的研究发现, 当两点间存在连接的概率在p<0.283时, 对网络中一个节点施加控制, 可以使得整个网络达到牵制同步、且各个节点状态相一致, 在0.283<p<0.562时, 尽管对一个节点施加控制, 可以使得整个网络所有状态都达到平衡点, 但平衡状态的同步性越来越差, 各个节点显示出状态分离的趋势。而在p>0.562时, 已经无法通过控制系统中一个节点, 达到整个系统状态同步的目的, 此时对于小世界网络而言, 控制一个节点达到牵制同步不可行。
3 结论
本文通过Matlab仿真算法, 将小世界网络的连接情况用蔡氏振荡器的实例模拟出来, 同时三个状态量又能很好的刻画系统的同步情况, 为同步判据的检验提供了一个直观的方法。
实验结果表明了由于小世界网络介于规则和随机网络之间的特性, 对于小世界网络的牵制控制效果, 与小世界网络中任意两点间存在连接的概率有很大关系, 从而找到两个概率临界值0.283和0.562。当任意两点存在连接的概率小于0.283时, 验证了判据的有效性, 整个系统可以通过控制一个节点达到同步。而当任意两点存在连接的概率大于0.562时, 整个系统无法通过控制一个节点达到同步, 牵制控制方法不可行。
参考文献
[1]Liu Z X, Chen Z Q, Yuan Z Z.Pinning control of weighted general complex dynamical networks with time delay[J].Journal of Physics A, 2007, (01) :345-354.
[2]Li K, Lai C H.Adaptive-impulsive synchronization of uncertain complex dynamical networks[J].Physics Letters A, 2008, 1601-1606.
耦合非线性系统 篇2
(1.东北大学机械工程与自动化学院， 辽宁 沈阳 110819；2.沈阳鼓风机集团安装检修配件有限公司技术工程部， 辽宁 沈阳 110869)
多跨转子系统耦合故障定量诊断方法
许 琦1， 吴 昊1， 赵立超2， 姚红良1， 闻邦椿1
(1.东北大学机械工程与自动化学院， 辽宁 沈阳 110819；2.沈阳鼓风机集团安装检修配件有限公司技术工程部， 辽宁 沈阳 110869)
以有限元理论和谐波平衡理论为基础分析多跨转子系统的动态特性并提出多跨转子系统耦合故障定量诊断方法。利用故障前后转子系统的振动响应,各次谐波分量与系统剩余量方程的频率响应矩阵之间的关系推导诊断理论并获得诊断方程，以转-定子碰摩和联轴器不对中耦合故障为例建立多跨转子系统和实验台转子系统的有限元模型，通过数值仿真和转子实验台实验准确地诊断出故障位置，验证方法的有效性和稳健性。通过同一转速下故障前后转子系统至少n+1个测点的振动响应数据可以在线确定具有n个耦合故障转子系统的故障位置。对于转子系统的主要振动故障，如碰摩、不对中、裂纹等，诊断方法具有良好的适用性。
故障诊断； 多跨转子系统； 耦合故障； 谐波分量； 剩余量方程
引 言
转子系统作为旋转机械的核心部件，常常出现各种单一或耦合的振动故障，如质量不平衡、联轴器不对中、转-定子碰摩、裂纹、轴承和支承故障等，不仅影响其正常工作，严重时会发生机毁人亡的事故，造成重大的损失[1-2]。及时发现并掌握故障信息，对于处理这些故障有着重要的意义。
目前，学者对于转子系统的故障诊断做了大量的研究，主要分为两类方法：一类是非模型的故障诊断方法，另一类是数学模型的诊断方法。
非模型的诊断方法如振动信号分析方法、模态和信息处理方法，人工智能、专家系统的诊断方法，神经网络的诊断方法等。将测量的振动信号进行滤波和去噪，经过各种信号处理技术，如快速傅里叶变换(FFT)、小波变换(WT)、经验模态分解(EEMD)和Hilbert-Huang变换等，进行转子系统的各种故障诊断[3]。如Hongkai Jiang等[4]提出多小波包改进的经验模态分解(EEMD)法的旋转机械多重故障诊断方法，提高了分解结果的精度和准确性；Jimeng Li等[5]利用随机共振(SR)法增强非线性系统中的微弱信号，诊断淹没在噪声中的微弱信号，提出了Morlet小波变换利用噪声控制的二阶增强随机共振法进行故障诊断；Long Zhang等[6]采用多尺度熵(MSE)和自适应的神经模糊推理系统(ANFIS)诊断轴承故障，此方法不仅能可靠地分离不同故障类别，而且能够识别故障严重程度；Yaguo Lei等[7-8]提出统计分析、改进的距离评估技术和自适应的神经模糊推理系统(ANFIS)的智能故障诊断方法，并应用在滚动轴承的故障诊断。
目前多数故障诊断方法是非模型的诊断方法，能够诊断出故障的存在、类型及估算出故障所在位置，在精度要求不高的前提下信号分析方法直接、有效，但需要精确诊断的时候，难以达到故障位置的定量诊断。
数学模型的诊断方法即利用转子系统动力学特性及振动信号相结合进行故障诊断的方法，能够诊断出故障的具体位置，以及时准确地处理严重故障[9]。如Arun Kr Jalan等[10]提出一种基于模型的转子系统不对中和质量不平衡故障诊断方法，采用剩余量生成法结合振动信号，得到故障状态特性和故障位置；Mohit Lal等[11]提出受迫响应信息的频域最小二乘法拟合技术诊断算法估计汽轮发电机系统模型的多重故障参数，并具有一定的抗噪声影响能力；A S Sekhar[12]提出一种基于有限元模型的在线诊断转子裂纹故障的方法，采用等效力和剩余振动量的数学模型，结合信号分析方法快速傅里叶变换(FFT)得出裂纹的位置和特性；G N D S Sudhakar等[13]采用等效力最小化方法和振动最小化方法相结合诊断转子系统不平衡故障位置和严重程度，解决了多故障参数下较少振动测量值的问题；姚红良等[14]根据碰摩故障转子系统中任意两节点之间高次谐波分量之比等于无故障转子系统频率响应矩阵的相关元素之比的关系，利用广义等效力、剩余振动量和谐波平衡理论提出了诊断转子系统单碰摩故障的方法。
数学模型的诊断方法精度高，能够将转子系统的故障诊断具体到“点”上，同时可以得到一些故障的严重程度，但根据作者的了解，目前或仅对转子系统单故障进行诊断，或加速度、速度和位移响应均需测量，或仅利用工频响应而忽视二倍频、三倍频等高频响应等，有一定的缺陷和不足。
本文的目的是建立基于模型的多跨转子系统耦合故障定量诊断方法。首先，利用位移振动信号的振幅、频率和相位信息以提高诊断精度；其次，结合有限元理论和谐波平衡理论，利用故障转子系统的振动响应各次谐波分量与系统剩余量方程的频率响应矩阵之间的关系，推导并建立诊断方程；第三，利用有限元模型中一些节点的位移响应诊断出故障位置节点；最后，通过数值仿真和实验验证诊断方法的有效性和稳健性。
1 转子系统动力学模型
典型的转子系统由一些离散的叶轮、具有分布质量及弹性的轴段和轴承座等部件组成，将系统离散成由刚性圆盘、轴段等单元连接成的模型，各个单元间在节点处联接，忽略转子系统的轴向变形。弹性轴段单元的广义坐标分别为两端节点的位移的转角，如图1所示，其复数表示为
(1)
图1 轴段单元的有限元模型Fig.1 Finite element model of shaft unit
将多跨转子系统离散成具有N个节点、N-1个单元组成的有限元动力学模型，其整体动力学方程为
(2)
式中M为整体质量矩阵;C=D+ωG，D为整体阻尼矩阵，G为陀螺力矩矩阵，ω为转子转速;K为整体刚度矩阵，矩阵为2N×2N阶对角方阵，具体形式参见文献[15];u为振动响应矢量;F为质量不平衡引起的外激励矢量。
(3)
2 定量诊断方法
式(3)-(2)，得
(4)
当系统稳定运行，故障激励，也就是各广义故障力包含谐波分量时，系统的响应也包含各次谐波分量，故可以将Δu展开成各阶谐波分量和的形式，即
(5)
式中 Δuiejiωt为Δu的第i阶谐波分量。
将广义故障力也展开成各阶谐波分量和的形式，即
(6)
式中Aiejiωt，Biejiωt，…，Niejiωt分别为广义故障力的第i阶谐波分量。
根据谐波平衡理论[16]，有
(7)
(8)
由式(8)就可以定量计算多跨故障转子系统的故障位置。设多跨转子系统任意两个节点k1和k2，有
E(2k1-1,2L2-1)(jlω)Bl+…+E(2k1-1,2Ln-1)(jlω)Nl]/
[E(2k2-1,2L1-1)(jlω)Al+E(2k2-1,2L2-1)(jlω)Bl+
…+E(2k2-1,2Ln-1)(jlω)Nl]
(9)
式中 相关参数Al，Bl，…，Nl由下式求解。
(10)
对于具有n个故障的转子系统，相关参数也有n个，即需要n个节点的响应。
将转子系统所有节点号Li分别代入式(9)中计算δL1L2…Ln，当δL1L2…Ln的绝对值最小时，此时的L1，L2，…，Ln的取值即为转子系统的故障节点位置。即定量诊断n个故障节点位置至少需要n+1个节点的振动响应。
3 诊断计算步骤
本方法的具体计算步骤如下：
当多跨转子系统有n个耦合故障时
(1)分别测量转子系统故障发生前n+1个点的振动响应；
(2)分别测量转子系统故障发生后同一n+1个点的振动响应；
(3)计算响应剩余量并展开成谐波分量和的形式；
(4)计算频率响应函数矩阵；
(6)求解相关参数Al，Bl，…，Nl；
(7)将转子系统所有节点号Li分别代入计算δL1L2…Ln；
4 数值仿真
4.1 简单多跨转子模型
建立多跨转子系统的有限元模型，如图2所示。
具体参数如下：将整体转子系统划分为52个单元，即有53个节点，弹性模量E=210 GPa；其中轴段直径d=10 mm，长度l=20 mm；圆盘直径d1=80 mm，长度l1=10 mm；将联轴器简化为当量轴段，直径d2=20 mm，长度l2=10 mm；将支承视为等刚度弹性支承，支承刚度为k=1×106N/m；质量偏心位于节点9和30处，不平衡质量me=5×10-6kg，偏心距e=2 mm，不平衡相位为0。
假设此转子系统的故障为定-转子碰摩和联轴器不对中耦合故障。碰摩故障模型如图3所示，位于节点26处，故障力的表达式见文献[17]，碰摩刚度krub=2×105N/m，碰摩间隙erub=2.5 μm，碰摩类型为局部碰摩，碰摩角度为0°≤θ≤45°。
图2 多跨转子系统的有限元模型Fig.2 Finite element model of multi-span rotor system
图3 碰摩力模型Fig.3 Model of rub force
联轴器不对中故障模型如图4所示，位于节点17处，故障力的表达式见文献[10]，初始不对中量为Δx0=5 mm，α0=2°。
和节点17
图4 不对中力模型Fig.4 Model of misalignment force
图5 诊断结果Fig.5 Diagnosis result
4.2 支承刚度和振动响应影响分析
假定实际支承刚度为k=1×106N/m，而诊断时采用的支承刚度为k=1×1010N/m，数值模拟实验得到的结果如图6所示，碰摩位置和不对中位置仍分别在节点26和节点17处，即支承刚度的变化对本诊断方法几乎没有影响。
假定实际情况测得的振动响应不准确，其数值服从方差为0.05的正态分布，分别取2，3，4和5阶谐波分量进行总计400次仿真计算，得出的结果如图7所示。其中碰摩故障位置所在的节点26的准确率在56%；不对中故障位置所在节点17的准确率在72%，均远高于其他节点所占的百分比，可见此时仍可以诊断出耦合故障的大致位置，表明此方法具有很好的稳健性。
图6 支承刚度不准确时的诊断结果Fig.6 Diagnosis result when support stiffness is not accurate
图7 振动响应数据不准确时的诊断结果Fig.7 Diagnosis result when vibration response data is not accurate
4.3 工程多跨转子实例
建立某大型空分设备压缩机-汽轮机转子系统的有限元模型，如图8所示。
图8 压缩机-汽轮机转子系统模型Fig.8 Finite element model of compressor-turbine rotor system
具体参数如下：将整体转子系统划分为136个单元，即有137个节点；转子总长14 880 mm；偏心me=16 kg·m位于节点28和98处；工作转速为4 000 r/min；滑动轴承支承，分别在节点4，65，72和120处，其动力学方程见文献[15]。
图9 诊断结果Fig.9 Diagnosis result
5 实 验
本方法的实验验证在Bently转子实验台上完成，利用B&K3560B信号采集器采集并处理信号。实验设备及其有限元模型如图10所示，转子模型由2个圆盘、38个轴段单元、2个支承及1个联轴器组成，尺寸参数如表1所示。
4个位移传感器分别放置在节点9，15，21和31处，碰摩和不对中故障分别加在节点16和联轴器(节点37～43)上，在转速为2 400 r/min时采集振动数据并处理，其中节点15的故障响应曲线如图11所示，信号中有明显了高倍频出现。
图10 实验设备及其有限元模型Fig.10 Experimental equipment and finite element model of Bently rotor test rig
参数名称数值单元数42节点数43弹性模量/GPa210单元的直径和长度/mmds=10，ls=15圆盘的直径和厚度mmdd=75，ld=25节点1与2之间的长度/mmls1=20节点11与12，14与15，22与23，25与26之间的长度/mmls2=17．5联轴器长度/mmlc=61
图11 节点15的故障响应曲线Fig.11 Response curve on node 15
诊断结果如图12所示，图中显示最小值的坐标为(16,37)，因此，碰摩和不对中故障分别在节点16和37处。诊断方法成功地诊断出了故障位置。
图12 诊断结果Fig.12 Diagnosis result
6 结 论
(1)本文提出多跨转子系统耦合故障定量诊断方法，并以碰摩和不对中耦合故障多跨转子为例进行故障诊断，通过数值仿真和实验成功地将故障定位在“点”上，验证了方法的有效性和稳健性；
(2)对于具有n个耦合故障的多跨转子系统，仅需同一转速下故障发生前后的n+1个点的振动响应数据即可诊断，方法简单；
(3) 诊断方法适用于各种具有局部故障的转子系统中，即在微分方程中体现为某点存在周期性的广义故障力，诊断方法对转子系统的主要故障如转子不平衡、碰摩、不对中、裂纹、轴承及支承等均具有良好的适用性。
[1] 闻邦椿，顾家柳，夏松波，等. 高等转子动力学[M]. 北京:机械工业出版社, 2000.
[2] 韩清凯，于涛，王德友，等. 故障转子系统的非线性振动分析与诊断方法[M]. 北京:科学出版社, 2010.
[3] Feng Z P, Liang M, Chu F L. Recent advances in time-frequency analysis methods for machinery fault diagnosis: A review with application examples[J]. Mechanical Systems and Signal Processing, 2013, 38(1):162—205.
[4] Jiang H K, Li C L, Li H X. An improved EEMD with multiwavelet packet for rotating machinery multi-fault diagnosis[J]. Mechanical Systems and Signal Processing, 2013, 36(2):225—239.
[5] Li J M, Chen X F, Du Z H, et al. A new noise-controlled second-order enhanced stochastic resonance method with its application in wind turbine drivetrain fault diagnosis[J]. Renewable Energy, 2013, 60:7—19.
[6] Zhang L, Xiong G L, Liu H S, et al. Bearing fault diagnosis using multi-scale entropy and adaptive neuro-fuzzy inference[J]. Expert Systems with Applications, 2010, 37(8):6 077—6 085.
[7] Lei Y G, He Z J, Zi Y Y. A new approach to intelligent fault diagnosis of rotating machinery[J]. Expert Systems with Applications, 2008, 35(4):1 593—1 600.
[8] Lei Y G, He Z J, Zi Y Y, et al. Fault diagnosis of rotating machinery based on multiple ANFIS combination with GAs[J]. Mechanical Systems and Signal Processing, 2007, 21(5):2 280—2 294.
[9] Lees A W, Sinha J K, Friswell M I. Model-based identi?cation of rotating machines[J]. Mechanical Systems and Signal Processing, 2009, 23(6):1 884—1 893.
[10]Jalan A K, Mohanty A R. Model based fault diagnosis of a rotor-bearing system for misalignment and unbalance under steady-state condition[J]. Journal of Sound and Vibration, 2009, 327(3-5):604—622.
[11]Lal M, Tiwari R. Multi-fault identification in simple rotor-bearing-coupling systems based on forced response measurements[J]. Mechanism and Machine Theory, 2012, 51:87—109.
[12]Sekhar A S. Crack identification in a rotor system: a model-based approach[J]. Journal of Sound and Vibration, 2004, 270(4-5):887—902.
[13]Sudhakar G N D S, Sekhar A S. Identification of unbalance in a rotor bearing system[J]. Journal of Sound and Vibration, 2011, 330(10):2 299—2 313.
[14]姚红良，韩清凯，李凌轩，等. 基于谐波分量的转子系统碰摩故障定量诊断方法[J]. 机械工程学报，2012, 48(5):43—48.
Yao Hongliang, Han Qingkai, Li Lingxuan, et al. Method for detecting rubbing fault in rotor system based on harmonic components[J]. Chinese Journal of Mechanical Engineering, 2012, 48(5):43—48.
[15]钟一谔，何衍宗，王正,等. 转子动力学[M]. 北京：清华大学出版社, 1987.
[16]Grollg V, Ewins D J. The harmonic balance method with arc-length continuation in rotor/stator contact problems[J]. Journal of Sound and Vibration, 2001, 241(2):223—233.
[17]Zhang W M, Meng G, Chen D, et al. Nonlinear dynamics of a rub-impact micro-rotor system with scale-dependent friction model[J]. Journal of Sound and Vibration, 2008, 309(3/5):756—777.
Quantitative coupling fault diagnosis method of multi-span rotor based on harmonic components
XUQi1,WUHao1,ZHAOLi-chao2,YAOHong-liang1,WENBang-chun1
(1.School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China; 2.The Engineering Department, Shenyang Blower Works Group Corporation Installation Accessories Co., LTD., Shenyang 110869, China)
A model based method for online quantitative multi-fault diagnosis of multi-span rotor system is presented. Based on finite element method and harmonic balance theory, the dynamic characteristics of the rotor system are analyzed, and the relationship between the harmonic components of vibration responses and the frequency response matrix of residual equation in the fault rotor system is used to derive diagnosis theory and equation. The finite element model of a multi-span rotor system with rub-impacts and coupling misalignment faults is established. The effectiveness and robustness of the proposed method is verified by numerical simulation and then by test rig experiment. The advantage of the proposed method is that n fault locations can be detected on line by using at least vibration response data onn+1 nodes of rotor system before and after faults occur at the same rotating speed, and also the method has a good applicability for detecting the main vibration faults in rotor system such as rub-impact, misalignment and cracks.
fault diagnosis; multi-span rotor system; coupling fault; harmonic component; residual equation
2013-11-12；
2015-01-04
国家重点基础研究发展计划项目(2011CB706504);国家自然科学基金资助项目(51005042);中央高校基本科研业务费专项资金资助项目(N100403005)
TH165+．3； TK267；
A
1004-4523(2015)03-0495-08
10.16385/j.cnki.issn.1004-4523.2015.03.021
耦合非线性系统 篇3
本文研究了非线性分数阶发展方程耦合系统的非局部柯西问题
这里0<p, q<1, cDp, cDq是两个分数阶导数, f, g:[0, ∞) ×E→E]和w:C ([0, ∞) ×E) →E是已知的满足假设条件的函数, u0, v0是Banach空间E中的元素.
Oldham和Spanier[1]中系统的陈述了分数阶微分方程的应用, 详细请参阅Miller和Ross[2]和Kilbas等人的[3]
分数阶微分方程耦合系统的研究是相当重要的, 很多人做了研究参阅参考[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14].最近, Fang[15]研究了非线性分数阶微分方程奇耦合系统正解的存在性.Su[16]讨论了分数阶微分方程耦合系统边界值问题.Ahmad和Nieto[17]研究了三点边界问题的分数阶微分方程耦合系统存在性结果.
在本文中, 假设E是范数为|·|的Banach空间.令JR, C (J, E) 是从J到E, 范数为的连续函数Banach空间, 这里x∈C (J, E) .
对于E上的任意强连续半群 (即C0半群) {T (t) }t≥0, 在E上定义算子:
其定义域D (A) 是所有E上极限存在的x集合, 且是稠密的, A是闭的, 详情请参阅[13].
1预备知识
在本节中, 我们将介绍文中涉及到的空间、基本定义及用到的引理 (详见[18])
设B (E) 是E到E范数为||Q||B (E) =sup{|Q (u) |:|u|=1}的所有有界线性算子构成的空间, 这里Q∈B (E) , u∈E.全文中, 设A是E中一致有界算子C0半群{T (t) }t≥0上的无穷小生成元.明显的,
定义1.1函数的次数为α, 极限为0的分数阶积分定义如下:
假设右边是定义在[0, ∞) 上的点态, 其中Γ (·) 是gamma函数.
定义1.2下界为0的阶数为α函数f∈AC[0, ∞) 的Riemann-Liouville导数能够被写为:
定义:1.3阶数为α函数f∈AC[0, ∞) 的Caputo导数表示为:
注记1.1 (1) 如果f (t) ∈C1[0, ∞) , 则:
(2) 常数的Caputo导数等于0;
(3) 如果f是值域在E的抽象函数, 则:1.1—1.3中积分定义是在Bochner意义下得到的.
假设JR, 1≤p≤∞, 对于可测函数m:J→R, 定义范数
其中上的Lebesgue测度.令Lp (J, R) 是所有范数||·||LpJ<∞的Lebesgue可测函数m:J→R构成的Babach空间.
引理1.1 (H觟lder不等式) 如果|H|是Lebesgue可积的, 则可测函数H:[0, a]→E是Bochner可积的.
引理1.2 (Bochner'定理) 如果|H|是Lebesgue可积的, 则可测函数H:[0, 1]→E是Bochner可积的.
引理1.3 (Schauder不动点定理) 如果B是Banach空间E中的有界闭凸子集, F:B→B完全连续, 那么F在B内有一个不动点.
2主要结果
定义空间X={u (t) |u (t) ∈C ([0, 1], E) }和Y={v (t) |v (t) ∈C ([0, 1], E) }.依据[15]中的结论, X和Y是Banach空间.
对 (u, v) ∈X×Y令
显然 (X×Y, ||·||X×Y) 是一个Banach空间.
基于以上的论证, 给出方程组 (1.1) mild解的定义.
定义2.1若非局部的柯西问题 (1.1) 的解 (u, v) ∈X×Y满足下式:
称 (u, v) 是方程组 (1.1) 的mild解.
定义算子F:X×Y→X×Y,
其中
对任意的常数k, 设:
显然Uk在Banach空间X×Y中是有界闭凸子集.
在证明主要结果之前, 先介绍下面的假设.
(H1) 对任意t>0, T (t) 是一个紧算子;
(H2) 对每个t∈[0, 1]函数f (t, ·) :X→X和g (t, ·) :Y→Y是连续的, 任意 (u, v) ∈X×Y函数f (·, u) :[0, 1]→E和g (·, v) :[0, 1]→E是强可测的;
(H3) 对所有的 (u, v) ∈X×Y和几乎所有t∈[0, 1], 存在常数p1∈[0, p) 和q1∈[0, q) , 使得|f (t, u) |≤m1 (t) 和|g (t, v) |≤m2 (t) 成立;
(H4) w:C ([0, 1], E) →E是一致连续, 以及对所有x∈C ([0, 1], E) , 存在正常数L1, L2使得|w (x) |≤L1||x||+L2.
以下非局部柯西问题 (1) 的存在性结果是以Schauder不动点定理为基础.
定理1
定理2.1如果 (H1) - (H4) 满足, ML<1, 那么方程组 (1) 有mild解.
证明:对任意 (u, v) ∈Uk, 由于, 于是:
直接计算得到, 当t∈[0, 1]和p1∈[0, p) , 令:
利用引理:1.1 (Hlder不等式) 和 (H3) , 当t∈[0, 1], 得到:
类似的, 有:
因此, 对所有t∈[0, 1], 当s∈[0, t], |0∫∞θ (t-s) p-1hp (θ) T ( (t-s) pθ) f (s, v (s) ) dθ和0∫∞θ (t-s) q-1hq (θ) T ( (t-s) qθ) f (s, v (s) ) dθ|是Lebesgue可积的.由引理1.2 (Bochner'定理) , 对所有的t∈[0, 1], s∈[0, t], 可知0∫∞ (t-s) p-1hp (θ) T ( (t-s) pθ) f (s, v (s) ) dθ和0∫∞θ (t-s) q-1hq (θ) T ( (t-s) qθ) f (s, v (s) ) dθ是Bochner可积的.
接下来用Schauder不动点定理, 证明结果.
定义
其中
观察, Uk1显然是Banach空间X×Y中的有界闭的凸子集.
下面分两部分证明F在Uk1内有一个不动点.
第一步:F:Uk1→Uk1.
对所有的 (u, v) ∈Uk1和t∈[0, 1], 有:
类似的, 有:
因此, 对所有的 (u, v) ∈Uk1, 有F:Uk1→Uk1.
第二步:F是全连续算子.
因此, 推断出:
另一方面, 当t∈[0, 1]
即F1是连续的.同样的, 得到F2也是连续的.也就是, 算子F:Uk1→Uk1是连续的.
其次, 证明{F (u, v) , (u, v) ∈Uk1}是相对紧的.这就可以证明函数族{F (u, v) , (u, v) ∈Uk1}是一致有界和同等连续的.
对任意 (u, v) ∈Uk1有||F1v||≤k1, ||F2u||≤k1, 从而||F (u, v) ||≤k1.因此{F (u, v) , (u, v) ∈Uk1}是一致有界的.在下文中, 将证{F (u, v) , (u, v) ∈Uk1}是同等连续函数族.
对每个 (u, v) ∈Uk, 0≤t1<t2≤1, 得到:
其中:
运用 (11) 和 (12) 式中类似的证明, 得到:
当t1=0, 0<t2≤1, 很容易得到I3=0.当t1>0, ε>0足够小, 当θ∈ (0, ∞) , 有:
由于 (H1) 表明T (t) (t>0) 连续, 推断出F1是同等连续的.类似的F2也是同等连续的, 因此, F (Uk1) 是同等连续的.当t2-t1→0, 与 (u, v) ∈Uk1无关, | (F1v) (t1) - (F1v) (t2) |趋近于零.这意味着{F (u, v) , (u, v) ∈Uk1}是同等连续的.
耦合非线性系统 篇4
近年来,由于脑电信号检测与分析技术的发展,使中枢疲劳的监测技术有了很大进步,脑电已成为最广泛的评定中枢神经系统变化的指标之一,被誉为监测疲劳的“金标准”[1]。Murata A利用事件相关电位P300的幅度和潜伏期大小对中枢疲劳进行分析,发现随着中枢疲劳程度的增加,潜伏期延长,幅度降低[2]。Jung TP用脑电功率谱对大脑警觉性水平进行分析,发现脑电功率谱可以反映大脑警觉状态的波动情况[3]。虽然先前研究已取得了一些有意义的成果,但中枢疲劳脑电的分析基本都是针对单一电极进行的,没有充分利用不同脑区采集得到的脑电信息,忽略了疲劳对不同脑区间功能联系的影响和研究。本文尝试采用非线性回归(NLR)系数方法对不同中枢疲劳状态下的不同脑区功能耦合程度的变化进行分析,将有助于更好地解释中枢疲劳现象及其机理。
1 非线性回归系数
Pijn等人最先建议使用NLR系数来量化两个EEG信号之间的幅度耦合程度[4]。假定x和y分别是两个导联上的脑电信号。给定信号x,信号y的期望值用表示µy|x,它是y对x的递归曲线
式中p(y|x)是给定信号x时信号y的条件概率。根据这条递归曲线,从x预测y的值得到的y的方差减小,就是由下式给出的关联测量:
式中从这条递归曲线估计的E[(y-µy|x(x))2]被称为解释的方差,即它是在信号x的基础上被解释或预测。通过从总方差var(y)中减去被解释的方差,可获得未被解释的方差。关联测量的基本思想是,如果把信号y的幅度看作是信号x的幅度的一个函数,给定某个x值,y的值可按照一条NLR曲线进行预测。这种测量的估计被称为NLR系数h2。
实际工作中,数据是离散的,p(y|x)未知,因此需要由离散数据估计p(y|x)。设有两组N点离散数据:xn=[x1~xN],yn=[y1~yN],把它们的对应点[x1,y1]~[xN,yN]画在直角坐标图上。为了获得这条递归曲线的一个近似,将x的幅度值等概率地分成M段,从而保证落在每一份中的样本点数大致相等。计算每段中点的x值pi(i=1~M)。求出各段中y(n)的均值qi(i=1~M),将它视为pi已知情况下y在这一段的估计值,并将各pi、qi标示于同一图上。最后,利用这些已知点,可以通过插值的方法得到在各己知点x处y的对应估计,也就是在x已知条件下y的回归曲线。通过回归曲线可以逐段求得2y|x。因此,NLR系数可以按照下面的表达式计算:η
式中是这条递归曲线的线性分段近似,而
2 中枢疲劳实验设计及脑电数据采集
实验对象是在校学生自愿者50名,年龄为20~27岁,均为男性,身体健康,无重大脑部疾病史、心理疾患及器质性疾病史。实验前要求受试者保证充足的睡眠,避免服用酒、茶和咖啡等任何使中枢兴奋或抑制的食物和药物,避免剧烈运动。实验开始时间为上午8点左右。实验前受试者完全了解整个实验程序和要求,均自愿参加,合作态度良好。
实验前,向受试者讲明实验步骤及测试方法,并进行预实验,以适应和熟悉测试程序。对测试方法也进行反复练习,达到能正确理解指导语,能熟练操作,以消除“学习效应”对实验的影响。
实验开始后,首先对受试者进行主观疲劳程度的心理评测,而后要求受试者安静闭目尽可能什么也不想,尽可能放松,记录5min的脑电信号,作为疲劳实验前的状态记录。然后,要求受试者完成3种不同类型的脑力劳动任务直至筋疲力尽,或两小时结束。而后记录受试者安静闭目状态5min的脑电信号,作为疲劳实验后的状态记录。同样,对受试者的主观疲劳程度的进行心理评测。
脑电数据的采集使用Neuroscan 32导脑电仪,采用10~20导联头皮电极系统,记录了Fp2、Fp1、F4、F3、C4、C3、P4、P3、Fz、Cz和Pz共11通道脑电信号,同时记录了水平和垂直眼电信号。参考电极置于双侧乳突连线上,放大器的通带频带为0.01~100 Hz,采样频率为250 Hz,各导联阻抗均小于5 KΩ。实验是在在室内正常照明、安静、尽量舒适的环境(室温25℃)下进行的。
采用SPSS 13.0软件对实验数据进行分析,连续脑力劳动任务前后受试者的主观疲劳量表得分以及NLR系数分别采用配对T检验。
3 实验结果
3.1 中枢疲劳的主观评价
对实验前后的中枢疲劳程度分别采用SSS、KSS、SFC、Samn-Perelli和CR10量表进行主观评价[6]~[10],得分越大,表明主观感觉疲劳程度越高。实验前后中枢疲劳程度的主观评价结果如图1所示。
从图1可以看出,实验后主观疲劳量表得分显著增加(P<0.005)。这表明受试者经过连续长时间脑力劳动任务后,主观感觉疲劳程度增加,说明连续长时间脑力劳动任务对受试者的心理产生了较大的影响,导致主观感觉疲劳程度增加。
3.2 中枢疲劳的幅度耦合分析
在计算NLR系数之前,首先选取连续1min的脑电数据作为研究对象,对各导联脑电信号进行分段,选取前10s采样数据为基本数据,每次移动1s的数据。对不同中枢疲劳状态下脑区内和脑区导联间的幅度耦合程度的变化进行分析,分别了计算了θ和β脑电节律的NLR系数值。为减少NLR系数值波动的影响,计算1min内所有数据段NLR系数值的平均值,代表脑电信号在1min内的NLR系数值。实验前后脑区内和脑区间导联对的NLR系数在θ和β脑电节律的平均值的统计分析,结果如图2所示。
从图2的比较结果可以看出,完成实验任务后,β脑电节律脑区内导联之间的NLR系数平均值在中央和顶区显著降低(P<0.05);而脑区间导联之间的NLR系数平均值在额-中央、额-顶以及中央-顶区间显著降低(P<0.005);而θ脑电节律脑区间导联之间的NLR系数平均值则在F4-C4和Fz-Cz间显著增加(P<0.05);脑区内导联之间的NLR系数平均值在F4-Fz和C3-Cz间也有
所增加。并达到了显著性意义的水平(P<0.05)。
4 讨论与结论
一般认为连续长时间脑力劳动对受试者的心理和行为会产生较大的影响,导致主观感觉疲劳程度增加,脑能力下降,警觉性和注意力水平降低,大脑处理信息的能力下降,中枢疲劳程度增加[11]。为了验证本文的实验任务对中枢疲劳的影响效果,采用5种中枢疲劳程度的主观评价量表对连续脑力劳动前后受试者的主观感觉疲劳程度进行评测,实际分析结果与理论分析完全一致。这说明受试者经过长时间脑力劳动后,中枢疲劳程度会增加。
NLR系数测量了两个信号之间幅度的非线性耦合程度,反映了它们之间的幅度同步性高低。本文的研究结果表明,随着中枢疲劳程度的加深,中央和顶脑区内导联之间β脑电节律的幅度耦合程度下降,额-中央、额-顶以及中央-顶区导联间β脑电节律的幅度耦合程度也降低;而右侧和中间脑区的额-中央导联间θ脑电节律的幅度耦合程度增加,额区右-中间和中央区左-中间导联之间θ脑电节律的幅度耦合程度也增加。以往中枢疲劳的研究表明,中枢疲劳状态与θ和β波密切相关,随着疲劳的加深,大脑的警觉性和兴奋性下降,脑电慢波成分增加和快波成分减少,β脑电节律的活动降低,θ脑电节律的活动增强[12]。中枢疲劳脑电的幅度耦合程度变化可以从另一个侧面反映脑电节律随中枢疲劳状态的变化情况。
通过本文分析可以看出,非线性回归(NLR)系数可以灵敏地反映出不同中枢疲劳状态下脑电节律的非线性幅度耦合程度的变化,为中枢疲劳研究开辟了一条新途径,其方法可成为研究中枢疲劳脑电的有力工具。
参考文献
[1]曹雪亮,苗丹民,刘练红.脑力疲劳评定方法现状[J].第四军医大学学报,2006,27(4):382-384.
[2]Murata A,Takasawa Y,Takasawa Y.Evaluation of mental fatigue using feature parameter extracted from event-related potential[J].International Journal of Industrial Ergonomics,2005,35(8):761-770.
[3]Jung T P,Makeig S,Stensmo M.Estimating alertness from the EEG power spectrum[J].IEEE Trans on Biomedical Engineering,1997,44(1):60-69.
[4]Pijn J P M,Vijn P C M,Lopes da Silva F H,et al.The use of signal-analysis for the localization of an epileptogenic focus:a new approach[J].Adv Epileptol,1989,17:272-276.
[5]Wei Q,Wang Y,Gao X,et al.Amplitude and phase coupling measures for feature extraction in an EEG-based brain-computer interface[J].Journal of Neural Engineering,2007,4(2):120-129.
[6]Li Z Y,Jiao K,Chen M,et al.Effect of magnitopuncture on sympathetic and parasympathetic nerve activities in healthy drivers--assessment by power spectrum analysis of heart rate variability[J].European Journal of Applied Physiology,2003,88(4-5):404-410.
[7]Hoddes E,Zarcone V,Smythe H,et al.Quantification of sleepiness:A new approach[J].Psychophysiology,1973,10:431-436.
[8]Akerstedt T,Gillberg M.Subjective and objective sleepiness in the active individual[J].International Journal of Neuroscience,1990,52:29-37.
[9]Samn S W,Perelli L P.Estimating aircrew fatigue:A technique with implications to airlift operations[M].Brooks AFB,TX:USAF School of Aerospace Medicine,Technical Report No.SAM-TR-82-21,1982.
[10]Borg G.Borg’s perceived exertion and pain scales[M].Champaign,IL:Human Kinetics,1998.
[11]Vidulich M A.The mental psychology of subjective mental workload[C].In Hancock PA,Meshkati N,eds.Human mental workload.North-Holland:Elsevier Science Publishers,1988,219-229.
耦合非线性系统 篇5
近年来,以形状记忆合金、电致伸缩材料、磁致伸缩材料和压电材料主导的智能材料获得了迅速发展,其中的压电材料成为学者研究的焦点[1],适应于各场合的各类微型驱动装置层出不穷。Dragan等[2]由千足虫的爬行得到启发,利用两个U形压电双晶片,研制了一台低频压电电机;Toyama[3]将设计的球形压电超声电机作为相机作动器用在管状探测机器人上;Tomoaki[4]研制了一台定子体积只有1mm3的微型超声电机,成为最小的压电电机之一;赵淳生团队研发的压电超声电机首次用于“嫦娥三号”探测器,实现了其在月球上的完美着陆[5]。
叠堆型压电驱动器在承受较大压力的同时能输出较大位移,一直受到学者的青睐。 陈维山等［６］利用20个压电叠堆驱动器设计和制造了一台利用径向弯曲模态的行波压电电机,其最大转速、最大转矩分别为146r/min、1.0N·m;Oli-ver［７］提出一种利用8个5mm×5mm×50mm的压电驱动器进行传动的谐波压电电机,该电机堵转转矩为0.75N·m;笔者利用2 个5 mm×5mm×20mm的压电驱动器设计了一种机电集成压电谐波传动系统［８］,该传动系统将压电驱动、谐波传动和活齿传动集成为一体,具有传动比大、输出转矩大和寿命长等优点。
压电驱动器的动力学特性将会对驱动力和输出位移产生重要影响,不少学者对叠堆压电驱动器进行了动力学分析。Vahid等[9]通过有限元方法对三自由度锥形压电驱动装置进行了动力学建模与分析;王光庆[10]对压电叠堆式发电装置进行了建模与仿真分析。然而研究者对压电的非线性动力学分析主要集中在压电片与梁或板的层叠结构上[11,12,13],对于使用率较高的叠堆型压电驱动器却未曾见到相关研究。因此,本文在压电非线性效应和位移非线性效应的基础上对压电驱动器进行机电耦合动力学建模,运用Linz Ted-Poincaré法(L-P法)对驱动器弱非线性自由振动、接近共振时受迫振动和亚谐波振动进行分析与求解,最后通过四阶Runge-Kutta数值法对文中的动力学推导进行验证。
1 非线性压电效应
由于压电陶瓷磁滞效应的存在,使其在通入激励信号后产生的应力和应变呈非线性变化,故考虑非线性时的压电应变方程为［１４］
式中,S33为压电材料弹性柔度系数,m2/N;d33为压电应变常数,m/V;T3为压电驱动器预应力,Pa;E3为电场强度,E3=U/lp;lp为压电片厚度,mm;U为驱动信号,U=Up-p(1+cosωt)/2;Up-p为驱动电压峰峰值,V;ω为驱动信号频率,rad/s;d333为二次非线性压电系数;K333为机电耦合矩阵中的元素。
当T３=0时,式(1)化简为
根据力和应力的关系式σｐ= Fｐ/Aｐ和广义胡克定律σｐ=c３３S３,可得压电驱动器末端非线性输出力为
式中,c３３为弹性刚度系数;Aｐ为驱动器横截面积。
2 非线性动力学方程
2.1 机电耦合动力学模型
考虑压电材料的位移非线性效应,引入量纲一小参数ε,则无激励时驱动器的非线性轴向应力为
其中,εｙ为无激励信号时压电堆内轴向应变,且εｙ=∂v/∂y,v和y分别为压电堆的轴向振动位移和纵坐标。
对压电驱动器施加激励信号时,设激励力幅值与小参数ε同数量级,则驱动器总的内力为
式中,lnp为压电驱动器总长度。
压电驱动器动力学模型如图1所示,假设各压电陶瓷片之间是理想黏结的,驱动器在整体上是一个连续杆,对图1中微元dy的受力在y向应用牛顿定律,可得压电驱动器非线性机电耦合动力学方程为
式中,ρｐ为压电叠堆的密度。
令v(y,t)=ф(y)q(t),其中,ф(y)和q(t)分别为压电堆轴向振动的模态函数和时间响应函数,代入式(6)化简得
式中,фｉ为第i阶模态函数。
2.2 弱非线性自由振动
当式(7)中激励信号Uｐ－ｐ=0且ε充分小时,系统为弱非线性自由振动系统,此时系统中只存在位移非线性,没有压电非线性。由此可得弱非线性自由振动方程为
式中,ω０为线性系统的固有频率。
采用L -P法求解非线性方程,将式(8)的解q(t,ε)和振动频率ω 展成ε 的幂级数:
引入新的自变量κ=ωt,将原微分改为对κ的微分,并将式(9)和式(10)代入式(8)。令ε的同次幂的每项系数都为0,导出前二次的线性方程组为
设压电驱动器的初始位移为δ０,初始速度为0,初始位移可根据线性系统受迫振动位移给出,这里不再列出具体过程。由初始条件和零次近似方程可得出零次近似解为
将式(12)代入一次近似方程,同时为避免出现久期项,令cosκ 的系数为0,得出 σ１=-3b３δ０２/(4ω０２),故可得一次近似解为
同理,将式(13)代入二次近似方程,令cosκ的系数为0可求得σ２,对二次近似方程求解可得
忽略二阶以上高阶项,驱动器非线性近似解为
压电驱动器非线性系统振动频率ω与初始位移δ０间的关系为
2.3 接近共振时非线性受迫振动
当Uｐ－ｐ≠0时,式(7)即为接近共振时非线性机电耦合动力学方程,此时压电激励力的幅值与小参数ε同数量级,激励频率ω接近线性固有频率ω０,引入频方差εσ=ω２-ω０２(σ为中间变量),将频方差和式(9)代入式(7),令ε的同次幂系数相等,导出下列近似方程组:
设零次近似方程的解为
式中,A０、B０为常系数。
由初始条件得A０=δ０,将式(18)代入一次方程,为避免久期项,令cosωt和sinωt的系数为0,可得
求解式(19),可得B０=0,同时得出
则压电驱动器接近共振时一次非线性近似解为
将式(21)代入二次近似方程,同时为避免出现久期项,令cosωt的系数为0,可得
故可求解出二次近似方程的解为
忽略高阶项,将式(18)、式(21)和式(23)代入式(15),可得压电驱动器接近共振时的非线性近似解。
2.4 亚谐波共振响应
当固有频率ω０接近激励频率ω 的1/3倍时,系统也会发生强烈的共振,这种现象为亚谐波响应。由于亚谐波响应是由较强的激励引起的,故激励力与ε不再是同数量级,其动力学方程变为
设(ω０-ω/3)与ε同数量级,令
将式(9)和式(25)代入式(24),令ε同次幂系数相等,导出零次和一次近似方程为
亚谐波零次近似方程的解为
将式(27)中二次谐波项去掉后代入一次方程,为避免出现久期项,令cos(ωt/3)的系数为零,得出
求解式(28)可得
压电驱动器亚谐波共振产生的条件是
由一次近似方程解出一次非线性近似解为
式中,C０、C１、C２为系数。
忽略一阶以上各项,亚谐波共振非线性近似解表示为q(t,ε)=q０(t)+εq１(t)。
3 算例求解与分析
3.1 压电驱动器非线性输出特性
本文采用5mm×5mm×20mm的压电驱动器作为研究对象,参数如表1所示。将参数代入式(2)和式(3),可得驱动器输出应变和输出力随电场强度E３变化曲线以及增压和减压时输出位移随电压的变化曲线,如图2所示。改变参数压电应变常数d３３和压电片厚度lｐ,得到压电驱动器非线性应变随参数变化曲线,如图3所示。由图2和图3可知:
(1)随着电场强度E３的增大,线性应变和非线性应变之间的差值增大,在E３=2V/μm时,非线性应变比线性应变小22.8%。同理,非线性输出力和线性输出力之差也随E３的增大而增大,在E３=2V/μm时,非线性输出力比线性输出力小25.8%。
(2)参数d３３和lｐ对压电驱动器非线性应变S３都有较大影响,非线性应变S３随d３３的增大而增大,随lｐ的增大而减小,且S３随d３３和lｐ的变化幅度较均衡。
(3)在增压过程中,非线性输出位移曲线向上凸起,减压过程中,非线性输出位移曲线向下凹。出现这种现象的原因是压电存在磁滞效应,使得电压下降时输出位移不能按照原路返回。
3.2 非线性幅频特性
取小参数ε=±0.2、初始位移δ０=0.6mm、激励信号峰峰值Uｐ－ｐ=150V,将ε和δ０代入式(16),得到前4阶弱非线性自由振动的固有频率,如表2所示。分别选取小参数ε、电压峰峰值Uｐ－ｐ、压电应变常数d３３和弹性刚度系数c３３作为研究对象,分析参数改变时接近共振时幅频曲线的变化情况,作出一阶幅频响应随参数变化图,见图4。
由表2和图4可得:
(1)当ε<0时,非线性固有频率ω小于线性固有频率ω０;当ε>0时,ω大于ω０;阶数相同时,ε<0或ε>0时|ω０-ω|的值恒定;随着频率阶数的增加,频率变化率|ω０-ω|/ω０变小,非线性现象减弱。
(2)与线性系统相比,非线性共振不出现在ω =ω０处及其附近,而是出现在偏离ω０较远处。当ω恒定时,对应的振幅|δ０|可以取到3个值,这种现象即是非线性中的跳跃现象。幅频响应中的骨架线主导了频响曲线的形状,反映了不同激励下振幅与激励频率的关系,且F１越大时频响曲线偏离骨架线越远。
(3)当ε增大时,频响曲线骨架线弯曲程度变大且向横轴靠近,相同频率对应的振幅值减小,共振曲线随骨架线变化趋势相同。
(4)随着Uｐ－ｐ改变,幅频骨架线没有发生变化,共振曲线偏离骨架线的程度发生变化,且Uｐ－ｐ越大,偏离程度越大。主要原因是Uｐ－ｐ通过改变F１的值来影响共振曲线的,而骨架线不受F１的影响。
(5)d３３对幅频响应影响较小,随着d３３的增大,骨架线无变化,共振曲线偏离骨架线程度增大。d３３和Uｐ－ｐ对幅频的影响都是通过改变F１的值实现的。
(6)当c３３增大时,频响骨架线连同共振曲线一同沿频率增大的方向平移;c３３取不同值时,各频响曲线形状没有发生变化。可见,c３３的改变使得线性固有频率发生了变化,进而使非线性频响曲线偏移。
一般情况下,压电驱动器通常与弹簧或质量块固连作为一个系统进行工作。以弹簧为例,当压电驱动器与一个线径、中径和长度分别为0.5mm、5mm和15mm的压缩弹簧串联时,系统的固有频率会迅速下降,如表3所示。故压电驱动器与弹性元件配合使用时在频率较低时也可能发生共振。随着技术的不断提高,压电驱动器的性能也会有所提升,压电驱动器在高频激励下的驱动有望实现。本文中高频研究的价值一方面会在层叠片数较少的压电叠堆中得到应用,另一方面可在未来的压电驱动器中得到体现。
Hz
3.3 非线性动态响应比较
取电压峰峰值Uｐ－ｐ=150V,分别对非线性自由振动、接近共振时受迫振动及亚谐波共振时一阶线性与非线性响应进行对比,三种振动对应的激励频率分别为33 958 Hz、33 958 Hz和101 874Hz。非线性自由振动不存在激励信号,故只存在位移非线性;接近共振时受迫振动和亚谐波受迫振动存在位移非线性和压电非线性,故分4种情况:同时存在两种非线性、只存在位移非线性、只存在压电非线性、线性,进行对比分析,如图5所示。表4所示是接近共振时4种情况振幅最大与最小值对比。
由图5和表4知:
(1)三种振动形式中,线性振动的幅值始终大于非线性振动的幅值,非线性对接近共振时受迫振动影响最大,对亚谐波受迫振动影响最小。
(2)在非线性自由振动中,非线性响应幅值比线性响应幅值小3.3%;在接近共振受迫振动中,两种非线性效应作用时非线性幅值与线性幅值差距最大,此时非线性幅值比线性幅值小10.8%;亚谐波受迫振动中,线性与非线性幅值变化不明显。
(3)由表4知,当只考虑位移非线性时,非线性响应幅值比线性幅值小3.2%;而当只考虑压电非线性时,非线性幅值比线性幅值小8.0%;压电非线性是位移非线性影响力的2.5倍,故在压电驱动器中压电非线性的影响起主要作用。
(4)在接近共振时受迫振动响应曲线中平衡位置不在零线处,这是由于激励信号是带偏置的余弦信号,偏置信号对接近共振受迫振动起了作用,而偏置对亚谐波受迫振动的作用却很弱。
4 数值验证
取电压峰峰值Uｐ－ｐ=150V,激励频率分别为33 958Hz、33 958Hz和101 874Hz。采用MATLAB的四阶Runge-Kutta指令对式(8)、式(7)及式(24)进行数值求解,并将一阶数值结果与解析解进行对比。图6所示是非线性响应一阶数值解,其中,图6a是三种振动的相图,1、2、3分别代表非线性自由振动、接近共振时受迫振动及亚谐波受迫振动。由图6得出规律:
(1)非线性振动相图曲线是由一系列椭圆曲线叠加而成,曲线呈闭合状,是外加激励频率与固有频率共同作用形成的周期运动。三种振动响应的相图收敛于闭合曲线,可见三种振动的动力学方程的解是收敛的,振幅是稳定的。
(2)弱非线性自由振动时,数值解与解析解的相位相同,幅值最大误差为8.2%。
(3)接近共振受迫振动时,数值解和解析解幅值和相位都存在误差。在一阶响应中,解析幅值比数值幅值小11.9%,相位随着时间的增加而增大,在0.1ms时相位差为0.15!。
(4)在亚谐波振动中,数值解和解析解相位同步,幅值误差最大为13%。
5 实验分析
采用德国OptoMET公司Vector类型的激光测振仪对压电驱动器-压缩弹簧系统进行振动测试,如图7所示。图中,数据采集器对测试数据进行处理然后通过软件在电脑输出频谱,压电驱动电源为XMT三通道具有反馈的信号放大装置。实验时,对压电驱动器通入50Hz具有正偏置的正弦激励信号,压电驱动器-弹簧系统在激励作用下产生振动,通过PicoScope软件对振动波形频谱分析可得到共振频率,表5所示为实验测试频率与理论频率比较。由表5可知:实验频率与理论频率最大误差控制在10%以内,且实验频率与非线性频率更接近。
6 结论
耦合非线性系统 篇6
目前,国内胡建华、王连华等学者通过对索结构单元进行柔性迭代分析,依据内力改变与位移增量的关系,建立了一种具有较高精度的索结构几何非线性分析方法———悬链线索单元法[2]; 西南交通大学的冉志红,李乔等人在考虑垂度和斜度的基础上,利用奇异摄动解法,推导出了斜拉索振动频率和振型的解析式[3]。国外,Kiisa Martti、Idnurm Juhan和Idnurm Siim三位学者在集中载荷的作用下,对弹性索进行了离散分析[4]; 而Vairo Giuseppe教授根据应变与弦长变化的关系,针对斜拉索的非线性响应提出了一个精致的闭式二阶模型[5]。但是,这些研究都是针对悬索桥或斜拉索桥的单一索结构,并且也没有将影响索结构几何非线性的因素进行完全考虑。悬索管桥的索结构作为一个体系,不但各索之间相互耦合作用,多种几何非线性行为也相互耦合,具有独特的结构特性和受力特点,笔者将在悬索管桥正常运营的现状态下,充分考虑几何非线性因素耦合作用的基础上,研究其索结构体系的应力分布规律,为以后悬索管桥的优化设计、可靠性分析以及风险评价提供参考依据。
1 索结构体系
一座悬索式管桥的完整索结构体系一般应该由1 对主索、1 对风索、1 对下稳定索、若干对吊索及系索组成。其中最重要的是主索,主索两端与塔架相连,整体几何形状在竖直平面内呈抛物线型; 它在整个悬索管桥的结构中起着至关重要的作用,是最主要的受力结构,基本上承载着整个管桥的所有内在载荷和外在载荷; 一旦主索发生失效破坏,将对悬索管桥产生致命的影响,会导致整个管桥的倒塌、断裂。风索的两端固定在水泥制成的锚固墩上,在水平面内呈抛物线型,在竖直平面内通常伴随着桥桁架的形状,可能呈拱形分布,也可能水平分布。由于悬索管桥的具有大柔性、大跨度的特点,桥身结构在横向极其不稳定,受到横向载荷作用时,容易出现大幅度的漂移,导致管桥毁坏,发生事故,风索的作用就是加强管桥在横向的约束,是管桥承受横向载荷的主要结构。下稳定索的两端也固定锚固墩上,索体与桥桁架底部相连,因此沿着桥桁架分布,在横向和纵向都起着稳定整个管桥的作用。吊索是连接主索和桥桁架的结构,是除主索以外的主要受力索。系索则连接着风索和桥桁架,正常情况下与吊索的对数相同,和风索一起对管桥进行横向约束。
2 理论分析
2. 1 几何非线性行为
在结构分析中,非线性问题包括材料非线性、几何非线性和状态非线性三类; 其实质是求解过程中结构的刚度矩阵不断变化。几何非线性是指结构经受大变形,几何形状的变化引起结构响应的非线性。悬索管桥的索结构属于大跨度、低阻尼、超柔性体系,变形特征为大位移、小应变,是典型的几何非线性结构[6]。虽然索结构的刚度与位移呈明显的非线性关系,但应力、应变仍符合虎克定律,即材料仍是线性的。索结构体系的几何非线性行为主要表现为:
2. 1. 1 结构大位移
这是影响索结构几何非线性最主要的因数。在外载荷作用下,主索、风索、下稳定索及系索会发生下挠,吊索将会出现伸长和倾斜,这些位移通常情况下比较大,不可避免地将对索结构的内力产生影响。由于力和变形是非线性的,因此在对索结构进行分析时,要依据变形后的几何位置来建立力的平衡方程[7]。
2. 1. 2 初始应力的影响
为了使柔性的索结构保持一定的几何形状,在悬索管桥成桥时,会对各个索施加一定的初始应力,这必然会使索结构具有初始刚度。在受到后续的载荷作用时,初始应力产生的初始刚度对后续的变形存在一定的抵抗力,而后续载荷建立新的平衡时必须打破原有平衡,这个过程中必然存在着能量消耗,也就是说初始应力可以提高结构刚度。
2. 1. 3 自重垂度的作用
进行有限元分析时,索单元通常选取直杆单元。但实际上索结构在自重的作用下将会产生一定的垂度,受到轴向力作用时,其弦长变化由材料的弹性伸长和结构的垂度变化两部分组成,其中后者使轴向力与弦长形成几何非线性关系[8]。
2. 2 耦合刚度矩阵
对悬索管桥索结构进行分析时,充分考虑以上三种非线性行为的耦合作用,采用有限位移理论[9],其耦合的刚度矩阵表达式如下:
式( 1) 中Ke为结构弹性刚度矩阵,KL为大位移刚度矩阵,R为恒活载内力,Kσ为结构几何刚度矩阵,δ为节点位移向量。
选取两节点组成的直杆单元来模拟索单元。取一个杆单元,长度为l ,横截面积为A ,弹性模量为E ,节点位移的列阵形式为qe= {ui,vi,uj,vj},由于索结构在实际状态下除了自重,只受轴向拉力,所以杆单元也只承受轴向拉力的作用,其线性应变根据经典的弹性理论得:
但是索结构是非线性的,横向位移v会引起单元的附加伸长,因此杆单元的应变可表达为:
式( 3) 中为应变的非线性项。不考虑杆单元的弯曲变形,那么位移便是x的线性函数:
式( 4) 中I = I2为二阶的单位矩阵,Ni、Nj为形函数,,由公式( 3) 和公式(4) 可以得到:
式( 5) 中T为体载荷向量,。因为杆在变形后,转角相对较小,故。
根据胡克定律,杆的拉伸应力为:
考虑到索结构还受自重垂度的非线性影响,所以上式的弹性模量E应依据Ernst公式[10]进行换算,采用等效弹性模量Eeq:
式( 7) 中W为单位长度索的重力,l' 为索弦长在水平面的投影长度,F为索内的张力。故杆的应力变为:
将各个参数值代入式( 8) 可得,结构弹性刚度矩阵为:
几何刚度矩阵为:
因为F是索的初始轴向力,即[K]σ也称为初始应力刚度矩阵,体现了成桥时初始轴向力对横向位移的作用。
大位移刚度矩阵为:
故最终索结构体系的耦合刚度矩阵为:
综上所述,所构建的刚度矩阵综合考虑了索结构的大位移、初始应力以及自重垂度三种几何非线性行为的影响,理论上具有足够的精度,可以用来进行下面的有限元分析。
3 有限元模型
3. 1 模型的建立
选取索结构体系比较完整的黄河悬索管桥为建模对象,管桥横跨黄河两岸,一端位于陕西省府谷县,另一端位于山西省保德县。桥跨长度270 m,塔架高44 m,桥面为桁架结构,两边高、中间低,呈微拱形,高度差为6 m,输气管道为 Φ660 × 14. 3 mm,最大输气压力6. 4 MPa,钢丝索包括1 对主索、1 对风索、1 对下稳定索、53 对吊索和53 对系索,每对吊、系索间相隔5 m。根据管桥各部分结构的受力特点和材料特性[11]选择单元类型,具体的单元材料特性如表1 所示。
结合工程实际情况,进行结构模型的简化及优化和边界条件的约束,依据施工图纸以及现场测量结果,运用ANSYS软件建立管桥的三维空间有限元仿真模型[12],为了方便研究,管桥桥跨中点为坐标原点,从府谷端到桥跨中点的X轴坐标为负值,从桥跨中点到保德端的X轴坐标为正值,并沿着管道内部天然气介质输送的方向逐渐增大,同时将每对吊索和系索进行编号,从小到大依次标号1 到53。图为1 全桥模型,图为2 索结构体系模型。
3. 2 索结构找形
以悬索管桥的成桥状态为模型初始平衡状态[13]进行索结构找形,因为此时管桥各部分结构的几何位置是已知的,所以模型在受到成桥恒载作用时应保持与管桥成桥状态相同的几何位置,这可以通过不断调整索结构的初始应变来近似地实现,假设主索的初始应变为 ε主0,风索的初始应变为 ε风0,下稳定索的初始应变为 ε下0,吊索的初始应变为ε吊0,系索的初始应变为 ε系0,则管桥模型成桥状态最大位移值f0max是 ε主0、ε风0、ε下0、ε吊0、ε系0的函数,记为f0max= f( ε主0,ε风0,ε下0,ε吊0,ε系0) ,那么当f0max等于零或接近于零时,所得到的 ε主0、ε风0、ε下0、ε吊0、ε系0值便是各个索结构在成桥状态下的初始应变值,所得模型便是后续分析的理想模型。考虑到主索、风索、下稳定索、吊索、系索的耦合作用,找形计算[14]通过循环程序实现,找形结果见表2。
4 模拟分析与验证
在管桥正常运营的工况下对索结构体系进行几何非线性有限元分析。该管桥已运行15 年,各部分结构都存在着一定程度的腐蚀损耗,因此必须依据现场实测的结构现尺寸数据对管桥模型进行修正[15]。正常运营工况是指只考虑输气管道内输送的天然气介质重力载荷和管桥结构自重,不考虑其他外载荷。通过加载计算,得到了主索、风索、下稳定索、吊索、系索五种索结构的应力分布规律,并将模拟结果与工程检测[16]数据进行对比验证。其中主索、风索、下稳定索的实测数据来自于该悬索管桥的安全监测预警系统[17],每索9 个测点。吊索和系索由于数量较多,首先通过SL - 30T型钢丝绳张力测量仪现场逐一测得每对索的轴向张力,然后换算成轴向应力。另外,由于主索、风索、下稳定索在整个索结构体系中起着主要作用,而且各只有1 对,所以模拟分析它们整条索不同位置的应力变化情况,图为3 为主索轴向应力曲线,图4 为风索轴向应力曲线,图5 为下稳定索轴向应力曲线; 对于吊索和系索,因为数量众多,长短不一,且每对索所处的位置也不同,不考虑每对索具体的应力变化情况,而是把它们作为一个整体来研究,图6 为吊索体系应力分布图,图7 为系索体系应力分布图。
从以上各图可以看出,仿真模拟结果与现场实测数据吻合良好。其中主索的轴向应力呈“U”型分布,在管桥中点处最小,距离管桥中点越远,其值越大,最大值为290 MPa左右; 风索的轴向应力值也呈现中间低、两边高的走势,但变化较为平缓,最大值约80 MPa,最小值约70 MPa; 而下稳定索的轴向应力基本保持不变,为一个恒定值,仿真结果为94. 4MPa,现场实测结果稍大,稳定在100 MPa左右; 对于吊索体系,管桥两端的几对应力值明显高于其它吊索,且越接近桥端,其值越大,其中第1 对和第53对的轴向应力值达到130 MPa,剩余吊索的应力变化比较平缓,保持在55 MPa左右; 系索体系,每对索的轴向应力越接近管桥中心越高,但峰值并不是出现在中间的第27 对系索上,而是出现在中间附近两边的系索上。
5 结论
( 1) 悬索管桥各个索的应力值都相对较小,远小于其标准破坏强度1 570 MPa,说明该索结构体系确实是小应变、大位移的几何非线性结构。
( 2) 由于管桥结构的对称性,每条主索、风索、下稳定索的轴向应力沿桥中心呈对称分布,吊索体系、系索体系的各对索的轴向应力以各自的第27 对索为中心,也呈对称分布。
( 3) 系索的最大应力没有出现在第27 对索上,而是在其两侧附近的系索上,可能是受管桥桥面桁架微拱形结构的影响。
耦合非线性系统 篇7
目前,电能传输方式主要由导线通过插头插座直接接触进行传送。由于这种方式存在物理接触和电气接触,在诸如潮湿、易燃、易爆等环境中的应用受到限制,而且可靠性差。松耦合感应电能传输技术是一种基于电磁感应耦合理论、现代电力电子能量变换技术及控制理论于一体的新型能量传输技术。通过采用一、二次侧可分离的松散耦合变压器将电能从电源侧经气隙传递给一个或多个负载,从而安全、可靠、高效、灵活地实现了供电线路和用电设备之间无物理连接下的能量传输,在交通运输、航空航天、医疗、照明、矿井和水下场合有着广泛的应用前景[1-3]。
由于开关器件及非线性控制方法的存在,松耦合感应电能传输系统必然是一个高阶多参数变化的非线性系统。其电路中也必然存在着分岔、混沌和共存吸引子等非线性现象[4],这些非线性行为将会对此类系统的设计带来一定的影响。但是,目前的研究总是采用各种线性化的方法来分析系统,而忽视其非线性行为,其非线性动力学行为至今为止都没有得到深入探索。近二十年对电力电子变换器的分岔和混沌现象的研究已经非常深入,发现了不同类型电力电子变换器中所出现的分岔和混沌现象[5-7],并形成了一套行之有效的研究方法,如用数值迭代建立低维DC/DC变换器的离散模型[8-11],基于雅可比矩阵稳定性判断[12]和分析了功率因数校正(PFC)变换器及逆变器的失稳现象[12-13]。为此, 利用电力电子变换器的非线性分析方法分析松散耦合感应电能传输系统的非线性行为将有助于设计开发出更为稳定的变换器用于无线电能传输系统。
本文以全桥谐振变换器为研究对象,基于松散耦合变压器互感模型建立了电压外环控制的系统状态方程,分析了系统谐振时的非线性行为。运用分岔图、频闪映射图以及频谱图,描述了全桥谐振变换器随着耦合系数的改变而产生的非线性分岔现象。
1基于松散耦合变压器的全桥谐振变换器及其控制方法
松散耦合变压器是非接触电能传输系统的关键部分,与常规变压器相比,其原、副边之间存在较大的气隙,空气磁路长度远远超过了常规变压器的空气磁路长度,使之相当一部分磁动势消耗在空气磁路部分,变压器漏感较大,耦合系数不高。松散耦合变压器属于疏松耦合磁结构,这大大影响了系统能量传输的功率和效率。基于松散耦合变压器的全桥谐振变换器如图1所示。系统主电路原边由全桥逆变电路组成,副边由整流滤波电路组成,本文系统采用原副边串联谐振。图1中:E为输入直流电压源; S1,S2,S3,S4为金属氧化物半导体场效应晶体管(MOSFET)构成的全桥逆变电路;Cp和Lp分别为原边谐振电容和电感;M为原副边的松耦合互感; Cs和Ls分别为副边谐振电容和电感;Rp和Rs分别为原副边电感的寄生电阻;D1,D2,D3,D4为二极管构成的整流电路;L和C组成副边滤波电路;R为负载。采样副边输出电压vC通过无线通信技术传输到原边进行控制。参考电压vref减去通过比例环G1的vC再通过比例环G2形成控制电压vcon与锯齿波信号vramp比较生成控制开关管的脉宽调制(PWM) 信号控制系统。具体的控制脉冲波形如图2所示。 当vcon>vramp,u为正驱动信号,S1和S4导通,此时S2和S3截止;当vcon
采用Cp,Cs和C的电容电压vCp,vCs和vC,以及Lp,Ls和L的电感电流iLp,iLs和iL为系统的状态变量,即X=[vCp,vCs,vC,iLp,iLs,iL]T。由于此系统中8个开关器件两两组合工作并且根据此系统的控制特点,可知本系统共有如表1所示的6种开关状态。
系统的状态方程如下(具体变量说明见附录A):
以锯齿波周期采样来构造频闪映射(即庞加莱截面法),解式(1)可得:
式中:ti为开关状态i的工作时间;I为单位矩阵。
如表1所示,在开关状态5与6时,D1,D2,D3和D4都处于关断状态,此时iLs=iL=0,说明系统处于电流断续模式(DCM)中。在开关状态1,2,3,4中,iLs和iL都不等于0,说明系统处于电流连续模式(CCM)中。令系统在一个开关周期中经历的开关状态数为j,如表1所示,各个开关状态可以通过vcon与vramp的关系以及iLs的大小区分计算。令各个开关模式的工作时间分别为t1,t2,…,tj,则可由式(2)得到一个从Xn到Xn+1的离散映射:
Xn+1=fj(fj-1(…f1(Xn,t1),…,tj-1),tj)(3) 自此,可利用式(3)对系统进行理论建模分析。
2基于松散耦合变压器的全桥谐振变换器谐振条件及其分岔行为
为了简化分析,假设原副边的磁芯一致、绕线方法一致,原副边的匝数相同,则原副边的自感为Lp=Ls,系统的耦合系数为:
本系统采用原边与副边同时串联谐振补偿,因此,原边及副边的谐振频率分别为:
根据文献[14],只有开关频率及原副边的谐振频率都到一致时才能达到系统谐振,取Cp=Cs。令系统的部分参数如下:E=30V,Cp=Cs=2.2μF, Rp=Rs=0.1Ω,Lp=Ls=116μH,L=220μH,C= 220μF,R=10Ω,G1=0.5,G2=0.4,vref=20V,vU=8V,vL=2V和k=0.85。则原副边的谐振频率为fp=fs=10 000 Hz。分别取开关频率f=1/T为5 000,10 000,20 000Hz,运用式(3)进行离散迭代映射,取系统进行稳态后的10 000点画入图3中。图中:n为迭代次数;vCn为第n次迭代后的副边输出电压。由图3的稳态迭代输出电压可知,3种开关频率下,系统都处于周期一轨道中,且系统达到全谐振f=10 000Hz时输出电压的增益最大。
利用以上的系统参数,并取系统进入全谐振的频率为开关频率,以耦合系数k∈(0.05,0.80)为分岔参数,得到图4所示的系统分岔图。图中,iLn为第n次迭代后的副边滤波电感电流。为了更清晰地分析,图5显示了系统典型k值下的采样点相图,即庞加莱截面图。
从图4和图5可知,当k属于0.65至0.80时, 系统处于周期一轨道中,因为迭代采样点为一个点(图5(a)),此时系统在一个开关周期中的开关状态序列为1→3→4→2。减小k至0.64时(此处可通过改变松散耦合变压器的相对位置减小k的值),系统发生了Hopf分岔进入低频振荡的准周期轨道,采样点为一环面(图5(b))。继续减小到k至0.52时,L进入DCM,此时系统发生边界碰撞导致环面不断增大且不再光滑(图5(c)),即出现环面破裂现象。减小k环面继续增大,当k为0.15时,环面达到最大(图5(d))。此后,随着k的减小,系统又进入CCM,并且环面不断减小(图5(e))。 图5(f)显示,当k减小至0.08时,系统又进入周期一轨道。
3基于松散耦合变压器的全桥谐振变换器电路仿真验证
根据图1所示的电路图和第2节的参数对全桥谐振变换器进行电路仿真。
典型的vC和iL波形、iL和vC相图及vC谐波频谱分析分别如图6、图7、图8、图9和图10所示。 当k=0.80时,电路仿真系统处于CCM周期一轨道(开关状态序列为1→3→4→2),此时的vC谐波频谱主频率集中在开关频率f=10 000 Hz处(图6)。随着耦合系数k的减小,系统进入准周期轨道,并且系统在一个开关周期中的开关状态序列依然为1→3→4→2,vC与iL的波形呈低频振荡(图7)。此时,vC谐波频谱出现了2 000 Hz的频点,且开关频率的倍频幅值也较大,如图7(d)所示。
继续减小k的值,随后系统在部分开关周期中进入DCM,其部分开关状态序列为1→5→3→4→ 2。随着k的减小,系统的低频振荡也越发严重,k= 0.11时振荡达到最大(图8),其部分开关状态序列为6→4→6→2→6和5→6→4→6→2→6,此时出现了严重的环面破裂。图8中vC与iL的谐波幅值都非常大,分别达到6V和8.5A,致使环面增大,并且谐波处于远低于开关频率的低频段。与上节的数值分析一样,当k继续减小时,系统又进入到CCM的准周期轨道,其开关状态序列在一个开关周期中重新回归到1→3→4→2。图9显示,k=0.05的系统准周期低频振荡与较大k值的准周期低频振荡最大的不同点在于,谐波频率集中于低频段且不存在开关频率的谐波。随着k的继续减小,系统又发生一系列的边界碰撞现象重新进入周期一轨道。在一个开关周期中,系统的开关状态序列先经过1→2→ 4→2,而后又重新回到1→3→4→2(图10)。图10 (d)所示的vC稳态时的谐波频谱显示其基频为f= 2 000kHz,说明此周期一具有2倍于开关频率的谐波,这是由于系统随着耦合系数的减小而使漏感增大及全桥变换器的对称拓扑结构所引起的。
此处分岔参数的具体状态轨道与数值分析的结果相比有所偏差,但得到的趋势一致。这是由于电路仿真中系统考虑了更多的寄生电阻、电感和电容等参数所致。
4结语
感应式无线电能传输系统在工程实践中得到了广泛的应用,以往研究总是采用各种线性化的方法来分析系统,而忽视其非线性行为,其非线性动力学行为至今为止都没有得到深入探索。
本文研究了电压外环控制下基于松散耦合变压器的全桥谐振变换器的非线性动力学行为,建立了系统的非线性模型;指出该系统随着耦合参数的变化存在Hopf分岔现象,分析了分岔过程系统发生的一系列的边界碰撞现象,并进一步分析了稳定状态及准周期状态下系统的时域和频域特性。本文将开关变换器的非线性研究拓展到无线电能传输系统中,为深入研究其非线性行为建立了平台,对无线电能传输系统的可靠设计也有重要指导意义。
耦合非线性系统 篇8
1 基本方程
1.1 激发方程
令第次本征模式轴向电场为
可得扰动电子注激励的电场为[5,6,7]:
式中:定义Kn=为n次谐波的耦合阻抗;(r⊥)为电子注横向分布函数ds。所有本征模中只有个别模式与电子注同步,且除了电子流i(z)激发的同步波之外,还有输入的“冷波”,即,则具有外加激励源的同步场为:
式中:K0,Γ0,β0分别为该同步模式的耦合阻抗、传播常数、相位常数。式(1)两边同时对z求导2次得到熟知的激发方程:
1.2 运动方程
相对论下电子的运动方程为:
能量守恒方程为:
式(4)代入式(3)可得一维电子运动方程:
又由:
所以式(5)可写为:
其中:Ez=Ecz(线路场)+Esz(空间电荷场);γ=(1-v2/c2)-1/2为相对论因子,c为真空中的光速,v为电子速度。
1.3 电子流复振幅方程
电子流是时间的非简谐周期函数,含有高次谐波,用傅氏分析。
2 慢变系统中归一化
上述是在实验室坐标系下得到的迅变方程,为了处理问题的方便和计算结果普遍性强,通常将其归一化到电子坐标系内,获得慢变方程。
为了表述方便,先引入迅变坐标系的归一化量:归一化距离为ξ=ω/v0z=βez;归一化时间为φ=ωt=2πt/T,归一化场为f=|e|E/mv0ω。则慢变系统中的归一化:归一化轴向距离为θ=Cξ=Cβez;归一化相位Υe=φ-ξ;归一化场幅值为Fcn(θ)=;归一化电流幅值为。
激发方程变为:
运动方程变为:
电子流复振幅方程为:
式(8)～(10)组成了行波管大信号注-波互作用基本工作方程组。其中=I0Kcn/4V0为n次谐波增益参量,非同步参量,dn=α0n/βeCn为衰减常数,rn=bn-idn。
3 空间电荷场的计算
由文献[8]可知z0处圆盘在z处圆盘平面上各点产生的平均空间电荷场为:
其中:Q为圆盘所带电量;b为电子注半径;a为漂移管半径,如图1所示,μ0n为零阶Bessel函数的第n个根。由此可知场是关于z的函数,可以表示为:
其中:B(|z-z0|)是以|z-z0|为变量的函数,由式(11)可以做出如图1所示曲线。
由图1可知,如果用近似式Esz=QB(|z-z0|)代替式(11),所引起的误差很小,而计算式却大大简化了。采用该式时,所有圆盘产生的空间电后为:
其中:ωp为电子等离子体频率;η为电子荷质比;v0为电子注初速度。仿真中用式(13)代替式(9)中的Fsz。
4 休斯结构耦合腔
休斯型耦合腔的剖面结构为图2(a)所示,图2(b)为其截面图,表1为设计S波段耦合腔行波管的结构参数。
5 仿真结果及讨论
在上述尺寸下,对2.87～3.35 GHz频率范围内的休斯结构耦合腔行波管进行了数值分析,仿真中电子注的注电压U0=21 kV,注电流I0=1 A,波的输入功率Pin=300 W。采用四阶龙格-库塔法[9]求解互作用方程组式(8)～(10),其结果如图3～图6所示。
图给出了在中心频率时电子的相位轨迹。由图可知,电子在归一化轴向距离Z=2.4附近获得了较好的群聚,此即为最佳互作用距离。由图4可以看出该频率的波在此处获得最大增益,此后电子注离开最佳互作用区,效率降低,增益下降。
图5给出了中心频率时电子注效率随轴向距离的变化曲线。图中实线为不考虑空间电荷力的情况,虚线为考虑空间电荷力的情况。由图可知,空间电荷力的作用使得饱和位置推后,增益下降,即空间电荷力起发散作用。图中效率刚开始时为负值,是因为电子刚进入互作用区时要得到部分能量,表现电子效率为负值。
图6给出考虑空间电荷场时,2.87～3.35 GHz内各频点的微波增益。由图可知,在给定电子注注电压,注电流,和波的输入功率等参量的情况下,频带内微波增益均大于18.5 dB。
6 结语
对S波段休斯结构耦合腔行波管非线性注-波互作用工作方程组进行了数值求解,求出考虑和不考虑空间电荷场时中心频率f=3.100 7 GHz处的效率,说明空间电荷场对互作用起散焦作用,与文献[2,10]中结论很吻合。求出工作在2.87～3.35 GHz频率范围内耦合腔行波瞬时带宽。仿真中所用管子已制作完成,实验数据对后期管子的热测试有很好的指导意义。
参考文献
[1]李文君,许州,林郁正,等.S波段宽带大功率连续波耦合腔行波管三维模拟设计[J].强激光与粒子束,2006,18(6):965-968.
[2]徐林,杨中海,莫元龙.行波管三维非线性计算机模拟的改进[J].强激光与粒子束,1999,11(3):342-346.
[3]雷文强,杨中海,廖平.耦合腔慢波特性及其注-波互作用的三维软件模拟[J].强激光与粒子束,2003,15(7):673-676.
[4]Kantrowitx F,Tammaru I.Three-dimensional Simulations ofFrequency-phase Measurements of Arbitrary Coupled-cavityRF Circuit[J].IEEE Trans.on Electron Devices,1988,35(11):2 018-2 026.
[5]Vaninstin L A.Waveguide Excited by the Electron Beam[J].Journal of Technology Physics,1953(4):654-657.
[6]李建清,江丽军,莫元龙.等离子体-腔混合模耦合腔行波管非线性注-波互作用分析[J].强激光与粒子束,2004,16(11):1 429-1 433.
[7]张昭洪,吴鸿适.耦合腔慢波结构特性研究[J].电子学报,1986,14(1):7-8.
[8]刘盛刚.微波电子学导论[M].北京:国防工业出版社,1985.
[9]John H Mathews,Kurtis D Fink.数值方法[M].周璐,译.北京:电子工业出版社,2005:399-413.