非线性动响应

非线性动响应（精选七篇）

非线性动响应 篇1

关键词：高压输电线路,导线跌落,非线性动响应,分析

绝缘子通常在高压输电线路中起到了导线与杆塔和横担之间的绝缘作用, 其要承受导线水平方向的拉力以及垂直方向的荷载力, 所以高压输电线路中的绝缘子必须要具有良好的电气性能同时要保证其的机械强度。而绝缘子在运行的过程中由于受到很多因素的影响发生断裂的不可避免的, 绝缘子一旦断裂就会导致导线跌落, 而在导线跌落的过程中会给当前塔和相邻的输电塔都会产生连续的动力冲击作用。截止目前为止, 对导线跌落给输电线路造成的影响方面研究的不多, 因为众多的原因, 很难从理论上来确定具体的影响程度。下文将以某375米高的大跨越输电塔为实例, 通过有限元的方法对高压输电线路导线跌落的过程进行了模拟, 分析了导线跌落过程中的非线性动响应情况, 并对导线在跌落的过程中和地面的接触状况进行分析, 最后对整个分析过程进行了论证, 以证明整个分析过程是正确的。

1 输电线路的有限元建模

某大跨越输电塔总高度为375米一共分为13层, 塔头的总高度为46米, 塔身是正方形。输电塔的主杆六层以下采用中空夹层钢管混凝土复合构件的形式, 其余的各层采用薄壁钢管, 横隔和斜撑采用薄壁钢管, 采用角钢来完成塔头部分的连杆。该输电塔的跨越距离为2750米, 垂直距离为266.9米, 耐张段的长度为925米。输电线总共分为三层, 上层为两根地线, 下层和中间分别为四相四分裂以及两相导线, 导线的挂线绝缘子的长度为11.5米, 地线的挂线绝缘子长度为1.5米。

进行有限元分析的时候, 通过空间梁单元来对中空夹层钢管以及薄壁钢管混凝土复合构件进行模拟, 通过空间杆单元对角钢进行模拟, 通过质量单元来对节点、爬梯以及其余的连接构件进行模拟, 通过杆单元来实现绝缘子的模拟, 通过单向受拉杆来对导线和地线进行模拟。通过两塔三线模型的建立来实现对高压输电线路的模拟, 其中将每一个杆塔都分为865个梁单元, 52个质量单元以及122个杆单元, 将每个绝缘子分为一个杆单元, 将每根导线分为五百个杆单元, 下图图1为高压输电线路的有限元模型图。

下面将阐述单塔和塔线体系的前六阶的频率以及主要振型, 前六阶单塔体系的频率分别为0.555Hz、0.557 Hz、0.956 Hz、0.973 Hz、0.969 Hz、1.31 Hz, 振型分别为方向1方向一阶弯曲、2方向一阶弯曲、1方向二阶弯曲、2方向二阶弯曲、3方向一阶扭转、3方向二阶扭转;前六阶塔线体系的频率分别为0.539 Hz、0.554 Hz、0.913Hz、0.967 Hz、1.535 Hz、1.732 Hz, 振型分别为2方向一阶弯曲、1方向一阶弯曲、1方向二阶弯曲、2方向二阶弯曲、1方向三阶弯曲、2方向三阶弯曲。

从以上的振型中可以看出, 单塔体系的主要振型为1方向和2方向的平动振型, 而且振动的频率很相似;塔线体系中由于导线在平面内的质量作用, 2方向第一节的弯曲振型的频率要稍微比单塔模型的完全频率小, 1方向二者的频率几乎是一致的, 这时因为平面外受到导线的振动作用的影响很小。

2 高压输电线路导线跌落的过程

文中主要是对绝缘子断裂而引起导线的跌落给输电线路所造成的影响进行分析, 所以先假设有绝缘子已经发生了断裂现象。首先假设导线跌落没有接触地面。在对绝缘子突然出现断裂导致导线跌落的模拟中, 首先对塔线体系在自重作用下的平衡计算出来, 然后哦假设绝缘子在很短的时间之内发生断裂, 绝缘子原本对导线的线性作用力已经消失, 最后对绝缘子断裂之后导线的跌落过程进行计算。

下图图2为导线跌落的过程中在的位移曲线图, 从图中可以看出, 导线在跌落的过程中, 出现了明显的回弹, 而在跨越短的中点位置出现了最大的回弹, 在原挂线点位置导线跌落的速度是最快的, 但之后跌落的速度减小了。

图3和图4分别为A塔8#绝缘子与和导线在跨越段中点的变化曲线, 这两个图中可以看出, 导线在跌落的过程中产生了极为明显的动力效应, 最大的应力值比和稳定状态先的应力相比, 差别比较明显, 有的杆件的最大应力值比稳定状态下的应力值要高出2-3倍。

通过对A塔典型杆件在跌落中的最大应力、自重作用应力、稳定状态的应力以及动力系数的分析中可以得知:第一, 当导线跌落到稳定的状态下, 杆件的内力和跌落之前的内力之间相差的比较小;第二, 在导线跌落的过程中, 所有的杆件的内力都出现了明显的波动, 产生的动力效应也比较明显;第三, 挂线杆中, 导线跌落所导致的最大应力表现的较为明显, 因此对挂线杆要引起重视。

下图图5和图6分别为斜杆和主力干的轴应力时程图, 从这两个图中可以看出, 每对杆件的稳定状态应力都偶反对称的现象, 相当于在外力作用的影响下, 横道线方向出现了左右摇摆, 这种现象实际上就是摇摆效应。

通常情况下, 导线跌落都会和地面发生接触, 通常接触面上的接触条件可以分为接触力和接触位移条件, 而接触位移条件有分为切向和法向两种。接触的状态可以分为滑动状态、分离状态以及黏结状态三种。所以该种问题属于边界的非线性问题。由于文中的输电塔处于山区, 假如导线跌落到地面被杂草绊住动能就会完全得打吸收, 导线和地面接触之后就不会分开。通过刚性解析表面来对地面进行模拟, 下图图7为导线和地面接触的模型图, 图8为和地面接触的导线跌落的过程图, 图9为跌落导线在跨越段的中点以及原挂点位置的竖向位移图。从这三个图中可以看出, 导线跌落的过程中, 依然会产生明显的回弹现象, 和导线的自由跌落相比, 回弹的幅度相对来说比较小;断裂点的导线在和地面接触之后没有出现回弹现象, 但依然出现了沿着2方向的运动;由于阻尼的影响, 导线在跌落之后的100秒已经逐渐稳定。

下图图10为导线跌落过程中A塔8#绝缘子的轴应力的时间和路程, 从此图中可以看出, 在整个跌落的过程出现了较为明显的动力效应, 但是和自由跌落过程相比, 振荡的幅度明显的减小了。

通过对A塔典型杆件在跌落过程中应力值进行计算分析可以得知, 导线在跌落的过程中所有的构件都出现了明显的内力波动, 出现了明显的动力效应, 到那时和和导线自由跌落过程相比, 计算的结果没有过大的差别, 所有当绝缘子断裂之后导线如果和地面的接触段不长的话可以按照自由跌落来处理。

4 分析导线跌落的能量

在对导线跌落能量的分析中, 根据能量守恒定律可以得知当导线自由跌落时体系的总能量基本上是一个常数, 如果导线和地面接触, 就会产生一定能量损失, 系体系的能量将不会是一个常数。而导线和地面接触的面积是有限的, 因此系统产生的能量损失是很小的, 由此可间上文中的分析是合理的。

结语

综上所述, 文中通过单塔和塔线体系有限元模型的建立, 经过分析可以得知单塔体系中主要的振型为1和2方向的平动, 而且这两个方向上的振动特性和相似, 他线体系由于受到质量的影响, 频率出现了明显的下降。导线在自由跌落的过程中, 出现明显的回弹现象, 输电线路中所有构件的内力会产生明显的波动, 产生的动力效应非常明显。中轴面对称的杆件有摇摆效应产生, 在导线跌落的过程中有些杆件会出现比较大的峰值应力;和地面接触的导线在跌落的过程中, 导线依然有回弹现象出现, 和自由跌落相比回弹的幅度较小, 通过对导线两种跌落情况的分析比较, 二者产生的内力差别不大, 因此和地面接触跌落可以按照导线自由跌落的情况来对待。文中最后根据能力守恒定律, 对系统的总能量进行了简单的分析, 对文中两种模拟分析的科学性进行了验证。

参考文献

[1]夏开全, 刘云, 钱振东.高压输电线路导线跌落非线性动响应分析[J].东南大学学报 (英文版) .2009 (1) .

[2]梁政平, 李黎, 王乘, 夏正春.大跨越输电塔在线路断线作用下的动力响应[J].华中科技大学学报 (城市科学版) .2009 (1) .

非线性动响应 篇2

摘要：提出了一种基于动力响应主分量的瞬时频率和幅值的非线性模型修正方法。首先通过解析模式分解和希尔伯特变换提取结构动力响应主分量的瞬时频率和瞬时幅值；然后选取瞬时幅值和瞬时频率慢变部分的有限多个时间点值来反映结构的非线性特性；最后采用响应面模型，并基于试验数据和有限元模型数据之间瞬时幅值和频率的残差建立目标函数，进行结构的非线性模型修正。通过三层框架的数值模拟分析，其结果表明该方法能精确有效地修正非线性结构模型。

关键词：解析模式分解；主分量；瞬时频率；瞬时幅值；非线性模型修正

引言

土木工程结构在其服役期限内会不可避免地受到较强的荷载激励，从而表现出非线性特性。例如，结构在地震、台风等极端荷载作用下，在一定程度上会表现出较强的非线性行为，如材料本身的非线性导致结构的非线性，结构发生大变形时产生的几何非线性，结构阻尼耗散的非线性，结构边界条件及状态的非线性。非线性结构系统往往又表现出复杂的动力学行为，因此，研究强荷载作用下的结构非线性模型修正，不仅可以对结构的安全状态、劣化行为进行准确评估，同时也可以利用修正的非线性模型，预测或者重现下一次强荷载作用下的结构响应，并利用相应的分析结果，提前做好安全加固工作，保障结构的安全性。

近三十年来，以提高有限元模型计算精度为目标的有限元模型修正技术得到了深入的发展。但是，以往研究的有限元模型修正大多数是线弹性结构关系，其应力应变关系和所建立的状态平衡方程都是线性的，修正中利用的也是固有频率等线性系统特征量。目前，国内外关于非线性有限元模型修正的研究文献较少。Hemez和Doebling等学者首先提出了结构非线性模型和实验测试数据相关性的理论，依靠实验数据处理非线性结构的实验一分析相关性和逆问题，并提出了非线性模型修正的概念，讨论了非线性模型修正的必要性和面临的一些挑战。Song等提出了一种基于低幅值环境振动数据的非线性结构模型修正方法。Schmidt等通过匹配数值时程响应和测试结果研究了局部非线性有限元模型的修正问题，如库伦摩擦力、问隙和局部塑性。Meyer等基于等效线性化运动方程使用谐波平衡法识别局部非线性刚度和阻尼参数来获得频域内模型的合理描述。silva等比较了一些非线性模型修正方法的优缺点，如谐波平衡法、恢复力曲面法和正交分解法等。费庆国等研究了基于神经网络的非线性结构模型修正方法，并讨论了应用神经网络进行模型修正的关键问题。最近，随着时频分析方法的兴起，基于时频分析的非线性模型修正方法开始引起一定的关注。Asgarieh等提出了基于时变模态参数的非线性模型修正方法，其时变参数通过确定的随机子空间方法提取。该方法利用识别的瞬时频率和振型作为优化目标函数来修正非线性模型的滞回材料参数，并基于Bouc-Wen模型假定利用振动台实验数据验证了该方法的可行性。

