机电耦合动力学

机电耦合动力学（共9篇）

机电耦合动力学 篇1

摘要：以弹性薄板小挠度弯曲理论为基础, 采用林滋泰德—庞加莱法对微板机电耦合动力学响应进行摄动, 求出满足精度要求的动力学响应解。在算例分析中, 以静电微泵泵膜为研究对象, 进行了实例计算。计算结果表明, 这种非线性分析方法得到的动力学响应在精度上更加具有优势, 能够有效增加系统设计的合理性, 提高设计精度。

关键词：微机电系统,微板,机电耦合,非线性动力学,林滋泰德—庞加莱法,动力学响应

0 引言

微板作为一种典型的微机电系统 (micro-electromechanical system, MEMS) 结构元件, 在现代工程技术中得到了广泛的应用, 微型阀、微型泵、微传感器等微机械器件设计中常采用各种微板结构[1]。微机械通常采用静电力驱动, 其构件受到静电力的作用而发生变形, 而结构的尺寸和各构件的空间位置又会反过来影响电场的分布, 这是一种很强的耦合作用, 静电力的变化和结构的变形之间并非简单、规则的线性关系, 因此这也是一个典型的非线性问题[2,3]。通常采用多能量域耦合方程进行仿真计算, 但该方法计算量大, 且大大降低了系统设计精度。

为解决此问题, 本文以泊松-克希霍夫 (Poisson-Kirchhoff G.) 建立的弹性薄板小挠度弯曲理论为基础, 建立了静电力作用下微板机电耦合结构模型, 采用林滋泰德—庞加莱法 (Lindstedt-Poincaré) 求解其动力学响应, 进行了模态分析。以静电微泵泵膜[4,5]为研究对象, 进行了实例计算。

1 微板机电耦合动力学方程的推导

图1所示为静电力作用下微板机电耦合结构模型, 包含机械系统、电系统和耦合部分, 机械部分是一个施加电场力的微板, 电场力沿微板横向均匀分布;电系统包括电源、电阻和电容。机械系统和电系统通过电场力达到机电耦合状态。微板遵循弹性小挠度理论的基本假定:①直法线假设;②忽略沿中面垂直方向的法向应力;③无沿中面内方向的变形;④只计入质量的移动惯性力, 略去其转动惯性力矩[6]。

基于以上基本假设, 以中面上的矩形微元代替微元体, 表示微元的受力情况, 如图2所示, q (x, y, t) 为分布在单位面积上的静电场力, -ρ h∂2w/∂t2为单位面积微板的惯性力, Mx、My为弯矩, Mxy为扭矩, Qx、Qy为剪力, 根据文献[6]可得微板运动方程:

$\frac{\partial Q{}_x}{\partial x}+\frac{\partial Q{}_y}{\partial y}+q-\rho h\frac{\partial {}^{2}w}{\partial t{}^2}=0$ (1)

剪力表达式为

$Q{}_x=-D\frac{\partial }{\partial x}(\frac{\partial {}^{2}w}{\partial x{}^2}+\frac{\partial {}^{2}w}{\partial y{}^2})=-D\frac{\partial }{\partial x}\nabla {}^{2}w$ (2)

$Q{}_y=-D\frac{\partial }{\partial y}(\frac{\partial {}^{2}w}{\partial x{}^2}+\frac{\partial {}^{2}w}{\partial y{}^2})=-D\frac{\partial }{\partial y}\nabla {}^{2}w$ (3)

$\begin{array}{l}\nabla {}^2=\frac{\partial {}^2}{\partial x{}^2}+\frac{\partial {}^2}{\partial y{}^2}\\ D=\frac{Eh{}^3}{12(1-\nu {}^2)}\end{array}$

式中, ∇2为拉普拉斯算子;ρ为微板质量密度;D为微板抗弯刚度;ν为泊松比;E为弹性模量;w (x, y, t) 为微板中面挠曲函数;h为上层微板厚度。

由式 (1) ～式 (3) 得到直角坐标系中微板横向振动基本微分方程:

$\frac{\partial {}^{4}w}{\partial x{}^4}+2\frac{\partial {}^{4}w}{\partial x{}^2\partial y{}^2}+\frac{\partial {}^{4}w}{\partial y{}^4}+\frac{\rho h}{D}\frac{\partial {}^{2}w}{\partial t{}^2}=$

∇2∇$2w+\frac{\rho h}{D}\frac{\partial {}^{2}w}{\partial t{}^2}=\frac{q(x,y,t)}{D}$ (4)

当在两板之间施加工作电压时, 固定电极与上层微板之间形成静电场, 忽略电场的边缘部分, 可以将微板机电耦合系统视为平行板电容器, 根据文献[7], 单位面积静电场力大小为

$q{}_{\text{e}}=\frac{1}{2}\frac{U{}^2\epsilon {}_0}{(t{}_0-w+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^2}$ (5)

式中, ε0为真空介电常数, ε0=8.85×10-12C2/ (N·m2) , εr为相对介电常数;t0为两微板之间的初始间隙大小;dc为上层微板内壁的涂层厚度;U为两板之间施加的工作电压。

将电容C、电压U、单位面积静电场力qe表示为静态和动态两部分, 电容对动态间隙求导, 忽略高阶无穷小量, 得到动态电场力密度Δq1:

$\text{Δ}q{}_1=\frac{U{}_0^2\epsilon {}_0}{(t{}_0-w{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3}\text{Δ}w$ (6)

式中, U0为平均电压值;Δw为上层微板的主振动振幅;w0为上层微板的初始静态位移。

将式 (6) 代入式 (4) 得到微板机电耦合的动态弯曲振动方程:

$\begin{array}{l}\frac{\partial {}^4\text{Δ}w}{\partial x{}^4}+2\frac{\partial {}^4\text{Δ}w}{\partial x{}^2\partial y{}^2}+\frac{\partial {}^4\text{Δ}w}{\partial y{}^4}+\\ \frac{\rho h}{D}\frac{\partial {}^2\text{Δ}w}{\partial t{}^2}=\frac{\text{Δ}q{}_1(x,y,t)}{D}(7)\end{array}$

考虑电场力的非线性因素, 将动态电场力密度在$w=\overline{w}{}_0$静态平均位移处进行泰勒级数 (Taylor) 展开:

$\begin{array}{l}\text{Δ}q{}_1=\frac{U{}_0^2\epsilon {}_0\text{Δ}w}{(t{}_0-\overline{w}{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3}+\\ \frac{3U{}_0^2\epsilon {}_0\text{Δ}w{}^2}{2(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3\overline{w}{}_0}(\frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}})[1+\\ \frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}}+(\frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}}){}^2+\cdots ]{}^4+\\ \frac{2U{}_0^2\epsilon {}_0\text{Δ}w{}^3}{(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^4\overline{w}{}_0}(\frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}})[1+\\ \frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}}+(\frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}}){}^2+\cdots ]{}^5+\cdots (8)\end{array}$

引进与变量无关的实数小参数$\epsilon ,\epsilon =\frac{\overline{w}{}_0}{t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}},0＜\epsilon ＜1$, 得

