CN103942374A - 一种刻画振动系统周期振幅与系统参数的关系的方法 - Google Patents

Abstract

本发明公开了一种刻画振动系统周期振幅与系统参数的关系的方法。许多自然现象都具有周期性，大到天体运行，小到电子旋转。社会现象也其周期特性，如社会经济的涨跌，甚至民族的兴衰。研究这些周期现象可以预测至控制其运动演化规律。本发明给出的方法可以用于含有任意可数维数参数变量的任意自由度振动系统，可以计算参数变量的任意给定区域，可以求出具有多个周期解的振动系统的所有周期解，为振动系统的研究提供了有效的手段，为周期现象的研究提供了可靠的方法。

Description

一种刻画振动系统周期振幅与系统参数的关系的方法

技术领域

本发明公开了一种刻画振动系统周期振幅与系统参数的关系的方法，属于振动控制领域。

背景技术

狭义的振动又称机械振动，是指物体在平衡位置附近作的往复运动。广义的振动是指描述系统的状态变量在其基准值上下交替变化的过程。振动可分为线性振动和非线性振动。线性振动通常是对非线性振动简化和近似，自然界中真实的振动都是非线性的。周期振动是常见的振动形式。许多自然现象都具有周期性，大到天体运行，小到电子旋转。社会现象也其周期特性，如社会经济的涨跌，甚至民族的兴衰。研究这些周期现象可以预测至控制其运动演化规律。

振动系统通常用微分方程来描述，形式如下：

$\stackrel{\·}{x}=g(t,x,\mathrm{\ϵ}).---\left(1\right)$

其中，t是标量，表示时间。x通常是向量，用来描述系统的状态变量，如物体的位移、电流的大小、声音的强弱，种群的数量，经济的总量等。表示状态变量对时间的导数，通常描述系统状态变量的变化速率。ε通常是向量，用来描述系统的参数变量，如惯性、弹性、激励的幅值或频率。g(t,x,ε)是关于时间、状态变量以及参数变量的非线性向量函数。

利用摄动法得到原系统的一阶近似解的调制方程具有如下形式：

X′＝G(X,ε).            (2)

其中X描述原系统状态变量的一阶近似解的幅值和相位，X′表示X对一阶时间尺度的导数，G(X,ε)是非线性向量函数。

方程(2)的右端等于零时，调制方程有平衡解，即原系统状态变量的一阶近似解的幅值和相位都是常量，这对应原系统出现周期解。求解方程(2)的平衡解，等价于求解

G(X,ε)＝0.             (3)

方程(3)称为调制方程的平衡解等价方程，是原系统的周期振动幅值和相位X关于系统参数变量ε的隐式函数。由于G(X,ε)通常是非线性函数，X很难表示成ε的显式函数。因此，数值求解变系数非线性方程(3)成为得到原系统周期振动幅值和参数变量的关系的重要手段。有些软件含有求解非线性方程组或带参数非线性方程组的函数，如MATLAB的fsolve函数。该函数需要事先给出初始值，而且只能得到由该初始值求出的单个参数点处的一个解，但是振动系统周期振幅是参数变量的多值函数，而且刻画周期振幅与参数变量的关系需要在参数变量的整个取值范围进行。

根据一些科技文献，国外有些机构研发了一些用于非线性分析的软件，可以用来计算微分振动系统的周期解，如AUTO-07P，AnT，XPP-AUT等。这些软件通常运行在Linux环境下，或者需要在Windows系统中虚拟出Linux环境，而且其安装运行都很繁杂。国内有些文献中出现过类似的计算，但是其中的计算方法都是针对简单的系统，例如单自由度系统，或者单参数系统。

本发明针对以上情况给出一种刻画振动系统周期振幅与系统参数的关系的方法。本发明所述的方法利用嵌套循环结构遍历所有参数变量的取值范围，在每个参数点处计算周期振幅，从而刻画振动系统周期振幅与系统参数的关系。本发明可以在参数变量的取值范围内求出振动系统的多个周期解，并假设每个解在各自参数点处的邻域是连续的，由此以数值解的收敛性作为退出条件向该邻域的上下边界扩展，从而提高计算效率。本发明可以应用于任意可数维参数变量和任意自由度的振动系统，为振动系统的研究提供了有效的手段，为周期现象的研究提供了可靠的方法。

