非线性行为（共7篇）
非线性行为 篇1
摘要：提出了一种利用Legendre正交多项式简化Volterra级数模型的建模方法(Legendre Volterra Series,LVS)。传统Volter-ra模型具有高精度和物理意义明确的优点,但其缺点在于计算量大,对硬件资源的消耗比较大,而且不利于对信号进行实时处理,该模型很大程度上减少了运算量,同时也保持了传统Volterra模型的优点。
关键词：LVS,功放非线性,Volterra,建模
0 引言
在广播数字电视的发射系统中,高功率放大器是其中的关键模块,其性能的优劣对通信的质量与效率有着重要的影响,然而高功率放大器的非线性失真难以避免[1]。目前一些新兴的传输格式,特别是地面数字电视体制,如国家标准DTMB,信号经OFDM调制后峰均比比较高,非线性失真愈加严重,因此需要通过一定的线性化技术校正。而建立一个有效而简洁的非线性数学模型正是功率放大器线性化技术的基础[2],具有重要的研究价值。
Volterra级数模型是目前最能体现功放记忆效应的模型[3],能够清晰准确地对系统物理意义进行描述,并能以任意精度逼近连续函数[4],但随着系统阶数的增加,功放的记忆效应和非线性较强时,Volterra级数的维数呈指数快速增长,使得辨识过程出现了维数灾难,以至计算量大大增加,影响模型的实际应用。因此目前对Volterra级数模型的研究大多集中在改进结构和克服维数灾难的问题上[5]。
本文从降低维数灾难的角度出发,提出了一种通过特殊函数(Legendre正交多项式),对通用Volterra级数模型结构做出改进,得到了基于Legendre正交函数的数学表达式和简化模型(LVS模型),从而有效降低了模型的辨识维数,达到了简化目的,并通过仿真说明该模型的优越性。
1 传统Volterra级数模型
高功率放大器非线性系统的输入输出关系可表示为
式中:yk[x(t)]是非线性系统的第k阶输出,hk(τ1τ2…τk)称为k阶Volterra核或k阶冲激响应[6]。从而可推出,Volterra级数的离散截断形式为
一般情况下,只有奇数阶互调频率分量可以落在通带以内,所以可以忽略偶数项的影响,从工程应用角度出发,考虑到高功率放大器非线性系统中通常采用复包络信号作为输入,因此输入输出信号的形式分别为
可以推导出离散有记忆奇数阶复数形式Volterra级数模型的表达形式为
式中:Q和K分别为非线性的记忆深度和截取的阶数。
2 LVS模型的建立与推导
为了避免应用Volterra级数建模中普遍存在的维数灾难及由之引发的繁重计算[7],可以采用通过Legendre正交多项式,以减少Volterra级数中核的系数数量方法,对通用Volterra级数模型结构做出近似修正。由式(5)可以推出,其中包含的每个相加项离散形式的多维分式为
对式(6)进行傅里叶变换可得
由式(7)可得,其中每一个相加项的傅里叶变换为
由式(8)可以得到
式(9)中引入的两类新的变量意义如下
根据Legendre多项式的数学形式,可以设其变量为ω/B,于是其完整代数表示形式为
得到其递推公式为
经过推导计算,可以得出ui(t)和ωi(t)两者之间的相互替代关系,进而可以得到LVS模型表达式为
由此可见,Volterra级数模型在Legendre多项式的辅助下得到了简化。在LVS模型中,系数数量得到有效减少,达到了预期提高计算效率的目的。
3 仿真过程与分析
为了验证模型的有效性,本文对LVS模型进行了仿真。系统框图如图1所示。
地面数字电视国标信号参数设置为:星座映射为16QAM,FEC码率为0.8,载波为3 780,帧头为420,作为输入信号x(k)(频谱如图2所示),分别通过传统Volterra级数模型和LVS模型,得到输出信号y1(k),y2(k)的频谱如图3和图4所示。
表1是两模型在记忆长度为3,阶数为6时,与实测数据d(k)的比较。
从频谱和表1中可以看出,新的LVS模型很好地实现了原有模型的功能,符合数字电视传输上对功放的基本要求,模拟出了实际功放的频谱特性,并且由测试得出,当记忆长度长,阶数多的时候,模型简化效果更为明显。
4 小结
本文通过Legendre正交多项式,对传统Volterra级数模型进行了近似和简化,得到改进后的LVS模型。仿真结果表明,LVS模型在误差允许范围内,保持了传统Volterra模型的优点,并且很好地描述功率放大器的非线性特性及记忆效应,有效地减少了计算量,降低了维数灾难,减少了硬件资源消耗,为高效的功放线性化提供了有力的保证。
参考文献
[1]ISAKSSON M,WISELL D,RONNOW D.A comparative analysis of be-havioral models for RF power amplifiers[J].IEEE Trans.MicrowaveTheory and Techniques,2006,54(1):348-359.
[2]郑百衡.一种基于Volterra级数的基带数字预失真[J].电视技术,2010,34(S1):28-29.
[3]张涛涛,唐世刚,潘长勇,等.基于多项式的功率放大器预失真技术[J].电视技术,2006,30(11):51-53.
[4]ZHU A,PEDRO J C,BRAZIL T J.Dynamic deviation reductionbasedVolterra behavioral modeling of RF power amplifiers[J].IEEE Trans.Microwave Theory and Techniques,2006,54(12):4323-4332.
[5]CRESPO-CADENAS C,REINA-TOSINA J,MADERO-AYORA M J,etal.A new approach to pruning Volterra models for power amplifiers[J].IEEE Trans.Signal Processing,2010,58(4):2113-2120.
[6]ZHU A,PEDRO J C,CUNHA T R.Pruning the volterra series for behav-ioral modeling of power amplifiers using physical knowledge[J].IEEETrans.Microwave Theory and Techniques,2007,55(5):813-821.
[7]ZHU A,BRAZIL T J.Behavioral modeling of RF power amplifiers basedon pruned volterra series[J].IEEE Trans.Microwave and WirelessComponents Letters,2004,14(12):563-565.
非线性行为 篇2
覆冰导线受风力作用产生的低频、大幅自激振动称为舞动, 舞动对线路安全运行有极大威胁, 有时威胁甚至是灾难性的, 舞动机理及防舞技术研究方兴未艾[1,2]。就舞动机理而言, 普遍为人所接受的是Den Hartog[3]的垂直舞动机理, Nigol等[4]的扭转舞动机理和Yu等[5]提出的偏心惯性覆冰导线舞动过程中产生的是交变张力, 就连续体模型等效为具有单自由度或多自由度耦合的集中参数模型, 但该模型无法考虑导线在重力作用下具有一定垂度和柔性及运动大位移、小应变的几何强非线性特性。
本文基于有关索的建模方法[6]建立了覆冰导线风致振动连续体非线性动力学偏微分方程, 在此基础上, 通过Galerkin积分得到了便于解析求解和数值模拟的常微分方程, 同时探讨了结构参数如导线跨度、水平张力等对各阶模态线性固有频率的影响。应用Mathematica程序数值模拟导线振动特性并与有限元分析结果进行比较, 比较结果表明本文连续体模型能准确地呈现系统的固有频率特性及振动幅频特性。
1 覆冰导线连续体动力学模型
考虑到受几何非线性和垂度影响的覆冰导线舞动问题很复杂, 本文在建立分布参数连续体非线性动力学模型时作了如下基本假设:基于准定常气动力假设, 不考虑空气湍流;导线在重力作用下的垂度曲线为抛物线;不考虑次档距振荡, 将分裂导线视为一根横向振动类似索, 其扭转振动类似具有一定刚度的梁的振动;不考虑导线安装初伸长, 本构关系服从胡克定律;不考虑导线轴向运动。
截取一微导线单元如图1所示, 其中, oxy面为垂直振动 (y向) 、oyz面为水平振动 (z向) , 绕导线轴线 (分裂导线绕分裂圆圆心) 为扭振。
导线微元段的动力学方程为
其中, T、t′为覆冰导线的初始张力和舞动过程中动张力增量;和l分别为导线在静平衡状态下的垂度曲线和档距;v (x, t) 、w (x, t) 和θ (x, t) 为导线舞动过程中相对于静平衡位置沿y、z和扭转方向的位移;fy (x, t) 、fz (x, t) 和fθ (x, t) 分别为单位长度覆冰四分裂导线在y、z和扭转方向所受到的气动力;m为单位长度覆冰导线质量;s和ds分别为弧长坐标和微元长度;Ic、kθ分别为分裂导线对旋转中心的转动惯量和扭转刚度;e、θ0分别为覆冰导线截面偏心距和覆冰角;cy=2mξyωy、cz=2mξzωz和cθ=2Icξθωθ分别为单位长度导线在y、z和扭转方向的线性结构阻尼系数;ξy、ξz、ξθ分别为各向阻尼比;ωy、ωz、ωθ分别为各向振动圆频率。
为简化式 (1) , 引入下列变换[6]:
由于分裂导线振动过程中扭转角很小, 本文取
将上述各式代入式 (1) , 得
其中, yx为y关于x的一阶导数, H、h为覆冰导线水平张力和水平张力增量。y、z和扭转方向的振型函数可分别取为φ1i (x) =sin (iπx/l) 、φ2i (x) =sin (iπx/l) 和φ3i (x) =sin (iπx/l) (i=1, 2, 3) 。由于舞动波形主要为低阶谐波[7], 本文以1阶模态 (i=1) 为例对振动方程进行模态截断, 令v (x, t) =φ1i (x) q1 (t) 、w (x, t) =φ2i (x) q2 (t) 和θ (x, t) =φ3i (x) q3 (t) , 其中q1 (t) 、q2 (t) 和q3 (t) 为描述导线上各点振动位移的时间函数, 经Galerkin积分可得非线性动力学常微分方程如下:
式中, a11, a12, …, a18;a21, a22, …, a27;a31, a32, …, a35分别为通过Galerkin积分得到的垂直、水平和扭转方向覆冰导线各项系数。
整档覆冰导线所受气动力为
其中, U为水平风速;ρ为空气密度;D为裸导线直径;Cy (α) 、Cz (α) 和Cθ (α) 分别为作用于三个方向上气动力的系数, 是气动攻角α的函数。取气动攻角表达式为
其中, θ0为由导线不均匀覆冰而具有的初始攻角 (覆冰角) , 分别为由扭转及垂直振动引起的攻角变化量。气动力参数取自文献[8], 平均升力、阻力和扭矩系数关于α的函数式为
将式 (7) 代入式 (5) 并积分, 结果代入式 (4) 可得耦合很强烈的动力学常微分方程组:
其中, b11, b12, …, b19;b21, b22, …, b29;b31, b32, …, b39分别为通过Galerkin积分得到的垂直、水平和扭转方向覆冰导线所受各项气动力的系数。前3阶模态振动方程形式相似, 此耦合动力学方程能够全面地反应连续体系统的振动特性。
2 覆冰分裂导线线性固有频率
覆冰导线低频、大幅舞动是典型的流固耦合弹性自激振动, 振动频率接近某一低阶频率[1]。从振动响应特性来看, 内共振与非内共振有本质区别。导线是具有一定垂度和柔性的三维弹性连续体且有无限多个固有频率及相应模态, 当两个或更多个频率的比值为有理数或接近有理数时就会发生内共振[8]。
对于式 (4) 来说, 横向 (垂直、水平) 和扭转振动各阶模态的线性固有频率均是分析振动原因和振动特性的重要依据。垂直、水平和扭转一阶振动的圆频率分别为
其中, a13、a23、a33为各向振动1阶模态刚度、a11、a21、a31为各向振动1阶模态质量 (扭转为转动惯量) , 是与结构参数如导线自重、跨度、水平张力等有关的积分常数 (积分式略) 。单根覆冰导线参数取自文献[9], 见表1。
2.1 横向振动线性固有频率
两支撑端等高线档前3阶模态示意图见图2, 实线为导线初始静平衡位形, 虚线为振型曲线。
覆冰导线舞动过程中产生的是交变张力, 就舞动特性而言, 水平张力及张力差起控制作用。
由文献[10]知, 振动过程中水平张力变化量为
其中, E为弹性模量, A为导线截面面积, n为阶数, 导线垂直振动前3阶模态中由于自重引起的水平张力变化量分别为
由式 (11) 可知, 奇数阶中重力引起水平张力增量与阶数成反比, 偶数阶中水平张力增量为零。结合图2可知, 1阶模态导线各点离开静平衡位置向同侧运动, 要克服重力恢复到静平衡位置势必会增大张力;2阶模态档距中点两侧导线反向运动水平张力增量抵消为零;3阶模态中间与两侧导线振动反向, 水平张力增量由1/3档导线重力引起。因此, 在计算垂直振动各阶模态线性刚度时, 偶数阶不考虑自重, 奇数阶考虑自重的1/n。
在跨度一定的工况下, 水平张力H对各阶模态线性刚度影响显著。应用Mathematica程序计算并拟合频率随水平张力变化关系, 结果如图3所示。
图3a中, 垂直振动1、3阶模态线性固有频率先随水平张力增大而减小, 经最低值后又随张力增大而增大, 2阶模态线性固有频率始终随张力增大而增大;图3b中, 水平振动前3阶模态线性固有频率均随水平张力增大而增大。
选取l=126m、H=30kN导线来计算其圆频率, 并建立有限元模型数值模拟以进行比较, 结果见表2、表3。
由表2、表3可见, 连续体模型计算导线线性固有频率与有限元分析结果非常吻合。
2.2 扭转振动线性固有频率
计算分裂导线中点等效扭转刚度的方法只有文献[1]中提到的张紧弦法和柔索法 (三分裂) 两种。本文推导了应用柔索法计算二分裂、四分裂导线中点的扭转刚度公式, 列于表4。
表中, D1为分裂圆直径。
对于l=126m, H=30kN四分裂导线, 其扭振各阶模态线性固有频率见表5。由表5可见, 对于小档距浅弧垂导线, 张紧弦法与有限元结果吻合。
3 数值模拟研究
受非线性动力系统理论发展及研究方法制约, 对于式 (8) 这样复杂的方程求解析解用于定量分析很困难, 可用数值模拟方法研究系统振动特性。跨度l=126m的覆冰导线, 当风速U=12m/s时档距中点的运动轨迹图和时程曲线如图4所示。
由Routh-Hurwitz判据可知, 系统属于Den Hartog垂直舞动激发模式。水平和扭转振动为垂直振动通过惯性耦合将动力传递到相应方向引起的受迫振动。图4a表明导线中点运动轨迹为椭圆, 图4b中q1、q2和q3分别为垂直、水平和扭转振动历程曲线。通过分析知, 三向振动一阶频率为, 小于同阶垂直振动固有频率fy1=0.608Hz, 且二者差5.48%。表明非线性因素使振幅对系统振动频率产生了影响, 即系统的频率失去“等时性”。
图5为有限元数值模拟得到的导线垂直振动 (y向) 频谱图[11], 所选参数与本文相同。垂直振动频率与连续体模型分析结果0.575Hz相对误差仅为0.87%, 且图5中振幅y值与图4中表示垂直振动位移的q1也十分吻合。
4 结论
本文建立了覆冰分裂导线垂直、水平和扭转振动耦合的连续体动力学方程, 在此基础上深入探讨了结构参数对各阶模态线性固有频率的影响, 结果表明, 导线水平张力和跨度对频率影响显著, 且在小张力范围内自重对垂直振动频率影响大。对耦合动力学方程进行了数值模拟, 模拟结果表明各向振动均收敛于稳定的极限环。1阶模态振动频率小于该阶模态的线性固有频率。上述分析结果与有限元结果吻合, 表明连续体模型能很好地反应振动系统结构和动力学特性。
参考文献
[1]郭应龙, 李国兴, 尤传永.输电线路舞动[M].北京:中国电力出版社, 2003.
[2]王吉岱, 连金玲.气动式除冰机器人及其运动学分析[J].中国机械工程, 2013, 24 (5) :610-613.Wang Jidai, Lian Jinling.Kinematics Analysis of a Pneumatic De-icing Robot[J].China Mechanical Engineering, 2013, 24 (5) :610-613.
[3]Den Hartog J P.Mechanical Vibrations[M].New York:McGraw-Hill, 1956.
[4]Nigol O, Buchan P G.Conductor Galloping:1.Torsion Mechanism[J].IEEE Transactions on Power Apparatus and Systems, 1981, 100 (2) :699-707.
[5]Yu P, Shah A H, Popplewell N.Inertially Coupled Galloping of Iced Conductors[J].Journal of Applied Mechanics, 1992, 59 (1) :140-145.
[6]李伟义.连续体斜拉索风雨激振的非线性研究[D].天津:天津大学, 2008.
[7]李碧辉, 田丰.浅谈线路覆冰的危害及防治[J].江西电力, 2008, 32 (1) :47-50.Li Bihui, Tian Feng.On the Damage and Control of Transmission Line Icing[J].Jiangxi Electric, 2008, 32 (1) :47-50.
[8]王洪礼, 张琪昌.非线性动力学理论与应用[M].天津:天津科学技术出版社, 2002.
[9]Zhang Q, Popplewell N, Shah A H.Galloping of Bundle Conductor[J].Journal of Sound and Vibration, 2000, 234 (1) :115-134.
[10]杨风利, 杨靖波, 付东杰, 等.输电线路导线舞动载荷分析[J].中国电机工程学报, 2011, 31 (16) :102-107.Yang Fengli, Yang Jingbo, Fu Dongjie, et al.Analysis on the Laods from Galloping Conductors of Transmission Lines[J].Proceedings of the CSEE, 2011, 31 (16) :102-107.
非线性行为 篇3
电液伺服系统执行机构在运行过程中，在油源压力和负载压力的作用下，可因油液的可压缩性而形成动态液压弹簧[1]。弹簧刚度的非线性会使运动过程中系统的固有频率不恒定、响应稳定区域变得复杂。液压弹簧与负载质量相互作用可构成一个液压弹簧-质量系统。该系统在一定条件下会引起伺服系统发生非线性振动。因此，液压弹簧力对电液伺服系统动态特征的影响值得关注。
目前，对液压系统动态特性的研究一般采用系统建模和数值仿真方法[2,3]。系统建模时一般对非线性因素进行线性化处理[4,5]，研究所得结论与实际情况有较大差异，很难解释实际动态测试中出现的时域波形复杂、频域尖峰繁多等异常现象[6]。所依据的理论多是线性动力学理论和经典控制理论[7,8]，而较少运用非线性动力学理论[9,10]。故针对非线性液压弹簧力作用下电液伺服系统动态特征的研究尚不多见。
本文以电液伺服系统为研究对象，重点探究非线性液压弹簧力对系统动态特征的影响规律。根据非线性动力学原理，建立系统的非线性动力学模型。通过数值分析，揭示系统内在的分岔现象及典型非线性动力学行为。用非线性动力学研究方法对实测动态数据进行深入分析，以揭示液压弹簧的软硬弹簧特性引起的“跳跃现象”。旨在揭示伺服系统非线性振动的机理及诱因，使综合分析系统的动态特征变得更接近实际。
1 执行机构的动力学模型
电液伺服系统的执行机构为伺服液压缸，本文以双作用单活塞杆液压缸为例进行分析，其工作原理如图1所示。
上述执行机构的动力学方程为
式中，m为活塞及惯性负载的折合质量；x为活塞位移；Fc为黏性力；Fs为弹性力；Ff为摩擦力；FL为负载力；p1、p2分别为无杆腔和有杆腔的压力；A1、A2分别为无杆腔和有杆腔的活塞有效作用面积。
2 非线性液压弹簧力
电液伺服系统执行机构的弹簧刚度由活塞杆刚度和液压油刚度串联合成。钢的体积弹性模量是液压油的100～150倍，故可以把活塞杆作为刚体处理。因此，系统的弹簧力主要由受控液压油所构成的液体弹簧产生[11]。
液压缸活塞的移动会导致其两侧液体弹簧长度的变化，进而引起液压弹簧刚度的变化，其变化规律为[12]
式中，K为油液体积弹性模量；L为液压缸总行程；L1为活塞初始位置，即无杆腔液柱长度；VL1为阀与无杆腔之间管道内油液体积；VL2为阀与有杆腔之间管道内油液体积。
令y为在工作点x附近的振动位移，即y=Δx。由泰勒级数可知，非线性弹簧刚度在工作点附近可表达为
记，并将其代入式（3），得
略去式（4）中的高阶无穷小项o(y2），则液压缸系统的液压弹簧力可以表示为
由于弹簧弹性势能U具有对称性，可以表示为
故液压弹簧力可以进一步表达为
式（7）中，k3<0表示软弹簧特性；k3>0表示硬弹簧特性；k3=0为线性弹簧特性[13]。
3 液压弹簧力非线性动态特征
3.1 非线性动力学模型
本文为集中研究非线性液压弹簧力对系统动态特征的影响，暂不考虑摩擦力、系统阻尼等非线性因素。则系统方程（式（1））在工作点x附近的特性可表达为
式中，c0为结构阻尼系数；c1为线性摩擦阻尼系数；Ff(v）为工作点处的摩擦力；v为活塞的移动速度。
进一步整理，得
式中，c为线性阻尼系数，c=c0+c1。
由于油源压力脉动、阀口流量-压力非线性等因素的影响，导致进油压力有微观波动，服从简谐振动规律，式（9）右边的输入项可近似表示为Fsin（ωt+φ0），是系统的激振源[14]。其中，F为激振力；ω为激振角频率；φ0为激振力的初相角。
据上述分析，系统动力学方程（式（9））可写为
由Duffing方程的结构形式可知，式（10）是含有阻尼的Duffing方程，可为研究电液伺服系统的非线性液压弹簧力的动态特征提供结构模型。
3.2 求解方程
把液压弹簧力的非线性动态特征归结为Duffing方程，就可以通过借助Duffing方程的特性来揭示系统内在的基本规律。
为便于求解计算，将式（10）化为如下形式：
式中，ξ为线性阻尼比；ω0为非线性项系数β=0时线性简谐振子的自然频率；F0为单位质量所受的激振力幅值。
下面采用非线性动力学中的定量分析法———谐波平衡法[15]求解式（11），令
取线性谐振子在谐波激励下的稳定解作为系统的一种形式解，即
式中，A为零次近似解的振幅。
由三角函数恒等变换公式知
将式（12）～式（14）代入式（11），得
略去高次谐波项，并使式（15）等号两端sinωt、cosωt项的系数分别相等，得
将式（16）、式（17）分别平方后相加，得幅频关系式：
进一步求解式（18），可得
将式（16）、式（17）相除，得相频关系式：
3.3 解分析
3.3.1 尾部弯曲曲线
由幅频关系式可得系统方程（式（10））的幅频特性曲线（图2），β>0时，幅频特性曲线为尾部右偏曲线；β<0时，幅频特性曲线为尾部左偏曲线；阻尼的作用限制了共振振幅的无限上升。
3.3.2 跳跃现象
对于线性系统的受迫振动来说，激振频率的连续变化只会导致响应幅值的连续变化。但是，对于非线性系统，即使激振频率连续变化，在某些特定点上也会发生振幅突跳现象。
从图2可以看出，系统方程（式（10））的幅频特性曲线并非单值曲线。在某些区间内，同一频率对应3个不同的振幅。当激励频率连续变化时，会发生振幅突然变化的“跳跃现象”[15]。随着系统参数的变化，系统的运动状态发生突变的现象称为“动态分岔”。“跳跃现象”是非线性系统所特有的现象之一，它是一种特殊的动态分岔现象。
3.3.3 多重定态
在没有外加周期力扰动时，系统方程（式（10））的状态方程可以表示为
在研究动力系统状态随时间变化的规律时，有一类状态———“定态”具有特殊重要的意义，它是指所有状态变量对时间的导数全都等于零时的状态，即
定态在相空间中的代表点称为“定点”或“平衡点”。由于非线性的存在，系统在运动过程中往往会出现“多重定态”或“多重定点”。
如图3所示，系统方程在相平面上有3个平衡点：鞍点S(0,0）、稳定定点。若初始条件不同，则系统在不同流域中的轨线将趋于不同的稳定定点，图3中阴影区中的轨线将流向F1，非阴影区中的轨线将流向F2。当外加周期力F≠0时，系统很可能穿越不同流域的分界线，在不同流域之间来回跳动，形成复杂的振荡状态，从而使系统可能在原来那些稳定定点的周围或不同定点之间做各种复杂的运动[13]。
4 数值试验
为了探索液压弹簧力非线性项系数β和外加激振力F0对系统动态特征的影响，以系统方程（式（10））的具体算例：
进行数值试验研究。
4.1 分岔特性研究
激振力F0取不同值时，以β为分岔参数绘制分岔图。图4中的横轴为液压弹簧力非线性项系数β，纵轴为振动位移y。由图4可知，参数β、F0取不同值时，系统发生了不同程度的分岔：(1)系统方程存在单根、多根和无穷多个根时，在分岔图上表现为单值曲线、多值曲线和涂黑区等不同区段（分别对应于单周期、多周期和混沌等不同运动状态）。(2)解曲线在某些点处会发生中断和跳跃，说明随着参数的变化，系统会发生振幅突然变化的“跳跃现象”。(3)随着参数的变化，系统会发生运动状态突然变化的动态分岔现象。由周期运动进入混沌运动主要是通过倍周期分岔途径实现的。
4.2 运动形态仿真
为了形象地体现系统在不同参数下的运动形态，在MATLAB中建立仿真模型，对系统典型的非线性动力学行为进行仿真。仿真中采用Runge-Kutta算法，采样频率选100Hz（远大于外控力频率fp=ω/（2π）=1/（2π）=0.16Hz），终了时间为1000s。
β=0.2N/（mm·kg),F0=0.2N/kg时，仿真结果如图5所示。由图5可知，时间历程呈周期重复；功率谱在基频fp及其倍频处出现尖峰；相轨迹在有限的区域内重复，呈封闭曲线，即有极限环存在；庞加莱图在一定的区域内只有1个孤立点存在。这是明显的单周期运动特征表现，说明此时系统处于极限环型振荡状态。
β=0.5N/（mm·kg),F0=0.2N/kg时，仿真结果如图6所示。由图6可知，时间历程呈周期重复；功率谱在分频fp/2和它的倍频处存在尖峰；相轨迹在有限的区域内重复，呈封闭曲线；庞加莱图在一定的区域上有2个孤立点存在，说明此时系统处于2倍周期运动状态。
β=0.4N/（mm·kg),F0=0.4N/kg时，仿真结果如图7所示。由图7可知，时间历程呈周期重复；功率谱在分频fp/3及其倍频处存在尖峰；相轨迹在有限的区域内重复，呈封闭曲线；庞加莱图在一定的区域上有3个孤立点存在，说明此时系统处于3倍周期运动状态。
β=1.8N/（mm·kg),F0=20N/kg时，仿真结果如图8所示。由图8可知，时间历程无规律；功率谱出现噪声背景和宽峰；相轨迹在有限的区域内不重复；庞加莱图有无限个孤立点存在，且分布在有限的区域内，说明此时系统处于混沌运动状态[16,17]。
由以上数值试验分析可知，当液压弹簧力非线性项系数β和外加激振力F0取不同值时，系统在运行过程中蕴含丰富的非线性动力学行为。系统可能做单周期运动、倍周期运动，进而通向混沌运动。
5 电液伺服系统动态实验
本节利用非线性动力学研究方法对实测的电液伺服系统的动态数据进行深入分析，以揭示非线性液压弹簧力软硬弹簧特性引起的“跳跃现象”。
5.1 实验系统组成
本文实验按图9所示的系统原理搭建电液伺服系统振动测试实验台。该实验系统可在不同供油压力和负载压力下采集电液伺服系统的状态数据。系统通过调节溢流阀阀口开度来改变系统供油压力；通过调节节流阀阀口开度来改变系统负载压力，以实现系统外加阻尼大小的调整；用精密压力表对系统进、回油路压力进行监测；用振动加速度传感器对执行机构轴向振动信号进行监测；用位移传感器对执行机构实时位置进行监测；用数据采集卡采集传感器输出信号，并传输至计算机系统进行分析处理。
5.2 振动信号的采集及处理
5.2.1 振动信号的采集
在液压缸活塞杆伸出运动状态下，按表1所示的不同工况对液压缸的不同工作状态进行动态测试。其中，输入信号为计算机控制系统给伺服放大器的电压，以控制伺服阀的阀口开度。根据液压缸无杆腔封闭液柱的相对受力情况将液压弹簧的工作特性分为3类：全程软弹簧、半程软弹簧半程硬弹簧、全程硬弹簧。根据供油压力和负载压力大小将外加阻尼大小界定为4类：无、较小、适中、较大。固定输入信号为0.2V，调整主溢流阀及节流阀的阀口开度，使系统分别在表1所示
的12种工况下运行。同时用振动加速度传感器对液压缸整个运行过程中的轴向振动信号进行采集，采样频率为10kHz。
5.2.2 振动信号的处理
采用图10所示的数据处理方案对采集的振动加速度信号进行预处理，并采用非线性动力学研究方法中的时间历程、频闪采样、功率谱图等有效方法对预处理数据进行分析研究[16,17]。
5.3 实验结果分析
供油压力为8MPa时，实验结果如图11～图15所示[18]。图11所示为供油压力为8MPa时，采用图10所示的数据处理方案对采集的振动加速度信号进行处理所得到的振动位移信号的时域波形。比较4种工况可以发现，在整个运行过程中，振动幅值随着活塞位移的变化而变化，其变化规律随工况不同而存在明显差异，这主要与液压弹簧刚度随位移变化有关。由此可以看出，在执行机构的运行过程中，系统动态性能随活塞位移的变化而变化。
图12～图15为供油压力为8MPa时4种工况的分段功率谱图。根据执行机构运行总时间长度，将其分成等分的4段：始段、中前段、中后段、终段。比较4种工况可以发现，振动能量值随负载压力的增大（系统阻尼增大）而逐渐降低，说明随系统阻尼的增大，振动幅值被抑制。不同工况下，功率谱图均由突跳部分和平缓波动段组成，说明均有“跳跃现象”的存在。由于弹簧力软、硬特性交替，波动区覆盖面较大，尖峰数量较多，间隔大小不均，较难分辨，说明“跳跃现象”发生在不同的频率点上。
图16为供油压力为8MPa时4种工况的全程频闪采样图。由图16可知，每种工况都有1个极限环，工况1、2这主要是由于摩擦力作用引起的极限环型振荡现象产生的。图16a、图16b（工况1、2）的轮廓边界由许多离散点构成，图16c、图16d（工况3、4）的轮廓边界比较清晰，这是因为工况1、2所受外加阻尼较小，发生了比较强烈的“跳跃现象”，工况3、4所受外加阻尼较大，“跳跃现象”受到了抑制。
为了验证上述所得结论的普遍性，采用与供油压力为8MPa时相同的数据处理方法，分别对供油压力为6MPa和4MPa时的实验数据进行了进一步的分析研究。通过比较分析，同样可以得到非线性液压弹簧力软、硬弹簧特性会引发“跳跃现象”的结论。
6 结论
(1）电液伺服系统执行机构在运动过程中，液压弹簧刚度随活塞位移的变化而变化，根据工况不同呈现出软弹簧特性或硬弹簧特性。
(2）液压弹簧力的非线性作用可以用含阻尼的Duffing方程来描述，其软硬弹簧特性决定了幅频特性曲线峰值尾部的弯曲方向。阻尼的作用限制了共振振幅的无限上升。激励频率连续变化时，会发生振幅突然变化的“跳跃现象”。
(3）液压弹簧力非线性项系数和外加激振力的大小影响系统的运动状态。当二者参数取不同值时，系统可能做单周期运动、倍周期运动，进而通向混沌运动。
(4）非线性液压弹簧力的软硬弹簧特性引发的“跳跃现象”会使系统响应稳定区域变得复杂，进而造成系统动态特性变得复杂和多变。因此在系统建模和动态特性研究时应该将液压弹簧力的非线性作用考虑在内。
摘要：探究了非线性液压弹簧力对电液伺服系统动态特征的影响。根据非线性动力学原理,建立了系统的动力学模型。通过理论研究指出,非线性液压弹簧力作用可以用Duffing方程描述。通过数值分析揭示了系统内在的分岔现象及典型非线性动力学行为。通过对实测数据进行深入的分析,揭示了液压弹簧的软硬弹簧特性引起的“跳跃现象”。发现液压弹簧力的非线性作用会引发非线性振动,在系统建模与动态特性研究时应该将其非线性作用考虑在内。
非线性行为 篇4
干气密封内部气膜平衡间隙的微小变化极有可能导致动静密封环间的干摩擦或泄漏量增大,因此保证气膜—密封环动态稳定是干气密封可靠运行的关键[1]。国内外学者对此问题也进行了较多的研究。Zirkelback[2]采用有限元法求解微扰雷诺方程,得到的几何参数能够保证密封在具有较大刚度和阻尼系数的同时,具有较低的泄漏率。2002年Miller等[3]对螺旋槽端面密封环在轴向和2个角向自由度上的运动加以分析,用直接数值频率响应法计算了气膜的刚度和阻尼特性;2003年Miller等[4]又利用半解析法求解瞬态雷诺方程,获得了气膜特性规律。Zhang等[5]建立了三自由度(1个轴向、2个角向)的微扰运动方程,并用正交分解法求得了密封环三维运动规律。国内学者大多采用小扰动法研究气体端面密封动力学特性,如刘雨川[6]利用有限元方法求解微扰雷诺方程,并与运动方程联解,迭代解出了密封气膜的动特性系数,并得到了有关密封气膜稳定性的判据;文献[7,8,9]分别用有限元法、近似解析法对轴向微扰或角向涡动下一些气膜动态特性参数进行了计算。目前未见在气膜和密封环流固耦合系统方面的非线性动力学行为研究。本文在文献[8]利用近似解析法求出气膜轴向线性刚度解析式的基础上,通过改变气膜厚度使其随振动位移发生变化,从而将线性刚度转换为非线性刚度,继而应用非线性振动混沌理论来研究干气密封润滑系统的稳定性问题,寻求控制系统稳定运行的结构参数区域,从而指导干气密封的优化设计和研制,保证螺旋槽干气密封系统安全可靠运行。
1 气膜—密封动环系统轴向振动方程的建立
1.1 干气密封工作原理和结构
干气密封结构主要由加载弹簧(波纹管)、O形圈、静环及动环组成(图1)。
1.动环 2.静环 3.弹簧 4.静环座 5,8.O形圈 6.转轴 7.轴套
干气密封的工作原理为:缓冲气体注入到密封装置中,动静环在流体静压力和弹簧力的作用下保持贴合,起到密封的作用。当动环旋转时,会将被密封气体周向吸入槽内(图2),导致气体的压力升高,即产生了流体动压力。当压力达到一定数值时,具有挠性支承的静环将从动环表面被推开,这样使得密封面之间始终保持一层极薄的气膜,厚度为3～5μm。
1.2 轴向振动下气膜—密封动环系统动力学模型
以图1中的动环为振子,气膜为弹簧,轴的激振力为F(t),建立气膜—密封动环系统的轴向振动模型(图3),其动力学数学模型建立的假设条件如下:①将气膜—密封动环系统视为单自由度受迫振动系统;②气膜假定为具有非线性刚度的弹簧;③瞬态激振力来源于转轴对系统的干扰力,且假定为简谐激振力。其振动方程为
式中,m为动环质量;c为系统阻尼系数;K为气膜非线性刚度;F(t)为瞬态激振力;F为稳态激振力;f为系统转子的转动频率;x为轴向振动位移。
2 气膜非线性刚度K(x)计算
应用PH线性化方法及变分运算瞬态微尺度流动场的非线性雷诺方程,得到了轴向微扰下气膜反作用力的增量,继而利用复数转换和迭代法对静态下气膜边值问题进行了求解,获得了气膜轴向刚度的近似解析解[8]。
量纲一气膜刚度a可表示为
β1=n2+β
α2=nχ/ε ω=nφ+β0ζ
式中,A、A1、A2、B、B1、B2、c10、c11、c20、c21均为积分常数;E为槽深的1/2,m;n为螺旋槽数;p0为内外介质压力之比;x为轴向位移,m;χ为中间变量;α为螺旋角,rad;β0为槽斜度系数;δ为两密封环间隙,m;ε为迭代摄动小参数;η为槽深度变化的相对幅度;η1、η2分别为量纲一气膜压力表达式的实部和虚部;ζ为量纲一极径;ζ0为量纲一外径;φ为量纲一极角;ω为当量螺旋角,rad;ω0为当ϕ=0时的当量螺旋角。
为了表示气膜刚度的非线性,已将文献[8]中的
气膜刚度K可表示为
式中,R为动环内圈半径;pi为环境压力。
3 轴向振动非线性动力学行为分析。
选取文献[10]中的实验参数:实验气体为空气,内径Ri=58.42mm,外径Ro=77.78mm,介质压力p01=4.5852MPa,环境压力pi=0.1013MPa,螺旋槽数n=10,螺旋角α=75°,转速nr=10 380r/min,黏度μ=1.8×10-5Pa·s,槽深2E=5μm,密封间隙(气膜厚度)δ=3.05μm,动环质量m=8.04×10-2kg,系统阻尼系数c=0.1N·s/m,轴的激振力F=0.05N。
3.1 动力学模型合理性分析
为了验证本文所建立的力学模型的合理性,在实验参数α=75°、槽深2E=5μm条件下,寻求气膜—密封动环系统振动响应规律,利用Maple程序求解气膜刚度(式(2)、式(3))和振动方程(式(1)),获得了振动响应的相轨(动环轴向速度
3.2 混沌现象
在实验参数α=75°邻域内分析螺旋角响应优化,利用Maple程序求解了气膜刚度(式(2)、式(3))和振动方程(式(1))。当螺旋角α=1.309rad=75°02′18″时出现了如图5所示的动环振动混沌现象,在主共振情况下,相轨图中存在混沌吸引子,Poincaré映射图中(映射周期为转轴旋转周期)出现了典型混沌现象的马蹄形形状。
3.3 混沌分析
非线性气膜刚度K导致了振动混沌现象。当螺旋角α=75°02′18″时,利用Maple程序求得其气膜刚度K的近似函数表达式为
K=2.370 807 656+1.275 846 670×106x+
2.288 228 232×1011x2+1.367 718 331×1016x3 (4)
3.3.1 激振力变化响应
其余条件不变,将激振力F分别由0.05N增加为0.10N和0.15N时,其混沌现象更加明显,振动位移和速度也增大,如图6、图7所示。
3.3.2 阻尼变化响应
其余条件不变,将系统阻尼系数c由0.10N·s/m减小为0.01N·s/m时,其混沌现象愈加明显,振动位移和速度增大,如图8所示。这与物理现象是相吻合的。
3.4 混沌控制
由于混沌现象增大了振动位移,如图6、图7、图8中的振幅均超过6μm,因密封间隙δ=3.05μm,动静环将发生碰摩现象,导致密封失效。为了提高干气密封运行的稳定性,必须控制混沌。本文通过改变干气密封螺旋槽结构参数来调整气膜刚度K,从而避免混沌现象的出现。将螺旋角α=75°02′18″改为α=74°或将槽深2E=5μm改为槽深2E=6μm,由图9、图10中的相轨图及时间历程图可以看出都能得到一准周期的振动图形,且振动幅值较小,最大位移2μm和最大速度1.2mm/s。
4 结束语
干气密封系统本身是一非线性系统,其动力学特性应具有非线性,通过特例验证了混沌现象的存在性,并分析了其非线性动力学行为。为了保证干气密封运行的可靠性,可通过调整干气密封螺旋槽结构参数来控制混沌现象,其机理是通过改变螺旋槽结构参数来增大气膜刚度K,使其值大于气膜失稳时的临界刚度Kcr,本例中通过改变螺旋角和槽深来达到此目的。本计算程序具有普遍适用性,今后可利用此程序对不同工况下的干气密封螺旋槽结构进行响应优化,本研究为干气密封优化设计奠定了理论基础。
摘要：建立了轴向振动下气膜-密封动环系统动力学模型,利用Maple程序求解了轴向振动方程,获得了螺旋槽结构参数响应的振动相轨图、Poincaré映射图和时间历程图,进而分析了螺旋槽干气密封系统非线性动力学行为。研究结果表明:在特定的螺旋槽结构参数范围内存在着振动混沌现象,通过选择合理的螺旋槽结构参数可以控制混沌,为干气密封动态优化设计提供了理论基础。
关键词：干气密封,螺旋槽,非线性动力学,混沌,结构优化
参考文献
[1]曹登峰,宋鹏云,李伟,等.螺旋槽气体端面密封动力学研究进展[J].润滑与密封,2006(5):178-182.
[2]Zirkelback N.Parametric Study of Spiral Groove Gas Face Seals[J].Tribology Transactions,2000,43(2):337-343.
[3]Mi1ler B A,Green I.Numerical Technique for Computing Rotor Dynamic Properties of Mechanical Gas Face Seals[J].Journal of Tribology,2002,124(4):755-761.
[4]Miller B,Green I.Semi-analytical Dynamic Anal-ysis of Spiral Grooved Mechanical Gas Face Seals[J].Journal of Tribology,2003,125(2):403-413.
[5]Zhang Haojiong,Miller B A,Landers R G.Nonlin-ear Modeling of Mechanical Gas Face Seal Systems Using Proper Orthogonal Decomposition[J].Journal of Tribology,2006,128(10):817-827.
[6]刘雨川.端面气膜密封特性研究[D].北京:北京航空航天大学,1999.
[7]李双喜,蔡纪宁,陈罕,等.高速螺旋槽气体密封轴向微扰的有限元分析[J].北京化工大学学报,2003,30(1):52-56.
[8]杜兆年,丁雪兴,俞树荣,等.轴向微扰下干气密封螺旋槽润滑气膜的稳定性分析[J].润滑与密封,2006(10):127-130.
[9]丁雪兴,王悦,张伟政.螺旋槽干气密封润滑气膜角向涡动的稳定性分析[J].北京化工大学学报,2008,35(2):82-86.
非线性行为 篇5
关键词：箍筋约束混凝土柱,率相关效应,本构模型,动力非线性行为
箍筋约束混凝土柱因其良好的抗压性能和延性而被广泛应用于建筑结构中。关于箍筋约束混凝土柱非线性行为研究的核心在于混凝土在箍筋约束作用下的本构关系, 对该本构关系的研究已经有100多年的历史。1903年, Considere发现了箍筋的约束效应[1]。随后, Kent和Park提出了著名的Kent-Park模型[2]。1980年, Sheikh提出了重要的“横向有效约束面积”和“有效约束系数”两个概念[3]。1988年, Mander提出了著名的全曲线本构方程[4]。随后关于箍筋约束混凝土本构的研究大多是基于以上研究成果。
虽然箍筋约束混凝土材料本构领域已经存在大量的研究, 但是仍然存在一些重要的不足:动力本构研究较少[5];经验模型居多;大部分本构模型无法定量表述混凝土的非线性演化过程。这些问题限制了箍筋约束混凝土本构模型在柱动力非线性分析时的应用, 亟待对本构模型进行进一步的动力化和物理化。
本文基于李杰等[6]的混凝土细观随机断裂理论提出物理模型, 利用已有的在不同配箍率、不同加载速率下的混凝土棱柱体试件单轴受压试验得到的试验结论建立模型参数表达式, 该模型同时考虑了箍筋约束效应和率相关效应, 定量的表达了混凝土的损伤演化和塑性演化, 模型结果与试验结果吻合良好。
1 箍筋约束混凝土动力本构模型
1.1 细观随机断裂理论
由于不考虑受拉塑性应变, 箍筋对于受拉混凝土的约束效应较小, 受拉损伤和动力提高系数较小, 因此混凝土受拉本构根据实测混凝土受拉强度采用规范[7]素混凝土本构, 而受压情况则采用细观随机断裂模型。混凝土的单轴性能是研究基础, 根据混凝土细观随机断裂理论, 将混凝土细化为多个混凝土微单元, 见图1。
每个微单元由弹性单元和塑性单元组成, 弹性单元类似于弹簧, 只产生弹性应变εe, 承担全部单元力, 塑性单元只产生塑性应变εp, 不承担单元力。每个微单元的应变发展是相同的, 其中总应变ε为:
受压塑性应变εp的演化法则根据经验确定, 为:
其中, d为受压损伤;ξp, ηp分别为塑性系数和塑性指数, 为模型的两个系数, 一般取为ξp=0.2, ηp=0.5。当弹性应变εe达到断裂应变Δ时, 微单元破坏, 混凝土产生损伤。每个微单元的断裂应变Δ不同, 服从对数正态分布, 该正态分布的均值和标准差定义为λ和ζ, 经推导得损伤d为:
其中, Φ为标准正态分布的累积分布函数。总应力σ为:
其中, E为约束混凝土的动力弹性模量, 在微观上代表弹性单元的弹性模量。
1.2 模型的箍筋约束效应和动力效应引入
以上是素混凝土的细观随机断裂模型体系, 对于动力作用下的箍筋约束混凝土, 本文认为箍筋和动力作用影响的是混凝土的材料属性, 包括混凝土的弹性模量E、断裂应变场均值λ和断裂应变场标准差ζ。箍筋约束效应采用箍筋约束系数cs表达, cs在数值上等于Mander (1988) 提出的“等效横向约束力 (MPa) ”[4];动力效应采用动力系数cr表达, cr采用传统的形式:
其中, 分别为动力、静力应变率, 取值为
经过大量的不同配箍率、不同加载应变率下的混凝土柱单轴受压试验数据分析, 可以得到如下关系式:
其中, x1~z3均为线性系数, 取决于混凝土强度等级。以上公式说明, 模型参数E, λ, ζ分别与箍筋约束系数cs和动力系数cr成线性关系, 根据箍筋配置和加载应变率可以得到以上模型参数, 进而确定整个箍筋约束混凝土动力本构模型。
2 有限元分析验证
为了验证本文所提模型的合理性, 选取相应的试验数据, 并将本文模型导入到Abaqus软件, 对试验数据进行对比分析。
2.1 选取试验数据
选取博士生王德斌 (2013) [8]所做的一组钢筋混凝土柱静力单推和动力推覆试验数据作为研究对象, 试验编号分别为S1和D2, 两种加载方式的对象为同一种柱, 柱的详细信息如图2和表1所示, 静力水平加载速率恒定为0.1 mm/s, 动力水平推覆速率恒定为40 mm/s, S1, D2的加载力—加载位移试验曲线分别如图3, 图4所示。
2.2 有限元建模
本文模型属于一维模型, 因此采用纤维梁单元进行建模, 截面分为三部分, 分别是受箍筋约束的混凝土、纵筋等效层、保护层素混凝土, 如图5所示, 根据试验情况设置约束和荷载, Abaqus模型如图6所示。
2.3 本构关系
约束混凝土采用本文模型。根据图2和表1中箍筋间距、箍筋强度和箍筋布置等信息, 利用Mander的等效横向约束力公式计算得出箍筋约束系数:
动力系数的计算可以根据几何运算得到。图7给出了1 s内柱子的位置变化和长度变化。
右侧粗实线为柱子的初始位置, 假设柱子向左运动, 1 s后柱子移动到左侧虚线处, 柱顶的水平位移为40 mm, 柱顶始终保持在同一水平面上, 则柱高在700 mm的初始值上会减小, 经过几何运算可以得到柱高降低了1 mm, 若从左向右运动则柱高会升高1 mm, 因此动力应变率:
根据公式可得动力系数为:
根据表1中混凝土强度, 引入根据试验确定的相应强度等级的线性系数x1~z3, 具体数值见表2。将箍筋约束系数cs和动力系数cr代入到公式中得到S1和D2的模型参数见表3。
将表3中模型参数数据代入到1.1小节中模型表达式即可得到引入箍筋约束效应和率相关效应的约束混凝土静力、动力本构关系, 通过VUMAT接口将其导入到Abaqus中。纵筋等效层采用两段式本构关系, 如图8所示, 保护层素混凝土的本构关系根据试验数据结合规范[7]得到。
2.4 结果分析
结合以上工作, 利用Abaqus分析即可得到与图3, 图4相对应的静力单推和动力推覆有限元模拟曲线, 将其与试验曲线放到一起进行比对得到图9和图10, 从中可以看出试验曲线与本文模型的Abaqus模拟结果吻合很好, 印证了模型的合理性与准确性。图11和图12分别给出了在最后一次最大水平位移处柱的静力单推和动力推覆损伤分布和塑性分布, 损伤取的是绝对值。静力单推和动力推覆的损伤和塑性都集中于柱底, 并从柱底到柱顶逐渐减小, 符合直观的认识。本文模型不仅仅是一条应力—应变曲线, 还包含了定量化的非线性演化过程。
3 结语
1) 以一种新型的方法提出了箍筋约束混凝土动力本构模型, 即以混凝土细观随机断裂模型体系为框架, 同时定量的引入箍筋约束效应和动力效应。
2) 本文模型不仅能够提供应力—应变全曲线, 还能给出非线性演化过程, 包括损伤演化和塑性演化。
3) 本文模型能够通过VUMAT子程序接口与Abaqus有限元软件结合分析箍筋约束混凝土柱的动力非线性行为, 有限元分析结果与试验结果吻合很好, 证实了模型的准确性和合理性。
参考文献
[1]Considere A.Experimental research on reinforced concrete[M].New York:Mc Graw, 1903:188.
[2]Kent DC, Park R.Flexural members with confined concrete[J].Journal of the Structural Division, 1971, 97 (7) :1969-1990.
[3]Sheikh S A, Uzumeri S M.Strength and ductility of tied concrete columns[J].1980, 106 (5) :1079-1102.
[4]Mander J, Priestley M, Park R.Theoretical stress-strain model for confined concrete[J].Journal of Structural Engineering, 1988, 114 (8) :1804-1826.
[5]魏公涛.箍筋约束混凝土动力本构试验和理论研究[D].上海:同济大学, 2014.
[6]李杰, 吴建营, 陈建兵.混凝土随机损伤力学[M].北京:科学出版社, 2014.
[7]GB 50010—2010, 混凝土结构设计规范[S].
非线性行为 篇6
控制理论和系统动力学方法[1]为探寻复杂系统提供了研究框架,研究者把此方法引入到库存控制系统的研究中,并取得了一定的成果[2]。Disney等[3]建立了供应链线性模型,采用改进的APVIOBPCS提高服务水平;D.Bijulal[4]考察了线性库存控制系统稳定区域内的成本变化情况;Disney[5]利用Jury判据分析供应链库存的准周期、牛鞭效应和稳定性;Cannella等[6]研究不同控制策略对供应链操作性能和服务水平的影响,发现保持系统稳定有时要以牺牲服务水平为代价;Ciancimino等[7]研究发现,供应链同步可以提高系统稳定性,弱化牛鞭效应;Lin等[8]以真实制造企业为原型建立了集成供应链系统动力学模型,发现在订货策略中考虑服务能力调整可以弱化牛鞭效应;Wang等[9]假设系统输入为确定的阶跃需求模式,建立了禁止退货的库存模型,给出了系统的稳定、周期、准周期以及混沌的边界;Garcia等[10]设计切换控制策略以减少缺货率,提高系统稳定性;Hwarng等[11]把对混沌行为的研究引入到不确定需求的供应链环境下。
通过文献梳理发现,提高系统稳定性、改善服务水平和降低库存成本是库存管理研究的主要目标,但综合考虑这三个目标的文献并不多,特别是在更符合实际的非线性系统模型下,大多研究者关注系统稳定性的研究,但系统稳定性与服务水平和成本之间的关系问题涉及较少。本文为探寻三个目标之间的关系,建立更符合实际的分段线性库存模型,分析其李雅普诺夫稳定边界,比较确定需求和不确定需求模式下,服务水平和库存成本随库存策略的变化以及与稳定边界的关联,给出同时满足三个目标的库存策略区域供管理者参考。
1 库存控制系统模型
1.1 基本模型
本文模型适用于处于供应链中的企业,其向上游采购原材料或成品,向下游开展销售活动,其与上下游企业的信息流是畅通的。采用APVIOBPCS订货策略,此策略利用库存调整和在途调整形成对系统的反馈控制,根据调整参数的不同,可以表示多种订货策略,常用的Order-up-to(OUT)订货策略即是此策略的特例。根据APVIOBPCS订货策略,订货活动以周期为单位展开,系统状态在每期末进行更新。每个周期的事件流程为:在第t期初,目标企业收到其在L+1期之前的订货(L为运输时间,本文称为提前期),到货量记为At;然后目标企业从实际库存St中发货(发货量为Ct)以满足下游的需求Dt,由于供货能力有限,目标企业可能不能完全满足其下游的订货需求,可表示为Ct≤Dt;接着,供应链成员决定本期的订货量Ot.每期的订货量由需求预测、库存调整和在途调整三部分组成。此订货策略表达式如式(1):
其中:
S*和WIP*分别为目标库存量和目标在途量,此处设两者都为^Dt.因此订货量公式表达如下:
库存表达式为:
注意,在实际企业运作中,目标企业由于库存限制并不能完全满足市场需求,即企业的服务水平并非都是100% .在此种情况下,St-1+At<0,即企业未能满足客户所有需求,本文即用此条件判断服务水平,具体的服务水平公式详见1.3节。
在途表达式为:
假设供货企业在接到订单后立即发货,则运输延迟为:
产品到货量公式为:
利用指数平滑法进行需求预测,需求预测公式为:
在以上公式中,式(4)加了非负约束条件,使供应链模型更贴近现实,但同时也提高了供应链模型的复杂程度,使整个供应链系统变为非线性系统。
若提前期L =1,设P1= {X|O ≥0},P2= {X|O<0},以上式(4)、式(5)、式(6)、式(7)、式(8)、式(9)可以表示为系统方程:
其中,
1.2 需求模式
为比较不同性质的需求分布对库存系统性能的影响,本文选取确定性需求和随机需求两种需求分布。确定需求为阶跃需求分布,随机需求为ARIMA需求。阶跃需求是考察供应链系统稳定性时常用的需求模式;ARIMA需求模式为平稳可逆的时间序列,更符合现实市场规律,是供应链研究中较常用的一种需求模式[2,12],但由于其随机性,一般研究者在研究供应链动态特性时很少使用。现实的市场需求在大多数情况下都具有绥机性,本文引入ARIMA需求模式,力图模拟更加真实的供应链环境。ARIMA需求公式如下:
其中,μ为常数,ρ为自回归系数,表示本期需求与上一期需求的相关程度(|ρ|<1),εt为市场需求的波动误差,为服从(0,σ2)独立分布的随机变量,β为误差系数(|β|<1)。本文设μ=12,ρ=0.5,β=0.5。
1.3 服务水平和成本
本文选取服务水平(SLR)和平均周期成本(ACOST)作为衡量库存系统性能的指标。服务水平是指在系统运行周期n中,能够100%满足客户需求的比例。根据本文所建立的模型,每期是否能够满足客户需求可以通过检查本期期初库存(即上期期末库存)和本期到货量之和是否大于本期需求量来获得,用公式表达如下:
库存成本假设由两部分组成:库存持有成本h和缺货成本b.每期末计算本周的总库存成本。平均成本公式可表示如下:
其中,St+和St-分别代表实际库存量和缺货量,用公式可以表示为:
2 系统稳定性分析
2.1 系统李雅普诺夫稳定条件
有关系统稳定性的定义有很多,本文采用较常用的李雅普诺夫稳定。由式(10)可以看出,系统在两种状态下切换,根据文献[13],为判断此切换系统的稳定性,分别求得A1,A2和A1×A2的特征值。
当λA1模小于1时,系统渐近稳定。结合以上公式,可得出当提前期为1时库存系统的渐近稳定条件:
对于库存管理,系统稳定是管理的最佳状态,但在达不到这种最佳状态时,若能够达到周期状态,使得管理有规律可循,也不失为一个比较好的方案。根据文献[13],系统有周期解的条件是矩阵A1iA2j渐近稳定,因此可以通过求解A1iA2j特征根,并令其模小于1求得系统有周期解的条件。由于A2j=A2,因此,只需考虑A1iA2.如|λA1A2|<1,可求得系统保持周期2的区域:
若,可求得系统保持周期3的区域:
由于i的取值可以有很多,因此系统的周期区域也有很多,本文只考虑i=1,2,3,即系统的渐近稳定区和周期2、周期3、周期4区域。通过后面实验发现,本文考察的最优库存策略区域都在此范围内,因此本文没有计算更大周期的情况。根据以上分析,可绘制出系统的李雅普诺夫稳定及周期2、周期3、周期4区域,如图1的白色部分。
2.2 最大李雅普诺夫指数
在上小节求解提前期为1的系统稳定边界时发现,随着提前期的增加,求解难度会逐渐增大,因此该方法适合提前期较小的情况,为了使研究可以更方便的考察较大提前期的情况,本节用求解最大李雅普诺夫指数法的方法确定系统的稳定边界。
本文采用wolf重构法计算李雅普诺夫指数,应用构建的供应链模型通过仿真,得到订货量时间序列O(t),以Takens的延迟坐标重构m维相空间中的一条轨道,公式如下:
设轨道O(t)的初始点为O(t0),其邻近轨道O0(t)的初始点为O0(t0),这两个点之间的距离记为L(t0)。设后面时刻为t1,距离L′(t1)满足以下条件:
找到另一个点O1(t1)使得与O(t1)的距离尽量小,即满足以下条件:
迭代上面的过程,直到O(t)为轨道的最后一个点。则最大李雅普诺夫指数值(以下简称LE值)可由式(24)求得。当LE值小于0时,系统处于稳定、总期或准周期态,若LE值大于0,说明系统不稳定。
Wolf法需要估计时间序列的平均轨道周期p,本文根据已有研究的成果[11],设p =15。另外,设嵌入维m =20,时间延迟τ=1。
如图2所示的白色区域为用wolf法求得的本文系统模型对应的LE值小于0的部分,即系统稳定、周期和准周期的区域。与图1比较,两者的渐近稳定及周期稳定区域基本相似,说明论文采用的求李雅普诺夫指数的方法是可行的-模型模拟使用的算法是正确的。比较图1和图2还发现,图2的稳定区域分支要比图1多,这是因为图1中只包括了周期4以内的解,而图3给出的稳定区域范围则包括了系统所有的周期解,这说明,用计算李雅普诺夫指数的方法求解稳定区域更加准确。基于以上分析,本文用计算李雅普诺夫指数的方法得到提前期为2和3时系统的稳定区域。图3(a)和图3(b)给出了当提前期分别为2和3时LE值小于0的区域。从图3可以看出,随着提前期的增加,库存策略的稳定域在减小。在管理实践和研究中都发现这种现象,即某些库存策略只适合提前期较小的情况,而当提前期增大时则不适用。本文的实验结果可以很好的解释这种现象:不同提前期下系统的稳定库存策略域是不同的,要区别加以对待。
3 最优库存策略域分析
3.1 最优库存策略域的含义
本文最终的目的是找到能够同时满足系统稳定、较高服务水平和较低成本这三个条件的库存策略,因此文中所指最优库存策略域为同时满足以下三个条件所对应的库存策略(αSL,αS)区域:(1) 李雅普诺夫指数LE<0;(2) 服务水平SL≥95%;(3) 平均库存成本ACOST≤50。若库存策略(αSL,αS)对应的数值模拟结果同时满足三个条件,即表示此策略属于最优库存策略域,简称最优域。
利用matlab在不同的需求模式下进行仿真实验,运行周期为2000即n=2000,若周期以天计算,则运行周期为5年,对于库存管理,数据量是充足的;根据前节推导的系统稳定域条件,设:αSL∈[-4,4],αS∈ [0,4];以步长0.1进行模拟,即分别计算3321个(αSL,αS)组合下对应的LE,SL和ACOST值,然后找出满足三个条件的最优库存策略域。本文以contour图的形式表达各种(αSL,αS)组合下的LE,SL和ACOST值,如图4所示,给出了确定需求模式下提前期为1时系统的服务水平和平均成本contour图。
3.2 确定需求模式下最优域分析
在库存管理中,确定需求模式并且提前期为1是较为理想的情况。从图4中可以发现,在此理想状况下,大部分的(αSL,αS)区域都能达到较高的服务水平和较低的库存成本;但是仍然有部分区域使得SL<70% 或ACOST>50。这说明即使在理想状态下,若库存策略不合理,仍然会导致性能低下的供应链系统。当αS=0时,若αSL<1,则服务水平SL<70%;反之,若αSL= 0,αS在很大范围内都能保证100%的服务水平。这说明库存调整为0 更容易导致较低的服务水平,考察文献中选取的库存模型也发现,研究者更加注重对库存的调整。通过图4与图2比较发现,稳定域与较高服务水平和较低成本的区域并不完全一致,但三者有重合的区域。
为进一步分析系统稳定性、服务水平和成本之间的关系,图5给出了在确定需求模式下,提前期分别为1、2、3时,库存系统的最优库存策略域(灰色部分)。分析图5发现,图5(a)、图5(b)、图5(c)所示的区域都在同条件下系统稳定域的范围,特别是当提前期为1时,如图5(a)所示的区域与图1所示的系统稳定以及2周期稳定域基本相同。随着提前期的增加,最优域在减小,但三者属于包含关系。表1给出了三种提前期下最优域的数量和占总策略组合的比例。表1显示,虽然提前期L增加幅度都为1,但L从1到2和从2到3最优域的变化幅度并不相同,提前期从2变化到3时,最优库存策略域的减小幅度远远大于提前期从1到2时的变化,说明系统受提前期的影响较大,应尽量缩短提前期。提前期对系统的影响已得到管理者的重视,因此在管理实践中企业会采取各种措施提高效率,尽量压缩提前期,但缩短提前期可能耗费更多管理成本,因此在没有能力缩短提前期时,可考虑把库存策略调整在三种情况的公共最优库存策略域以减少不确定性,提高库存系统性能。
3.3 ARIMA随机需求模式下最优域分析
确定需求是较理想的情况,但大多数情况下,需求是变化的、未知的,因此讨论随机需求下的最优库存策略域更具现实意义。本节选取ARIMA需求模型作为系统的输入,考察当输入变为随机函数时最优库存策略域的变化。图6给出了ARIMA需求模式下,提前期分别为1、2和3时系统的最优域。当提前期为1时,ARIMA需求驱动下的最优域与确定需求时差距并不大,说明若提前期为1,即使需求存在随机因素,对系统的影响也并不大;但随着提前期增加到2和3,最优域的范围大幅度减小,表2给出了ARIMA需求三种提前期下最优域的数量和占总策略组合的比例。
虽然随着提前期的增加,最优域急剧缩小,但各区域仍然存在包含关系,因此在选择库存策略时可以尽量选取这三个区域的交集,即图6(c)显示的区域,值得注意的是常用的order-up-to策略,即(1,1)并不在此交集内。统计(0,0)点附近的最优库存策略(αSL,αS)有七组,分别为:(-0.1,0.1),(-0.1,0.2),(0,0.1),(0,0.2),(0,0.3),(0.1,0.3),(0.1,0.4)。当选择这些库存策略时,可以在ARIMA需求模式下,提前期小于3的情况下,保证系统的性能,即较高的服务水平和较低的成本。
4 结论
提前期的作用及影响一直是供应链管理研究的焦点问题之一,本文的研究不但得出与以往研究相似的结论:提前期对系统稳定性、服务水平和成本都有不同程度的影响,而且给出了不同需求模式下提前期对系统性能的影响程度及不同提前期下都可以保证系统稳定的最优库存策略域。实验结果显示,各提前期的最优域存在包含关系,管理者应尽量选择它们的交集作为库存策略;虽然随机需求是系统不稳定的主要因素,但在ARIMA随机需求模式下,仍然存在最优库存策略域,本文给出了具体的最优库存策略供管理者参考。
非线性行为 篇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分岔现象,分析了分岔过程系统发生的一系列的边界碰撞现象,并进一步分析了稳定状态及准周期状态下系统的时域和频域特性。本文将开关变换器的非线性研究拓展到无线电能传输系统中,为深入研究其非线性行为建立了平台,对无线电能传输系统的可靠设计也有重要指导意义。