$\text{Δ}q{}_1=\frac{U{}_0^2\epsilon {}_0\text{Δ}w}{(t{}_0-\overline{w}{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3}+\frac{3U{}_0^2\epsilon {}_0\text{Δ}w{}^2}{2(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3\overline{w}{}_0}\epsilon +$

$\frac{2U{}_0^2\epsilon {}_0\text{Δ}w{}^3}{(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^4\overline{w}{}_0}\epsilon +\cdots$ (9)

将式 (9) 代入式 (7) , 忽略高阶无穷小量, 得到系统振动方程:

$\begin{array}{l}\frac{\partial {}^4\text{Δ}w}{\partial x{}^4}+2\frac{\partial {}^4\text{Δ}w}{\partial x{}^2\partial y{}^2}+\frac{\partial {}^4\text{Δ}w}{\partial y{}^4}+\frac{\rho h}{D}\frac{\partial {}^2\text{Δ}w}{\partial t{}^2}=\\ \frac{U{}_0^2\epsilon {}_0\text{Δ}w}{D(t{}_0-\overline{w}{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3}+\frac{3U{}_0^2\epsilon {}_0\text{Δ}w{}^2}{2D(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3\overline{w}{}_0}\epsilon +\end{array}$

$\frac{2U{}_0^2\epsilon {}_0\text{Δ}w{}^3}{D(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^4\overline{w}{}_0}\epsilon$ (10)

采用分离变量法求解, 将系统的主振动写成模态函数ϕ (x, y) 与广义坐标q′ (t) 乘积的形式:

Δw (x, y, t) =ϕ (x, y) q′ (t) (11)

将式 (11) 代入式 (10) 并化简得到关于时间函数的动力学方程:

${\ddot{q}}^{\prime }(t)+q^{\prime }(t)(\omega {}_0^2-\frac{b{}_1}{\rho h}\epsilon q^{\prime }(t)-\frac{b{}_2}{\rho h}\epsilon q^{\prime }{}^2(t))=0(12)$

我们认为作用力的非线性因素不影响模态函数ϕ (x, y) , 设ω0为不计非线性因素时系统的固有圆频率, 令:

$\begin{array}{l}{\rm P}=\frac{U{}_0^2\epsilon {}_0}{(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}-\overline{w}{}_0){}^3}b{}_1=\frac{3U{}_0^2\epsilon {}_0\overline{\text{ϕ}}}{2(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^3\overline{w}{}_0}\\ b{}_2=\frac{2U{}_0^2\epsilon {}_0\overline{\text{ϕ}}{}^2}{(t{}_0+d{}_{\text{c}}/\epsilon {}_{\text{r}}){}^4\overline{w}{}_0}\overline{\text{ϕ}}\approx \frac{1}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^b{\displaystyle \sum \text{ϕ}}}}{}_i(x,y)\text{d}x\text{d}y\\ \overline{\text{ϕ}}{}^2\approx \frac{1}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^b{\displaystyle \sum \text{ϕ}}}}{}_i^2(x,y)\text{d}x\text{d}y\\ \frac{{\rm P}\text{ϕ}(x,y)-D\nabla {}^2[CX]\nabla {}^2[CX]\text{ϕ}(x,y)}{\rho h\text{ϕ}(x,y)}=-\omega {}_0^2\end{array}$

当小参数ε=0时, 式 (12) 退化为原系统的派生系统:

${\ddot{q}}^{\prime }(t)+\omega {}_0^2{q}^{\prime }(t)=0$ (13)

2 微板机电耦合动力响应计算

根据文献[8], 采用林滋泰德—庞加莱法进行求解, 将原系统的解展开成ε的幂级数:

q′ (t, ε) =q′0 (t) +ε q′1 (t) +ε2q′2 (t) +… (14)

将原系统的自由振动圆频率ω也展开成ε的幂级数:

ω2=ω${}_0^2$ (1+ε σ1+ε2σ2+…) (15)

系统的初始位移激励为Q0, 各方程的初始条件为

将式 (14) 、式 (15) 代入式 (12) , 引入新的自变量φ=ω t, 将原来的微分符号改定义为对φ的微分, 并令ε的同次幂的每一项系数为零, 导出各阶近似的线性方程组, 利用各方程的初始条件式 (16) 进行迭代计算, 最终推导出能够满足所需精确度的时间函数周期解:

得到非线性系统自由振动圆频率与振幅的关系式:

由于ε2项对ω数值的影响很小, 将式 (18) 简化为

$\omega {}^2=\omega {}_0^2(1-\epsilon \frac{3b{}_{2}Q{}_0^2}{4\rho h\omega {}^2})$ (19)

3 算例分析

以上述理论为基础, 以静电微泵泵膜为研究对象, 求解其动力学响应。为了增加泵腔容积, 提高微泵的输出效率和运动灵敏度, 泵膜采用柔性铰链结构, 泵膜可以简化为图1所示的四边简支微板机电耦合系统, 设计参数如表1所示, 由式 (11) 、式 (17) 可得微板机电耦合自由振动的各阶位移的时域动态响应 (图3) 以及幅频响应特性关系曲线 (图4) , 在图5中, 将非线性分析方法求解的位移时间响应与传统的线性分析结果进行了对比。

以静电微泵泵膜为研究对象, 对四边简支边界条件下微板机电耦合的非线性自由振动进行了定量分析, 计算结果如图3～图5所示, 结果表明:①阶次不同, 微板振动的峰值位置和位移函数的变化周期是不同的, 高阶模态固有圆频率值越大, 达到振动峰值的点越多, 振动周期越小;②非线性解的动态响应出现了不对称现象, 这是由于非线性解中的非线性因素使自由振动出现了直流和谐波成分, 使得响应不再对称于系统的平衡位置, 这种动态响应更加符合微机械系统实际的工作状态;③考虑了非线性因素的幅频特性曲线并非直线而是朝频率减小的方向弯曲, 呈现出软弹簧特性, 这种特性与非线性因素ε有关;④振动系统的自由振动圆频率不仅与系统本身属性有关, 还与振幅的变化相关, 此点也不同于线性系统, 在进行微系统动态性能设计时应当给予特别的注意;⑤由时间函数的周期解还可以看出, 周期解中除基频为ω的谐波以外, 还有频率为3ω、5ω的高次谐波存在, 在进行微系统动态性能设计时还要注意避免倍频共振的发生。

4 结论

(1) 采用非线性分析方法求解微板机电耦合系统的动力学响应是合理的, 非线性方法更加符合微泵等微机械器件的实际工作状态, 比采用传统的线性分析方法更为准确。

(2) 得到微板系统的机电耦合非线性动力学方程, 为静电驱动MEMS系统典型元件的动力学设计提供了理论基础。

(3) 采用林滋泰德—庞加莱法求得的非线性多项式解可以校正线性分析方法的计算误差, 提高设计精度。对于静电场较强的MEMS系统, 采用非线性分析方法具有更加重要的意义。

(4) 本文针对四边简支微板的非线性动力响应进行了分析, 对于其他条件下的微板, 只要选取适当的边界条件, 同样可以得到对应的非线性动力响应解。

(5) 这种非线性分析方法不仅适用于微泵的系统设计, 也适合微传感器等具有其他边界条件微板结构的微机械器件的设计。

参考文献

[1]Kazanc Z, Mecitoglu Z.Nonlinear Dynamic Behaviorof Simply Supported Laminated Composite PlatesSubjected to Blast Load[J].J.Sound and Vibration, 2008, 19 (3) :883-897.

[2]杨卫.微纳米尺度的力学行为[J].世界科技研究与发展, 2004, 26 (4) :2-6.

[3]王超, 郭早阳.微机电系统应用中的非线性力学问题分析[J].机电工程技术, 2005, 34 (8) :18-20.

[4]杜立群, 杨季玲, 徐征.无阀微泵流体特性的数值计算和结构优化设计[J].中国机械工程, 2005, 16 (7) :116-118.

[5]沙菁丰女刀, 侯丽雅, 章维一.数字化无阀微泵内微流体泵送运动的理论解析[J].中国机械工程, 2007, 18 (20) :2472-2476.

[6]曹志远.板壳振动理论[M].北京:中国铁道出版社, 1989:13-25.

[7]Herbert H W, James R M.机电动力学[M].华中机电教研室, 译.北京:机械工业出版社, 1982.

[8]刘延柱, 陈文良, 陈立群.振动力学[M].北京:高等教育出版社, 1998.

机电耦合动力学 篇2

摘要： 利用改进小参数平均法，推导出了平面运动双质体振动系统中两偏心转子的无量纲耦合方程，并依据其零解的存在及稳定性，得出了实现同步与同步稳定性运行的条件，给出了系统负载系数及同步能力系数定义。通过数值分析，讨论了系统动力学参数对耦合同步特性的影响，得到了双偏心转子自同步的稳定运行的参数区间。结果表明：随着系统动力学参数的变化，同步极值点既可能是负载系数的极小值点，也可能是负载系数的极大值点。其原因是由于双质体的耦合作用，来自于机体运动的同步力矩，在一定参数范围内驱动两个偏心转子相位差向负载系数最大值点趋近；而在此区域以外，驱动两偏心转子相位差向负载系数极最小值点趋近。通过数值仿真，验证了耦合分析结果的正确。

关键词： 双质体振动系统； 自同步； 稳定性； 耦合动力学

中图分类号： TH113.1； O347.6 文献标志码： A 文章编号： 1004-4523（2016）03-0521-11

DOI：10.16385/j.cnki.issn.10044523.2016.03.019

引 言

振动系统中偏心转子的自同步现象为振动技术提供了一个新的应用领域[12]，并导致了振动利用工程建立[35]。

最早研究同一刚体上偏心转子同步理论的是前苏联学者Blekhman，他利用小参数PoincareLyapunov方法与运动稳定性理论，从物理上解释了两偏心转子的自同步现象，由此发展成振动系统中偏心转子自同步分析的基于直接运动分离的小参数平均法[12]。在20世纪70年代始，中国学者闻邦椿院士将双偏心转子的运动微分方程合并为相位差扰动参数的微分方程，简化了系统的稳定性分析，发明了一系列的振动机械，创建了振动利用工程新学科[35]。作者通过设置两个偏心转子平均转速与相位差的扰动量，推导出偏心转子的无量纲耦合方程，改进了小参数平均法，将偏心转子的同步问题转换为无量纲耦合方程的零解存在及稳定性问题[67]。但以上内容主要是针对过共振系统中两个或多个偏心转子的同步问题进行研究的。在此类振动系统响应中，系统负载总是与阻尼相位角正弦量项有关，阻尼相位角余弦量项总是影响系统的同步力矩与稳定性。其原因是阻尼相位角接近于0，与阻尼相位角正弦量相关项对称稳定性影响可以忽略，进而在系统稳定性分析过程中，可以忽略阻尼相位角正弦量的耦合力矩对系统稳定性的影响[610]。

机电耦合动力学 篇3

近年来,以形状记忆合金、电致伸缩材料、磁致伸缩材料和压电材料主导的智能材料获得了迅速发展,其中的压电材料成为学者研究的焦点[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 结论

机电耦合动力学 篇4

（重庆大学资源及环境科学学院，重庆 400030）

引 言

旋转结构的动特性和动力响应分析是进行旋转机械和旋转结构可靠性设计和安全性评价的理论基础。进行风力发电机风轮、直升机螺旋桨、离心机叶轮及其叶片设计时，了解旋转结构的特征频率变化既是材料选取和结构优化的依据，也为避免结构共振实施控制策略提供关键参数。旋转结构的动力分析涉及结构的变形以及应力状态等，为结构变形控制和强度校核、确定动态测试方案提供基础的分析数据。由于旋转结构运行在离心环境下，结构自身的弹性变形与旋转运动耦合在一起，造成离心环境严重影响结构的动特性和动力学响应。譬如风力机风轮的坎贝尔图是关于风轮的各阶频率与风轮转速的关系［1］，可为风轮共振分析和风电机组控制提供设计和分析参数。因为旋转结构动力学在实际工程中的巨大应用价值，几十年来引起了学术界和工程界的持久兴趣。直到最近，大量的文献依然关注旋转的杆梁结构、旋转环、旋转盘和板、旋转壳的振动特性和响应［2～23］，其中旋转结构的科氏力效应，动力刚化问题，旋转角速度对结构特征频率的影响，复值的特征矢量和复模态处理等是研究重点，特别是旋转和变形的耦合机制还不完全清楚，解决这些问题有赖于揭示刚体旋转与弹性变形的耦合机理，建立精确可靠的动力学分析模型，发展合适的算法以及实验或工程验证。

旋转结构较之非旋转结构的复杂性在于结构进一步承受由于运动耦合带来的附加惯性力，这种附加惯性力影响结构动力学，同时结构变形影响结构的刚性旋转运动。为揭示运动变形耦合以及对结构动力学造成的影响，本文引入双自由度的质量弹簧振动系统作已知的刚体旋转，系统可考虑刚体旋转形成的初始离心力、质量点振动造成的离心力改变以及科氏力等。这种简化模型既不失去结构变形与刚体旋转耦合的典型特性，又可通过考察简单物理模型来探讨动力学规律，为复杂结构的动力学分析提供清晰的参考。

1 双自由度离心振动系统的动力学模型

考虑刚体绕定轴作匀速圆周运动，转动的角速度为ω，刚体上任意一点r的线速度v为

二阶反对称张量Ω可表示旋转轴和绕轴旋转角速度，即有［24］

式中e为笛卡尔坐标系下的置换张量，角速度ω为反对称张量Ω的反偶矢量或轴矢量。二阶反对称张量Ω对任意矢量u的映射为

二阶反对称张量只关联着转动角速度，正交张量能描述刚体作定轴有限转动的旋转角速度和旋转角。现令单位矢量n代表旋转轴方向，θ为旋转角，引进反对称张量A＝－e·n，对正交张量R有［25］

上式为Euler－Rodrigues旋转公式。对匀速转动，正交张量R可进一步改写为

ω为角速度矢量ω的模。

现在关注一个质量点m与4个相同弹簧组成的双自由度振动系统作平面刚体旋转运动（忽略振动造成的离心振动系统转动惯量的变化），振动系统固定在刚体上以垂直于旋转平面的恒角速度ω作逆时针旋转。建立固定和转动两个坐标系，见图1，质量点的初始位置r0为转动坐标系的原点，两个坐标基分别沿转动的径向和切向，质量点在转动坐标系的振动位移为r，则质量点在转动坐标系的位置为r0＋r。由于正交张量描述了刚体有限转动，则任意时刻质量点在固定坐标系的位置为

图1 双自由度振动系统作定轴平面转动Fig．1 The vibrating system with two degrees of freedom rotates around with a fixed axis

这里正交张量R代表转动坐标系与固定坐标系构成的两点张量。对式（6）取时间导数，利用式（3）和（5）得质量点在固定坐标系的绝对速度为

再对式（7）取时间导数并再次应用式（3）和（5），得固定坐标系下质量点的绝对加速度为

式（7）和（8）中˙r和¨r表示质量点在转动坐标系的振动速度和加速度。利用正交张量的转置张量RT立即得到在转动坐标系下质量点的绝对加速度

式中ω×（ω×r0）为振动系统的初始向心加速度、ω×（ω×r）是由于质量点振动导致的向心加速度增量、2ω×˙r为科氏（Coriolis）加速度。由此可知，如果质量点在垂直于旋转平面的方向上运动，不产生科氏加速度和向心加速度的增量。

应用达朗贝尔原理，忽略质量点重力作用，可建立质量点m在转动坐标系的动力学方程

k为弹簧的弹性系数。由式（10）可知，质量点m受到弹簧回复力、相对惯性力、初始离心力、离心力增量和科氏力的共同作用，注意到质点动力学方程（10）与旋转角θ无关。由于质量弹簧系统的刚体转动和振动都在旋转平面上（图1），在转动坐标系中有

u和v表示在转动坐标系质量点的两个位移分量，进一步令r0的模为r0。利用式（11）以及矢量的二重叉积公式，改写式（10）得到沿转动坐标系两个坐标轴e1和e2的动力学方程

把式（12）和（13）联立成为矩阵形式

改写式（14）得离心振动系统动力学方程的一般形式

式（15）中质量矩阵M，科氏矩阵C，振动刚度矩阵K，离心矩阵KS分别为：

科氏矩阵C为反对称矩阵，尽管在动力学方程出现一阶项但不造成系统的阻尼衰减。式（15）也是线弹性情况下旋转结构动力学方程的典型形式。观察式（14）可知在旋转平面内，由于刚体旋转与弹性振动的运动耦合产生科氏力，导致u和v两个方向的运动也是耦合的，而运动耦合也反之影响科氏力，这正是离心振动系统复杂性的本质。此外，离心振动系统刚度阵K－KS的系数应为正，表明为保持系统稳定应存在最大刚体转速，下面将进一步说明。

2 双自由度离心振动系统的动频和动模态

为研究离心振动系统的动力学性质，将动力学方程（15）改写成下述形式

记zT＝｛u˙uv˙v｝T，则式（17）在形式上成为一阶非齐次线性微分方程组

式（18）中的T和g取与式（17）右端对应的方阵和列阵。方程组（18）的齐次方程解的形式为z＝Cφeλt，对应的特征方程为｜T－λI｜＝0，容易得到系统的特征值为

ω1＝分别表示振动系统两阶模态的圆频率，而f1和f2分别为随转速变化的两阶动频。作为例子，令m＝0．1kg，k＝10 N／m。系统两阶特征频率随转速变化的曲线见图2。当转速为零时刚体旋转对振动的影响消失，系统的两阶动频退化为系统的固有频率（基频）。随刚体转速增加系统的第一阶动频f1线性减小而第二阶动频f2线性增加，当转速达到固有频率时，第一阶动频趋于零，表示刚体转动造成了振动系统失稳，所以，刚体的最大转速等于质点振动的固有频率。

图2 离心振动系统的两阶动频与旋转运动（转速）的关系Fig．2 Two order dynamic frequencies of the centrifugal vibrating system vary with rotational velocity

对式（18）的齐次形式以及式（19）的特征值，可以确定特征值对应的广义特征矢量。令λ1＝－ω1i，λ2＝ω1i，λ3＝－ω2i，λ4＝ω2i，所对应的广义特征矢量分别为

特征值λ1，λ2，λ3，λ4所对应的复值解分别为

利用线性微分方程组解的可叠加性，对每组复值解作简单的代数运算可得到每个特征值所对应的线性独立的实值解，即为两阶模态，它们分别为

两阶模态都是刚体转速ω的函数。κ1和κ2表示系统的第1阶动模态，κ3和κ4表示系统的第2阶动模态。注意到这些实值解列阵中的第1行和第3行分别为质量点的两个模态位移分量，显然，两个模态位移分量的平方和均为1，表明系统的两阶动模态都是半径为1的极化圆周运动，模态圆周运动以及模态位移的初始位置（相位）见图3。系统第1阶动模态的转动方向为逆时针，圆周运动的角速度即为系统第1阶圆频率ω1。第2阶动模态的转动方向为顺时针，圆周运动的角速度即为系统第2阶圆频率ω2。第1阶动模态的转动速度随系统转速ω的增加变得越来越慢，直至趋于第1阶动频的最小值，而第2阶动模态的转动速度随转速ω的增加会变得越来越快。

图3 离心振动系统两阶动模态的极化圆周运动及其相位Fig．3 The polarized circular motion with the two order dynamic modal and the phase of the centrifugal vibrating system

3 双自由度离心振动系统的质点运动轨迹

根据式（23）可写出方程组（18）的齐次形式的基解矩阵为

由于Z（0）≠I，可以选取基解矩阵为Z′（t）＝Z（t）Z－1（0），则非齐次微分方程组（18）的通解形式为［26，27］

将式（24）带入到式（25）可得

令振动系统的初始条件为t＝0时，u（0）＝u0，˙u（0）＝˙u0；v（0）＝v0，˙v（0）＝˙v0。将初始条件代入式（26）和（27）中可得振动系统响应解的系数

质点的动力学响应即为质点的平面运动轨迹，由式（26），（27）以及（28）可知，质点的运动轨迹取决于振动初始条件、振动刚度、刚体转速、质点初始位置r0。由于质点振动的同时作平面刚体转动，所以质点初始位置可视为离心振动系统的初始偏心位置。

为考察质点的运动轨迹，进一步令振动的初始条件t＝0．0s时，u（0）＝0．1m，v（0）＝0．05m，˙u（0）＝˙v（0）＝0．0m／s，初始偏心位置r0＝0．4m。当刚体转速分别为ω＝0．25r／s和ω＝1．0r／s时，质点的位移时程曲线见图4和5。根据质点随时间变化可以得到质点在转动坐标系下的平面运动轨迹，图6为刚体转速ω＝0．25r／s和ω＝1．0r／s时的质点运动轨迹，对给定的初始条件和振动刚度，质点的运动轨迹和振幅强烈依赖于刚体转速。图7表示质点振幅与转速的关系曲线，在一定的振动刚度和振动初始条件下，存在一个临界转速（在本例中为1．0r／s），超过临界转速后，质点振动幅值会迅速增加，超越线性振动的范围。前面关于系统动频的分析提到过刚体的最大转速，但对质点动力响应的分析表明，为保持系统的线性振动，刚体的转动应该小于刚体临界转速。

观察式（26）和（27）右端可以发现，质点的运动实际上由非偏心运动和偏心运动两部分构成，非偏心运动取决于振动初始条件、刚体转速和振动刚度，而偏心运动依赖于初始偏心位置、刚体转速和振动刚度。

图4 刚体转速为0．25r／s时的位移时程曲线Fig．4 Mass displacement history atω＝0．25r／s

图5 刚体转速为1．0r／s时的位移时程曲线Fig．5 Mass displacement history atω＝1．0r／s

图6 转速为0．25r／s和1．0r／s时质点的运动轨迹Fig．6 Mass trajectories atω＝0．25r／s andω＝1．0r／s

图7 质点振幅随刚体转速的变化Fig．7 Mass displacement amplitude vary with the rotational velocityω

图8和9分别为刚体转速为0．2 5r／s和1．0 r／s时分解的非偏心运动和偏心运动轨迹，注意到非偏心运动的幅值在不同刚体转速下是相同的，而偏心运动的幅值强烈依赖于刚体转速。当刚体转速为0．25r／s时，质点的运动轨迹主要是非偏心运动轨迹，偏心运动的影响很小。当刚体转速为1．0r／s时，偏心运动的影响显著影响质点总的运动轨迹。当初始偏心位置和振动刚度一定，偏心运动幅值随刚体转速增加而增大，考虑到非偏心运动和偏心运动的相位因素，这解释了图6和7的质点振幅随刚体转速先降后升的变化，特别是刚体转速超过某一临界值时，偏心运动逐渐主导了质点的运动轨迹。

图8 转速为0．25r／s时的非偏心运动和偏心运动轨迹Fig．8 Split non－eccentric and eccentric motional trajectories of the mass atω＝0．25r／s under rotating

图9 转速为1．0r／s时的非偏心运动和偏心运动轨迹Fig．9 Split non－eccentric and eccentric motional trajectories of the mass atω＝1．0r／s under rotating

4 结 论

本文借助双自由度离心振动系统的简单模型，讨论了大范围刚体转动和双自由度振动的耦合动力学问题。在已知刚体转动中，质点受到弹簧回复力、相对惯性力、初始离心力、离心力增量和科氏力。随着刚体转速的增加，第1阶动频线性降低而第2阶动频线性增加，刚体的最大转速就是振动系统的固有频率。双自由度离心振动系统的两阶模态都是极化的圆周运动，圆周运动的角速度对应于系统的两阶圆频率。

离心振动系统的质点运动轨迹与振动初始条件、振动刚度、刚体转速和初始偏心距离有关。质点运动轨迹由非偏心运动和偏心运动两部分构成，非偏心运动取决于振动初始条件、振动刚度和刚体转速，偏心运动取决于初始偏心距离、振动刚度和刚体转速，当初始偏心距离和振动刚度一定，偏心运动的振动幅值随刚体转速增加而增加。在有初始偏心情况下，离心振动系统存在临界刚体转速，超过临界转速后质点振动幅值迅速升高。对于大范围刚体运动下的线弹性振动系统，为保持线弹性范围内的振动，应考虑临界的刚体转速。

［1］ Park J H，Park H Y，Jeong S Y，et al．Linear vibration analysis of rotating wind－turbine blade［J］．Current Applied Physics，2010，10（2）：332—334．

［2］ Shum W S，Zhang L．Dynamic motion of whirling rods with Coriolis effect［J］．Applied Mathematical Modelling，2010，34：1 203—1 216．

［3］ Huang C L，Lin W Y，Hsiao K M．Free vibration analysis of rotating Euler beams at high angular velocity［J］．Computers and Structures，2010，88（17－18）：991—1 001．

［4］ Turhan O，Bulut G．On nonlinear vibrations of a rotating beam［J］．J．Sound and Vibration，2009，322（1－2）：314—335．

［5］ WEBER H I．A note on the stability of a rod subjected to compression by centrifugal force［J］．J．Sound and Vibration，1976，46（1）：105—111．

［6］ Lin S C，Hsiao K M．Vibration analysis of a rotating Timoshenko beam［J］．J．Sound and Vibration，2001，240（2）：303—322．

［7］ Cai G P，Hong J Z，Yang S X．Dynamic analysis of a flexible hub－beam system with tip mass［J］．Mechanics Research Communications，2005，32（2）：173—190．

［8］ 蒋丽忠，洪嘉振．大范围运动对弹性梁振动频率及模态的影响［J］．振动与冲击，1999，18（1）：12—16．

Jiang Lizhong，Hong Jiazhen．The frequency and mode of elastic beams in large overall motions［J］．Journal of Vibration and Shock，1999，18（1）：12—16．

［9］ 刘锦阳，洪嘉振．大范围运动空间梁的耦合动力学模型［J］．上海交通大学学报，2003，37（4）：532—534．

Liu Jinyang，Hong Jiazhen．Coupling dynamic modeling of the spatial beam under going large overall motions［J］．Journal of Shanghai Jiao Tong University，2003，37（4）：532—534．

［10］刘锦阳，李彬，洪嘉振．作大范围空间运动柔性梁的刚－柔耦合动力学［J］．力学学报，2006，38（2）：276—282．

Liu Jinyang，Li Bin，Hong Jiazhen．Rigid－flexible coupling dynamics of a flexible beam with three－dimensional large overall motion［J］．Chinese Journal of Theoretical and Applied Mechanics，2006，38（2）：276—282．

［11］章定国，余纪邦．作大范围运动的柔性梁的动力学分析［J］．振动工程学报，2006，19（4）：475—480．

Zhang Dingguo，Yu Jibang．Dynamical analysis of a flexible cantilever beam with large overall motions［J］．Journal of Vibration Engineering，2006，19（4）：475—480．

［12］Lin J L，Soedel W．On general in－plane vibrations of rotating thick and thin rings［J］．J．Sound and Vibration，1988，122（3）：547—570．

［13］Kim W，Chung J．Free non－linear vibration of a rotating thin ring with the in－plane and out－of－plane motions［J］．J．Sound and Vibration 2002，258（1）：167—178．

［14］Eley R，Fox CHJ，McWilliam S．Coriolis coupling effects on the vibration of rotating rings［J］．J．Sound and Vibration，2000，238（3）：459—480．

［15］Chen J S．On the linearization of the equations of motion of a rotating disk［J］．Applied Mathematical Modelling，2011，35（1）：392—397．

［16］Baddour N，Zu J W．A revisit of spinning disk models．Part I：Derivation of equations of motion［J］．Applied Mathematical Modelling，2001，25（7）：541—559．

［17］Mignolet M P，Eick C D，Harish M V．Free vibration of flexible rotating disk［J］．J．Sound and Vibration，1996，196（5）：537—577．

［18］Hashemi S H，Farhadi S，Carra S．Free vibration analysis of rotating thick plates［J］．J．of Sound and Vibration，2009，323（1－2）：366—384．

［19］Maretic R．Transverse vibration and stability of an eccentric rotating circular plate［J］．J．of Sound and Vibration 2005，280（3－5）：467—478．

［20］刘锦阳，洪嘉振．作大范围运动矩形薄板的建模理论和有限元离散方法［J］．振动工程学报，2003，16（2）：175—179．

Liu Jinyang，Hong Jiazhen．Dynamic modeling theory and finite element method for a rectangular plate un－dergoing large overall motion［J］．Journal of Vibration Engineering，2003，16（2）：175—179．

［21］赵飞云，洪嘉振．作大范围运动矩形板的动力学建模理论研究［J］．计算力学学报，2008，25（6）：868—873．

Zhao Feiyun，Hong Jiazhen．Study on dynamic modeling of rectangular plates undergoing large overall motion［J］．Chinese Journal of Compu－tational Mechanics，2008，25（6）：868—873．

［22］Wang Y Q，Guo X H，Chang H H，et al．Nonlinear dynamic response of rotating circular cylindrical shells with precession of vibrating shape－Part I：Numerical solution［J］．Inter．J．Mechanical Sciences，2010，52（9）：1 217—1 224．

［23］Liew K M，Ng T Y，Zhao X，et al．Harmonic reproducing kernel particle method for free vibration analysis of rotating cylindrical shells［J］．Computer Methods in Applied Mechanics and Engineering，2002，191（37－38）：4 141—4 157．

［24］黄克智，薛明德，陆明万．张量分析［M］．北京：清华大学出版社，2003．

Huang Kezhi，Xue Mingde，Lu Mingwan．Tensor A－nalysis［M］．Beijing：Tsinghua University Press，2003．

［25］Palais B，Palais R．Euler’s fixed point theorem：the axis of a rotation［J］．J．Fixed Point Theory and Applications，2007，2：215—220．

［26］蔡燧林，盛骤．常微分方程组与稳定性理论［M］．北京：高等教育出版社，1988．

Cai Suilin，Cheng Zhou．Ordinary Differential Equations and Stability Theory［M］．Beijing：Higher Education Press，1988．

［27］马知恩，王绵森．工科数学分析基础（下册）［M］．北京：高等教育出版社，1999．

机电耦合动力学 篇5

伺服系统用来控制被控对象的某种状态, 属于自动控制系统的一种。机电相结合的伺服系统主要是由机械传动部分和控制部分组成的, 除此之外往往还有润滑、液压部分, 系统通过各部分的耦合实现控制。对于传动系统, 若在非平稳的状态下发生机电耦合振动, 将对系统的安全运行造成威胁。因此, 分析伺服系统中精密传动系统的耦合机理就显得十分重要, 这也对工业生产有重大意义。本文介绍了永磁交流伺服精密驱动系统机电耦合及其相关模型的建立, 同时结合具体实例, 分析了该系统在实际运作中的特点以及一些对应参数的变化规律, 以便理解机电耦合模型, 同时也给我们的建模提供一些建设性的意见。

1 永磁交流伺服精密驱动系统机电耦合概述

作为现代交流伺服系统的主流, 永磁同步电动机具有结构简单、功率高、易于散热、转动惯量低等许多优点。精密传动是一种非常重要的基础性零件, 也是双向高精度传递运动机械传动形式的总称。这类机械传动形式是做高精度机械传递运动, 具有极高的承载能力和运动精度, 并且具有刚度高、体积小、响应快等特点。精密传动包括很多种, 例如常见的齿轮传动和螺旋传动, 以及不常见的少齿差行星传动、谐波传动和轮系传动等。同样, 永磁交流伺服驱动系统也包含很多种子系统, 如驱动系统、传动系统、冷却系统以及负载系统等, 这些子系统相结合才组成了完整的主系统。各子系统间存在大量的物理过程和多参量耦合关系, 这些耦合关系以及关系的输入和输出决定了相应的永磁交流伺服精密驱动系统的动态性能, 而各个子系统的动态性能又受到其他子系统动态性能的影响。因此, 在分析永磁交流伺服精密驱动系统的动态性能时要先建立局部的耦合模型, 通过对局部耦合模型、耦合方式及规律的分析得出整体的规律, 才能更加快捷准确地解决问题, 找准原始的规律, 这对更好地分析永磁交流伺服精密驱动系统的运作规律及方式有很大帮助[1]。

2 伺服系统中耦合的建模分析

高精度伺服系统具有精密控制系统运行、保持过程高度稳定性的特点, 并且对制造过程有精密的自动调节能力。在对伺服系统进行耦合分析时要采取耦合事实提取、耦合问题建模、耦合参数模型解耦等步骤, 下面我们就通过建模的方法对伺服系统中的局部耦合和整体耦合作相应分析[2]。

2.1 局部耦合建模分析

要想对局部耦合问题进行建模分析, 就得先把问题从整体中分离出来, 利用现有的数学、物理公式建立模型, 并作相应分析。通常, 一个伺服系统都是由许多基本元件组成的, 如电感、电阻、弹簧、电容等, 可以对应地写出其各自的特性规律方程, 即:

(1) 电感的物理特性方程:u=Ldi/dt, 即流经电感L的电流i和电感端电压的关系;

(2) 电阻的物理特性方程:u=Ri, 即电阻两端电压与流经电阻的电流i的关系;

(3) 弹簧的物理特性方程:F=K∫vdt+F0, 即在外力F作用下的刚度为K的弹簧与其速度v的关系;

(4) 电容的物理特性方程:i=Cdu/dt, 即电容在充电或放电的条件下与电流的关系。

机电回路是由机电元件连接组成的, 因而机电系统中的输入和输出关系满足基尔霍夫电流和电压定律, 即在任意一个闭合回路中有u=0或经过任意一个节点的i=0。同时, 在机电系统中仍然遵循达朗贝尔空间连接定律, 即在机械网络中, 绕任一回路的位移和速度的和均为0。通过对每一个元件的特性方程的列举, 并结合对应满足的定律就可以建立起相应的数学模型, 而通过对建立的模型进行试验和分析就可以解决部分耦合问题。

2.2 全局耦合建模分析

通过对耦合参数基本规律的了解, 可以对全局耦合结构进行相应的分析, 从而对系统的全局耦合进行建模。伺服系统是由很多个子系统组成的, 其通过子系统之间的耦合来实现整体的联系, 因此必须根据各部分之间的关系, 并结合对应的方程才能建立起关于整体耦合的模型。具体要解决以下几个问题: (1) 分析耦合参数与主运动的关系及全局耦合模型如何解耦; (2) 分析奇异点的条件和动态规律, 进行解耦单元设计; (3) 对部分进行分析解离, 然后找到各部分之间的耦合方式, 从而实现数学模型的建立。

3 机电耦合模型的算例分析

永磁精密传动装置的机电耦合动力学方程组是多变量及非线性的, 因此在求解时不宜用解析法, 最好用数值计算法来对其动力学规律进行研究。假设存在一个传动系统, 相关的参数如下:定值电阻5.6Ω, dq轴电感11.57 mH, 转子永磁体磁通量0.125 Wb, 极对数为4, 转子转动惯量为0.384×10-4 kg·m2, 额定转矩为1.47N·m, 额定转速4 000r/min, 电机粘滞摩擦系数为0, 减速器转动惯量为1.25×10-6 kg·m2, 减速器粘滞摩擦系数为0.001kg·m2/s。系统空载启动时电机转速设定为800r/min, 空载三相定子电流波形在0.015s后为一条平滑的水平线, 在开始到0.015s之间为波动的波形, 且平稳后的电流显示为0。对应的空载转速波形, 同样是在0.015s之前波动极大, 0.015s后稳定在大约800r/min。同样, 空载的转矩波形也是在0.015s前后呈现不同的变化规律, 随后稳定在0值。为了进一步分析, 在系统运作的0.06s时突加负载转矩2N·m, 且初始转速不变, 此时波形图出现明显变化, 并且三相定子电流在0.06s后出现极大的波动, 波动区间为-5~5A。三相定子电流在持续0.1s后仍是一种波动的形态, 此时, 对应的转速波形在0.06~0.065s之间波动, 在0.065s后稳定在与前一阶段相同的转速值上。此时, 转矩波形则不同, 即在0.06~0.07s的整个区间内都出现波动, 而在0.07s后的稳定值明显大于前一阶段的稳定值。通过对初始转速相同, 系统启动后不加载和启动后加载这2种情况的分析可看出, 加载后系统的转速波形和转矩波形在0.06s都有变化, 但是最终保持稳定。这种现象正好符合PMSM特性, 说明对应的PMSM机械特性比较显著。同时, 上述分析也成功模仿了系统空载和突加负载的情况, 也正好验证了所建立数学模型的正确性[3]。

4 结语

永磁交流伺服精密驱动系统是一种机电耦合系统, 它的主要作用是实现机电能量的转换, 因此在分析其耦合特性时, 必须从复杂机电系统的角度对其进行局部耦合和全局耦合的分析, 从而达到建立相应的数学模型的目的, 以方便分析和实验。只有对该系统的每一部分都进行仔细的分析, 并结合适用的方程及原理, 才能建立耦合模型, 从而通过对模型的研究实验来解决机电耦合中的问题。同时, 这对现代工业化生产也是很有利的, 不但可以提升工作效率, 还可以降低机电耦合问题的发生率, 减少机械生产的损失。如今, 永磁交流伺服精密驱动系统在我们的工业化生产中运用较多, 其具有运作完美流畅的明显优点, 不过仍需要我们进行仔细的分析, 从根本上了解其运作规律和方法, 以便更好地运用。

参考文献

[1]温熙森, 邱静, 陶俊勇.机电系统分析动力学及其应用[M].北京:科学出版社, 2003

[2]周超群, 陈小安, 合烨.伺服系统中精密传动装置精度分析[J].现代制造工程, 2007 (4)

单轴微加速度器机电耦合分析 篇6

1 惯性系统与机电耦合模型的建立

加速度测量系统是由感应单元、弹簧结构、阻尼系统的微结构组成。当加速度器基底感应到加速度作用时,感应结构就会受到与相反方向的惯性力作用,导致感应单元发生位移[3,4]。文中通过建立动态系统模型(系统模型如图1),便于计算在静电力作用下的结构位移。

根据动态系统原理,加速度器系统阻尼力和刚度分别与感应质量块的速度和位移是线性关系。可得如下系统方程:

m是感应单元质量,c是黏性阻尼系数,k是弹簧刚度,F是加速度器系统所受外力,Fa是惯性力,Fθ是静电力[5]。

加速度器由一个感应单元和一组弹簧结构组成[6]。加速度器之所以可以测量加速度作用及感应结构位移,是因为加速度器系统结构中有电压被施加在梳齿(极板)上,这样可以将输入的机械能转化电能输出。先将加速度作用转变为位移变化转,然后进一步转变为电容变化。

梳齿(电极板)是带电的,所以静电力是库仑力作用的。静电作用力有效地驱动的平行板电容器结构。

微观结构之间的静电势可以由势能公式计算静电场存储装置。梳子(电极)的存储的能量可以表示为:

2 加速度器结构设计与计算

一般加速度器机电耦合分析较少考虑静电力对微加速度器的影响。但当系统中测量电压被施加在梳齿(电极板)间,静电力作用对加速度器极板间位移有较大影响,经理论推导发现两者之间近似线性关系。因此,当感应质量块发生位移时,可以将静电力作用近似为静电弹簧Kθ作用。可以将公式(1)转化为:

其中,k’=k-kθ,Fa=a·m

加速度器加速度量变为电信号测量的最大难点在于微弱信号提取。电路寄生电容包括测量电路的寄生电容110f F和路由寄生电容10f F,加速度输入引起系统结构电容变化比系统固有寄生电容小得多[7,8]。此外,将临近梳齿当作电容器的电极板,作为一个平行板电容器的感应电容器,电气信号的灵敏度可以表达为

其中,Vm调制电压为2V,d0梳齿间距为1.5μm,经计算灵敏度约为2.3m V/g。

本文利用COVENTORWARE软件对感应质量块进行动态分析,为惯性系统施加1g加速度作用和2V检测电压。根据静电力和动态系统理论分析,计算出质量块位移是3.282nm,与计算机仿真结果仅有1.86%差异。

若感应质量块沿感应(X轴)方向运动时,固定梳齿与活动梳齿间距发生变化,差分电容也随之变化。当有1g加速度作用时,感应质量位移约为3.343nm,产生差分电容0.1557f F与软件模拟结果电容变化结果相差1.29%。加速度作用与结构位移间的关系见表1。

以上所提出的加速度计的设计,当加速度被施加在Y、Z轴的方向,电路输出的输出电压为零,通过仿真软件,其结果从电容变化上分析。利用此差分电路输出,既能提高性能,也能抑制交叉轴(非测量方向)敏感性。

3 结论

本文对单轴MEMS梳齿加速度器进行动态分析研究,对比软件仿真结果,证明了在低g值加速度作用下,该系统中感应质量块位移与差分电容输出是线性的。同时,该系统设计在建模与结构方面还可以通过动态分析进行进一步优化,为后续微弱信号提取电路设计提供参考。

参考文献

[1]樊建勋.MEMS传感器及其在航空领域的应用[J].电子技术与软件工程,2014,04.

[2]李瑞养,林强.基于MEMS技术的惯性传感器的发展及在组合导航系统的应用研究[J].中国科技全景,2010,7.

[3]刘恒,等.一种硅微机械结构振动幅度的电学测量方法及实验研究[J].传感器技术学报,2013,12.

[4]谌贵辉,等.地震勘探用三维MEMS加速度传感器的研究[J].传感器与微系统,2012,04.

[5]高奂文,董慧颖.一种振梁电容式硅微加速度计的设计与分析[J].沈阳工业学院学报,2004,4.

[6]张莉,等.基于微电子系统(MEMS)加速度地震勘探三分量数字检波器简介[J].中国煤田地质,2006,02.

[7]赵楷,等.基于Verilog-A的电容式MEMS加速度计模型[J].传感器技术学报,2004,4.

机电耦合动力学 篇7

双Stewart平台推力向量控制系统是自主创新成果, 将应用于我国新一代先进常规战术导弹II级固体火箭发动机轻质化、小力矩、高精度的先进推力向量控制系统。利用直线电机驱动的双Stewart平台, 通过调节喷管可动扩张段的直径大小和方位实现推力向量系统控制, 克服了目前国内火箭全周摆动推力向量系统存在的驱动力矩大、运动部件接触应力、变形大、摆心漂移等缺陷, 为机械类固体火箭发动机推力向量控制系统未来型号研制奠定理论基础[2]。双Stewart平台推力向量控制系统具有同时调节喷管喉径和推力矢量的功能, 保证了在调节推力大小和方向的过程中载荷集中于喷管扩张段的固定体, 少量载荷分布于由直线电机驱动的喷管扩张片, 使机构的运动件从高载荷、高温高压状态转变为相对的低载荷、低温区域。由于喷管离火箭重心的距离较远, 一个很小的偏转角, 便会产生足够大的偏转力矩, 一般偏转角只需要6°~8°, 此时, 对轴向推力的影响只有10%左右。因此, 即使是推力很大的火箭发动机, 采用矢量喷管方案所需驱动力也非常有限, 这为利用直线电机进行直接驱动提供了条件。减小系统驱动力就意味着减轻系统重量 (发动机每减少1kg, 可以增加射程20~30km) 。直线电机驱动与传统液压驱动相比, 在减少火箭重量的同时提高了喷管的控制精度[2,11]。直线电机驱动的双Stewart平台推力向量系统组成结构如图1所示。该装置的基本结构由喷管扩散段 (固定体) 、A2转向调节环、A1调节环、收敛调节片、扩张调节片、收敛密封片、扩张密封片、十字转接头、拉杆等组成。

从机构学的角度讲, 该装置为双Stewart平台的复杂空间机构。并联于双Stewart平台之间的十几组空间铰链机构, 导引收敛调节片和扩张调节片围成的空间成为收扩式喷管, 收敛调节片的位姿由A1调节环确定:扩张调节片的位姿由A2转向调节环和A1调节环联合控制。A1、A2调节环由直线电机驱动, 以实现喷管摆动运动。该装置由前部、中部和后部三个组成部分。前部为与喷管扩散段相连的部分, 用于控制A1的面积。中部为扩张调节片驱动机构, 后部为扩张调节片部分。收敛调节片、扩张调节片和收敛密封片共同组成喷管, 承受火箭燃气射流的作用。在推力矢量调节过程中, 调节片和密封片始终受到燃气的压力作用, 该压力随喷管的变化而发生变化。在工作过程中, 在直线电机驱动力和燃气射流共同作用下, 如何保证收敛 (扩张) 调节片协调运动, 可操控地完成推力矢量调节, 是该新型喷管研制的核心问题。

2 Stewart平台的机构模型

2.1 Stewart平台的姿势

并联机构的Stewart平台一般由负载和基台上下平台组成, 而且其由铰副连接6个移动副构成。运动分析是要建立其末端构件的位姿与输入关节变量之间的关系。Stewart平台的运动学设计有两个最基本的问题:正解和反解。在本文中, 对于固体火箭发动机喷管的双Stewart平台来说, 其运动学反解就是给定动平台的位置参数即3个线性坐标参量 (X, Y, Z) 和姿态即3个旋转坐标参量 (α, β, γ) , 确定Stewart平台6个分支杆的长度[15]。对于串联机构而言, 求解运动学正解比较容易, 而对于并联机构来说, 求解运动学位置反解比较容易。并联机构的位置反解是位移分析的一部分, 研究高速、方便的位置反解方法对并联机构的控制具有实际的应用意义[15]。

描述动平台在固定坐标系O-XYZ的姿态, 采用的是欧拉角描述法。动平台的姿态可通过3次连续的绕坐标轴旋转来实现, 在这里采用绕Z-Y-X轴旋转角度的变换矩阵方式, 即:首先将动平台绕自身动坐标系O1-X1Y1Z1的Z1轴旋转α角, 再绕当前动坐标系的Y1轴旋转β角, 最后绕当前动坐标系的X1轴旋转γ角[15], 三个基本旋转矩阵为:

式中的“s”、“c”分别表示“sin”、“cos”, 以下各章中均采用此种标记方法。则动坐标系与静坐标系之间对应的旋转矩阵为3个基本旋转矩阵之积。即:

2.2 Stewart平台的运动规则

并联机构的位置反解是给定动平台的位置和姿态, 求能够实现该位置和姿态的驱动关节的位移。对于本文分析的喷管的Stewart平台是在理想的运动状态下。因此, 在实际的工程应用中, 固体火箭发动机喷管的Stewart平台的运动轨迹可分解为直线轨迹和圆弧轨迹。

3 喷管的热分析

在喷管的动力学分析基础上, 对固体火箭发动机喷管的耦合场进行分析。固体火箭发动机一般通过燃烧三氧化二铝产生燃气射流, 提供可控轴向和侧向飞行推力, 以保证飞行器稳定性, 按规定弹道飞行或作必要机动飞行, 故首先对喷管进行热分析。火箭发动机喷管热分析的相关参数导热系数为1.2, 单元类型定义Thermal Solid的Tet 10node 87。由于模型过于复杂, 从Solid Works2012中导出X_T格式, 再导入ANSYS中, 模型只保存与流体接触的1/4。采用自由划分, 勾选上smart size。定义载荷:第一部分:喷射管喉管部分, 载荷大小为927K, 温度最高;第二部分:收敛调节片和收敛密封片部分, 载荷为800K;第三部分:十字转接头, 载荷大小为700K;第四部分:扩散调节片和扩散密封片部分, 载荷大小为500K。进行运算, 计算结果。

通过上面的热分析可知, 辐射换热对内壁面温度分布的影响较大。内壁表面温度开始很低, 当然由此可知它的烧蚀率也不高。后来随着化学反应的加剧, 温度逐渐升高, 而对流换热系数随着内壁温度的升高却降低。在开始时, 在进口处以及收敛段, 由于燃气温度较高且受外流低温气体影响较小, 考虑辐射换热计算所得内壁面温度明显高于未考虑辐射时计算所得温度;在扩张段, 由于燃气温度下降以及外流低温气体介质的影响, 考虑辐射换热与否对内壁面温度数值大小的影响逐步减弱;在接近出口段, 外流场低温气体介质对喷管内壁面热辐射的吸收起到了主导作用, 以致考虑辐射换热所得壁面温度低于未考虑辐射时所得温度[12]。而且温度场从喉管到扩张片的温度分布是递减的, 实现了高温气体向低温气体的转化[1,3]。

4 火箭发动机矢量喷管的热-结构耦合分析

本文的固体火箭发动机喷管热结构计算的方法是由一维等熵流和巴兹经验公式得到喷管内壁面的表面传热系数, 利用ANSYS软件的热结构耦合分析模块将表面传热系数作为喷管的热边界条件, 进行热结构计算[4,5]。

(1) 火箭发动机矢量喷管的热-结构耦合分析的参数:结构分析的材料为30Cr Mn Si A钢, 弹性模量E=210GPa, μ=0.3, ρ=7800kg/m3, 导热系数k=43.2, 热膨胀系数为1.62×10-5。单元类型为SOLID 187这些参数在火箭发动机矢量喷管的热-结构耦合分析时需要设定。

(2) 矢量喷管的热-结构耦合分析的步骤

(1) 清除之前热分析的物理环境。 (2) 将温度单元转换为结构单元。这一步非常重要, 这是实现热-结构耦合的桥梁, 通过转换温度单元转换为结构单元, 这样才能定义结构单元的相关参数。 (3) 施加温度载荷。将之前热分析结果作为结构分析的面载荷, 这就是顺序耦合法, 将前一个分析的结果当作后一个分析载荷, 从而实现耦合。 (4) 定义固定端。 (5) 应力场求解。

从应力场分布我们可以看出, 最大应力为11.9GPa, 最小应力为1712Pa, 压力由进口到出口递减, 跟温度场变化趋势相同。从x向应力分布上看压力, 以x轴为中心向四周逐渐增加。作用在喉衬上的机械载荷推挤扩张段后, 移在扩张段和扩张段背衬, 所以喷管内壁面受热严重最高的应力区出现在喉衬和扩张段入口处。除此之外, 可知对于整个结构在金属支撑件存在倒角区域。

5 矢量喷管流体-热-结构耦合分析

流固热耦合问题是指在流体流动、固体变形和热量传递系统中三者之间的相互作用。固体火箭发动机喷管内流动的工质是高温燃气, 燃气热量产生的热应力与其流动所产生的流体压力导致固体变形, 固体变形导致流场变化, 进而使温度场发生变化。流、固、热三个物理场场间的耦合作用是同时发生、双向相互的。在耦合过程中, 流体问题涉及大量的非线性问题、固体问题涉及复杂的几何和弹塑性非线性问题以及复杂的温度场问题。因此, 整个流固热耦合问题是一个十分复杂的非线性问题[14]。

首先对本文涉及的流体力学、传热学和弹性固体力学基本方程作如下假设: (1) 针对二维问题; (2) 不考虑质量力; (3) 物性参数不随温度变化; (4) 流体为粘性不可压缩流体; (5) 弹性固体为各向同性体。

6 火箭发动机矢量喷管的流固热场耦合分析

6.1 流体-热-结构耦合分析的相关参数

流体单元的参数:流体的单元类型为FLOTRAN CFD中的3D FLOTRAN 142单元;流体物性为AI-SI;进口压力为1MPa, 落压比NPR=2~8。出口压力为0.075MPa。

温度单元的参数:温度单元类型为FLOTRAN CFD中的3D FLOTRAN 142单元;导热系数为1.2。

结构单元的参数:结构分析的材料为钢材为30Cr Mn Si A, 弹性模量E=210GPa, μ=0.3, ρ=7800kg/m3, 导热系数k=43.2, 热膨胀系数为1.62×10-5。单元类型为SOLID 187。

6.2 流体热耦合的相关步骤

选择FLOTRAN CFD模块, 设置分析选项: (1) 定义工作文件名和工作标题。 (2) 定义单元类型。单元类型为FLOTRAN CFD中的3D FLOTRAN 142。 (3) 导入流体分析区域的模型。 (4) 划分有限元网格。 (5) 施加边界条件。施加温度、流体速度、进出口压力载荷。 (6) 设置求解选项:将“Adiabatic or thermal?”项设置为thermal。 (7) 稳态控制设置:将载荷步设置为100, 其他选项默认。 (8) 流体物性设置:将流体设置为标准空气。 (9) 设置重力加速度。 (10) 选择CFD压力求解器:选择TDMA压力求解器, 压力取100 (默认) 。 (11) 求解。

流体热耦合结束后得到结果文件, 这个结果文件将会是下一个分析的载荷。

6.3 流体-热-结构耦合分析步骤

(1) 设置分析选项, 选择结构模块。 (2) 定义工作文件名和工作标题。 (3) 定义单元类型, 单元类型为10node 187单元。 (4) 导入固体火发动机轴对称矢量喷管的UG模型。 (5) 定义材料属性。结构分析的材料为钢材为30Cr Mn Si A, 弹性模量E=210GPa, μ=0.3, ρ=7800kg/m3, 导热系数k=43.2, 热膨胀系数为1.62×10-5。单元类型为SOLID 187。 (6) 划分有限元网格。 (7) 施加载荷:来自FLOTRAN分析的温度和压力载荷。 (8) 进行求解后处理。

机电耦合动力学 篇8

混合动力汽车 (HEV) 与传统汽车最大的区别在于它有多个动力转化装置, 能根据行驶工况的变化随时改变动力输出源, 以达到最佳的燃油经济性和排放性能。如何将多个动力源有效地耦合, 使HEV的各项功能得到充分发挥是一项关键技术[1,2]。目前国外HEV动力耦合装置的功能已经十分完备, 可以实现HEV所有的动力传递要求。相比之下, 国内关于动力耦合装置的研究较少, 当前HEV采用的动力耦合装置结构过于简单, 难以实现HEV的全部功能, 已经严重阻碍了我国HEV产业化的进程。

1 行星机构运动特性分析

本文研究的行星机构动力耦合装置由1组双行星齿轮机构、2组湿式多片离合器 (C1、C2) 和1组制动器 (B1) 构成, 如图1所示。行星机构有2个自由度, 通过离合器和制动器控制1个自由度, 最终实现整车动力的稳定传递。

1.1 动力耦合装置的机械连接

行星机构的太阳轮与发动机曲轴直接相连;行星架一端与电机转子相连, 另一端与离合器C1的主动部分相连;齿圈与制动器B1的可动部分相连, 制动器B1的固定部分与变速器壳体固连, 一旦制动器B1接合齿圈被制动, 转速即变为0;齿圈还与离合器C2的主动部分连接;变速器的输入轴与离合器C1、C2的从动部分相连, 只要离合器C1或C2接合, 动力耦合装置的动力就会传递到变速器输入轴上。行星机构动力耦合装置的机械连接情况如图1所示。

1.2 行星机构的动力特性

本文所研究的行星机构是双行星轮机构, 其转速特性和力矩传递特性如下。

1.2.1 转速特性

双行星轮行星机构各构件的转速满足下式:

ωR=ρ ωS+ (1-ρ) ωH (1)

ρ=zS/zR

式中, ωS为太阳轮的绝对转速;ωR为齿圈的绝对转速;ωH为行星架的绝对转速;ρ为行星机构的特征参数;zS为太阳轮齿数;zR为齿圈齿数。

为了更直观地表示行星机构3个构件的转速关系, 本文引入了虚拟杠杆作为分析工具, 如图1所示。行星机构的3个构件分别由3个纵轴表示, 纵轴上的坐标值表示3个构件的转速。无论行星机构的转速怎样变化, 代表3个构件转速的节点都在一条直线上, 这样, 已知行星机构任意2个构件的转速, 通过这两点作直线与第3坐标轴的交点即为第3构件的转速, 行星机构3个构件的转速点仿佛在一条“杠杆”上。

1.2.2 力矩传递特性

在忽略摩擦的情况下, 根据内力矩平衡和功率的平衡原理可知行星机构的3个构件的力矩满足下式:

$\frac{{\rm T}{}_{\text{S}}}{\rho }=\frac{{\rm T}{}_{\text{Η}}}{1-\rho }=-{\rm T}{}_{\text{R}}$ (2)

式中, TS为太阳轮的内力矩;TR为齿圈的内力矩;TH为行星架的内力矩。

利用虚拟杠杆进行力矩分析时, 作用在行星机构各构件的内力矩相当于作用在虚拟杠杆上的“力”, “力臂”就是3个纵轴间的距离 (由行星机构的结构参数ρ决定) , 如图1所示。已知任意1个构件的内力矩, 就可以通过杠杆平衡原理求出其余2个构件的内力矩[3]。

1.3 行星机构动力耦合装置的运动方程

为了更方便地分析HEV的各种工作模式, 将行星机构动力耦合装置与相关部件的运动和力矩传递关系总结成如下方程:

太阳轮运动方程

${\rm T}{}_{\text{e}}-{\rm T}{}_{\text{S}}=\dot{\omega }{}_{\text{e}}{\rm I}{}_{\text{e}}$ (3)

行星架运动方程

${\rm T}{}_{\text{m}}-{\rm T}{}_{\text{Η}}-{\rm T}{}_{\text{C}1}=\dot{\omega }{}_{\text{m}}{\rm I}{}_{\text{m}}$ (4)

齿圈运动方程

$-{\rm T}{}_{\text{R}}-{\rm T}{}_{\text{C}2}-{\rm T}{}_{\text{B}1}=\dot{\omega }{}_{\text{R}}{\rm I}{}_{\text{R}}$ (5)

变速器运动方程

${\rm T}{}_{\text{C}1}+{\rm T}{}_{\text{C}2}-{\rm T}{}_{\text{V}}=\dot{\omega }{}_{\text{V}}{\rm I}{}_{\text{V}}$ (6)

角加速度约束方程

$\dot{\omega }{}_{\text{R}}=\rho \dot{\omega }{}_{\text{e}}+(1-\rho )\dot{\omega }{}_{\text{m}}$ (7)

式中, Te为发动机输出转矩;Tm为电机输出转矩;TV为变速器输入轴上的等效阻力矩;TC1为离合器C1传递的转矩;TC2为离合器C2传递的转矩;TB1为制动器B1的制动力矩;$\dot{\omega }{}_{\text{R}}$为齿圈角加速度;$\dot{\omega }{}_{\text{e}}$为发动机 (太阳轮) 角加速度;$\dot{\omega }{}_{\text{m}}$为电机 (行星架) 角加速度;$\dot{\omega }{}_{\text{V}}$为变速器输入轴角加速度;Ie为发动机曲轴+太阳轮转动惯量;Im为电机转子+行星架的转动惯量;IR为齿圈的转动惯量;IV为变速器输入轴等效转动惯量。

上述方程与式 (2) 构成了行星机构动力耦合装置的运动方程。当HEV处于不同工作模式时, 方程将会蜕化成各种简洁形式。

2 HEV的各种工作模式

行星机构动力耦合装置能实现HEV的10种工作模式, 这里仅以最具代表性的电力变矩模式为例说明动力耦合装置的工作状态。

2.1 电力变矩模式分析

当车速较低且驱动力矩要求较高, 电机驱动力矩达不到要求 (如低速爬坡工况) 时, 或荷电状态SOC较低不能实现电机起步时, HEV进入电力变矩模式[4]。该模式是行星机构动力耦合装置最具特色的一种工作模式。

汽车处于起步阶段, 发动机转速一般在800r/min以上, 而此时车速为0, 变速器输入轴转速也要由0逐渐提升。传统汽车通过液力变矩器 (或其他形式的起步离合器) 来消除这一阶段的转速差, 但是液力变矩器在起步阶段要产生大量的热以耗散发动机的盈余功率, 这对于降低整车油耗和起步离合器的磨损很不利。电力变矩模式通过行星机构动力耦合装置使电机反转发电, 吸收发动机起步过程中的盈余功率, 这样既减少了整车能耗, 又有效地避免了产生热量带来的不利影响。

电力变矩模式要求动力耦合装置的离合器C2接合, 离合器C1和制动器B1分离, 如图2所示。动力耦合装置将发动机的动力分成两部分:一部分拖动电机反转发电, 另一部分驱动车辆行驶。发动机对太阳轮输出动力, 太阳轮带动行星架和齿圈转动。行星架带动电机转子反转发电, 并调节齿圈的转速, 齿圈以合适的转速和转矩通过离合器C2带动变速器输入轴驱动车辆以低转速、大扭矩行驶。

整个电力变矩过程中, 发动机的转速稳定在最低经济转速, 开始时齿圈与变速器输入轴的转速为0, 对应车辆静止状态;电力变矩模式结束时行星架转速变为0, 齿圈与变速器输入轴的转速得到提升, 这时汽车已达到一定的行驶速度, 如图2所示。

电力变矩模式的约束条件为${\rm T}{}_{\text{C}1}=0‚{\rm T}{}_{\text{B}1}=0‚\dot{\omega }{}_{\text{R}}=\dot{\omega }{}_{\text{V}}$, 把这些约束条件代入行星机构动力耦合装置的运动方程得蜕化方程:

${\rm T}{}_{\text{e}}-{\rm T}{}_{\text{S}}=\dot{\omega }{}_{\text{e}}{\rm I}{}_{\text{e}}$

${\rm T}{}_{\text{m}}-{\rm T}{}_{\text{Η}}=\dot{\omega }{}_{\text{m}}{\rm I}{}_{\text{m}}$

$-{\rm T}{}_{\text{R}}-{\rm T}{}_{\text{C}2}=\dot{\omega }{}_{\text{V}}{\rm I}{}_{\text{R}}$

${\rm T}{}_{\text{C}2}-{\rm T}{}_{\text{V}}=\dot{\omega }{}_{\text{V}}{\rm I}{}_{\text{V}}$

$\dot{\omega }{}_{\text{R}}=\rho \dot{\omega }{}_{\text{e}}+(1-\rho )\dot{\omega }{}_{\text{m}}$

$\frac{{\rm T}{}_{\text{S}}}{\rho }=\frac{{\rm T}{}_{\text{Η}}}{1-\rho }=-{\rm T}{}_{\text{R}}$

变速器输入轴的力矩来自离合器C2, 由上述方程得TC2=Te/ρ (ρ<1) , 即驱动力矩相当于把发动机转矩放大了1/ρ倍, 电力变矩的名称也由此而来[5]。

2.2 HEV的10种工作模式分析

用同样的方法可以得到另外9种工作模式各构件的运动状态。行星机构动力耦合装置通过离合器C1、C2和制动器B1的接合/分离控制实现动力传递的不同路径;动力源的各种工作状态实现了能量的不同流向, 两者配合共同实现了HEV的10种工作模式[6], 如表1所示。

注:○分离, ●接合。

3 动力耦合装置的模式切换控制

由表1可知, 当HEV处于某种驱动模式时, 行星机构动力耦合装置的离合器和制动器首先进入相应的接合状态, 以实现该模式动力传递路径的畅通, 动力源随即进入相应的工作状态, 最终整车实现了与行驶工况相适应的工作模式。对于行星机构动力耦合装置来说, 模式切换控制就是对离合器C1、C2和制动器B1进行控制。

3.1 动力耦合装置的控制要求

由表1可知:要求离合器C1接合的模式有7种, 当电力变矩模式向发动机驱动模式切换时需要控制C1的接合速度, 以满足整车舒适性的要求;要求离合器C2接合的模式有4种, 当电机驱动模式向发动机驱动模式切换时需要控制C2的接合速度;要求制动器B1接合的模式有2种, 当发动机倒车模式时需要控制B1的接合速度。

HEV工作模式切换时, 既要满足动力性要求又要满足舒适性要求, 行星机构动力耦合装置依靠离合器接合控制满足这两方面要求。

3.2 离合器的接合控制

湿式多片离合器采用液压远程控制, 接合过程中离合器传递的力矩与控制压力成正比, 文献[7]详细阐述了湿式多片离合器的控制压力与其传递力矩的关系。

离合器的接合过程如图3所示, 大致可以分为三个阶段[8]:第一阶段, 根据驾驶员加速踏板开度α和开度变化率$\dot{\alpha }‚$确定一个合适的初始压力ps;第二阶段, 根据驾驶员加速踏板开度变化率$\dot{\alpha }$和离合器主从动片的转速差Δω, 确定一个合适的压力变化率$\dot{p}{}^{[9]}$;第三阶段, 离合器主从动片的转速差小于一定的门限值时, 离合器控制压力迅速上升至最大压力, 实现离合器的锁止。

3.3 离合器的模糊控制策略

为了缩短离合器的接合时间, 同时兼顾驾驶员的操纵要求, 需要控制离合器接合第一、第二阶段的压力。本文采用模糊控制算法, 以驾驶员的加速踏板开度和离合器的主从动片转速等作为输入量, 计算这两个阶段的离合器控制压力。图4为离合器C1在第一阶段的初始压力模糊规则MAP图, 图5为离合器C1在第二阶段的压力变化率模糊规则MAP图。

当HEV模式切换时, 由加速踏板信号和主从片转速信号根据图4、图5得初始压力ps和压力变化率$\dot{p}‚$离合器C1的控制压力随时间变化关系为

$p{}_{\text{c}\text{l}}({\rm T})=p{}_{\text{s}}+{\displaystyle \int }{}_0^{{\rm T}}t\dot{p}(t)\text{d}t$ (8)

式中, pcl为离合器的控制压力;T为离合器接合过程的时间历程。

离合器C2和制动器B1的控制压力的计算方法与离合器C1相似。

4 行星机构动力耦合装置的试验研究

在现有车型CFA6470E的基础上搭建了混合动力试验样车, 样车上搭载了行星机构动力耦合装置并将其模糊控制策略下载到整车控制器中, 完成了10种工作模式的测试, 图6所示为电力变矩模式向发动机驱动模式切换的试验结果。

试验中, 离合器的接合状态用电磁阀的占空比 (PWM) 表示, 对于离合器C1, 0表示完全分离, 60%表示完全接合;对于离合器C2, 0表示完全接合, 50%表示完全分离。HEV由电力变矩模式起步, 当车速上升至3m/s时, 切换至发动机驱动模式。

模式切换过程中, 离合器C1先根据加速踏板开度和踏板开度变化率查MAP图 (图4) 得到一个合适的初始压力ps (转化成占空比) , 然后根据加速踏板开度变化率和离合器的主从动片转速差查MAP图 (图5) 得到一个压力变化率$\dot{p}‚$离合器接合过程中, 通过式 (8) 计算离合器的控制压力, 当主从动片转速差达到阈值时, 离合器控制压力立即上升到锁止压力, 完成模式切换过程。

电力变矩模式完全符合前面的分析;模式切换过程离合器C1、C2的控制也按照预先制定的控制策略接合;由于发动机驱动模式下离合器C1、C2均接合, 动力耦合装置蜕化成一个刚体, 行星机构3个构件连同变速器输入轴以同样的转速转动。进入发动机驱动模式后车速随加速踏板开度的增大而迅速提升。

5 结语

动力耦合装置是混合动力汽车的关键部件之一, 关系到整车动力性和乘坐舒适性。行星机构动力耦合装置巧妙地利用了行星机构的两自由度结构, 并通过离合器和制动器限制了一个自由度, 实现了混合动力汽车10种工作模式的动力传递要求。试验测试表明, 行星机构动力耦合装置不但实现了混合动力汽车的所有功能, 在模式切换过程中, 在保证汽车动力性的前提下还有效地避免了冲击, 提高了混合动力汽车的舒适性。

参考文献

[1]黄贤广, 何洪文.混合动力车辆动力耦合装置特性研究[J].上海汽车, 2008 (14) :4-7.

[2]黄贤广, 林逸, 何洪文, 等.混合动力汽车机电动力耦合系统现状及发展趋势[J].上海汽车, 2006 (7) :2-5.

[3]李兴华, 何国旗.等效杠杆法分析行星齿轮传动[J].机械设计, 2004, 21 (1) :44-45.

[4]Endo Hiroatsu, Ito Masatoshi, Ozeki Tatsuya.De-velopment of Toyota’s Transaxle for Mini-van Hybrid Vehicles[J].SAE Paper, 2002-01-0931.

[5]Oba Hidehiro, Yamanaka Akihiro, Katsuta Hiroshi, et al.Development of a Hybrid Powertrain System Using CVTin a Minivan[J].SAE Paper, 2002-01-0991.

[6]邹乃威, 刘金刚, 周云山, 等.混合动力汽车行星机构动力耦合器控制策略仿真[J].农业机械学报, 2008, 29 (3) :5-9.

[7]葛安林, 吴锦秋, 郭万富.离合器最佳接合规律的研讨[J].汽车工程, 1988 (2) :54-65.

[8]邹乃威.无级变速混合动力汽车动力耦合及速比控制研究[D].长春:吉林大学, 2009.

机电耦合动力学 篇9

随着现代机械向高速、重载、高精度方向发展, 机械系统的结构逐渐多样化。采用构件刚性假设的传统多刚体系统动力学模型逐渐不适用于越来越复杂的机械系统, 其分析结果会产生一定的偏差以至于不能满足精度要求, 因此考虑系统中构件的弹性变形已成为多体系统动力学的重要研究方向之一[1,2]。

ANSYS是一款大型通用有限元分析软件, ADAMS主要用于多刚体运动学与动力学仿真, 将两款软件结合使用, 即ANSYS-ADAMS联合仿真分析是解决复杂机械系统动力学问题的一种有效的方法[3,4]。

本文基于柔性体理论, 以3-PUS并联机构为研究对象, 采用ANSYS-ADAMS联合仿真的分析方法, 对3-PUS并联机构的刚柔耦合模型进行了仿真分析, 以了解这一机构的动力学特性, 为进一步优化提供理论依据。

1 柔性体理论

1.1 柔性体系统建模理论基础

ADAMS用模态来表示物体的弹性, 这种弹性变形近似地用模态的线性叠加来表示, 刚体的位姿用反映刚体位置的笛卡尔坐标x= (x, y, z) 和反映刚体姿态的欧拉角u= (u, v, w) 表示, 模态坐标用q={q1, q2…qi} (i为模态坐标数) 表示。在完全约束条件下, 唯一能够确定质点系位置的独立变量可以用作广义坐标[5], 则柔性体的广义坐标表示为

根据选取的坐标系的不同, 柔性体建模可以分为3种方法, 包括浮动坐标系法、惯性坐标系法和随转坐标系法[6]。本文采用应用最为广泛的浮动坐标系法建立柔性体的动力学模型。

柔性体系统中建立静坐标系 (er) 和浮动坐标系 (eb) , 如图1所示。静坐标系是固定不变的, 作为参考坐标系, 浮动坐标系建立在柔性体上, 以其相对于静坐标系的平移和转动来表述柔性体的位姿。柔性体相对于静坐标系的运动可以分解为:刚体平移、刚体转动和柔性体的变形运动, 即柔性体上任意一点的位置的向量表示为

式中:ri为柔性体中一点在静坐标系中的向量, A为动坐标系到静坐标系的变换矩阵, rr为浮动坐标系原点在静坐标系中的向量, rb为柔性体中一点在浮动坐标系中的向量, up为相对变形量, up=φpζ, 其中φp为点p满足里兹基向量要求的假设变形模态矩阵, ζ为变形的广义坐标。

1.2 柔性体系统的动力学方程

柔性体的动能表达式为

式中, M (ξ) 为质量矩阵, 是一个3×3维的矩阵。柔性体的势能包括重力势能和弹性势能两部分, 其表达式为

式中:Wg (ξ) 为重力势能;2/1ξTξ为弹性势能;K为模态刚度矩阵, 是常值矩阵。

能量损耗函数为

式中:为模态坐标对时间的一阶导数;矩阵D为模态阻尼矩阵, 为常值对称阵。

为得到柔性体系统动力学方程, 将式 (3) ~式 (5) 代入下面的拉格朗日方程:

式中:L为拉格朗日项, L=T-W;Ψ为约束方程。得到柔性体系统的动力学方程为

2 3-PUS并联机构刚柔耦合模型仿真[7]

2.1 构建3-PUS并联机构的多刚体模型

本文研究的3-PUS并联机构包括动平台、静平台、3组连杆滑块和机架, 滑块和机架之间为滑动副, 滑块和连杆之间通过虎克铰连接, 连杆和动平台之间通过球铰相连, 通过电机传动使3个滑块沿导轨直线运动, 进而带动动平台做3个平移运动。将UG NX7.5中建立的3-PUS并联机构模型以Parasilid (*.x_t) 文件格式导入ADAMS中, 在ADAMS中将3-PUS并联机构的材料定义为45钢, 添加正确的运动副, 完成后并联机构的多刚体模型如图2。

2.2 构建3-PUS并联机构的刚柔耦合模型

在本文研究的3-PUS并联机构中, 其连杆的长度和直径的比值较大, 在工作中易发生弹性变形, 因此不能忽略连杆的柔性, 将ADAMS中3-PUS并联机构中原来的刚性连杆替换为柔性连杆并建立3-PUS并联机构的刚柔耦合模型再进行动力学仿真, 可得出较为精确的仿真结果。构建3-PUS并联机构刚柔耦合模型的方法如下:

1) 将UG NX 7.5中的3-PUS并联机构的三维装配模型中的连杆模型以Parasolid (*.x_t) 文件格式输出, ANSYS通过专门的接口读入Parasolid模型文件;

2) 连杆的材料为45钢, 在ANSYS中定义材料的弹性模量E=2.07×1011Pa, 泊松比PRXY=0.27, 密度ρ=7 801 kg/m3;单元类型选择Solid45和mass21;

3) 选用Solid45单元, 使用智能网格划分功能对连杆模型进行网格划分;

4) 根据ADAMS中3-PUS并联机构中构件的连接关系, 为ANSYS中的连杆模型添加Keypoints, 并使用mass21对Keypoints进行网格划分, 生成外部连接点;在ANSYS主菜单中选择“Preprocessor→Coupling/Ceqn→Rigid Region”创建刚性连接区域。

5) 在ANSYS的主菜单中选择“Solution→ADAMS Connection→Export to ADAMS”, 输出的连杆模型的模态中性文件, 如图3。