总体来讲，结构在遭遇强烈地震、强风等极端荷载作用时会呈现出明显的非线性特征，结构的非线性行为可以通过结构恢复力随结构变形变化的曲线滞回环来表现。滞回环包含的面积即为结构在一个周期内由摩擦所消耗的能量，滞回系统是包含了非线性刚度及非线性阻尼的典型非线性系统。近年来，Caughy双线性模型、Neilsen退化双线性模型、Clough退化双线性模型、Bouc-Wen模型等各种滞回曲线被用来表征结构的非线性特征。因此，非线性模型修正的整体过程是采用合适的非线性模型模拟实际结构，通过优化目标函数修正非线性结构模型的参数，使计算结果与实验数据更好地吻合，最终建立符合实际情况的非线性结构模型。

水下系泊支撑平台非线性动力响应 篇3

水下系泊支撑平台的工作环境为一定深度的水域中, 是一种新型的装置, 在工作过程中会受到来自于外界流场的不确定的流体载荷, 工作环境复杂。为了使这类新型的支撑装置安全、可靠, 需要详细研究它受到的外界荷载。这类新型的装置国内外对它的研究很少, 考虑到它的结构和水下悬浮隧道及各类海洋平台在系泊方式上的类似之处, 在课题研究过程中借鉴了它们的研究成果。

Gault[1]采用悬链线方法对系泊缆的恢复刚度进行了分析, 发现系泊缆的恢复系数具有很强的非线性, 必须考虑几何非线性的影响。但是悬链线方法忽略了系泊缆索所受的流体作用力等部分外力, 是对实际问题简化后的近似解。若要对系泊缆索的静力特性进行精确模拟, 必须考虑所有外部载荷的影响。唐友刚, 等[2]采用数值模拟方法研究了考虑流场作用时深海平台系缆形状和张力, 探讨了张力分布和上部浮体运动的影响。最近几年, 国内学者唐友刚等[3]和王丽元[4]发展了对系泊缆索进行动力分析的二阶摄动理论。此方法可以考虑系泊系统非线性和动态耦合的影响, 能够模拟系泊浮体的慢漂响应。本文在研究过程中考虑了系泊支撑平台具有长细比大、约束方式、平台结构特征的特殊性, 利用水力学理论和海洋工程结构力学理论, 建立了水下系泊支撑平台的数值分析模型, 采用伽辽金法和综合数值解法, 计算了支撑平台浮体某阶固有频率与漩涡泄放频率相近发生谐振时, 系泊浮体的非线性动力响应。

1 流场载荷

水下系泊支撑平台由系泊缆、浮体以及其他安装在上面的附件组成。在工作水域中它的动力学特性除受到环境载荷的影响外, 还受到自身结构形状的影响, 环境载荷中影响最大的是流场速度、方向、浮体的布放深度。处于一定深度水域的流场环境中浮体受到的主要作用力是流体力和绕流阻力, 流速垂直于浮体轴线方向的流体力和绕流阻力最大, 流速平行于浮体轴线方向的绕流阻力最小, 因此在进行受力分析时主要考虑浮体轴线垂直于流场速度的方向, 如图2所示。以浮体轴线的中点为坐标原点, 浮体轴线为X轴, Y轴垂直于浮体轴线指向水面, 流场速度的正向为Z轴方向。

浮体在水域中受到的垂直于流速方向的作用力主要有水流作用在速度方向的惯性力fi (x, t) 和流体拖曳力fd (x, t) , 以及垂直于流速指向流场表面的水流涡激升力fl (x, t) , 它是由涡激振动产生的。力的平衡方程为

根据流体动力学理论, 流场作用于浮体表面的惯性力fi (x, t) 和拖曳力fd (x, t) 可以根据莫里森公式计算[5]。

式中, CM为惯性力计算系数, D为浮体截面平均直径, CD为拖曳力计算系数;设江水密度为ρ。

涡激振动产生的流体升力fl (x, t) 可根据流体动力学理论“卡门涡旋”来进行求解。[6]

设ωs为流场在浮体表面产生漩涡脱落频率, , St为0.2[7], CL为升力计算系数, 设流场流速为v。浮体单位长度所受到的流场作用力为

2 分析模型

考虑到水下系泊支撑平台浮体在流体力的作用下, 在流速方向和垂直于流速方向会产生较大位移, 但是在流速方向位移很快趋于平稳, 垂直于流速方向的位移由于水流涡激升力的影响会产生振动[8]。在流场中流体作用力的影响下, 位于两个系泊缆的支撑之间的浮体的中部会产生较大的位移, 整个浮体会产生弯曲振动现象, 这种弯曲振动会对浮体的稳定性造成一定的影响, 模型的建立基础是采用在受到外界载荷作用变形之后的浮体微元模型, 考虑了变形对稳定性的影响, 此时, 考虑了几何关系的非线性因素建立的结构平衡方程也具有非线性的特征, 适用于进行结构非线性动力响应分析。

浮体在涡激振动过程中, 假设截面始终满足经典梁理论[9], 以y (x, t) 表示垂直于流速方向的振动位移, 取长度为dx微元进行分析, 如图3所示。

根据浮体微元受力分析图, 设浮体受到的轴向力为N, 浮体绕轴线的转角为θ, 浮体在垂直于流速方向受到的流体作用力为f, M为浮体在流场作用下产生的弯矩, fi为结构在振动过程中产生的惯性力和结构阻尼力。

涡激振动时, 考虑结构自身的振幅很小, 即θ变化很小, 有sinθ≈θ, cosθ≈1, 外力作用下浮体受到的力和力矩的变化为

建立浮体微元在涡激振动方向的平衡方程:

浮体受到的轴向力N为

式 (8) 中, A为浮体平均截面积;l为浮体整体长度;E为浮体材料弹性模量, 可得

为了求解方程, 将浮体在涡激振动方向的振动位移表达成振型的级数形式[10], 采用分离变量的方法, 即

把式 (10) 代入式 (9) , 根据伽辽金法, 可得

式 (1) 中, R (x, t) 称为留数。

式 (11) 是一个常微分方程组, 根据微积分理论, 它的个数与公式 (10) 中的系数n有关, 阻尼力在非线性方程组中的处理很困难, 因此, 式 (10) 中的系数n的取值不能过大, 否则会影响求解, 随着振动频率的提高, 前几阶振型对振动响应的影响较大, 随着频率的逐渐提高, 振动响应依次减小, 根据数值积分的理论, 引入以下符号:

式 (13) 中,

式 (13) 中, sgn与时间t及空间x有关的变化量, 可以采取数值分析的方法来进行求解, 通过查阅参考文献可知[11], 当项数n=3其动力响应的计算值就能满足实际需要, 对式 (11) 进行积分, 把方程组改写为

式 (15) 中, Dn可由式 (13) 进行积分求出, 它是yn (t) (n=1, 2, 3) 的二次函数。因此式 (15) 可以看成三个耦合的具有强迫振动的非线性方程组, 可以采用Runge-Kutta法求解。

3 设计数据

按照结构设计的具体要求, 根据理论分析及相关参考文献的推荐, 拟定数据见表1。

4 动力响应分析

水流与系泊支撑平台浮体相互作用的一个重要现象是同步效应, 即流体经过系泊支撑平台浮体时, 会出现分离产生涡旋, 当浮体的固有频率范围包含在漩涡释放的频率范围时, 涡旋的释放将由浮体的振动控制, 此时浮体振动频率与涡旋释放频率相一致, 和结构振动的自然频率相接近, 产生谐振[12]。求解式 (15) , 以系泊缆支撑为边界条件, 可得模态响应图, 图4~图6分别为前三阶谐振时三个模态的动力响应值Y1、Y2、Y3。

从图4~图6可以看出, 二阶谐振模态响应Y2和三阶谐振模态响应Y3比一阶谐振模态响应Y1小很多, 表明第一阶振型对振动位移的贡献较大, 高阶次振型对振动位移贡献较小, 由图可以看出, 随着振动的持续, 各阶模态的非线性特征表现十分明显;对于系泊支撑平台浮体, 由水流涡激作用引起的结构的共振要经历一段时间后才能形成幅值一定的周期振动, 随着时间增大其非线性特征逐渐增强。

求解式 (15) , 根据材料力学的理论可以求出弯矩M (x, t) 和剪切力Q (x, t) 的表达式为

求解上式计算得出系泊支撑平台浮体前三阶模态对应的动弯矩M1、M2、M3和动剪切力Q1、Q2、Q3的分布如图7、图8所示。

模态对位移、动弯矩、动切力影响的幅值百分比见表2。

从图7、图8可以看出, 在系泊支撑平台浮体的中点处出现最大动弯矩, 而在系泊支撑平台浮体的约束处出现最大动切力。由此可以看出, 进行系统的整体动力分析时, 虽然高次模态振型对振动位移产生的影响较小, 但是对动弯矩、动切力的影响要远远大于对振动位移的影响, 尤其是在很大程度上影响了动切力的数值, 因此在计算内力响应时, 应该多取几级振型来保证计算精度。

5 结论

通过分析水下系泊支撑平台在涡激振动下所受到的流体载荷, 建立了其非线性流固耦合动力学方程, 采用伽辽金法和综合数值解法, 计算了当支撑平台浮体某阶固有频率与漩涡释放频率相近发生共振时, 系泊浮体的非线性动力响应, 通过分析可以得到如下结论:

(1) 第一阶振型对振动位移的贡献最大, 其余高阶振型对位移的贡献依次减小, 各阶模态随时间变化呈现明显的非线性特征;

(2) 对于系泊支撑平台浮体, 由水流涡激作用诱发结构的共振要经历一段时间后才能形成幅值一定的周期振动, 其非线性特征随着时间增大逐渐增强;

(3) 在系泊平台浮体的中点处出现最大动弯矩, 而在系泊支撑平台浮体的约束处出现最大动切力;

(4) 虽然高次模态振型对振动位移产生的影响较小, 但是对动弯矩、动切力的影响要远远大于对振动位移的影响, 尤其是在很大程度上影响了动切力的数值, 因此在计算内力响应时, 应该多取几级振型来保证计算精度。

摘要：针对水下系泊支撑平台在水流涡激作用下的非线性动力响应, 基于水力学理论和海洋工程结构力学理论, 建立流固耦合非线性动力学方程。采用伽辽金法和综合数值解法, 计算了支撑平台浮体固有频率与漩涡泄放频率相近发生谐振时, 系泊浮体的非线性动力响应。结果表明, 第一阶振型对振动位移的贡献最大;高阶振型对振动位移的贡献依次减小, 且各阶模态随时间变化呈现明显的非线性特征。高阶模态对动弯矩、动切力的贡献要远远大于对振动位移的贡献, 尤其是对动切力的影响非常大。在计算动力响应时, 应该多取几阶振型来保证计算精度。