发明内容

本发明的目的在于，为克服现有方法的不足，给出一种刻画振动系统周期振幅与系统参数的关系的方法。针对具体振动系统的调制方程X′＝G(X,ε)，X＝[X_m]表示m维状态变量，ε＝[ε_n]表示n维参数变量，本发明可以应用于任意可数维参数变量和任意可数维状态变量的调制方程。

1.本发明首先录入具体振动系统的状态变量的一阶近似解的调制方程的平衡解等价方程,即方程(3)。为了控制程序的运行，本发明定义了一些计算所需要的全局变量：记录解的最大个数的变量，限定最大求解数的变量，存储计算结果的变量，记录平衡解数据大小的变量，判断相同解的容差变量，限定振动系统维数的变量，划分参数变量网格数的变量，限定参数变量范围的变量。本发明给出嵌套循环结构遍历参数变量的整个取值范围，此嵌套循环结构的每一层循环结构对应一维参数变量ε_i(i＝1...n)。在嵌套循环结构的最底层，记ε＝ε⁰，即每一维参数变量都取一个固定数值。在嵌套循环结构的最底层定义局部变量：记录所给参数点处平衡解的变量，记录所给参数点处平衡解个数的变量，记录试探次数的变量，限定最大试探次数的变量。在嵌套循环结构的最底层，再以试探次数不超过最大限定试探次数和求得的解的个数不超过限定最大求解数作为退出条件给出一个循环结构。在此循环结构的内部，在求解范围内随机给出初值X＝X⁰，并用迭代方法X＝X(X⁰)求出ξ使得调制方程的平衡解等价方程满足G(ξ,ε⁰)＝0。如果所求得的平衡解不与已经求解并记录的平衡解重复，利用局部变量记录此求得的平衡解。假设ξ在ε⁰处连续，即ξ+Δξ＝ξ(ε⁰+Δε)，以收敛性为退出条件向该平衡解所在的参数点的邻域的上下边界扩展即令ε＝ε⁰+Δε，以X＝ξ为初值迭代求解方程(3),若能得到收敛的数值解ξ+Δξ，则继续向边界扩展，否则就停止扩展。利用全局变量记录所有扩展后的收敛的平衡解，画出参数变量与平衡解的关系图，刻画出振动系统周期振幅与系统参数的关系X＝X(ε)。

其中，状态变量可以是任意随时间往复变化的特征量，例如机械的振动、电磁的振荡、种群数目的变化、经济总量的起伏。系统参数是指可以人为调节的固有的系统性质，如惯性、弹性、外激励幅值及频率。

本发明所给出的方法可用于任意可数维状态变量的调制方程，即用于任意自由度的振动系统。

本发明所给出的方法可以用于含有任意可数维数参数变量的振动系统。

所定义的限定最大求解数的变量，本发明所给出的方法可以求出具有多个周期解的振动系统的所有周期解。

本发明所给出的方法可以计算参数变量的任意取值区域。

所述的假设平衡解在其参数点的邻域内连续，以收敛性为退出条件向该平衡解所在的参数点的邻域的上下边界扩展，该方法提高了程序的计算效率。

本发明的有益效果是：1)本发明给出的方法可以用于任意自由度的振动系统；2)本发明给出的方法可以用于含有任意可数维数的参数变量的振动系统；3)本发明给出的方法可以计算参数变量的任意给定区域；4)本发明所给出的方法可以求出具有多个周期解的振动系统的所有周期解。

附图说明

图1是本发明所述方法的计算程序设计流程图。

图2是本发明所述方法实施所得一维参数变量四维调制方程的平衡解，描述二自由度振动系统的周期振幅与一维参数变量外激励频率的关系以及周期振幅与一维参数变量外激励幅值的关系。

图3是本发明所述方法实施所得二维参数变量四维调制变量的平衡解,描述二自由度振动系统的周期振幅与二维参数变量外激励幅值和频率的关系。

具体实施方式

本发明可以用任意编程语言实现。MATLAB是常用的数学计算软件，本发明利用MATLAB软件通过具体实施例并结合附图作进一步详细的描述。

实施例