6) 在ADAMS中, 将3-PUS并联机构中的6个连杆替换为柔性体连杆, 替换后柔性体连杆会继承原来刚性连杆上的运动副、载荷和Marker点等特征, 图4为3-PUS并联机构的刚柔耦合模型。

2.3 3-PUS并联机构刚柔耦合模型的仿真

在ADAMS环境中为图4左侧滑块移动副添加滑移驱动, 使该滑块在前1 s以50 mm/s2的加速度做匀加速运动, 在1~5 s以50 mm/s的速度做匀速运动, 其他2个滑块保持静止, 驱动函数为IF (time-1:50*time, 50, 50) 。

表1为样机模型在ADAMS中的仿真参数设置, 仿真I的仿真时间为5s, 以观察仿真结果的整体趋势, 便于进行全面的分析;仿真II的仿真时间为0.1s, 以观察仿真模型在启动瞬间的运动特性, 是对仿真I的分析结果的补充和深化。

仿真运动中的并联机构如图5所示, 可以看出运动中连杆上Von Mises应力最大值为76.2 MPa, 而45钢的屈服极限为300 MPa, 可知该机构运行安全。

3 仿真结果与分析

多刚体模型与刚柔耦合模型仿真均可得到速度、位移和应力随时间的变化曲线图, 但刚柔耦合模型仿真的结果, 考虑了部件的弹性变形等实际物理特性, 因此反映的运动学和动力学特性更为准确。为说明刚柔耦合这一特点, 本文将3-PUS并联机构的刚柔耦合模型 (部件中含柔性体) 仿真结果和多刚体模型 (部件中均为刚体) 的仿真结果进行对比说明。