关键词：系泊支撑平台,非线性,流固耦合,谐振,动力响应

参考文献

[1] Gault, Lee M Y, Devlin P, et al.Analytical studies of resonance in taut-moored system.OTC, 2012;2:401—416

[2] 唐友刚, 樊娟娟.高雷诺数下圆柱顺流向和横向涡激振动分析.振动与冲击, 2013;32 (13) :88—92Tang Yougang, Fan Juanjuan.In line and transverse vortex-induced vibration analysis for a circular cylinder under high Reynolds number.Journal of Vibration and Shock, 2013;32 (13) :88—92

[3] 曹菡, 唐友刚.半潜型风机浮式基础设计及幅频运动特性研究.海洋工程, 2013;31 (2) :61—67Cao Han, Tang Yougang.Design and frequency domain analysis of semi-submersible floating foundation for offshore wind turbine.The Ocean Engineering, 2013, 31 (2) :61—67

[4] 王丽元, 唐友刚.球形波浪浮标运动性能分析.海洋技术, 2012;31(4) :10—14Wang Liyuan, Tang Yougang.Analysis of the motion performance of a spherical wave buoy.Ocean Technology, 2012;31 (4) :10—14

[5] Tariverdilo S, Mirzapour J, Shahmardani M, et al.Vibration of submerged floating tunnels due to moving loads.Applied Mathematical Modeling, 2011;35:5413—5425

[6] 孙胜男, 陈健云, 苏志彬.悬浮隧道锚索黏弹性阻尼器的最优阻尼系数.中南大学学报, 2011;42 (6) :1791—1796Sun Shengnan, Chen Jianyun, Su Zhibin.Optimal damping coefficient of visco-elastic damper for submerged floating tunnel tether.Journal of Central South University, 2011;42 (6) :1791—1796

[7] 何孔德, 李友荣, 方子帆, 等.水下系泊支撑平台非线性振动稳定性分析.海洋工程, 2012;30 (3) :119—124He Kongde, Li Yourong, Fang Zifan, Yang Weihua.Stability analysis of nonlinear vibration for submerged mooring monitoring platforms.The Ocean Engineering, 2012:30 (3) :119—124

[8] 麦继婷.波流作用下悬浮隧道的响应研究.成都:西南交通大学, 2005Mai Juting.The study on responses of a submerged floating tunnel subjected to the wave and current.Chengdu:Southwest Jiao Tong University, 2005

[9] 唐友刚.海洋工程结构力学.天津:天津大学出版社, 2008Tang Yougang.The structural dynamics of ocean engineering.Tianjin:Tianjin University Press, 2008

[10] 麦继婷, 罗忠贤, 关宝树.流作用下悬浮隧道张力腿的涡激动力响应.西南交通大学学报, 2004;39 (5) :600—604Mai Juting, Luo Zhongxian, Guan Baoshu.The dynamical response for submerged tunnel by current and wave.Southwest Jiaotong University Press, 2004;39 (5) :600—604

[11] Gobat J I, Grosenbaugh M A.Time domain numerical simulation of ocean cable structures.Ocean Engineering, 2006; (3) :373—400

非线性动响应 篇4

电荷耦合器件(CCD)是一种广泛应用的光电检测器件,它能将入射到其光敏面上的光信号转换成电荷信号,并且在一定的时序驱动下,将这些电荷信号按照像元的顺序依次输出。这样,每个像元的输出信号就代表了入射至光敏元上的光信号。在CCD的工作范围内,它有着线性的光电响应特性[1],但是当光照很强时,CCD像元中存储的电荷就会饱和,过多的电荷会溢出到相邻的像元和CCD基底中,CCD的光电响应就不再是线性,从而影响光电测量。

为了准确地应用CCD进行光电检测,就必须要弄清楚它的光电响应特性。到目前为止,测量CCD光电响应特性的方法主要有尖劈法、双缝衍射法、小孔衍射法。尖劈法是将入射至尖劈上的光,经过两个表面的多次反射形成一列光强分布相似能量逐级递减的子光束[2],并且递减的倍数就是两个表面的反射率相乘。这种方法可利用的点太少,比如递减的倍数为0.9,反射20次之后也就只剩下原来的0.12倍了;双缝衍射法是利用双缝夫琅禾费衍射的原理,将衍射后的光强用CCD接收,由衍射的光强公式结合像元位置可以测得光电响应[3]。这种方法虽然可利用的像元点多,但是出现了CCD的像元与衍射条纹对不准的问题, 从而影响测量。小孔衍射法是对双峰衍射法的改进,很好的解决了像元对不准的问题[4],但是由于衍射之后的图像是个圆斑,圆斑的中心位置不容易找到,只能通过算法来寻找中心点。目前利用CCD进行光电测量的研究很多[5,6],用CCD进行光电测量是建立在CCD的输出是与积分时间成线性关系的基础之上的。 但是在曝光量很大的时候,CCD的输出不再与积分时间成一个线性关系[7]。 本文提出一种利用LED的发光来研究CCD的光电响应特性曲线,并且做了CCD的输出与积分时间的相关实验。对比实验结果,从理论上分析了实验现象,为更好地利用CCD进行光电检测提供了参考。

1LED发光法测量CCD的光强响应曲线

当CCD的电源电压一定的时候,有光照射在CCD的像元上,CCD的有效输出电压就是曝光量的函数, 可用响应函数f描述为

其中:V表示CCD输出的有效电压,E表示CCD像元接收的照度,t表示积分时间,E t表示曝光量[8]。 可见,如果保持曝光量中一个因数不变,CCD的有效输出电压就是另一个因数的函数。我们所做的这两个实验就是基于这种原理,分别来探究CCD的输出与光强和积分时间的函数关系。实验中用自己设计的CCD数据采集系统,CCD是采用东芝公司的TCD1304DG线阵CCD,它有3 648个有效像元,暗电流低,灵敏度高,光谱响应范围为350800 nm,已被广泛应用于光谱仪器中。并且对CCD的输出信号进行处理,最后用12位的模数转换器(ADC)采集,并用通用串行总线(USB 2.0)将数据发送到上位机中。

一般来说,评价LED性能的好坏,光通量是一个很重要的因数,在其额定工作范围内,LED的光通量与前向电流成正比,随着电流增加,LED光通量随之增大[9]。 实验中所采用的LED是Ga N材质的,它的典型的光通量电流(L-I)输出特性如图1中曲线所示,图1中虚线表示线性关系。加在Ga N LED上的前向电流越大,光通量就越大,在一定的范围内,光通量与前向电流成正比。但是,当前向电流增加到足够大之后,光通量逐渐出现饱和,这种现象称为droop效应[10]。图中可以看到在前向电流小于10 m A时,L-I的线性很好。

1.1实验

在实验中,采用的是 Ф5的单色Ga N LED,正常工作电流20 m A,而实验使用的电流范围为0到5 m A, 在这个范围内,LED的L-I有很好的线性。实验就是基于这个原理来对CCD的光强响应进行测量。

实验装置示意图如图2所示,实验采用一个500 的电阻串联一个LED灯,通过改变电阻和LED灯的电压值就可以改变流过LED的电流,同时用我们设计的CCD数据采集系统采集光信号,最后得出电流与光强的关系。为了保证测量结果的正确性,先测量在没有光照时CCD的输出,结果在423左右,这个背景在CCD输出时应该减掉。实验都是在暗室中进行的,以防止外部光对测量结果的影响。由于CCD属于高灵敏的光电检测器件,如果将LED发出的光直接照射到CCD的光敏面上,CCD的输出很可能就饱和, 所以在LED与CCD之间加了一个光阑,这样就可以使光线不至于直接照射在CCD上。由于所用的LED正向导通电压在2 V左右,实验就将供电电压从2 V开始,每次增加0.1 V,直到电压达到5 V。在每个电压值下记录CCD的读数,同时测量流过LED的电流。

1.2实验结果

由于LED的光斑是圆形的,光强是从圆心向外逐渐递减,所以实验取固定的10个像元数据,将这10个像元的平均值作为真值。得到的实验结果如图3所示,从图中可以看出,在LED的电流小于4 m A时, CCD的输出随电流的增加变化较大,而当电流大于4 m A时CCD的输出随电流的增加变化相对较小。如果将CCD的输出与电流进行线性拟合,得到的结果为

相关性系数R2为0.992。如果去掉4 m A以上的数据再进行拟合,R2为0.997,可见4 m A以上的数据影响了CCD响应的线性。对比与其他测量CCD光电响应特性的方法,我们的实验装置简单,不需要来回调节光路,操作方便。由于使用的像元固定,也避免了由于像元差异引起的误差。同时由于LED发光成像是一个圆斑,圆斑上在小范围内发光强度差别很小,就可以取相邻的几个像素点的均值作为拟合真值,减小了实验误差。本实验所采用的是线阵CCD,但是分析方法同样适用于面阵CCD,并且由于面阵CCD的单个像元面积较小,相对于线阵CCD来说,每个像元接收到的光信号差别更小,所以实验效果可能更好。

2CCD积分时间响应曲线

2.1实验

CCD的光电响应特性还与另一个因素有关,也就是光积分时间,这个实验就来探讨CCD的响应与积分时间的关系。实验示意图如图4所示,主要部件有标准灯,光栅分光系统,CCD数据采集系统以及上位机。实验中采用晶飞公司的D1001标准灯,由光纤将光信号导入自制的光栅分光系统中,用CCD接收色散之后的光信号,这样CCD的像元上就会依次响应不同波长的光信号。由于D1001标准灯是汞灯,有多条特征谱线,综合考虑CCD的波长响应率和汞灯的光谱强度,为了使测量的范围大而且容易观察,实验选择它在404.66 nm处的特征谱线作为实验的数据采集点。实验观察到,在积分时间为8 ms时,在像素点2 377(对应波长404.66 nm)刚有输出,而在积分时间为22 ms时,输出接近饱和,所以实验的积分时间选择在8 ms到22 ms之间。实验中分别在每个积分时间采集5组数据,用均值作为CCD输出有效值。

2.2实验结果与分析

图5所示为实验实际测得CCD的积分时间响应曲线,横坐标表示积分时间,纵坐标是CCD的输出信号经AD转换后的值。很明显可以看出,在积分时间8 ms到17 ms之间CCD的输出增加得比较快,而积分时间在18 ms到22 ms之间由于CCD的输出接近饱和,输出增加的较慢,这个现象与第一个实验结果相似。从实验结果可以看出,CCD的光电响应可以分为两个部分来考虑[11],即积分时间18 ms以下和18 ms

以上两个部分,第一部分CCD的输出随曝光量的增加快而第二部分相对较慢。

CCD的成像原理是内光电效应,随着曝光量增加内光电效应越明显,为了防止像元饱和后过多的电荷溢出到相邻像元中,目前CCD广泛采用垂直溢漏技术,其光电响应可用一个NPN型耗尽基底的晶体管模型来分析[12],N层相当于发射极,P层相当于基极和N型基底相当于集电极,在基底电压不变的时候,击穿电流IPT是晶体管N层电势的函数:

