带有非线性联轴器轴系稳态响应计算方法的研究

5-1带有非线性联轴器轴系力学和数学模型的建立
在第三章，根据试验的基础建立了具有非线性迟滞特性联袖器的恢复力模型。这一章，将研究一个带有这种联轴器的轴系如图5-1所示，该轴系有n个圆盘，由钢丝绳联轴器与主机相连接，按以下原则建立力学模型：
1.每个圆盘均视为刚性匀质，所有圆盘的质量mi都有不同程度的偏心距ei。转轴轴线垂直通过各圆盘的几何中心；
2.设轴承、轴承座以及联轴器各向同性，静坐标系如图5-1所示，支座处理成简支；
3.联轴器从动端处理成一集中质量m_b，由于联轴器主动端与主机轴相连接，主机袖相对轴系轴来说，较短较粗，刚性较大，变形较小，故将联轴器主动端、主机轴视为一体，位移为零，联轴器从动端与主动端之间由非线性弹簧及非线性阻尼器联结；
4.轴系为小位移振动，且忽略回转效应。

当主机以角速度ω转动时，轴系在各偏心力的作用下产生稳态振动，其运动微分方程可表示为：

式中 分别为轴系的质量和左阵、阻尼矩阵和刚度矩阵。其中阻尼矩阵，a，b为比例常数。

其中：y_i,z_i分别为轴系中第i个圆盘在y,z方向的位移分量； =dy_i/dt， =d²y_i/dt²， =dz_i/dt， =d²z_i/dt²分别为对应于位移分量y_i,z_i的速度和加速度，y_b,z_b为联轴器从动端集中质量块m_b在y,z方向的位移； =dy_b/dt， =d²y_b/dt²， =dz_i/dt， =d²z_d/dt²分别为对应于位移分量yb,zb的速度和加速度；,为y,z方向的激励力向量，其中m_ie_iω²cos（ωt+ ），m_ie_iω²sin（ωt+ ）分别为y,z方向的偏心力分量； （ ，y_b， ，ωt), （ ，z_b， ，ωt）分别为非线性弹性联轴器恢复力 在y,z方向的分量。 为第i个圆盘质量偏心的初相位角，即轴系静止时，第i个圆盘质心和几何中心连线与水平轴的夹角。
由于假设轴承、轴承座以及联轴器各向同性，式（5-1〕、式（5-2）解法相同，只需要讨论二者之一即可。

式（5-l）表示一个局部非线性弹性和阻尼元件的振动系统，对这种系统，按现有常规方法来求解是非常困难的，为此，本文以GILM为基础发展了一种称为SSGILM(Separate System-Gal-erkin and Improved Levenbery-Marquar-dt）的方法来求解此类微分方程组。

5-2  SSGILM法

一.振动微分方程组的改写和解耦
首先把有局部非线性系统振动微分方程组（5-1）式改写成只有线性常系数的微分方程组和一个具有非线性变参数的微分方程两部分：

式中：

式（5-4）的P_y中含有y_b和 而（5-5）的F_j中含有y_j，（j=1，2，… n）。这样，轴系分成运动微分方程耦联的两个子系统一线性轴系子系统和非线性联轴器子系统。

式（5-4）为线性方程，按常规方法可求出无阻尼的各阶固有频率p_i以及对应的主振型Y_i向量（i=1，2，…，n），Y_i分别除以相应广义质量的平方根 得到正则振型Y_Ni向量。引入正则振型坐标W_Ni，对（5-4）式进行坐标变换；

式中，W_Ni为正则振型向量W_N中的第i个分量， 和 ，分别为它对时间t的一次导数和二次导数；I为单位矩阵； s_i=（a+ ）/2P_i=c_ii/2p_im_i为振型比例阻尼比；P_Ni为P_N= P_y激励力向量中的第i个分量。
经正则坐标变换，虽然得到互不耦合的线性微分方程组（5-9)，但是（5-9）式仍然不能像单自由度振动系统那样求解，因为（5-9）式中的激励力包含未知的振动位移y_b和振动速度 。如果y_b已知，就可以从（5-9）式中解得W_Ni，进而可以由（5-7）式得到各y_i。因此需要求得y_b。
图5-1所示轴系在主机带动下转动时，各圆盘质量偏心将产生周期偏心激励力，根据第二章的试验结果，联轴器非线性恢复力Q是时间的周期函数，因此当轴系中有联轴器这种局部非线性元件时，可以设它的位移响应是周期性的，即假设y_b有下面的形式：