3.1 运动学仿真结果分析

仿真模型建好后, 按规划的路径驱动进行仿真, 仿真结束后, 在ADAMS/Postprocessor模块中观察仿真结果图 (图6、图7) , 图中Part2为并联机构的动平台。

图6 (a) 、图6 (b) 分别为仿真I中刚柔耦合模型和多刚体模型中动平台的速度和位移图。

图7 (a) 、图7 (b) 分别为仿真I中刚柔耦合模型和多刚体模型中动平台的X向速度和位移图 (动平台的Y向和Z向的速度和位移对比图与X向的速度和位移对比图与图7类似, 这里不再给出) 。

对比分析图6和图7, 可得出以下结论:

1) 比较两模型速度曲线可以看出, 速度曲线基本一致, 说明刚柔耦合模型的动平台的速度与多刚体模型的动平台的速度均值及整体变化趋势一致;但在仿真运动的开始前0.5 s内, 刚柔耦合模型中动平台的速度曲线有较为剧烈的波动现象。

2) 比较两模型位移曲线可以看出, 图6、图7中刚柔耦合模型和多刚体模型的位移曲线基本相同。

这里需要分析的问题是速度波动, 应该反映在位移曲线上 (位移是速度的积分) , 但仿真的结果两种模型的位移曲线却基本相同。原因可能是刚柔耦合模型在速度上的短时波动变化, 由于仿真步数的限制, 未能反映在更小时间单位上位移曲线存在波动的真实情况。