式中:β=q/k T,I0取决于基底电压工艺参数,η 为非理想性因数,VN表示晶体管N层电势。VN的变化:

其中:I是总电流,CPN是晶体管等效电容。当VN很小时,击穿电流IPT几乎为零。由式(4)可推导出:

其中T为积分时间。而当VN很大的时候,IPT约等于I,则可由式(3)推导出:

可以看出,CCD的输出前一部分与曝光量成线性关系,后一部分与曝光量成一个对数关系,将其分别称之为线性区和非线性区。将积分时间为18 ms以下的数据和18 ms以上的数据分别进行拟合得到图6。

图6(a)横坐标是积分时间,图6(b)横坐标是积分时间的对数,纵坐标都是CCD输出电压,拟合曲线为

相关性系数R2分别为0.993和0.995,可以看出拟合的效果较好。在实际测量中,由于采用的是12位的AD,电压范围为03.3 V,所以采样误差最大为±LSB,也就是在0.8 m V左右;同时由于采用的是标准汞灯,谱线准确度和强度有了保证。

根据以上分析可见,当曝光量过大时,CCD输出与曝光量不再是线性关系,所以建议在使用CCD时尽量不要使曝光量过大。当然也可以使用补偿的方法来实现精确的测量,用线性区的延长线作为校准线,来补偿所测到的结果。同时为了验证实验效果,本文把校准后的结果与Maya2000Pro光谱仪作比较。 对比图形如图7所示,横坐标为积分时间,纵坐标为归一化的CCD输出,可以看见经过校准之后,CCD的输出在18 ms以上有了很大的改善。经过拟合,在积分时间8 ms到19 ms之间实验测量结果线性度为0.998 1(校准之后积分时间22 ms已经饱和),已经非常接近Maya2000Pro的0.999 9,效果非常好。

3结论

非线性动响应 篇5

近些年来, 国内外学者对超大跨径混凝土拱桥的研究相当活跃, 不断提出修建超大跨径混凝土拱桥的设想, 并开展了相关研究[1,2,3]。混凝土拱桥向超大跨径方向发展有两个关键问题需要研究: 一是其结构自重大, 施工困难, 因此如何有效地减轻其自重是研究方向之一; 二是施工架设方法, 拱桥的施工方法是影响拱桥方案能否成立、能否被采用的最关键技术问题, 对于超大跨径拱桥更是如此[4]。

为此, 文献[5]创新性地提出了波形钢腹板- 混凝土组合拱桥新桥型的构想, 即采用轻质高强的波形钢板取代传统钢筋混凝土拱圈中厚重的混凝土腹板, 通过连接件与顶底混凝土板联结形成新型的拱圈结构, 从而形成波形钢腹板- 混凝土组合拱。相关的试设计研究结果表明, 这种新型的组合拱桥, 拱圈自重轻, 免除了在混凝土腹板内架设模板、绑扎钢筋等施工工序, 方便施工并缩短了工期, 具有相当的应用前景[6,7,8]。

在抗震性能研究方面, 文献[9]分析了波形钢腹板- 混凝土组合拱桥的线性地震响应特性, 但是分析时是将钢- 混凝土组合截面简化为单一的材料并假定其为弹性材料。显然, 对于这种适用于大跨径甚至是超大跨径的新型钢- 混凝土组合拱桥来说, 这样的抗震计算方法过于简单和粗糙了。大量的研究也表明, 在大跨度拱桥的抗震计算过程中材料非线性与几何非线性的影响是不能忽略的[10]。因此, 为了进一步论证这种新型组合拱结构在实际工程中的应用可行性, 笔者在文献[11]试设计的基础上, 建立波形钢腹板- 混凝土组合拱桥的空间有限元模型, 对其基本动力性能和非线性地震响应特性进行分析, 并与福建宁德岭兜大桥 ( 混凝土箱拱) 进行对比, 探讨这种新型组合拱桥的抗震性能, 为该桥型的实际应用提供参考。

1 有限元计算模型与计算方法

岭兜大桥主跨为净跨径160m的悬链线钢筋混凝土箱拱, 净矢跨比为1 /4, 拱轴系数为2. 114。该桥的总体布置如图1 所示。试设计的波形钢腹板- 混凝土组合拱桥采用与原桥相同的悬链线作为拱轴线, 主拱圈内横隔板的位置与原桥保持一致[11]。

为了更准确地考虑空间杆系结构的几何非线性和材料非线性, 本文采用基于梁单元的双重非线性有限元程序—NL_Beam3D[12], 建立了岭兜大桥和波形钢腹板- 混凝土组合拱桥的空间杆系有限元计算模型。

采用能自动同时考虑轴力和两方向弯矩非线性的纤维单元来考虑钢筋混凝土拱肋、波形钢腹板-混凝土组合拱肋以及钢筋混凝土拱上立柱的材料非线性。将桥面系简化为单主梁的形式, 桥面板和盖梁采用普通的梁单元来模拟, 二者之间的连接采用弹簧单元来模拟, 拱肋与立柱之间、立柱与盖梁之间均采用刚臂连接。在有限元分析中, 钢材的本构关系采用了双折线计算模型;混凝土采用不考虑抗拉的二次抛物线模型, 不考虑压应力下降。

有限元计算中, 自由振动计算采用子空间迭代法求解特征值;采用时程分析法进行非线性地震响应分析, 数值计算采用常加速度法来求解。用P-δ效应和动坐标法考虑大位移的影响。采用Rayleigh阻尼, 其两个阻尼系数均设为0.02, 两个频率根据地震作用方向的不同而设为不同值。

输入地震动为文献[13]中规定的板块边界型地震 (Type I, 如关东大地震) , Type I有三个波形, 本文在计算中采用T111地震波 (如图2所示) 。由于大跨度拱桥应同时考虑3个方向地震动的影响[14,15], 因此, 本文把地震波同时施加于桥梁三个方向 (横桥向、纵桥向、竖向) , 三个方向的地震动分量比为1:1:0.5。

2 结构自振特性

岭兜大桥 ( 混凝土箱拱) 和波形钢腹板- 混凝土组合拱桥前5 阶的固有振动频率、有效质量及振型分别见表1 和表2。 ( 为方便表述, 下文图表中统一采用RC表示混凝土箱拱, RC - CSW表示波形钢腹板- 混凝土组合拱。) 从表中可以看出, 两种拱桥的面外一阶固有振动均为对称模态, 且均早于面内一阶振动出现; 面外一阶反对称侧弯模态的基频分别为0. 903Hz ( RC) 和0. 407Hz ( RC - CSW) , 振型贡献率都非常小。同时, 波形钢腹板- 混凝土组合拱桥的面外一阶侧弯基频较之混凝土箱拱下降较多, 分别为26% ( 对称模态) 和55% ( 反对称模态) 。

对于面内振型, 二者的一阶竖弯振型均表现为拱桥特有的反对称振型, 自振频率分别为1. 003 Hz ( RC) 和1. 078 Hz ( RC - CSW) ; 二者的面内一阶对称模态出现的较早, 均在第4 阶时就出现了, 自振频率分别为1. 354 Hz ( RC) 和1. 421 Hz ( RC - CSW) 。由此可见, 二者的面内反对称和对称竖弯基频相差不到7% , 这说明波形钢腹板- 混凝土组合拱桥试设计中保证了桥梁的纵向刚度, 使得面内自振频率基本保持不变。

3 非线性地震响应特性

3. 1 横桥向地震作用时

当地震动T111 作用于横桥向时, 混凝土箱拱和波形钢腹板- 混凝土组合拱的主要截面 ( 拱脚和拱顶) 内力时程图见图3 ~ 图4。从图中可以看出, 两种拱桥横桥向受力均是以轴力和面外弯矩为主, 这与一般拱桥的抗震特性相符。同时由于拱肋截面产生的面外弯矩Mz比面内弯矩My大得多, 因此在横桥向地震动作用下不必考虑两方向弯矩相关的问题。也就是说, 可以采用轴力N与面外弯矩Mz的相关曲线对组合拱桥的横桥向非线性地震响应特征进行评价。

图5 ~ 图6 分别为横桥向地震动作用下两种拱桥主要截面的轴力N与面外弯矩Mz的相关曲线图。图中内圈曲线为截面中钢筋最外缘应变达到屈服应变时轴力与面外弯矩的包络线, 外圈曲线为截面中混凝土最外缘应变达到屈服应变时轴力与面外弯矩的包络线。从图中可清楚地看到, 在T111 地震动作用下, 两种拱桥的拱脚截面由于产生了较大的轴力和面外弯矩, 致使钢筋和混凝土最外缘应变均超过屈服强度, 从而使得拱脚截面进入塑性状态; 而拱顶截面在横向地震动作用下产生的轴力和面外弯矩响应都较小, 因此拱顶截面仍处于弹性阶段。

表3 为拱脚与拱顶截面的内力峰值对比情况。由于波形钢腹板- 混凝土组合拱的自重比混凝土箱拱小, 因此其在横桥向地震作用下产生的内力较后者有明显的降低。拱脚截面处, 波形钢腹板- 混凝土组合拱的最大面外弯矩、最大横向剪力和最大轴力分别为混凝土箱拱的95. 9% 、49. 4% 和81. 0% ; 而拱顶截面的面外弯矩, 波形钢腹板- 混凝土组合拱仅为混凝土箱拱的15. 7% 。

3. 2 纵桥向地震作用时

当地震波T111 作用于纵桥方向时, 由于地震动作用于纵桥向时产生的面外弯矩Mz较小, 两种拱桥纵桥向受力均是以轴力和面内弯矩为主, 这也符合一般拱桥的地震特性。

表4 为纵桥向地震作用下两种拱桥主要截面的内力峰值对比情况。从表中可以看出, 最大弯矩发生在拱脚处, 波形钢腹板- 混凝土组合拱的最大面内弯矩、最大竖向剪力和最大轴力分别为混凝土箱拱的95. 5% 、66. 7% 和85. 6% 。这说明尽管与横桥向地震动作用的影响相比, 在纵桥向地震作用下, 波形- 混凝土组合拱主要截面的内力响应较之混凝土箱拱减小的不多, 但采用波形钢腹板替代混凝土腹板仍然可以在一定程度上改善了拱桥的纵向抗震性能。

两种拱桥拱圈纵桥向和竖向的位移峰值对比情况见图7。从图中可以看出, 在三向地震作用下, 波形钢腹板- 混凝土组合拱的拱圈面内位移峰值与混凝土箱拱接近, 这也说明纵桥向地震动对拱圈的影响不如横桥向地震动来的大。

3. 3 纵、横桥向同时作用

当T111 地震动同时作用于纵、横桥向时, 两种拱桥拱脚截面产生的轴力N、面外弯矩Mz和面内弯矩My的时程曲线见图8 ~图9。为了便于比较, 图中也给出了单向地震动作用时的内力。从这些图中可以看出双向地震动同时作用时两种拱桥拱脚截面的轴力N和面外弯矩Mz与横桥向单独地震作用时的数值基本相等; 而面内弯矩My则是在纵桥向地震动作用下的响应比双向地震动同时作用时来得大。