本发明利用MATLAB编写程序计算文献《Three-to-One Internal Resonances inParametrically Excited Hinged-Clamped Beams》中调制方程(31-34)：

$\begin{array}{c}{p}_1^{\′}=-{\mathrm{\μ}}_1{p}_1-v_1{q}_1+{\mathrm{\γ}}_{11}{q}_1(p_1^2+q_1^2)+{\mathrm{\γ}}_{12}{q}_1(p_2^2+q_2^2)\\ -{\mathrm{\δ}}_1[2p_1{q}_1{p}_2-q_2(p_1^2-q_1^2)]-s_{11}F{q}_1+s_{12}F{q}_2,\end{array}---\left(4a\right)$

$\begin{array}{c}{q}_1^{\′}=-{\mathrm{\μ}}_1{q}_1+v_1{p}_1-{\mathrm{\γ}}_{11}{p}_1(p_1^2+q_1^2)-{\mathrm{\γ}}_{12}{p}_1(p_2^2+q_2^2)\\ -{\mathrm{\δ}}_1[2p_1{q}_1{q}_2+p_2(p_1^2-q_1^2)]-s_{11}{\mathrm{Fq}}_1-s_{12}F{q}_2,\end{array}---\left(4b\right)$

$p_2^{\′}=-{\mathrm{\μ}}_2{p}_2-v_2{q}_2+{\mathrm{\γ}}_{21}{q}_2(p_1^2+q_1^2)+{\mathrm{\γ}}_{22}{q}_2(p_2^2+q_2^2)+{\mathrm{\δ}}_2[3p_1^2{q}_1-q_1^3]+s_{21}{\mathrm{Fq}}_1,---\left(4c\right)$

$q_2^{\′}=-{\mathrm{\μ}}_2{q}_2+v_2{p}_2-{\mathrm{\γ}}_{21}{p}_2(p_1^2+q_1^2)-{\mathrm{\γ}}_{22}{p}_2(p_2^2+q_2^2)+{\mathrm{\δ}}_2[3p_1{q}_1^2-p_1^3]-s_{21}{\mathrm{Fq}}_1.---\left(4d\right)$

令[q′₁,p′₁,q′₂,p′₂]＝0，方程(4)有平衡解，其对应的平衡解等价方程具有方程(3)的形式，

G(X,ε)＝0.           (5)

其中

X＝[q₁,p₁,q₂,p₂],              (6a)

ε＝[F,σ₂],                      (6b)

$\begin{array}{c}{G}_1=-{\mathrm{\μ}}_1{p}_1-v_1{q}_1+{\mathrm{\γ}}_{11}{q}_1(p_1^2+q_1^2)+{\mathrm{\γ}}_{12}{q}_1(p_2^2+q_2^2)\\ -{\mathrm{\δ}}_1[2p_1{q}_1{p}_2-q_2(p_1^2-q_1^2)]-s_{11}{\mathrm{Fq}}_1+s_{12}{\mathrm{Fq}}_2,\end{array}---\left(6c\right)$

$\begin{array}{c}{G}_2=-{\mathrm{\μ}}_1{q}_1+v_1{p}_1-{\mathrm{\γ}}_{11}{p}_1(p_1^2+q_1^2)-{\mathrm{\γ}}_{12}{p}_1(p_2^2+q_2^2)\\ -{\mathrm{\δ}}_1[2p_1{q}_1{q}_2+p_2(p_1^2-q_1^2)]-s_{11}{\mathrm{Fq}}_1-s_{12}{\mathrm{Fq}}_2,\end{array}---\left(6d\right)$

$G_3=-{\mathrm{\μ}}_2{p}_2-v_2{q}_2+{\mathrm{\γ}}_{21}{q}_2(p_1^2+q_1^2)+{\mathrm{\γ}}_{22}{q}_2(p_2^2+q_2^2)+{\mathrm{\δ}}_2[3p_1^2{q}_1-q_1^3]+s_{21}{\mathrm{Fq}}_1,---\left(6e\right)$

$G_4=-{\mathrm{\μ}}_2{q}_2+v_2{p}_2-{\mathrm{\γ}}_{21}{p}_2(p_1^2+q_1^2)-{\mathrm{\γ}}_{22}{p}_2(p_2^2+q_2^2)+{\mathrm{\δ}}_2[3p_1{q}_1^2-p_1^3]-s_{21}{\mathrm{Fp}}_{1.}---\left(6f\right)$