为了得到刚柔耦合模型中动平台的真实位移曲线和更清晰的速度曲线, 以确认仿真结果的正确性, 将仿真时间缩小为0.1 s、步长为20 000即仿真II, 以进一步观察和分析两者的动平台的速度和位移曲线的对比图 (如图8和图9) 。

对比分析图8和图9, 可得出以下结论:

1) 速度曲线表明:在仿真运动的初始时刻, 刚柔耦合模型中动平台的速度曲线呈周期约为0.01 s的简谐波, 并且其振幅逐渐衰减, 在t=0时其速度为0, 在t=0.001 8 s时刻速度就达到了54.022 5 mm/s, 而多刚体模型中动平台的速度t=0.001 8 s时速度仅为0.116 3 mm/s。

刚柔耦合模型仿真结果说明在实际运动的初始时刻, 动平台的速度以较高的频率和较大的振幅振荡, 并联机构动平台由静止到运动瞬间存在一定的振动, 其原因是杆件存在弹性变形。依据振动的特性与波形, 可以进一步地优化杆件结构尺寸以减小速度曲线的频率和振幅, 从而使并联机构的启动阶段更加平稳。

2) 位移曲线表明:在图8 (b) 中可以看出刚柔耦合模型动平台的位移曲线的斜率随速度的波动出现变化, 图10 (a) 中表示在前0.1 s内, 动平台的X向位移曲线的斜率随X向速度的波动也出现变化, 动平台的位移曲线也呈现一定频率的波形。