表5 给出了双向与横桥向地震动作用下进入屈服的拱脚截面最外边缘应变最大值的对比情况。从表中可以看出, 混凝土箱拱在双向地震动同时作用时拱脚截面产生的最大应变是横桥向地震动单独作用时产生的应变的1. 17 倍。而与之不同的是, 波形钢腹板- 混凝土组合拱的拱脚截面在横桥向单独地震作用时产生的最大应变与双向地震动同时作用时产生的应变基本相同。这说明, 对波形钢腹板- 混凝土组合拱桥进行抗震设计时, 如果输入的地震波为板块边界型地震动, 则可直接对其施加横桥向地震动, 而不需要考虑纵、横向地震动同时作用的情况。

4 结论

( 1) 波形钢腹板- 混凝土组合拱的面外一阶对称和反对称侧弯基频较之混凝土箱拱分别下降了26% 和55%, 二者的面内一阶对称和反对称竖弯基频相差不到7%。这说明采用波形钢腹板替代混凝土腹板, 使得波形钢腹板- 混凝土组合拱的总体刚度较混凝土箱拱的刚度小。

( 2) 当T111 地震动单独作用于横桥向时, 波形钢腹板- 混凝土组合拱拱脚截面的最大面外弯矩、最大横向剪力和最大轴力分别为混凝土箱拱的95.9%、49.4%和81.0%, 而拱顶截面的面外弯矩, 前者仅为后者的15.7%。当T111地震动单独作用于纵桥向时, 最大弯矩也是发生在拱脚处, 波形钢腹板-混凝土组合拱的最大面内弯矩、最大竖向剪力和最大轴力分别为混凝土箱拱的95.5%、66.7%和85.6%。这表明在板块边界型地震动作用下, 采用波形钢腹板替代混凝土腹板可以改善拱桥的抗震性能。

( 3) 与混凝土箱拱不同的是, 波形钢腹板- 混凝土组合拱的拱脚截面在横桥向单独地震作用时产生的最大应变与纵、横向地震动同时作用时产生的应变基本相同。这说明在波形钢腹板-混凝土组合拱桥的非线性抗震计算中, 如果输入的地震波为板块边界型地震动, 则可直接对其施加横桥向地震动, 而不需要考虑纵、横向地震动同时作用的情况。

非线性动响应 篇6

齿轮系统的振动和噪声一直以来都是学界研究的重点。近年来,国内外学者以非线性振动理论为基础,以齿轮啮合过程中的时变刚度和齿侧间隙等非线性因素为核心,对齿轮系统的非线性振动进行了研究[1,2,3]。研究结果表明:齿轮摩擦也是齿轮非线性振动的重要的影响因素之一。Velex等[4]分析了直齿和斜齿的摩擦力。Vaishya等[5]从能量的角度证明了摩擦对齿轮振动性能的影响。He等[6]研究了实测刚度影响下齿轮摩擦对系统动力学的影响。

本文建立了周期双参变激励齿轮系统非线性振动模型,运用Floquet-Liapunov理论,用傅里叶级数将周期齿间摩擦力展开,用解析和数值相结合的方法计算得出:齿轮摩擦对系统的非线性振动有重要的影响。

1 齿轮摩擦非线性振动模型

齿轮摩擦非线性振动模型如图1所示。

1.1 系统参数假设

(1)载荷分配

$\begin{array}{l}{\rm N}{}_1(t)=\alpha {\rm N},{\rm N}{}_2(t)=(1-\alpha ){\rm N}t\in [0,t{}_{\text{a}})\\ {\rm N}{}_1(t)=0,{\rm N}{}_2(t)={\rm N}t\in [t{}_{\text{a}},t{}_{\text{c}})\end{array}\}(1)$

式中,N为法向载荷,N=Tp/Rbp,Tp为主动轮上的力矩,Rbp为主动齿轮的基圆半径;N1、N2分别表示啮合齿对1号和2号的法向载荷;ta为单齿啮合最低点的时间;tc为单齿啮合最高点的时间;α为载荷分配因子,0≤α≤1。

(2)啮合刚度k(t)

$\begin{array}{l}k(t)=k{}_0(1+\beta (t))t\in [0,t{}_{\text{a}})\\ \beta (t)=k{}_{\max }/k{}_{\min }-1t\in [0,t{}_{\text{a}})\\ \beta (t)=k{}_{\min }/k{}_{\max }-1t\in [t{}_{\text{a}},t{}_{\text{c}})\\ k{}_0=\frac{1}{\overline{{\rm T}}}{\displaystyle {\int }_0^{\overline{{\rm T}}}k}(t)\text{d}t=k{}_{\max }(\Gamma -1)+k{}_{\min }(2-\Gamma )\end{array}\}(2)$

式中,k0为平均啮合刚度;$\overline{{\rm T}}$为齿轮的啮合周期;kmax为双齿啮合刚度;kmin为单齿啮合刚度;Γ为齿轮的重合度。

(3)摩擦因子μ0

$\begin{array}{l}\mu {}_{01}=\mu {}_0\mathrm{sgn}(t-t{}_{\text{a}})t\in [0,t{}_{\text{a}})\\ \mu {}_{02}=-\mu {}_0\mathrm{sgn}(t-t{}_{\text{b}})t\in [t{}_{\text{a}},t{}_{\text{b}})\end{array}\}(3)$

(4)摩擦力臂ρ

$\begin{array}{l}\rho {}_{10}=R{}_{\text{b}\text{p}}\tan \phi {}_{\text{p}}\\ \rho {}_{1\text{A}}=R{}_{\text{b}\text{p}}\tan \phi {}_{\text{p}}+2\text{π}(\Gamma -1)R{}_{\text{b}\text{p}}\\ \rho {}_{20}=R{}_{\text{b}\text{p}}\tan \phi {}_{\text{p}}-{\rm P}\end{array}\}(4)$

式中,ρ10、ρ20分别为啮合齿对1号和2号的摩擦力臂;ρ1A为主动齿轮在啮合点处的摩擦力臂。φp为节圆压力角;P为法向齿距。

(5)摩擦力矩Tf

$\begin{array}{l}{\rm T}{}_{\text{f},1}=\alpha {\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t+{\rm P})t\in [0,t{}_{\text{a}})\\ {\rm T}{}_{\text{f},2}=-(1-\alpha ){\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [0,t{}_{\text{a}})\\ {\rm T}{}_{\text{f},1}=0,{\rm T}{}_{\text{f}‚2}=-{\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [t{}_{\text{a}},t{}_{\text{b}})\\ {\rm T}{}_{\text{f},1}=0,{\rm T}{}_{\text{f},2}={\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [t{}_{\text{b}},t{}_{\text{c}})\end{array}\}(5)$

式中,ωp为主动轮的角速度;tb为节点啮合时间;Tf,j为主动齿轮的摩擦力矩,j取1和2时分别表示啮合的齿对1号和2号。

(6)等效啮合阻尼c[7]

c=2ζ(kmR2bpR2bgJpJgR2bpJp+R2bgJg)12 (6)

式中,ζ为相对阻尼系数;Jp、Jg分别为主动齿轮、被动齿轮的等效转动惯量;Rbg为被动齿轮的基圆半径;km为轮齿的啮合刚度。

1.2 系统的数学模型

$\begin{array}{l}\ddot{x}+2\zeta (1+\mu {}_{1}p{}_1(\overline{\omega }\tau ))\dot{x}+\\ (1+\mu {}_{2}p{}_2(\overline{\omega }\tau ))x=f(\tau (\omega ))(7)\end{array}$

其中,x为振动位移,$\tau =\Omega t,\zeta =c{}_0/(2m\Omega ),\dot{x}=\text{d}x/\text{d}\tau ,\Omega {}^2=k{}_0/m,m$、c0分别为系统质量、等效平均啮合阻尼。Ω为系统的固有频率。ω、$\overline{\omega }$分别为系统参数的变化频率、系统的激励频率。$p{}_1(\overline{\omega }\tau )、p{}_2(\overline{\omega }\tau )$分别为系统参数的阻尼、刚度周期控制函数;μ1、μ2分别为系统参数的周期函数的幅值控制因子。f(τ(ω))为外界周期的激励。

引入状态向量$\mathit{X}(\tau )=(x(\tau ),\dot{x}(\tau )){}^{\text{Τ}}$,则式(7)可写成

$\dot{\mathit{X}}(\tau )=\mathit{G}(\tau )\mathit{X}(\tau )+\mathit{F}(\tau )(8)$

$\mathit{G}(\tau )=\left[\begin{array}{cc}0& 1\\ -(1+\mu {}_{2}p{}_2(\overline{\omega }\tau ))& -2\zeta (1+\mu {}_{1}p{}_1(\overline{\omega }\tau ))\end{array}\right](9)$

F(τ)=[0 f(τ)]T (10)

G为周期函数,$\mathit{G}(\tau +\overline{{\rm T}})=\mathit{G}(\tau ),\overline{{\rm T}}=2\text{π}/\overline{\omega }$。当f(τ)=0时,式(8)可转化为相应的齐次方程:

$\dot{\mathit{X}}(\tau )=\mathit{G}(\tau )\mathit{X}(\tau )(11)$

根据Floquet-Liapunov理论[8],如果xi(τ)是式(11)线性无关的基解(i=1,2),则系统朗斯基矩阵W(τ)为

$\mathit{W}(\tau )=[\mathit{X}{}_1(\tau )\mathit{X}{}_2(\tau )]=\left[\begin{array}{cc}x{}_1(\tau )& x{}_2(\tau )\\ \dot{x}{}_1(\tau )& \dot{x}{}_2(\tau )\end{array}\right](12)$

Xf0(t)=W(t)W-1(0)X(0)+∫t0W(t)W-1(τ)F(τ)dτ (13)

式中,Xf0(t)为在计及系统初始条件时,系统在外界激励下的响应。

X(0)是系统的初始条件,W-1(τ)是W(τ)的逆矩阵。在式(13)中常令

W(t)W-1(0)=Φ(t,0) (14)

Φ(t,0)表示系统(0,t)的状态转移矩阵。将式(14)代入式(13)得

Xf0(t)=Φ(t,0)X(0)+∫t0Φ(t,τ)F(τ)dτ (15)

如果不计及系统初始条件的瞬态影响,X(0)=0,系统在外界激励下的响应为

Xf(t)=∫t0Φ(t,τ)F(τ)dτ (16)

1.3 系统响应的计算

设周期双参变激励系统的两个线性无关的基解为

$\begin{array}{l}x{}_1(\tau )=\exp (-c{}_1^2\zeta \tau )\cos (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})\\ x{}_2(\tau )=\exp (-c{}_1^2\zeta \tau )\sin (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})\end{array}\}(17)$

$c{}_1=\sqrt{1+\mu {}_{1}p{}_1(\overline{\omega }\tau )},c{}_2=\sqrt{1+\mu {}_{2}p{}_2(\overline{\omega }\tau )}$

经证明:x1(τ)、x2(τ)为式(7)的两个线性无关的基解。

根据式(12),系统的朗斯基矩阵W(τ)为

$\mathit{W}(\tau )=\exp (-c{}_1^2\zeta \tau )\left[\begin{array}{cc}W{}_{11}& W{}_{12}\\ W{}_{21}& W{}_{22}\end{array}\right](18)$

