陳 寧, 張紅兵, 李萬祥
(蘭州交通大學 機電工程學院,蘭州 730070)
在非線性振動系統(tǒng)的定性研究中,龐加萊映射是極其重要的思想工具,可以將連續(xù)微分動力學系統(tǒng)轉(zhuǎn)化離散差分動力學系統(tǒng),使之深化對分叉行為和過程的研究。在非線性動力學研究中,定量分析是定性分析的基礎,定量分析的結(jié)果是定性分析中龐加萊映射得以構(gòu)建的基礎。但是在當前的研究中,對復雜的、難以解析或不可解析的振動系統(tǒng),定性研究的深度遠不及對其的定量分析。這是由于對這類系統(tǒng)無法得到精確的解析解,而在現(xiàn)行的定性分析方法中,一個可以被分析的龐加萊映射必須依賴精確的解析解,而非數(shù)值解。
非線性振動系統(tǒng)的定量分析方法可劃分為解析法和數(shù)值法兩個范疇,與之相對應的在定性分析過程中,理論上也具有依靠解析解和數(shù)值解兩種構(gòu)建龐加萊映射的基本方式。為便于下文闡述,將映射的可視化結(jié)果簡稱映像,將映射的迭代不動點簡稱為平衡點,分別將通過精確解析解、近似解析解和數(shù)值解構(gòu)架起來的龐加萊映射簡稱為解析映射、近似解析映射和數(shù)值映射。
龐加萊映射只是一個特殊的運算,從本質(zhì)上來說依靠精確的解都可以構(gòu)建起系統(tǒng)的龐加萊映射。Luo等[1-3]考慮了分段線性系統(tǒng)可以分段解析的特點,利用解析解以及周期運動的假設構(gòu)建系統(tǒng)的周期映射,深刻地揭示了含間隙振動系統(tǒng)的分叉過程和機理。其中得出擦邊可能是導致分叉奇異的發(fā)現(xiàn)和結(jié)論,這既是非線性動力學領域持續(xù)關注的熱點之一[4-6],也是機械振動系統(tǒng)中普遍需要解決的問題之一。但是大多數(shù)的系統(tǒng)都不能利用解析法進行定量分析,如文獻[7-9]中針對含庫倫摩擦振動系統(tǒng)做出的研究,雖然可以分段解析,但是系統(tǒng)中由于間隙和摩擦力的兩個相互不干涉的因素存在,無法完全以解析的方法構(gòu)建系統(tǒng)的龐加萊映射,必須采用分支程序以輔助完成對系統(tǒng)分叉的分析。類似這種結(jié)合解析解和數(shù)值方法的定量分析方法可以統(tǒng)稱為半解析法,是針對復雜分段線性系統(tǒng)切實有效的方法,如李萬祥等[10]就單自由度非光滑系統(tǒng)采用半解析法進行研究,得出單自由度碰撞振動系統(tǒng)也可以存在Hopf分叉的發(fā)現(xiàn)。
面對更加復雜的系統(tǒng),如分段光滑系統(tǒng)中,是難以分段精確解析的,陳洪月等[11]在對采煤機這類復雜設備的動力學響應研究中,就只能以數(shù)值模擬的方法求得系統(tǒng)的龐加萊截面圖??梢哉f,基于精確的解析解、半解析解和數(shù)值解都是可以構(gòu)建起系統(tǒng)的龐加萊映射。
對于解析映射可以方便地利用不動點處雅可比矩陣的特征值信息,對系統(tǒng)中發(fā)生的分叉等做出精確的判別,如文獻[12]中對內(nèi)伊馬克沙克-音叉分岔點處的兩參數(shù)開折問題的研究。但是對于數(shù)值映射,由于缺少后續(xù)對應的分析方法,對其研究只能止步于觀察數(shù)值映像的分叉圖和結(jié)合經(jīng)驗做出分析,然后利用經(jīng)驗對系統(tǒng)中的分叉行為等做出解釋,如文獻[13]中關于采煤機和復雜齒輪系統(tǒng)的動力學特性研究,就以數(shù)值模擬的方法結(jié)合龐加萊截面法求得系統(tǒng)的數(shù)值映像及分叉圖,但是并未做到文獻[14]中對分叉分析的深刻成度,僅依靠觀察分岔圖從而得出部分規(guī)律。上述現(xiàn)象反映出數(shù)值映射的可分析性不足,不像解析映射那樣可以對龐加萊映射平衡點做出分析,這也是數(shù)值映射不被廣大學者所公認的原因。
數(shù)值求解適合于所有振動系統(tǒng)的求解,這是解析分析所不具備的優(yōu)勢。雖然數(shù)值映射的概念尚未被認可,但是已有大部分學者在變相地使用,如文獻[15-19]中針對含摩擦系統(tǒng)、分段光滑系統(tǒng)、含間隙系統(tǒng)等復雜多體動力學系統(tǒng)的研究,通過大量龐加萊截面法來分析問題,其過程的本質(zhì)是龐加萊映射的數(shù)值化實現(xiàn)方法,其結(jié)果和解析映射的結(jié)果一樣為龐加萊映像圖。
數(shù)值計算的適用性更強,應用面更廣,但是數(shù)值映射在定性研究中發(fā)揮的作用卻相對更小,這其中也包含學者對數(shù)值解精確性有所懷疑的因素。大量的理論分析和試驗證明,只要嚴格控制算法穩(wěn)定性,必定可以保證系統(tǒng)的積分精度。鑒于其求解精度,數(shù)值解是目前大多數(shù)近似解析方法精度驗證的標準。如劉汝逾等[20-22]在解決非線性振動問題時,就以數(shù)值解為標準對各自所提出的近似解析法的求解精度進行驗證。大量拓展近似解析法的論文,都以數(shù)值解驗證近似解析解精度,這說明數(shù)值解具有比近似解析解更高的可靠性。
針對數(shù)值映射尚不能被分析研究的問題,部分學者嘗試采用近似解析解替代精確解析解建立系統(tǒng)的近似解析映射,但是目前還沒有實質(zhì)性的突破以便在實際中廣泛應用。雖然一個成熟的近似解析方法足以對系統(tǒng)做出有效的定量研究,但是函數(shù)映射的特點和近似解析解的本質(zhì)特征也決定利用近似解析映射可能導致定性分析結(jié)果失真的必然性。
因此,考慮數(shù)值映射的可靠性和廣泛的適用性,如果能夠?qū)?shù)值映射找到與之匹配的分析過程和方法,必定可以深化對不可解析甚至難以解析振動系統(tǒng)的定性分析。
分叉點處龐加萊映射迭代不動點的雅可比矩陣特征值是研究高余維分叉的關鍵,目前成熟的方法中缺乏對數(shù)值映射迭代不動點極其雅可比矩陣特征值的求解方法,這是導致對復雜系統(tǒng)定性分析深入度不足的主要原因,也將是本課題研究的主要工作。
系統(tǒng)龐加萊映射的迭代不動點是分叉行為分析的關鍵,因此利用數(shù)值映射的本質(zhì)特點尋找分叉點附近系統(tǒng)的迭代不動點是首要工作。
在系統(tǒng)分叉點的某個去心鄰域內(nèi)無論迭代過程是漸近收斂還是發(fā)散狀態(tài),其在時間進程中向穩(wěn)定狀態(tài)過渡的過程都十分緩慢。由于這種過程的緩慢性和穩(wěn)定性,通過自然迭代的方法求解不動點是不現(xiàn)實的,但是也為和其他方法的結(jié)合提供可能性和前提條件。
在此定義某種基于數(shù)值解建立的離散動力系統(tǒng),由映射Py描述,如式(1)所示。
ynyn+1=Py(μ,yn)yn
(1)
式中:yn+1為系統(tǒng)參數(shù)yn在映射Py作用下的象; 映射Py由原像yn和參數(shù)μ共同決定。
設映射Py的迭代不動點為Y,根據(jù)不動點的映射條件應該滿足如式(2)所示的代數(shù)方程,則Y同樣的和參數(shù)μ相關。
Y=Py(μ,Y)Y
(2)
參數(shù)yn可以表示為關于不動點Y和不動點的誤差εn之間的關系,如式(3)所示。將式(3)代入式(1),并將Py(μ,yn)按Taylor級數(shù)在不動點Y處展開,代入式(2)可以重構(gòu)與映射Py等價的映射Pε如式(4)所示, 映射Pε可以描述誤差參數(shù)εn和εn+1之間的映射關系,如式(5)所示。
yn=Y+εn
(3)
(4)
εnεn+1=Pε(μ,εn)εn
(5)
式中:εn為參數(shù)yn相對于不動點Y的誤差; 原點O即映射Pε的不動點,滿足式(6)所示的平衡關系。
O=Pε(μ,O)O
(6)
當系統(tǒng)的不動點穩(wěn)定時,映射Pε在原點O處的譜半徑ρ[Pε(μ,O)]<1。假設系統(tǒng)在參數(shù)μ=μ0處發(fā)生某種基于不動點的經(jīng)典分叉,此時必然有ρ[Pε(μ0,O)]=1,即存在Pε的雅可比矩陣的特征值λ滿足如式(7)所示的形式,和分叉參數(shù)μ0有關。
λj(μ0)=cosθ±sinθi∈diag[Pε(μ0,O)]
(7)
利用映射關系式(5)分析擾動εn的傳遞情況
(8)
假設存在轉(zhuǎn)化矩陣Qε(μ0,O)使Pε(μ0,O)可以對角化為Λε(μ0,O),滿足式(9)的形式。
Qε(μ0,O)εn+1=Λε(μ0,O)Qε(μ0,O)εn
(9)
令Qε(μ0,O)εn=Zn,則映射Λε可以表示誤差εn正則化后的誤差參數(shù)Zn和Zn+1之間的映射關系,如式(10)所示。
Zn+1=Λε(μ0,O)Zn
(10)
當條件式(7)成立時,初始擾動εn在某個空間方向或者空間曲面內(nèi)臨界穩(wěn)定,既不會收斂,也不會發(fā)散,保持相對初值的穩(wěn)定性,則式(10)在迭代過程中必然滿足式(11)的關系。
‖Zn+1‖=‖Zn‖
(11)
但是在一般情況下,由于計算的精度以及其他原因,只能夠確定分叉參數(shù)所在的某個鄰域,即僅僅可以找到某個μ∈U(μ0),可以表示為
μ=μ0+μΔ
(12)
則特征值式(7)可表示為含有誤差的形式
λj(μ)=[cos(θ+θΔ)±sin(θ+θΔ)i][1+Δλ]∈
diag[Pε(μ,εn)]
(13)
式中: Δλ為由于參數(shù)μ的誤差引起的譜半徑的變化,Δλ=ρ[Pε(μ,εn)]-1;θΔ為特征值輻角的誤差; 此時若有μΔ→0, 則Δλ→0,(1+Δλ)→1,且使條件式(14)成立。
‖Zn+m-Zn‖=
‖Zn‖[1+Δλ]m-‖Zn‖?‖Zn‖
(14)
說明在有限次數(shù)的迭代之后,由于系統(tǒng)的發(fā)散和收斂造成的攝動的變化量要遠遠小于初始攝動量。在系統(tǒng)的某個維度面上,如系統(tǒng)的正則范式所對應的基平面,狀態(tài)參數(shù)的某些分量繞特定的點發(fā)生旋轉(zhuǎn),此時其平衡點近似在旋轉(zhuǎn)中心位置所在,但是由于其衰減指數(shù)十分接近1,故而其旋轉(zhuǎn)半徑的攝動量十分小。利用映射關系式(5)可知即便是收斂情況,通過持續(xù)迭代求取迭代不動點也是極其困難的,當?shù)l(fā)散時,必然是一種不可取的方案。
在自然迭代過程中,如果初始狀態(tài)的攝動量‖Zn‖變化速度較為緩慢,則在有限的時間內(nèi),‖Zn‖的值會保持相對的穩(wěn)定,可以借助這種相對穩(wěn)定特征來加速迭代收斂的過程。
在此,可以借助分叉點處映射的特征進行迭代不動點的求解,如圖1所示。圖1中:O為其理論平衡點在系統(tǒng)特定低維度子空間中的投影;Zn,Zn+1,Zn+2為連續(xù)相鄰三次狀態(tài)構(gòu)成的龐加萊映像,近似滿足式(15)、式(16)的關系;Oμ為其外接圓心。
(15)
∠ZnOZn+1=∠Zn+1OZn+2=θ+θΔ
(16)
可設龐加萊映像中狀態(tài)參數(shù)的坐標如式(17)所示。
Zn[‖εn‖cos 0,‖εn‖sin 0],
Zn+1[‖εn+1‖cos(θ+θΔ),‖εn+1‖sin(θ+θΔ)],
Zn+2[‖εn+2‖cos(2θ+2θΔ),
‖εn+2‖sin(2θ+2θΔ)]
(17)
利用幾何關系可求得外接圓圓心Oμ坐標
(18)
計算的外接圓心的理論攝動半徑和初始狀態(tài)攝動半徑之比如式(19)所示。
(19)
當sinθ→0時,代入式(19)可得
(20)
由式(20)可以看出,當|sinθ|>|Δλ|時,經(jīng)過該方法計算后的攝動量是絕對小于自然迭代過程的,對平衡點失穩(wěn)的狀態(tài)同樣適用。初始狀態(tài)Zn經(jīng)過兩次自然迭代后的攝動半徑為‖Zn+2‖=‖Zn‖(1+Δλ)2,和自然迭代法相比,其收斂速度之比如式(21)所示。
(21)
當sinθ=0時,
(22)
同樣滿足當|sinθ|>|Δλ|時,該迭代計算過程是絕對收斂的,且效率絕對優(yōu)越于自然迭代。
式(22)也反映出對于sinθ=0的情況該算法效果不太理想,極有可能導致算法失穩(wěn)失效。對于這種情況,可以給出算法失效的判決方式和替代計算方法。
在此以運算前的初始狀態(tài)和運算后的初始狀態(tài),分別經(jīng)過m次的自然迭代,則可以輕易求得以原初始狀態(tài)Zn為起點的迭代序列有限集合s1,以及以運算后的結(jié)果Zoμ為初始條件的迭代序列s2。
s1: {Zn,Zn+1,Zn+2,Zn+3,…,Zn+m};
s2: {Zoμ,Zoμ+1,Zoμ+2,Zoμ+3,…,Zoμ+m},Zoμ=Oμ。
其集合S1和S2的構(gòu)建以外接矩形為例,如圖2所示。圖2和圖3中:黑點表示有限的龐加萊映像;矩形
表示按照特定方向最小外接矩形構(gòu)造的凸集;圓點表示以原初始條件映射序列的有限元素集合s1;三角形表示運算后的值為初始狀態(tài)的映射序列有限元素集合s2;菱形表示最小凸集的幾何中心。
以迭代計算前后的結(jié)果為起始點所構(gòu)造的兩個凸集的集合關系,如圖3所示。此時:若有S2∈S1,則算法必然穩(wěn)定;若有S1?S2,則可以判定為算法失穩(wěn)。當結(jié)果如圖3(a)所示時,運算可以連續(xù)壓縮平衡點所在區(qū)間,多次迭代后可得到有效的近似平衡點;當集合包含關系如圖3(b)所示時,算法處于失穩(wěn)狀態(tài);圖3(c)和圖3(d)所示情況說明,已經(jīng)達到計算機的或者程序求解的精度極限,或者原初始狀態(tài)已經(jīng)超出原平衡點的吸引域,并呈現(xiàn)指數(shù)發(fā)散狀態(tài)。對高維度的系統(tǒng),可以利用迭代映像在不同相平面上的投影來構(gòu)建所需要的兩個凸集,并結(jié)合圖3所對應的集合關系判斷算法的有效性。
基于上述的特點,無論龐加萊映射發(fā)散與否只要轉(zhuǎn)換矩陣的特征值不存在大于1的正實數(shù),可以在判斷原算法失效的條件下,改用多次迭代的凸集幾何中心或者若干次龐加萊映射的重心替代原算法中的圓心Oμ,作為迭代運算的結(jié)果,如式(23)或式(24)所示。
(23)
(24)
當離散樣本數(shù)量m較大時,理論上對各種基于周期運動發(fā)生的分叉都具有較好的收斂性和較為普遍的適應性。
對平衡點的分析主要分析其雅可比矩陣的特征值,確定近似平衡點特征值也是十分重要的工作。可以采用第1章提供的迭代方法求得在系統(tǒng)Λε中的近似平衡點Zn。
此時,可以保證‖Zn‖足夠小,假設存在若干初始狀態(tài)Zkn屬于列向量組Zsn,Zsn及其Zkn滿足式(25)所示的條件。
(25)
式中,E0=e0I,e0為對每個初始狀態(tài)Zn預設的初始誤差量,I為單位矩陣,且有rank(I)=rank(Zsn)。
對初始狀態(tài)組Zsn和Zn取相同次數(shù)的龐加萊映射,可以得到關于近似不動點和初始狀態(tài)組的映射結(jié)果,Zsn+m及Zn+m,誤差關系可表示為
(26)
式中,Em為m次復合映射后的誤差矩陣,根據(jù)離散系統(tǒng)的穩(wěn)定性理論,可以利用式(27)確定雅克比矩陣的所有特征值。
(27)
方程γm=1在復數(shù)域內(nèi)具有m個解,其解集γ的表達式為
(28)
假設近似平衡點Zn處的特征值λ的真解如式(29)所示。
λ=ρcosΦ±jρsinΦ
(29)
式中:ρ為特征值的模長,ρ=‖λ‖;Φ為特征值在復平面內(nèi)的角度。根據(jù)復數(shù)的乘積運算,特征值的m次冪可以表示為
λm=ρm{cos[mod(mΦ,2π)]+
jsin[mod(mΦ,2π)]}
(30)
通過對復數(shù)取根式可以得出的模接近真值,但是其角度信息并不一定正確。其中,當m=1時可以近似確定輻角位置,但是當m=1時的精確性較低,可以利用復數(shù)的指數(shù)運算規(guī)律,通過迭代法,分別運算m=[1,2,3,…]時的情況,利用漸近性的原理逐步逼近系統(tǒng)雅可比矩陣特征值的輻角真解,根據(jù)式(31)的迭代過程算法可以求得特征值輻角逼近解。
(31)
式中:Φm為m次誤差迭代對應的特征值輻角;Φ(m)為逐步逼近的不動點雅可比舉證的特征值輻角。
以下以一種含間隙兩自由度彈簧非線性柔性碰撞振動系統(tǒng)在含三倍諧波的情況下的動力學研究分析為例,對其中某個類似Flip-NS型余維二的分叉過程進行研究。
選取算例的力學模型如圖4所示。圖4中:M1和M2為振子的質(zhì)量;C1,C2和C為阻尼器的阻尼系數(shù);K1,K2和K為彈簧的剛度系數(shù);D為靜止狀態(tài)時的振子間的間隙;g為振子之間的實時動態(tài)間隙,負值表示壓縮量;F為外激勵;FD和-FD為由于彈簧的含間隙接觸緣故造成的作用與反作用力。由此可以建立系統(tǒng)的微分方程式(32)。
(32)
式中:彈簧剛度系數(shù)Kn具有非線性表達式(33);外激勵F含有高次諧波成分如式(34)所示;FD表達式如式(35)所示。
(33)
式中,χ為非線性剛度比例系數(shù)。
(34)
式中:P為無量綱的比例系數(shù);Pn為諧波強度;φn為n次諧波的相位。
(35)
由于間隙造成系統(tǒng)中含有一個非光滑切變界面ΣC,如式(36)所示,在其兩側(cè)系統(tǒng)的微分方程不連續(xù)。為保證系統(tǒng)狀態(tài)穿越界面ΣC的必然性,參考線性系統(tǒng)對間隙D做出約束,在此將彈簧力學性能曲線進行近似線性替代,并假設外激勵只存在主諧波,即理想化假設K1,2,3(Δ)=kn,F(xiàn)=P1cos(ω1t+φ1)。
ΣC:X1-X2=D
(36)
假設系統(tǒng)位移的穩(wěn)定待定解如式(37)所示。
(37)
式中,A1,A2,B1和B1為解的待定系數(shù),將其代入方程并分離三角函數(shù)的相關項可得系數(shù)的平衡表達式(38)。
(38)
其中,
在線性近似條件下系統(tǒng)振動時可穿越切變截面的充要條件如式(39)所示。
(A1-A2)2+(B1-B2)2>D2
(39)
根據(jù)表達式(38)并定義間隙系數(shù)δ表征靜態(tài)間隙量D,滿足如式(40)所示的關系。
(40)
取模型參數(shù):M1=50;M2=4.5;k1=500;k2=1 244;k=2 244;χn=500;C1=0.3;C2=20.03;C=5.4;δ=0.007;ω=15.23;P1=0.128 0;P2=0;P3=0.154 4;φ1=0.24;φ2=0;φ3=0.63。
以RK-4數(shù)值積分法進行仿真,采用穩(wěn)定的二分法進行求解穿過非光滑轉(zhuǎn)折界面的時間間斷點,求解精度取1×10-13。以系統(tǒng)主諧波外激勵周期定義龐加萊截面Σ如式(41)所示,在數(shù)值積分的過程中,利用龐加萊截面法構(gòu)建離散映射。
Σ:mod(ωt,2π)=φ1
(41)
通過數(shù)值仿真可以得到如圖5所示的系統(tǒng)分叉圖,具有雙穩(wěn)態(tài)的特征。圖5中:灰色分支的分叉過程比較明確,是單周期運動失穩(wěn)以激變的方式演化為擬周期運動,在擬周期區(qū)間發(fā)生的鎖相等行為也極其明確;黑色分支在P=12.5附近出現(xiàn)從擬周期演化為二倍周期運動的分叉過程,從分叉過程推斷,該分叉具有Flip-NS分叉的趨勢,但是由于本系統(tǒng)得不到系統(tǒng)的解析映射,故而無法利用解析映射的特征值來定義該分叉過程的類型。
利用第1章提供的方法,可求得圖5中所對應兩支分叉過程的平衡點變化規(guī)律,如圖6所示。通過對比圖5和圖6可以看出,第1章所提供的方法對基于單周期運動的分叉具有良好的適應性,平衡點的激變過程和圖5中的激變過程相一致。圖中求解結(jié)果顯示,對基于多周期運動的分叉過程適應性不足,這是由于方法機制所導致,對于多周期運動,可以通過構(gòu)建復合映射的方式將其轉(zhuǎn)化為單周期運動,以便于研究。但是對于混沌運動,由于混沌運動的瀝遍性,該方法已經(jīng)失效。
同時可以看出,對于圖5中存在疑點的分叉區(qū)間,圖6中對平衡點的求解是收斂且穩(wěn)定的,可以繼續(xù)利用第2章提出的迭代方法求得系統(tǒng)雅可比矩陣特征值的近似值。對平衡點的求解精度取1×10-13,然后對該近似平衡點的特征值進行37次迭代運算,求得系統(tǒng)的雅可比矩陣特征值變化如圖7所示。
圖7中,當P=12.5時,不動點的雅可比矩陣其中一對共軛特征值處于單位圓外部,具有明顯擬周期運動的特點,和理論分析的結(jié)果相一致。隨著參數(shù)的遞增,其雅可比矩陣特征值的虛部分量逐漸減小為0,直接導致系統(tǒng)的龐加萊映像由擬周期狀態(tài)遷變?yōu)槎吨芷跔顟B(tài)??梢钥闯?,始終僅僅有一對特征值在單位圓附近變化,并不符合發(fā)生Flip-NS組合分叉時兩對特征值同時穿越單位圓的特點。
上述分叉過程從分叉圖的變化規(guī)律分析,基本符合Flip-NS分叉的特征,但是從圖7內(nèi)特征值的變化規(guī)律中可以看出,其中一對復共軛特征值在單位圓外側(cè)退化為重實特征值,然后其中一個特征值從-1處穿越單位圓,該變化規(guī)律符合Hopf分叉周期二鎖相狀態(tài)的特征,其演化過程如圖8所示。
本文總結(jié)和分析了振動系統(tǒng)中不能利用數(shù)值解進行定性分析的現(xiàn)象及其原因。本質(zhì)原因是現(xiàn)行的定性分析方法對解析解的依賴;關鍵原因是無法利用數(shù)值解求得龐加萊映射的平衡點。文中對上述問題做出了針對性研究,總結(jié)了基于數(shù)值解求解平衡點及其龐加萊映射特征值的方法,并針對一個存在疑點的分叉行為成功地進行了分析。
文中提到的方法在算例模型上的實踐被證明是可行有效的,非線性系統(tǒng)的數(shù)值解和非線性定性理論之間可以實現(xiàn)對接,利用差分和迭代原理基于數(shù)值解可以完成非線性系統(tǒng)的定性分析,擺脫在定性分析時對解析解的依賴。