图6、图7中, 在0~5s时间内速度发生一定频率的波动, 但由于时间步长的原因, 不能在位移曲线发生相应的波动。而图8 (b) 、图9 (a) 在更小时间单位上, 刚柔耦合模型中动平台位移曲线证明了波动的存在, 从而反映了动平台的实际运动特性。

3.2动力学仿真结果分析

图10为仿真I中刚柔耦合模型和多刚体模型中球铰5处 (图4右前位置处) 沿X向、Y向、Z向受到的力及合力的对比图;图11为仿真II中球铰5处沿X向、Y向、Z向受到的力及合力对比图。

对比分析图10得知, 在仿真运动开始的0.5s, 刚柔耦合模型中球铰5处受到的力均以较高的频率和较大的振幅往复振荡, 之后球铰5处的力的变化趋于稳定, 总体上刚柔耦合模型和多刚体模型中球铰5处的力的变化趋势一致, 但是多刚体模型中球铰5的力学曲线并未反映出在运动的初始时刻球铰5处的力的振动现象, 在图11中可以看出刚柔耦合模型中球铰5处的受力曲线呈振幅逐渐衰减的简谐波形, 周期约为0.01 s, 在t=0.015 s处的力达到了最大值为2 886.735 4 N, 而刚体模型中球铰5在仿真运动的5s内最大力为1 531.419 7 N, 其原因在于在仿真运动初始时刻, 因惯性的存在, 刚柔耦合模型中的柔性连杆会发生弹性变形, 因此在运动初始时刻球铰5的应力变化更大。