$\begin{array}{l}W{}_{11}=\cos (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2}),W{}_{12}=\sin (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})\\ W{}_{21}=-c{}_1^2\zeta \cos (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})-\\ \sqrt{c{}_2^2-c{}_1^4\zeta {}^2}\sin (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})\\ W{}_{22}=-c{}_1^2\zeta \sin (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})+\\ \sqrt{c{}_2^2-c{}_1^4\zeta {}^2}\cos (\tau \sqrt{c{}_2^2-c{}_1^4\zeta {}^2})\end{array}$

状态转移矩阵Φ(τ2,τ1)为

$\begin{array}{l}\mathit{\Phi }(\tau {}_2,\tau {}_1)=\mathit{W}(\tau {}_2)\mathit{W}{}^{-1}(\tau {}_1)=\\ \exp (-c{}_1^2\zeta (\tau {}_2-\tau {}_1))\left[\begin{array}{cc}\text{ϕ}{}_{11}& \text{ϕ}{}_{12}\\ \text{ϕ}{}_{21}& \text{ϕ}{}_{22}\end{array}\right](19)\end{array}$

$\begin{array}{l}\text{ϕ}{}_{11}=\frac{c{}_1^2\zeta }{\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}}\sin (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))+\\ \cos (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))\\ \text{ϕ}{}_{12}=\frac{1}{\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}}\sin (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))\\ \text{ϕ}{}_{21}=-\frac{c{}_1^4\zeta {}^2}{\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}}\sin (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))-\\ \sqrt{c{}_2^2-c{}_1^4\zeta {}^2}\sin (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))\\ \text{ϕ}{}_{22}=-\frac{c{}_1^2\zeta }{\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}}\sin (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))+\\ \cos (\sqrt{c{}_2^2-c{}_1^4\zeta {}^2}(\tau {}_2-\tau {}_1))\end{array}$

对于分段周期系统[9]:

Φ(T,0)=Φ(T,tb)Φ(tb,ta)Φ(ta,0) 0≤t<tc (20)

Φ(t1,t2)=W(t2)W-1(t1) t1≤t<t2 (21)

当t=mT时,由Floquet-Liapunov理论[9]有:

Φ-1i(τ,0)=(Φ-1(T,0))i-1Φ${}_1^{-1}$(τ,0) (22)

Φ(mT,0)=(Φ(T,0))m(m为整数) (23)

将式(22)和式(23)代入式(16)得

$\begin{array}{l}\mathit{X}{}_{\text{f}}(\tau )=\mathit{\Phi }(m{\rm T},0){\displaystyle \sum _{i=1}^m{\displaystyle {\int }_{(i-1){\rm T}}^{i{\rm T}}\mathit{\Phi }}}{}_i^{-1}(\tau ,0)\mathit{F}(\tau )\text{d}\tau =\\ \mathit{\Phi }({\rm T},0){}^m{\displaystyle \sum _{i=1}^m(}\mathit{\Phi }{}^{-1}({\rm T},0)){}^{i-1}{\displaystyle {\int }_{(i-1){\rm T}}^{i{\rm T}}\mathit{\Phi }}{}_1^{-1}(\tau ,0)\mathit{F}(\tau )\text{d}\tau =\\ (\mathit{\Phi }(t{}_{\text{c}},0)){}^m{\displaystyle \sum _{i=1}^m(}\mathit{\Phi }{}^{-1}(t{}_{\text{c}},0)){}^{i-1}\mathit{{\rm H}}(\tau )(24)\end{array}$

$\begin{array}{l}\mathit{{\rm H}}(\tau )={\displaystyle {\int }_0^{t{}_{\text{a}}}\mathit{\Phi }}(\tau ,0)[B{}_0({\displaystyle \sum _{n=1}^{\infty }(}b{}_{n,1}+b{}_{n,2}+\\ b{}_{n,3})\sin \frac{n\text{π}\tau }{t{}_{\text{c}}}+{\displaystyle \sum _{n=1}^{\infty }c}{}_{n,1}\sin \frac{n\text{π}\tau }{t{}_{\text{c}}})+\\ C{}_0({\displaystyle \sum _{n=1}^{\infty }(}b{}_{n,1}^{*}+b{}_{n,2}^{*}+b{}_{n,3}^{*})\sin \frac{n\text{π}\tau }{t{}_{\text{c}}}+{\displaystyle \sum _{n=1}^{\infty }c}{}_{n,1}^{*}\sin \frac{n\text{π}\tau }{t{}_{\text{c}}})]\text{d}\tau +\\ \mathit{\Phi }(t{}_{\text{a}},0){\displaystyle {\int }_{t{}_{\text{a}}}^{t{}_{\text{b}}}\mathit{\Phi }}(\tau ,t{}_{\text{a}})[B{}_0{\displaystyle \sum _{n=1}^{\infty }c}{}_{n,2}\sin \frac{n\text{π}\tau }{t{}_{\text{c}}}+\\ C{}_0{\displaystyle \sum _{n=1}^{\infty }c}{}_{n,2}^{*}\sin \frac{n\text{π}\tau }{t{}_{\text{c}}}]\text{d}\tau +\mathit{\Phi }(t{}_{\text{a}},0)\mathit{\Phi }(t{}_{\text{b}},t{}_{\text{a}})\cdot \\ {\displaystyle {\int }_{t{}_{\text{b}}}^{t{}_{\text{c}}}\mathit{\Phi }}(\tau ,t{}_{\text{b}})[B{}_0{\displaystyle \sum _{n=1}^{\infty }c}{}_{n,3}\sin \frac{n\text{π}\tau }{t{}_{\text{c}}}+C{}_0{\displaystyle \sum _{n=1}^{\infty }c}{}_{n,3}^{*}\sin \frac{n\text{π}\tau }{t{}_{\text{c}}}]\text{d}\tau \end{array}$

式中,B0、C0分别为系统的摩擦响应系数,由齿轮系统的参数决定。

为了计算方便,将积分运算化为级数运算:

$\begin{array}{l}\mathit{{\rm H}}(\tau )={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }[}}B{}_0(b{}_{n,1}+b{}_{n,2}+b{}_{n,3}+c{}_{n,1})+\\ C{}_0(b{}_{n,1}^{*}+b{}_{n,2}^{*}+b{}_{n,3}^{*}+c{}_{n,1}^{*})]\mathit{\Phi }(\tau {}_i,0)\text{Δ}\tau {}_i\sin \frac{n\text{π}\tau {}_i}{t{}_{\text{c}}}+\\ \mathit{\Phi }(t{}_{\text{a}},0){\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }(}}B{}_{0}c{}_{n,2}+C{}_{0}c{}_{n,2}^{*})\mathit{\Phi }(\tau {}_i,t{}_{\text{a}})\text{Δ}\tau {}_i\sin \frac{n\text{π}\tau {}_i}{t{}_{\text{c}}}+\\ \mathit{\Phi }(t{}_{\text{a}},0)\mathit{\Phi }(t{}_{\text{b}},t{}_{\text{a}}){\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }(}}B{}_{0}c{}_{n,3}+\\ C{}_{0}c{}_{n,3}^{*})\mathit{\Phi }(\tau {}_i,t{}_{\text{b}})\text{Δ}\tau {}_i\sin \frac{n\text{π}\tau {}_i}{t{}_{\text{c}}}\end{array}$

齿轮系统的啮合周期与其他参数的关系为

ta=(Γ-1)tc (25)

$t{}_{\text{c}}=\frac{60}{\Pi {}_{\text{p}}{\rm Z}{}_{\text{p}}}(26)$

$t{}_{\text{b}}=\frac{30(\phi {}_{\text{p}}-\phi {}_0)}{\Pi {}_{\text{p}}\pi }(27)$

式中,Πp为主动轮的输入转速,r/min;Zp为主动轮的齿数;φ0为主动齿轮开始进入啮合时转过的角度。

齿轮系统的摩擦力矩为