方程(5)的解对应原系统的周期解，即文献中方程(30)中A_k为常数。方程(30)如下

$A_k\left(T_1\right)=0.5[p_k\left(T_1\right)-{\mathrm{iq}}_k\left(T_1\right)]e^{{\mathrm{i\λ}}_k},k=\mathrm{1,2}.---\left(7\right)$

A₁和A₂分别是文献中方程(23)的第一阶和第二阶周期振幅。方程(23)如下

$w_1(T_0,T_1,x)=A_1\left(T_1\right){\mathrm{\φ}}_1\left(x\right)e^{{\mathrm{i\ω}}_1{T}_0}+A_2\left(T_1\right){\mathrm{\φ}}_2\left(x\right)e^{i{\mathrm{\ω}}_2{T}_0}+\mathrm{cc}.---\left(8\right)$

其中，w₁是文献中梁结构的横向振动位移。

方程(5)是含有两个可变参数的四维调制方程，两个可变参数分别是外激励的幅值F，外激励频率与结构第二阶自然频率的差值σ₂，四维状态变量是[q₁,p₁,q₂,p₂]。方程(5)中其他参数均为已知。

1.图1给出N维参数变量任意可数维调制方程的平衡解的计算程序设计流程。为与文献中的计算结果进行对比，本实施例限定N＝2，以F和σ₂作为参数变量，即本实施例计算二维参数变量四维调制方程的平衡解。程序首先录入方程(5),然后定义全局变量：记录解的最大个数的变量“NumOfSol”，限定最大求解数的变量“MaxMark=6”，存储计算结果的变量“Result”，记录平衡解数据大小的变量“TotalCounter”，判断相同解的容差变量“ValTol=0.01”，限定振动系统维数的变量“Dim=4”，划分参数变量网格数的变量“NumOfF=20”和“NumOfSigma=20”，限定参数变量范围的变量“FLim=[0,50]”和“SigmaLim=[0,200]”。

2.如图1所示，给出嵌套循环结构，即第1层至第N层循环，每一层循环对应一维参数变量。内层循环遍历完本层对应参数范围后返回其最近外层继续遍历外层所对应参数范围，直到遍历完参数变量所给定范围，退出嵌套循环结构，结束计算。本实施例给出两层嵌套循环结构，第一层遍历参数F的取值范围，第二层遍历参数σ₂的取值范围。

3.如图1所示，在嵌套循环结构的最底层，即第N＝2层，定义局部变量：记录所给参数点处平衡解的变量“MyResult”，记录所给参数点处平衡解个数的变量“MyMark”，记录试探次数的变量，限定最大试探次数的变量“MyNumTol=30”。

4.如图1所示，嵌套循环结构的最底层再以试探次数“MyMark”不超过最大限定试探次数“MyNumTol”和求得的解的个数“MyMark”不超过限定最大求解数“MaxMark”作为退出条件给出一个循环结构。在此循环结构的内部，在求解范围内随机给出初始值“Ini”，并利用fsolve函数求解方程。如果所求得的平衡解符合收敛条件且不重复，即满足收敛条件且在全参数范围和该参数点处都不与已经求得的解相同，利用局部变量“MyResult”记录此求得的平衡解，直到试探次数“MyCounter”大于限定最大试探次数的变量“MyNumTol”或求得的解的个数“MyMark”大于解的最大个数“NumOfSol”。

5.如图1所示，以收敛性为退出条件向“MyResult”记录的平衡解所在的参数点的邻域的上下边界扩展，利用全局变量“Result”记录所有扩展后的收敛的平衡解。

6.文献中的结论与本实施例的计算结果对比如图二所示。左边的曲线图分别是文献《Three-to-One Internal Resonances in Parametrically ExcitedHinged-Clamped Beams》中的Figure2和Figure3，α₁与σ₂的曲线图表示周期振动一阶分量的幅值与外激励频率的关系，α₂与σ₂的曲线图表示周期振动二阶分量的幅值与外激励频率的关系，α₁与F的曲线图表示周期振动一阶分量与外激励幅值的关系，α₂与F的曲线图表示周期振动二阶分量与外激励幅值的关系。右边曲线图分别与其左边曲线图对应，是本实施例的计算结果。由图二可以看出，本发明所述方法能计算出一维参数变量四维调制方程的平衡解,刻画出了二自由度振动系统周期振幅与一维参数变量的关系。