动力学仿真结果进一步说明刚柔耦合模型可以更准确地反映部件受力情况, 从而为部件产品设计、优化和校核提供依据。

4 结论

1) 应用ADAMS等工程软件构建复杂机械系统的多刚体模型和刚柔耦合模型进行运动学和动力学仿真和求解, 为复杂机械系统的设计分析和优化提供了一种实用的方法。

2) 在建立虚拟样机时, 对含有长细比较大或在实际工况中易发生变形的部件的机械系统, 应考虑构建刚柔耦合的虚拟样机, 以获得更为真实有效的仿真结果。

3) 采用ANSYS和ADAMS联合仿真的分析方法, 对3-PUS的并联机构的刚柔耦合模型和多刚体模型进行了仿真, 分别将两者的结果进行了整体和局部的对比, 发现该并联机构在启动瞬间, 速度、位移、应力等均存在一定的波动, 刚柔耦合模型仿真结果反映了该并联机构真实的动态特性, 为并联机构的设计和优化提供了理论依据和方法指导。

摘要：以柔性体理论为基础, 使用多体动力学仿真软件A D A M S和有限元分析软件A N SY S建立了3-PU S并联机构的刚柔耦合模型, 对3-PU S并联机构的刚柔耦合模型和多刚体模型进行了仿真和分析。结果表明, 在仿真运动开始的瞬间, 刚柔耦合模型中动平台的速度、位移和铰接处的应力均有一定频率的波动现象, 反映了该并联机构实际的动态特性, 为并联机构的进一步的设计、优化和校核提供了依据。

关键词：并联机构,ADAMS,ANSYS,刚柔耦合模型

参考文献

[1]Shabana A A.Flexible multibody dynamics:review of past and recent developments[J].Journal of Multibodt System Dynamics, 1997 (2) :189-222.

[2]刘铸永, 洪嘉振.柔性多体系统动力学研究现状与展望[J].计算机力学学报, 2008, 25 (4) :411-415.

[3]张雷, 余克龙, 陈少钟, 等.基于ANSYS和ADAMS空间四连杆引纬机构的柔性动力学仿真[J].纺织学报, 2013, 34 (5) :116-120.

[4]张永德, 汪洋涛, 王沫楠, 等.基于ANSYS与ADAMS的柔性体联合仿真[J].系统仿真学报, 2008, 20 (17) :4501-4504.

[5]王爱勤, 冯振宇.也谈虚位移原理中的广义坐标的选择[J].大学物理, 2003, 22 (8) :13-14.

[6]Wasfy T M, Noor A K.Computational strategies for flexible multibody system[J].Appl Mech Rev, 2003, 56 (6) :553-613.