$\begin{array}{l}{\rm T}{}_{\text{f},1}=\{\begin{array}{l}\alpha {\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+{\rm P}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [0,t{}_{\text{a}})\\ 0t\in [t{}_{\text{a}},t{}_{\text{b}})\\ 0t\in [t{}_{\text{b}},t{}_{\text{c}})\end{array}\\ {\rm T}{}_{\text{f},2}=\{\begin{array}{l}-(1-\alpha ){\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [0‚t{}_{\text{a}})\\ -\mu {}_0{\rm N}(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [t{}_{\text{a}},t{}_{\text{b}})\\ {\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [t{}_{\text{b}},t{}_{\text{c}})\end{array}\\ {\rm T}{}_{\text{f},1}^{*}=\{\begin{array}{l}\alpha {\rm N}(-\mu {}_0)(D-\phi {}_{0}R{}_{\text{b}\text{p}}-{\rm P}-\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [0,t{}_{\text{a}})\\ 0t\in [t{}_{\text{a}},t{}_{\text{b}})\\ 0t\in [t{}_{\text{b}},t{}_{\text{c}})\end{array}\\ {\rm T}{}_{\text{f},2}^{*}=\{\begin{array}{l}(1-\alpha ){\rm N}\mu {}_0(D-\phi {}_{0}R{}_{\text{b}\text{p}}-\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [0,t{}_{\text{a}})\\ \mu {}_0{\rm N}(D-\phi {}_{0}R{}_{\text{b}\text{p}}-\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [t{}_{\text{a}},t{}_{\text{b}})\\ {\rm N}\mu {}_0(D-\phi {}_{0}R{}_{\text{b}\text{p}}+\omega {}_{\text{p}}R{}_{\text{b}\text{p}}t)t\in [t{}_{\text{b}},t{}_{\text{c}})\end{array}(28)\end{array}$

式中,D为啮合轮齿理论啮合线的长度;T*f,j为被动齿轮的摩擦力矩,j取1和2时分别表示啮合的齿对1号与2号。

摩擦力矩的傅里叶级数展开为

${\rm T}{}_{\text{f},1}={\displaystyle \sum _{n=1}^{\infty }(}b{}_{n,1}+b{}_{n,2}+b{}_{n,3})\sin \frac{n\text{π}t}{t{}_{\text{c}}}(29)$

${\rm T}{}_{\text{f}‚2}={\displaystyle \sum _{n=1}^{\infty }(}c{}_{n,1}+c{}_{n,2}+c{}_{n,3})\sin \frac{n\text{π}t}{t{}_{\text{c}}}(30)$

$\begin{array}{l}b{}_{n,1}=\alpha {\rm N}\mu {}_0(\phi {}_{0}R{}_{\text{b}\text{p}}+p)\frac{2}{n\text{π}}(1-\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}})\\ b{}_{n,2}=-2\alpha {\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\frac{t{}_{\text{a}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}\\ b{}_{n,2}=-2\alpha {\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\frac{t{}_{\text{a}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}\end{array}$

$\begin{array}{l}c{}_{n,1}=(\alpha -1){\rm N}\mu {}_0\phi {}_{0}R{}_{\text{b}\text{p}}\frac{2}{n\text{π}}(1-\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}})+\\ (\alpha -1){\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\frac{2t{}_{\text{a}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}+\\ (\alpha -1)\frac{2}{t{}_{\text{c}}}{\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\sin \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}\\ c{}_{n,2}=\mu {}_0{\rm N}\phi {}_{0}R{}_{\text{b}\text{p}}\frac{2}{n\text{π}}(\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}-\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}})+\\ \frac{2\mu {}_0{\rm N}\omega {}_{\text{p}}R{}_{\text{b}\text{p}}}{n\text{π}}(t{}_{\text{b}}\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}-t{}_{\text{a}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}})+\\ \frac{2}{t{}_{\text{c}}}\mu {}_0{\rm N}\omega {}_{\text{p}}R{}_{\text{b}\text{p}}(\sin \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}-\sin \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}})\\ c{}_{n,3}=\mu {}_0{\rm N}\phi {}_{0}R{}_{\text{b}\text{p}}\frac{2}{n\text{π}}(\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}-\cos n\text{π})+\\ \mu {}_0{\rm N}\omega {}_{\text{p}}R{}_{\text{b}\text{p}}(\frac{2t{}_{\text{b}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}-\frac{2t{}_{\text{c}}}{n\text{π}}\cos n\text{π}-\frac{2}{t{}_{\text{c}}}\sin \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}})\end{array}$

${\rm T}{}_{\text{f},1}^{*}={\displaystyle \sum _{n=1}^{\infty }(}b{}_{n,1}^{*}+b{}_{n,2}^{*}+b{}_{n,3}^{*})\sin \frac{n\text{π}t}{t{}_{\text{c}}}(31)$

${\rm T}{}_{\text{f},2}^{*}={\displaystyle \sum _{n=1}^{\infty }(}c{}_{n,1}^{*}+c{}_{n,2}^{*}+c{}_{n,3}^{*})\sin \frac{n\text{π}t}{t{}_{\text{c}}}(32)$

$\begin{array}{l}b{}_{n,1}^{*}=\alpha {\rm N}\mu {}_0(D-\phi {}_{0}R{}_{\text{b}\text{p}}-p)\frac{2}{n\text{π}}(\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}-1)\\ b{}_{n,2}^{*}=2\alpha {\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\frac{t{}_{\text{a}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}\\ b{}_{n,3}^{*}=-\frac{2}{t{}_{\text{c}}}\alpha {\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\sin \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}\\ c{}_{n,1}^{*}=(1-\alpha ){\rm N}\mu {}_0(D-\phi {}_{0}R{}_{\text{b}\text{p}})\frac{2}{n\text{π}}(1-\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}})+\\ 2(1-\alpha ){\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\frac{t{}_{\text{a}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}+\\ \frac{2}{t{}_{\text{c}}}(\alpha -1){\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}\sin \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}\\ c{}_{n,2}^{*}={\rm N}\mu {}_0(D-\phi {}_{0}R{}_{\text{b}\text{p}})\frac{2}{n\text{π}}(\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}-\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}})+\\ 2{\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}(\frac{t{}_{\text{b}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}-\frac{t{}_{\text{a}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}}+\\ \frac{1}{t{}_{\text{c}}}\sin \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}-\frac{1}{t{}_{\text{c}}}\sin \frac{n\text{π}t{}_{\text{a}}}{t{}_{\text{c}}})\\ c{}_{n,3}^{*}={\rm N}\mu {}_0(D-\phi {}_{0}R{}_{\text{b}\text{p}})\frac{2}{n\text{π}}(\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}-\cos n\text{π})+\\ {\rm N}\mu {}_0\omega {}_{\text{p}}R{}_{\text{b}\text{p}}(\frac{2t{}_{\text{c}}}{n\text{π}}\cos n\text{π}-\frac{2t{}_{\text{b}}}{n\text{π}}\cos \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}}+\frac{2}{t{}_{\text{c}}}\sin \frac{n\text{π}t{}_{\text{b}}}{t{}_{\text{c}}})\\ p=\sqrt{1+\mu {}_2}\end{array}$

系统在啮合周期内的状态转移矩阵Φ(tc,0)为

$\begin{array}{l}\mathit{\Phi }(t{}_{\text{c}},0)=\mathit{\Phi }(t{}_{\text{c}},t{}_{\text{a}})\mathit{\Phi }(t{}_{\text{a}},0)=\\ \exp (-(d{}^2\zeta (t{}_{\text{c}}-t{}_{\text{a}})+b{}^2\zeta t{}_{\text{a}}))\left[\begin{array}{cc}\text{ϕ}{}_{11}^{(d)}& \text{ϕ}{}_{12}^{(d)}\\ \text{ϕ}{}_{21}^{(d)}& \text{ϕ}{}_{22}^{(d)}\end{array}\right](33)\end{array}$

$\begin{array}{l}\text{ϕ}{}_{11}^{(d)}=\cos (A(t{}_{\text{c}}-t{}_{\text{a}}))\cos Bt{}_{\text{a}}+\frac{d{}^2\zeta }{A}\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\cos Bt{}_{\text{a}}+\\ \frac{b{}^2\zeta }{B}\cos (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}+(\frac{d{}^{2}b{}^2\zeta }{AB}-\frac{b{}^4\zeta {}^2}{A\sqrt{p{}^2-d{}^4\zeta {}^2}}-\\ \frac{B}{A})\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}\\ \phi {}_{12}^{(d)}=\frac{1}{B}\cos (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}+\frac{1}{A}\sin (A(t{}_{\text{c}}-\\ t{}_{\text{a}}))\cos Bt{}_{\text{a}}+\frac{(d{}^2-b{}^2)\zeta }{AB}\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}\\ \text{ϕ}{}_{21}^{(d)}=-(A+\frac{b{}^2\zeta A}{B}+\frac{d{}^2\xi B}{A})\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}-\end{array}$

$(B+\frac{d{}^4\zeta {}^2}{B})\cos (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}-\frac{d{}^4\zeta {}^2}{A}\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\cos Bt{}_{\text{a}}$

$\begin{array}{l}\text{ϕ}{}_{22}^{(d)}=\cos (A(t{}_{\text{c}}-t{}_{\text{a}}))\cos Bt{}_{\text{a}}-\frac{b{}^2\zeta }{B}\cos (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}-\\ \frac{d{}^2\zeta B}{A}\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\cos Bt{}_{\text{a}}+\\ (\frac{d{}^{2}b{}^2\zeta {}^2}{A\sqrt{p{}^2-d{}^4\zeta {}^2}}-\frac{d{}^4\zeta {}^2}{AB}-\frac{A}{B})\sin (A(t{}_{\text{c}}-t{}_{\text{a}}))\sin Bt{}_{\text{a}}\end{array}$

$\begin{array}{l}b=\sqrt{1+\mu {}_1},d=\sqrt{1-\mu {}_1},\\ q=\sqrt{1-\mu {}_2},A=\sqrt{q{}^2-d{}^4\xi {}^2},B=\sqrt{p{}^2-b{}^4\xi {}^2}\end{array}$

2 算例分析

齿轮系统参数:输入转矩Tp=1.0×105N·m,k0=8.75×108N/m,Γ=1.65,Πp=20r/min,主动齿轮的齿数Zp=23,被动齿轮的齿数Zg=27,模数m=14mm,压力角α0=20°,齿宽B=300mm,ta=0.084s,tb=0.111s,tc=0.130s。在齿轮系统中,$\omega =\overline{\omega }$,为了便于分析,将系统的摩擦力作用下的振动响应区域分为双齿啮合区和单齿啮合区。

根据上述推导,下面分析系统各参数对摩擦响应的影响。为保证计算精度,摩擦力矩傅里叶级数n取10项。

2.1 系统阻尼因子μ1

定义系统惯性比η=Jp/Jg。由图2可看出:①系统最大振动幅值的绝对值发生在单齿啮合区,且发生在节点附近,当阻尼因子μ1,1=0.1时,系统的最大振幅的绝对值为1.0μm。②系统阻尼因子的绝对值对单齿啮合区最大振幅的绝对值的影响与双齿啮合区最大振幅的绝对值的影响机理并不相同。

(b)单齿啮合区

在双齿啮合区,增大阻尼因子的绝对值,并不能使响应的最大幅值的绝对值减小;而在单齿啮合区,增大系统的阻尼因子的绝对值,可以有效地抑制系统最大幅值的绝对值。图2a中,当阻尼因子μ1,1=0.1时,系统的最大振幅的绝对值为0.7μm;当阻尼因子μ1,5=0.5时,系统的最大振幅的绝对值为0.85μm,差值Δ=0.15μm。图2b中,当阻尼因子μ1,5=0.5时,系统的最大振幅的绝对值为0.25μm;当阻尼因子μ1,1=0.1时,系统最大振幅的绝对值为1μm,差值Δ=-0.75μm。由此可以看出,在周期双参变激励齿轮系统中,增大系统的阻尼因子可以有效地减小摩擦激励引起的最大振幅的绝对值。

2.2 系统刚度因子μ2

由图3可看出:①系统的最大振幅的绝对值发生在单齿啮合区,且发生在节点附近。②图3b中,减小系统的刚度因子并不能使系统的最大振幅减小。当刚度因子μ2,1=0.1时,系统的最大振幅的绝对值为1.0μm;当刚度因子μ2,5=0.5时,系统的最大振幅的绝对值为0.8μm。③在图3a中,随着系统刚度因子的增大,系统在双齿区的次峰值得到有效的抑制。当刚度因子μ2,1=0.1时,系统的次峰值为0.7μm;当刚度因子μ2,5=0.5时,系统的次峰值为0.2μm。由此可以看出,在周期双参变激励齿轮系统中,减小系统的刚度因子并不能有效地减小摩擦激励引起的最大振幅的绝对值。

(b)单齿啮合区

2.3 载荷分配因子α

由图4可看出:①系统响应的最大振幅的绝对值发生在双齿啮合区。②单齿啮合区的最大幅值的绝对值发生在节点附近。在单齿啮合区,载荷分配因子的变化对最大幅值的绝对值无影响。③在双齿啮合区,载荷分配因子的变化可以明显地影响系统响应的最大振幅的绝对值。由图4a可以看出,在双齿啮合时,当载荷分配因子α=0.5时,系统响应的最大幅值的绝对值最小。载荷分配因子α为0.1和0.3时,系统有较小的系统响应峰值。

(b)单齿啮合区

2.4 系统重合度Γ

由图5可看出:①系统响应的最大幅值的绝对值发生在单齿啮合区,且发生在节点附近。但是系统重合度的增大对系统最大幅值的绝对值的影响并不明显(图5b)。②在双齿啮合区,当系统的重合度Γ=1.3时,系统的峰值响应最大。③在双齿啮合区,当系统的重合度Γ分别为1.0,1.9,2.0时,系统的峰值响应较小。由机械原理知,重合度Γ≥1时可以通过对节圆至齿顶的齿面修形、增加齿轮的名义重合度来减小系统的峰值响应。

(b)单齿啮合区

2.5 系统惯量比η

由图6可看出:①系统响应的最大幅值发生在双齿啮合区,且两处响应的峰值的大小比较接近。②随着系统的惯量比接近于1,系统在双齿啮合区的峰值响应迅速减小。③单齿啮合区的最大响应发生在节点附近。当系统的惯量比η为0.9,1.0时,系统的峰值响应较小。