7.本实施例的一个计算结果如图三所示。图三是文献《Three-to-One InternalResonances in Parametrically Excited Hinged-Clamped Beams》中的Figure2和Figure3综合。图三给出了周期振动一阶分量的幅值α₁以及周期振动二阶分量的幅值α₂分别在二维参数变量σ₂和F中的取值情况。文献《Three-to-One InternalResonances in Parametrically Excited Hinged-Clamped Beams》没有给出与图三对应的情况。基于本发明人所阅读文献，没有文献给出过二维以及高于二维参数变量四维调制方程的平衡解。由图三可以看出，本发明所述方法能计算出二维参数变量四维调制方程的平衡解，并且能够刻画出含有任意可数维参数变量的任意自由度振动系统的周期振幅与参数变量的关系。

Claims (2)

1.一种刻画振动系统周期振幅与系统参数的关系的方法，其特征在于：本发明首先录入具体振动系统的状态变量的一阶近似解的调制方程的平衡解等价方程；为了控制程序的运行，本发明定义了一些计算所需要的全局变量：记录解的最大个数的变量，限定最大求解数的变量，存储计算结果的变量，记录平衡解数据大小的变量，判断相同解的容差变量，限定振动系统维数的变量，划分参数变量网格数的变量，限定参数变量范围的变量；本发明给出嵌套循环结构遍历参数变量的整个取值范围，此嵌套循环结构的每一层循环结构对应一维参数变量ε_i(i＝1...n)；在嵌套循环结构的最底层，记ε＝ε⁰，即每一维参数变量都取一个固定数值；在嵌套循环结构的最底层定义局部变量：记录所给参数点处平衡解的变量，记录所给参数点处平衡解个数的变量，记录试探次数的变量，限定最大试探次数的变量；在嵌套循环结构的最底层，再以试探次数不超过最大限定试探次数和求得的解的个数不超过限定最大求解数作为退出条件给出一个循环结构；在此循环结构的内部，在求解范围内随机给出初值X＝X⁰，并用迭代方法X＝X(X⁰)求出ξ使得调制方程的平衡解等价方程满足G(ξ,ε⁰)＝0；如果所求得的平衡解不与已经求解并记录的平衡解重复，利用局部变量记录此求得的平衡解；假设ξ在ε⁰处连续，即ξ+Δξ＝ξ(ε⁰+Δε)，以收敛性为退出条件向该平衡解所在的参数点的邻域的上下边界扩展，即令ε＝ε⁰+Δε，以X＝ξ为初值迭代求解方程(3),若能得到收敛的数值解ξ+Δξ，则继续向边界扩展，否则就停止扩展；利用全局变量记录所有扩展后的收敛的平衡解，画出参数变量与平衡解的关系图，刻画出振动系统周期振幅与系统参数的关系X＝X(ε)；

所定义的限定调制方程维数的变量，本发明所给出的方法可用于任意可数维状态变量的调制方程，即用于任意自由度的振动系统；

所述的嵌套循环结构，本发明所给出的方法可以用于含有任意可数维数参数变量的振动系统；

所定义的限定最大求解数的变量，本发明所给出的方法可以求出具有多个周期解的振动系统的所有周期解；

限定参数变量取值范围的变量，本发明所给出的方法可以计算参数变量的任意取值区域；

所述的假设平衡解在其参数点的邻域内连续，以收敛性为退出条件向该平衡解所在的参数点的邻域的上下边界扩展，该方法提高了程序的计算效率。

2.根据权利1所述一种刻画振动系统周期振幅与系统参数的关系方法，其特征在于：所述状态变量的一阶近似解的调制方程描述的是周期振幅与系统参数的关系；其中，状态变量可以是任意随时间往复变化的特征量，例如机械的振动、电磁的振荡、种群数目的变化、经济总量的起伏；系统参数是指可以人为调节的固有的系统性质，如惯性、弹性、外激励幅值及频率。