由式（5-6）、（5-8）和（5-9）可知，激励力由二类力构成，一类是质量偏心力m_ie_iω²cos（ωt+ ，另一类是恢复力k_i（n+1）y_b-c_i（n+1） ，为求解（5-8）式方便，将（5-10）代入（5-8）式，并将正则激励力分解成：

P_N= P_y=P_N1+P_N2                                        （5-11）

式中：

由（5-11）和（5-12）式可知，若{a}已知，就可以由（5-12）式求出各W_Ni，代W_Ni入（5-7）式，可以得到各y_i。由此可以把求y_b的问题转化为求{a}=[a₀  a₁  的问题。

二、Galerkin法的应用

为了求出{a}，引入无量纲时间τ= ，则（5-10）式可以写成：

式中[q₀  q₁   … q_N  ]^T的各分量为联轴器恢复力 的傅立叶级数各简谐分量的系数分量。[f₁₀  f₁₁    …  f_1N    ]^T，…[f_jo  f_j1   … f_jN  ]^T，…[f_n0  f_n1   … f_nN  ]^T分别为弹性恢复力和阻尼恢复力的F₁=k_（n+1） +c_（n+1） …，F_j=k_（n+1）+ c_（n+1） …，F_n=k_（n+1）+_（n+1） 的傅立叶级数各简谐分量的系数向量。

由以上推导可知，求解带有非线性联轴器轴系的稳态响应，可以归结为联合求解2N+1个非线性代数方程组（5-15）和求解n个微分方程组（5-4)而重点又在于由（5-15）式求得2N