(b)单齿啮合区

2.6 摩擦因子μ0

由图7可看出:①系统响应的最大峰值的绝对值在单齿啮合区,且在节点附近。减小摩擦因子μ0可以有效地抑制系统响应的最大峰值。当μ06=0.10时,系统响应的最大峰值的绝对值为

(b)单齿啮合区

1.0μm(图7b);当μ01=0.01时,系统响应的最大峰值的绝对值为0.05μm。②在双齿啮合区,减小摩擦因子可以有效地降低系统的峰值响应。当μ06=0.1时,系统响应的最大峰值的绝对值为0.85μm;当μ01=0.01时,系统响应的最大峰值的绝对值为0.1μm。③在双齿啮合区,减小摩擦因子能较好地抑制次峰值响应的大小。

2.7 长期振动响应

运用递归法计算系统长期的振动响应。由图8可看出:在周期摩擦激励下,周期双参变激励系统的响应为周期响应(t为$10~13\overline{{\rm T}})$。

3 结论

本文基于Floquet-Liapunov 理论,为周期双参变激励非线性系统寻找一基本解系,将系统的周期外激励用傅里叶级数展开。通过计算研究发现此类非线性齿轮系统振动具有以下特征:

(1)在摩擦力的作用下,系统的最大响应峰值大都位于节点附近,而不是系统参数发生突变的地方。

(2)在摩擦力的作用下,系统的设计参数对系统的最大响应峰值有重要的影响。减小系统的摩擦因子,通过修形增大齿轮系统的重合度,增大系统的啮合阻尼,使齿轮副的转动惯量大致相等,可以有效地抑制由摩擦引起的振动。可以调整这些设计参数,从而有效地控制此类系统由摩擦激励引起的振幅响应。

(3)在周期摩擦力的作用下,周期双参变激励系统的响应为周期振动。

参考文献

[1]Ozguven H N,Houser D R.Mathematical ModelsUsed in Gear Dynamics-a Review[J].Journal ofSound and Vibration,1988,121(3):383-411.

[2]李润方,王建军.齿轮系统动力学—振动、冲击、噪声[M].北京:科学出版社,1997.

[3]王建军,李其汉,李润方.齿轮系统非线性研究进展[J].力学进展,2005,35(1):37-51.

[4]Velex P,Cahouet V.Experimental and NumericalInvestigations on the Influence of Tooth Friction inSpur and Helical Gear Dynamics[C]//Proc.ASME2000DET Conf..New York,2000:591-600.

[5]Vaishya M,Houser D R.Modeling and Analysis ofSliding Frictionin Gear Dynamics[C]//Proceedingsof the 2000 ASME Design Engineering TechnicalConference.Baltimore,2000:601-610.

[6]He S,Gunda R,Singh R.Effect of Sliding Frictionon the Dynamics of Spur Gear Pair with RealisticTime-varying Stiffness[J].Journal of Sound andVibration,2007,301(3/5):927-949.

[7]Lin H H,Hoston R L,Coy J J.On Dynamic Loadsin Parallel Shaft Transmissions I——Modeling andAnalysis[J].Transmissions of the ASME Journal ofMechanisms Transmissions and Automation in De-sin,1988,110(2):221-225.

[8]Richards J A.Analysis of Periodically Time-var-ying Systems[M].New York:Springer-Verlag,1983.

通用桥式起重机结构动响应研究 篇7

起重机结构动力学分析为起重机结构动态设计首先考虑的问题。目前研究起重机械动力学问题的方法有动载系数法、集中质量法[1,2]和有限元法[3,4,5]。其中有限元法借助现代计算机技术的发展, 已经成为模拟复杂工程问题强有力工具。于兰峰、王金诺[3,4]在频率域内建立基于有限元法的塔式起重机结构动态优化方法, 给出结构动态特性 (如固有频率) 对设计参数的灵敏度, 并据此选取结构优化设计变量。曾春[5]研究了QD 16—10.5桥式起重机的动态特性, 对起升机构采用两自由度简化模型, 估算上升启动和上升制动工况下起升钢丝绳内的动态力变化曲线。刘哲民[6]进行了MDG 30—42门式起重机的模态和抗倾覆性能分析。龙靖宇, 张轶蔚[7]采用多体动力学软件进行桥式起重机虚拟样机建模, 对起重机作业过程进行了动态仿真。然而, 对由钢丝绳、卷筒、小车架、桥架等组成的复杂系统的瞬态响应及稳态响应规律及随结构设计参数的变化规律, 相关系统的研究文献较少。

本文采用有限元法建立计入小车架质量和钢丝绳弹性的由吊重、钢丝绳、小车架、桥架组成复杂的桥式起重机动力学模型;参考文献[7]的工作, 以吊钩处的载荷时程曲线模拟吊重启动、起升、下降、制动工况;通过瞬态动力学分析和谐响应分析, 讨论桥式起重机结构参数和吊重下降速度对主梁跨中的瞬态响应的影响规律;并分析起重机金属结构的稳态频谱特性及随结构参数变化的规律。

1有限元模型

基于有限元的弹性系统的瞬态动力学方程为:

式 (1) 中[M]、[C]和[K]分别为系统的质量阵、阻尼阵和刚度阵;{x}为系统节点位移列向量;{F (t) }为系统的节点载荷列向量;特别地, 若节点位移和载荷向量有如下时间依赖关系

以某通用桥式起重机 (其主要参数见表1) 为研究对象进行结构动力学分析;其金属结构有限元模型如图1所示。结构参数按照表1的数据定义;不考虑小车架的弹性, 仅考虑小车架及其上面的卷筒、电机等的质量, 并按照质心近似的原则在图1中布置4个等效集中质量点, 程文明邓斌, 王金诺[1]指出通常将小车架假想为刚性, 在大多数情况下比较合理, 结果误差也较小。工程中, 主梁和端梁的连接出由于有高强度螺栓连接, 假定其各自的连接面板是刚性连接, 自由度完全耦合;三个水平方向的约束施加在大车轮的轴承座上。龙靖宇, 张轶蔚[7]进行了起重机载荷历程仿真, 给出吊钩拉板拉力时程曲线;本文在模拟起重机动力学问题时, 参照其工作, 在跨中起吊32t额定载荷, 假定吊钩处的载荷历程如图2所示;0—0.5s内将吊重从地面完全提起 (启动) ;0.5s—24.5s内匀速起升 (或匀速下降) ;24.5s—25s将吊重完全放到地面 (制动) 。

2结构参数对主梁瞬态响应的影响

本节采用有限元瞬态动力学模块分析该桥式起重机在图2所示载荷下主梁跨中的瞬态动响应, 分四种情况: (1) 模型尺寸按表1, (2) 端梁所有面板厚度增加到20mm, (3) 在 (2) 的基础上主梁面板厚度加倍, (4) 在 (2) 的基础上仅主梁下盖板厚度加倍;以此分析桥架结构参数改变对起重机结构瞬态响应的影响规律。

由图3, 在 (0—24.5) s内, 主梁跨中的动响应振幅较小, 而在 (24.5—25) s内时主梁跨中动响应最大;可见同样在0.5s时间间隔内, 吊钩处从零加载到满载[ (0—0.5) s内]与从满载卸载到零[ (24.5—25) s内]两种情况, 主梁动响应振幅差别较大, 卸载对主梁的动响应影响更大。在图3中, 各曲线对应不同的起重机金属结构参数;从中看出, 端梁加强可以降低主梁结构的动响应, 但不明显;在端梁加强基础上, 主梁下盖板厚度加倍, 可以明显降低匀速起升或下降阶段[ (0.5—24.5) s]的动响应, 由34.3mm降至23.6mm;而在端梁加强基础上, 主梁面板厚度加倍可以显著降低各个阶段的动响应, 如 (0.5—24.5) s阶段由34.3mm降至17.1mm。进一步观察图3, 在 (24.5—25s) 内, 结构加强 (如面板加厚、端梁加强) 可以降低主梁的动响应振幅。

图4为吊重下降速度对起重机主梁动响应的影响曲线;本节假定初始时刻系统内各节点的位移、速度、加速度都为零, 吊重下降速度在0.5s时间内分别由零线性提高到0.1m/s, 0.2m/s和0.4m/s, 计算主梁跨中的动响应。如图4, 随吊重下降速度的降低, 结构的动响应显著降低, 而且衰减的时间也明显缩短。

3结构参数对桥架频谱特性的影响

前节分析起重机在给定工况下, 结构参数对起重机主梁瞬态响应的影响规律;但不能得出哪些频率处, 结构振动更容易被激发, 而这是通常结构灵敏度分析和结构动态优化设计中比较重要的信息, 如于兰峰, 王金诺[3,4]在进行塔式起重机结构动态设计时, 通过频谱分析确定第二阶固有频率随设计变量的灵敏度。

在吊钩施加幅值为32t力的简谐激励, 频率范围为20Hz, 步长为0.2Hz, 计算起重机主梁跨中 (导电侧1 228点和非导电侧1 224点) 频率响应。如图5—图6, 可见无论何种形式的局部结构修改, 均不能明显峰值出现的频率 (首个共振峰在1.8Hz) ;但幅值却有明显变化。响应最小的两种情况分别是:主端梁加强连接;端梁所有面板加倍, 同时主梁下盖板厚度加倍;而其他情况, 无论侧重加强端梁或主梁, 其频谱幅值都比前者大得多;特别地若其他尺寸不变, 单纯加强下盖板, 对响应的影响远不如在主端梁连接处焊接加强板。由图5—图6还可以发现:1.8Hz处的共振峰峰值明显高于7.5Hz处的峰值, 建议在后续结构动态优化设计时应针对此频率进行灵敏度分析并确定优化目标。

4结语

本文采用有限元法建立计入小车架质量和钢丝绳弹性的由吊重、钢丝绳、小车架、桥架组成的桥式起重机动力学模型, 吊重启动、起升、下降、制动时吊钩处的载荷时程曲线参考文献[7]确定;通过瞬态动力学分析和频谱分析, 讨论了启动、起升、下降、制动额定载荷工况下主梁跨中的动响应规律及各结构参数的影响规律;吊重下降速度对结构动响应的影响规律;结构参数对起重机主粱频响特性的影响规律;可为起重机相关研究提供参考。

参考文献

[1]程文明, 邓斌, 王金诺.小车架为弹性结构时门式起重机的动态特性研究.西南交通大学学报, 2001;36 (2) :144—148

[2] Lightfoot E, Clarkson B L.Dynamic stresses in E.O.T.cranes dueto the hoisting and lowering of loads.Material Handling Mechanics, 1955;169 (10) :223—232

[3]于兰峰, 王金诺.塔式起重机结构系统动态优化设计.西南交通大学学报, 2007;42 (2) :206—210

[4]于兰峰, 塔式起重机结构刚性及动态优化研究.西南交通大学博士论文, 2006

[5]曾春.基于ANSYS的桥式起重机的桥架结构有限元动态结构分析研究.武汉:武汉理工大学硕士论文, 2006

[6]刘哲民.30 t_42 m门型起重机力学分析及结构参数优化.长沙:中南大学硕士论文, 2007