+1个未知量［a₀  a₁  ］^T。由(5-15）式可知，它不同于一般的非线性方程组，在一般的非线性方程组中，未知量是显含的，而（5-15）式的未知量{a}是隐含的，而且{q｝是{a｝的非线件函数。不能通过将y_b的傅立叶展开式代人恢复力O的表达式中来直接求得，这是因为Q表达式中的振幅 不能表示成{a}的显函数的缘故。因此要想由（5-1）式解得｛a}，还需要解决两个问题：一是如何求 ,二是用什么来解隐含{a}的方程组（5-15)。下面就决这些问题。

三、ILM法

将非线性方程组简写成：

由于向量函数{r（{a}）}中非线性项{q（{a}）}的复杂性，可行的求解方法是采用数值法。

对于式（5-17）的数值法求解问题，我们将其转化为与之等效的目标函数：

的极小值问题。即非线笥方程组（5-17）的解{a}可以通过求此极小值问题的最小二乘解得到。在此采用一种改进的LM算法（ILM）来求解这一问题。ILM法的基本思想简述如下：

当LM法用于（5-18）问题时，迭代公式为：

式中μ_k是LM参数，I为单位矩阵，{p^（k）（μ_k）}为搜索逼近（5-17）式解的校正向量，它依赖于参数μ_k值，控制搜索方向和长度。

LM算法的优点是，在计算（5-20）式时，可以不考虑[J（{a^（k）}）]是否满秩，μ_k可以起到 （{a}）下降参数的作用，{p^（k）（μ_k）}是 （{a}）在{a^（k）}处的下降方向，而且总存在这样的μ≥0，使 （{a^（k）}+ p^（k））< {a^（k）}）（对所有k），随着μ值的增大，{p^（k）（μ_k）}的方向越来越接近 （{a^（k）}的负梯度方向，此时，收敛速度变慢，但放宽了对初始近似{a^（0）}的要求。缺点是选择μ满足 （{a^（k+1）}）< （{a^（k）}）时，在一个迭代步骤中需要多次求解方程组（5-20），大大增加了计算量，因此有必要对LM法进行改进，以使选择μ_k时即不增加太多的计算量，又能保证满足下降性质，即有一定的收敛速度，又能改善（5-20）方程组的病态性。改进后的迭代公式为：

式中才 为矩阵 ^（k）的元素。 ， 分别为矩阵L^（k）， ^（k）的元素。公式（5-23）和（5-27）和（5-28）中不含参数声μ_k，因此改变μ_k时，只需重新计算(5-24),而不需要重新进行三角分解。显然，改进后的算法减少了由于选择合适的μ_k而带来的计算量。这种算法与LM法的不同之处在于用正定矩阵L_k 代替了单位矩阵I，它不仅调整了矩阵 ^（k）的主对角元素，而且对整个矩阵进行了调整，它比LM法有更好的收敛速度。
为了得到合适的μ_k值，使得搜索向量{p^（k）（μ_k）}更快更好地满足 （{a^（k）}+ p^（k））< {a^（k）}），加速收敛，改善LM法对初始值{a^（0）}要求过高的缺点，本文采用一种选择合适μ_k值的新算法，对LM法进行了改进，选择计算步骤如下：

①.取定初始近似{a^（0）}和初始值 ，为了调整μ_k值，引人调整系数v, 可以取得小一些，取 =10^-2，v=5。

②.按（5-21) 式首先算得［J（{a^（0）}）]，然后按（5-26）式算得0⁽⁰⁾ ，再根据（5-27）、（5-28）式算得矩阵L^(0）和 ^（0）中的各元素，按(5-25）、(5-23）和（5-24）算得{P^（0）（ ）}；按（5-17）算得{r（{a^（0）}）}；按（5-18）算 （{a^（0）}）。并计算{a^（1）}={a^（0）}+{P^（0）（ ）}和 （{a^（1）}）
③.(l）如果少 （{a^（1）}）≤ （{a^（0）}），则以 /v代替 重新计算{P^（0）（ /v）}和新的{ ^（1）}， （{ ^（1）}）。（a）如果有 （{ ^（1）}）≤ （{a^（0）}），则取定μ₀= /v；（b）如果 （{a^（1）}）> （{a^（0）}），则取定μ₀= ，（2） （{a^（1）}）> （{a^（0）}），则需增大μ值，取 =v^m ，以 代替 ，逐次对m=1，2，…时的 代入（5-24）进行计算，最后取定使不等式 （{a^（0）}）+{ P^（0）（ ）}）< （{a^（0）}）成立的最小正整数 对应的 作为μ₀以及{a^（1）}+{P（ ）}。

④.假定已求得第k次近似{ a^（k）}和相应的 =μ_k-1，算出[J（{ a^（k）}）]， ^（k），L^（k）， ^（k），r（{ a^（k）}）}， （{ a^（k）}），然后算得{P（ ）}，并计算{ a^（k+1）}={ a^（k）}+{P（ ）}和 （{ a^（k+1）}）。

⑤.（1）如果 （{ a^（k+1）}）≤ （{ a^（k）}），则减少 以/V代替 重新计算{P（/V ）}和对应的{ ^（k+1）}及 （{ ^（k+1）}）（a）如果 （{ ^（k+1）}）≤ （{ a^（k）}），则取定μ_k=/V ；（b）如果 （{ ^（k+1）}）> （{ a^（k）}），则取定μ_k= （2）如果 （{ ^（k+1）}）> （{ a^（k）}），则增大 值，可取 =v^m ,以 代替 ,逐次时m=1，2，…时的 代入（5-24）式进行计算，最后取定使不定式 ({ a^（k）}+{P（ ）}< （{ a^（k）}）。

⑥.目标函数 （{a}）的k次迭代和k+1次迭代如果满足绝对误差和相对误差精度，即可将{ a^（k+1）}作为近似解。如果不满足，则以{ a^（k+1）}代替{ a^（k）}，μ_k代替μ_k-1去执行步骤⑤，直到满足要求为止。

至此，（5-17）式的求解方法已经建立起来，日（5-17）式求解{a}还需要求出 ，{q}和（f_j）（j=1，2，…，n）。另外，在用（5-21）式求[J（{a}）]时遇到r_i（{a}）对隐含向量{a}求偏导数的麻烦，因此在计算[J（{a}）]时，需要做些处理。由（5-15）式可知，它是频域内的2N+1个非线性代数方程组，而非线性迟滞恢复力 是时域函数，因此在求解（5-15）时需要进行频域与时域相互变换，并做一些数学上的处理，下面逐一叙述。

5-3一步用SSGILM法计算带有非线笥联轴器轴系稳态响应

为了求得未知向量{a}，还需做以下工作：

一、向量函数{r（{ a^（k）}）}的计算

要计算{r（{ a^（k）}）}为联轴器恢复力 的傅立叶级数展开式各简谐分量的系数向量。对于给定的{ a^（k）}，不能由时域表达式求得频域内的{q（{ a^（k）}）}，而必须按以下步骤计算。

首先，由给定的频域向量{q（{ a^（k）}）}，根据 的N个时域离散值，记为Y（1）～Y（N）和DY（1）～DY（N），其次，由下式求得 幅值

然后，将 以及 时域离散值代入第三章中得到的恢复力 的表达式（3-8）、（3-12）、（3-27）和（3-41）式，算得对应于{ a^（k）}的恢复力Q在一个周期内N个Q的时域离散值进行FFT变换，即可得到非线性恢复力Q在频域内的各谐波分量的系数向量{q（{ a^（k）}）}

2. {f_j（{ a^（k）}）}的计算

{f_j（{ a^（k）}）}是式（5-6）中恢复力F_j的傅立叶级数展开式各谐波分量的系数向量。当向量{ a^（k）}给定时，首先由5-2节中所述方法算得各个 和 （j=1，2，…，n），其次由（5-6）式中恢复力F_j的各谐波分量的系数向量{f_j（{ a^（k）}）}。

算得{q（{ a^（k）}）}和{f_j（{ a^（k）}）}后，代入（5-15）式即可得到{r（{ a^（k）}）}。

二、雅可比矩阵[J（{ a^（k）}）]的计算

为计算方便起见，将雅可比矩阵式（5-21）分解成三个矩阵：

[J（{ a^（k）}）]=[J₁]+[J₂]+[J₃]                  （5-30）

1.[J₁]的计算

由（5-15）和（5-21）式可得：

[J₁]=

2. [J₂]的计算

[J₂]的计算实际上就是对 {q（{ a^（k）}）}/ { a^（k）}的计算。由于恢复力Q的函数表达式是{ a^（k）}的隐函数，因此不可能由Q的表达式直接进行FFT变换来求得{q（{ a^（k）}）}，也可能由{q（{ a^（k）}）}直接对{ a^（k）}；求偏导数来得到[J₂]：

由第三章中的式（3-8）、（3-12）、（3-41）和本章中的式（5-13）可得：

式中， ，， ，， 可以直接由Q的表达式和（5-13）进行解析计算得到， / { a^（k）}， / { a^（k）}也可以由（5-13）式直接进行解{ a^（k）}，选取增量△{ a^（k）}，然后计算一个周期内对应的 （{ a^（k）}+△{ a^（k）}）和 （{ a^（k）}）值，得：

由（5-32）式算得 / { a^（k）}的时域离散值后，进行FFT变换即可得到频域中的 {q（{ a^（k）}）}/ { a^（k）}。

3. [J₃]的计算

[J₃]的计算实际上就是偏导数 {f_j（{ a^（k）}）}/ { a^（k）}的计算。计算方法与计算 {q（{ a^（k）}）}的方法类似，故不再赘述。对于[J₂]，[J₃]、分别有式（5-34）、（5-35）。

在计算出{r（{ a^（k）}）}和{J（{ a^（k）}）}之后，就可以按5-2节中发展的ILM法进行搜索迭代计算，由初始向量{a^（0）}不断搜索逼近满足精度要求的{ a^（k+1）}，然后代入（5-10）、（5-11）、（5-12）和（5-7）即可求得联轴器从动动端质理块m_b的稳态响应y_b以及轴系中各圆盘转子稳态应Y。同时还可以得到对应的幅频图。

5-4小结
本章对带有非线性联抽器的n+1个自由度轴系在四盘质量偏心产生的激励力作用下稳态响应计算方法进行了研究，以GLM法为基础，提出了一种称为SSG

ILM (Separate System-Ga1erkin and Imp-r0ved Leenberg-Marquardt的方法，用于计算局部具有非线性特性轴系的稳态响应，
主要工作如下：
1.建立带有非线性联铀器n+1个自由度轴系的止）学模型和数学模型，为了求解这种局部含有非线性元件的轴系的稳态响应，将轴系分成两个子系统，即含有n个自由度的线性轴系和含有非线性变参数的联轴器子系统。
2.对线性轴系子系统运动微分方程组用正则坐标变换进行解藕处理；对非线性变参数的联轴器子系统运动微分方程用Galerkin法，变成一组隐含系数向量{a}的非线性代数方程组。
3.针对LM法的缺点，采用已有的选择合适μ值的方法和加快收敛速度的算法对LM法进行了改进，用改进算法对上面得到的非线生方程组中的未知向量{a}进行搜索计算，最后解出局都具有非线性特性轴系的稳态响应。