外文翻譯--流體動壓軸承-撓性轉子系統(tǒng)的非線性動態(tài)特性 中文版

1、<p>　　流體動壓軸承-撓性轉子系統(tǒng)的非線性動態(tài)特性</p><p><b>　　呂延軍 虞烈 劉恒</b></p><p>　　(西安交通大學潤滑理論及軸承研究所 中國西安 710049)</p><p>　　摘要：分析了帶液壓軸承座的撓性轉子系統(tǒng)的非線性動力性。軸通過使用能考慮物體的慣性和剪切力影響的有限元法仿效出來。根據(jù)非

2、線性流體動壓軸承-轉子系統(tǒng)，一種被修改的帶有自由接口的模態(tài)綜合技術是被用來降低柔性轉子系統(tǒng)的模型的自由度。根據(jù)油膜的物理特性，改變約束的方法是通過引入持續(xù)改變的每一階段的動態(tài)的集成與重復的雷諾茲方程式。采用等周圖形學的有限元方法解決雷諾茲邊值的液體潤滑問題而沒有增加計算。非線性油膜力及其Jacobians矩陣的數(shù)值解必須要具有協(xié)調(diào)一致的精度。周期通過使用Poincare-Newton-Floquet（PNF）方法而獲得，一種方法，把軌跡

3、預測追蹤的延續(xù)算法和PNF方法結合起來提出計算周期運動的分岔點是由于系統(tǒng)參數(shù)的改變而受到影響。局部的穩(wěn)定和周期運動的分岔現(xiàn)象是通過Floquet 理論獲得的。轉子系統(tǒng)的混亂運轉需進行能量譜的檢查。許多的例子顯示這項研究能節(jié)省計算量而且具有很好的精密度。</p><p>　　關鍵詞：非線性動力學 軸承轉子系統(tǒng) 穩(wěn)定性 分岔 混亂 有限元方法</p><p><b>

4、　　前言</b></p><p>　　旋轉機已經(jīng)應用于能量站、飛行器、機床夾具、汽車、以及家庭應用等各個領域。軸承轉子系統(tǒng)，就像一種轉子機器一樣，是一種典型非線性機械系統(tǒng)。這種非線性分析方法已應用于非線性的軸承轉子系統(tǒng)上就如同線性分析方法不能應用于分析的那樣。很多種方法成為重要的合適于分析的多自由度軸承轉子系統(tǒng)的穩(wěn)定性和分岔問題，因為此種方法不能提高工作效率。</p><p>

5、　　目前世界上有許多研究者致力于非線性動態(tài)轉子軸承領域的研究。不幸的是他們塑造的轉子軸承有一些不利條件被直接利用于指導那些力學產(chǎn)品的設計。因為非線性分析的復雜性，軸承轉子系統(tǒng)的非線性模型通常被當作有很少自由度的和分析形式上軸承力的非線性模型。例如，一個勻稱的硬的轉子[1，2]或者Jefrcott轉子模型[3，7]，多項式模型[8，9]，以及或長或短的軸承模型，都不能準確地描述實際的系統(tǒng)。在這些研究過程中的軸承的非線性的油膜力的分析形式而

6、軸承轉子系統(tǒng)不能在實踐中得以分析。然而，軸承轉子在本質上是非線性的，轉子支撐的非線性運動是由軸承引起的。非線性的油膜力按照轉子的幾個交點運動。軸承轉子系統(tǒng)的局部非線性和組分是連接在一起的。因此，來自 于非線性的影響是全球的。</p><p>　　Ref.[10]描述一個局部非線性的高位動力系統(tǒng)被運用于模態(tài)削減方法上，基于模態(tài)的固定接口的綜合技術，在Ref.[11]中提出了一種為解決局部非線性的動力系統(tǒng)的周期性和穩(wěn)

7、定性的方法。在這種理論上，軸被描述為多自由度有限元素使用了2個結點Timoshenko 梁軸有限元模型。根據(jù)非線性的轉子系統(tǒng)，一種被修改過的帶自由接口的模態(tài)綜合技術是被用來降低柔性轉子系統(tǒng)的自由度的有限元模型。在減少以后，系統(tǒng)仍舊保留它的非線性和保存系統(tǒng)的非線性分析力。根據(jù)改變約束方法的八個交點的等周圖形學的有限元方法是用于解決液體潤滑發(fā)生在雷諾茲邊界的橢圓不等式。一個擾動等式能夠得到直接的有限元等式。因此，非線性油膜力及其Jacobi

8、ans矩陣的數(shù)值解必須要具有協(xié)調(diào)一致的精度。這樣不能引起軸的摩擦和應力集中。采用PNF方法，F(xiàn)loquet 理論和軌跡預測追蹤的延續(xù)算法研究了不平衡響應，T周期運動和隨軸承系統(tǒng)設計參數(shù)的改變的分岔現(xiàn)象。軸承轉子系統(tǒng)的這種混亂狀態(tài)是從能量譜中調(diào)查出來的。</p><p><b>　　1 動力系統(tǒng)方程式</b></p><p>　　圖1所示的轉子-軸承系統(tǒng)是一個典型的非線

9、性動力系統(tǒng)。該系統(tǒng)由線性部分(柔性轉軸)和局部非線性部分(徑向軸承)組成。圖2所示的2節(jié)點具有8自由度的Timoshenko梁軸單元模型，由于其可以計及轉動慣量與剪切變形的影響，故更接近實際運行的轉子系統(tǒng)。因此采用有限元方法建立如下形式的柔性轉軸橫向振動方程</p><p><b>　　（1）</b></p><p>　　式中分別為轉軸的質量矩陣、陀螺矩陣、剛度矩陣、

10、周期為 的外力向量(包括重力和不平衡力)和軸承施加于轉軸的非線性力向量。對于具有P個節(jié)點的轉軸，節(jié)點位移向量可表示為</p><p>　　式中 分別為第j個節(jié)點沿水平和鉛垂方向的橫向位移與彎曲轉角。非線性力向量可表示為</p><p>　　式中Fxj，F(xiàn)yj分別為軸承作用在軸第j個節(jié)點上的水平和鉛錘方向的油膜力。由于軸承的非線性油膜力孤立地作用于轉子的個別節(jié)點上，因此對具有m個軸承支承的轉

11、子系統(tǒng)，軸承力</p><p><b>　　具有如下局部性質</b></p><p>　　式中xsb∈R4m, Fsb(xsb,xsb)∈R4m可被寫成：</p><p>　　圖1 轉子-軸承系統(tǒng)示意圖</p><p>　　圖2 轉軸有限單元模型</p><p>　　為了簡化符號，將式(1)中各元

12、素重新排序且表示為如下分塊形式</p><p>　　由于需要花很多的時間計算多自由度的轉子系統(tǒng)，在維持系統(tǒng)響應準確的情況下，減少系統(tǒng)的自由度是非常重要的。由于系統(tǒng)是局部非線性的，僅非線性自由度受控于非線性方程，而線性自由度又依賴于非線性自由度，因此可對線性自由度進行減縮，以使該系統(tǒng)降階。為了避免縮減自由度時，坐標轉換給系統(tǒng)的非線性因素帶來的數(shù)值誤差，僅將線性自由度轉換為模態(tài)坐標，而非線性自由度和決定系統(tǒng)動力特性的

13、非線性力仍保留在物理空間中，使降階后的系統(tǒng)仍具有局部非線性特征。</p><p>　　為了降低線性組合的自由度，將XS表示為列的線性組合</p><p>　　式中 </p><p>　　因此，矩陣的保留彈性特征模態(tài)的列是ωk∈(0，ωcut)的無阻尼特征值問題(ks-ω2jms)ψj=0(j=1,

14、…，nk)的質量正則解。矩陣的剩余柔性模態(tài)的列可表示為</p><p>　　其中對角矩陣是角頻率小于或等于時的譜矩陣。因此</p><p>　　從（11）式開始，可被寫成：</p><p>　　這樣就有如下整個變換</p><p>　　在以上等式中，矩陣變換T=T1T2T2，運用式(13)，可得縮減的系統(tǒng)方程</p><p

15、>　　通過縮減把n(n=nb +nc)階方程組減縮為s(s=nb +nc)階方程組，由式(11)- (I4)，可見轉軸的不平衡力及非線性項的影響全部保留在縮減的方程組(14)中?？紤]到圓盤的不平衡力的影響，可得動力系統(tǒng)的運動方程</p><p><b>　　有 </b></p><p>　　式中分md，Gd, Kd, Fdex別為圓盤的質量矩陣、阻

16、尼矩陣、剛度矩陣、不平衡力向量。</p><p>　　當引入狀態(tài)變量后，其相應的系統(tǒng)方程在狀態(tài)空間中為</p><p>　　2.非線性力的計算及其流體動壓軸承的Jacobians 矩陣</p><p>　　對于實際軸承，不具有解析形式的油膜力，而在計算系統(tǒng)非線性響應時，每一時刻動力積分均需要非線性力及其Jacobians矩陣的求解。非線性油膜力及Jacobians矩

17、陣協(xié)調(diào)一致的精度不僅影響到求周期解的PNF是否收斂，而且對周期解的穩(wěn)定性及分岔的分析有著極其重要的影響。同時任一時刻油膜力的Jacobians矩陣的準確性又影響著判斷周期解穩(wěn)定性的Floquet乘法的求解?；谝陨蠁栴}，運用有限元法求解具有變分不等式形式的流體潤滑Reynolds邊值條件(油膜區(qū)域下游邊值條件)問題。將油膜力視為某時刻軸頸中心位移及速度的函數(shù)，</p><p>　　由此可以得到一組微分方程，根據(jù)該

18、方程組的特點，在求出油膜力的同時，可很快求得Jacobians矩陣。</p><p>　　對于有限長軸承流體潤滑的Reynolds邊值問題:</p><p>　　式中 — —油膜壓力函數(shù)(表壓)</p><p>　　— — 潤滑油動力粘度</p><p>　　——軸承長徑比的倒數(shù)</p><p><b>

19、;　　— — 油膜厚度</b></p><p>　　— — 從Y軸負方向到油膜位置的角度</p><p><b>　　— — 偏位角</b></p><p>　　— — 偏位線與軸承中心連線至油膜位置的角度，如圖3和圖4所示</p><p>　　式(18)可等價于如下離散的橢圓型變分不等式[18]</p

20、><p>　　式中Φ(p,q)= 是H(Ω)× H(Ω) 上的強制、對稱，橢圓雙線性泛函,ψ(q)=,k={p∈H(Ω)，p≥0 in Ω}是線性泛函，H(Ω) )是Sobolev空間。是唯一的一層油，是的微分。是偏位線與軸承中心連線至油膜完整區(qū)和油膜破裂區(qū)交界線（此交界線是隨位移及速度擾動變化的曲線)的夾角。由以上等式，油膜厚度h和變角可寫成如下方程式：</p><p>　　圖3

21、單塊瓦計算坐標 圖4 橢圓軸承及其計算坐標</p><p>　　由于非線性油膜力是軸頸中心動態(tài)位置的函數(shù)，所以油膜力可被寫成如下等式：</p><p>　　油膜力Fx ,Fy過對動態(tài)壓力分布分別積分后得到:</p><p>　　運用8節(jié)點等參有限元法可求出油膜區(qū)域各節(jié)點的壓力分布Pi。壓力函數(shù)P可表示為：</p><p>　

22、　式中 Li——有限元插值函數(shù)，把式(23)代入不等式(19)可得如下n階離散不等</p><p><b>　　式方程：</b></p><p><b>　　式中 </b></p><p>　　為了求解等式（24），組成矩陣和矢量是必要的。把有限插值函數(shù)Li，Lj代入，得。然后等式（24）等價為約束迭代方程</p

23、><p>　　式中矩陣 和矢量分別為滿足橢圓邊值問題第一類及第二類約束條件的稀疏、帶狀、對稱矩陣和列矢量。解出P后，將式(23)代入式(22)可得：</p><p><b>　　式中 </b></p><p>　　均為常數(shù)列矢量。油膜力， 相對于x ，y 和 ，的Jacobians矩陣可表示為:</p><p><

24、b>　　式中 </b></p><p>　　將式(25)分別對x，y，，求偏導數(shù)可得如下擾動方程</p><p>　　將矩陣(k= x，y， ，)和矢量(x, y， ， )組合在一起可以得到如下等式：</p><p>　　把PK (k=x，Y， ，) 代入式(28),求得油膜力的Jacobians矩陣。由于式(31)與式(27)具有類似的形式，

25、不需要再次求解。所以在求解PK (k=x，Y， ，)時可由式(25)、(27)和式(30)的PK很快得到式(29)，而式(29)與式(25)具有相同的系數(shù)矩陣。因為式(29)與式(25)具有相同的系數(shù)矩陣，所以這樣就使得非線性分析所需要的油膜力和油膜Jacobian矩陣能夠同時計算完成。然而,為了保證非線性分析的準確性和可靠性，花費在計算每一步重復動態(tài)積分的Jacobians矩陣與油膜力自己本身的重復動態(tài)積分差不多。</p>

26、<p>　　3.非線性周期響應及其解法</p><p>　　運轉著的軸承轉子系統(tǒng)總是引起擾動。假定系統(tǒng)外在裝載根據(jù)周期T可得：</p><p>　　系統(tǒng)的穩(wěn)定現(xiàn)象，i，e 達到最大之后瞬間就會減弱，這種現(xiàn)象可能是由周期、半周期或是混亂引起的。系統(tǒng)的周期性問題可能在某些時間間隔變得不穩(wěn)定與系統(tǒng)的參數(shù)有關。例如，轉子的角速度，軸承的直徑與寬度比d/B，距離半徑比，凹槽周長與寬度比

27、（位于墊片上的凹槽周長），墊片圓弧，橢圓度，質心偏差 ，軸承的非線性是由不穩(wěn)定引起的，一般說來，系統(tǒng)的最大響應增加且線性系統(tǒng)的響應頻率也將作出響應。在實踐中，這將可能導致不利的軸摩擦。實際應用上，最重要的是決定這些不穩(wěn)定性間隔的設計和導致這些間隔的現(xiàn)象：周期、半周期或混亂。</p><p>　　3.1 PNF 方法</p><p>　　采用PNF(Poincare-Newton—Floqu

28、et)方法即求解兩點邊值問題的打靶法并結合Floquet分岔理論的周期解求解及穩(wěn)定性分析的方法，在選定的Poinca截面上，給出周期解不動點的初始迭代值Xs(t0)，結合方程（34）通過數(shù)字時間綜合化方法從一個周期T到另一個周期T間發(fā)現(xiàn)的問題，如果能滿足如下方程：</p><p>　　當系統(tǒng)控制參數(shù)u= 時，用Newton—Raphson迭代方法可求得 。雅可比矩陣:</p><p>　　

29、式中的 J= 可對系統(tǒng)方程式(17) 以(X5(t0),I)</p><p>　　為初始值積分一個Poincare映射周期而得到，其中δS(t0)=I， δS(t0+T)=J，式(35)中F=F(t，u，X )= 。 </p><p>　　3.2 周期解預測追蹤的延續(xù)算法和PNF方法</p><p>　　將延續(xù)算法和PNF法相結合的方法稱為連續(xù)周期方法，即計算

30、整個周期解的分支。當</p><p>　　為某一確定值時，就可以得到。連續(xù)周期方法就是通過觀察改變u在一周期內(nèi)的變化。開始從一個已知的，，得到周期后n步的方程式：</p><p>　　然后運用打靶法校正系統(tǒng)參數(shù)時的周期解。式中 可由系統(tǒng)方程(17)關于軌跡Xs(t0 +T)和</p><p><b>　　式中。</b></p>&

31、lt;p><b>　　4 應用算例及結果</b></p><p>　　帶有圓盤剛度（點D1）和橢圓軸承（點B1和點B2，墊弧為150o，潤滑油動力粘度取為0.0287Pa.s，寬度與直徑比B/d=0.8）的軸承轉子系統(tǒng)的不平衡響應已被許多例子分析了（如圖1和圖4）。軸（直徑為0.45m，長8.4m軸向力為2x1011Pa，剪切力為7.6923x1010Pa，質量密度為7800Kg/m3

32、）等同于引入六個有限元素。因此，系統(tǒng)的有限元素模型有七個點（28個自由度）和非線性軸承的兩個支撐。一個支撐的軸元素有四個自由度，橢圓軸承使用B1，B2兩個支撐點。在這種情況下，系統(tǒng)的四個自由度直接影響非線性力。因此轉子系統(tǒng)是一種典型的非線性動力系統(tǒng)。軸的質量偏心（）和圓盤的質量偏心（）有著相同的轉動角。有許多的例子采用了八固有模式方法。模態(tài)的影響使分析結果的準確性降低了。由于ψ=O.003， μ=0.028 7 Pa·s， =

33、0.556，B’／B=O，ω=l000 r／min,系統(tǒng)的周期問題分別用八固有模式方法和全固有模式方法解決。在選定的Poinca截面上，給出周期解不動點的初始迭代值Xs(t0)，用PNF方法可求得。Floquet乘子模的最大值為|fmax|=O.938</p><p>　　表 周期計算用PNF方法與不同算法（八固有模式方法或全固有模式方法）的對比</p><p>　　ω=1 000 r／

34、rim，δ=0.556，ψ=O.003，B’／B=0</p><p>　　運用PNF方法速度集中是非常迅速的。當使用八固有模式時，由于ψ=O.003，μ= 0.028 7 Pa·s，6=O．556，B'／B=O．4，0<ω<1 526 r／min（i．e．ω =1 526 r／min，|fmax|=1.0通過使用周期連續(xù)解決方法），因此周期是穩(wěn)定的。當環(huán)行值的注意時，在B1，D1點的

35、不平衡響應是一個半周期。此時，ψ=0.003, μ=0.0287 Pa. δ=0.556，B’／B=0.4，ω=1 580 r／min，i．e．選定Poincare截面投影在x-y軸的如圖5所示。當ψ=0.003，μ=0.0287 Pa·s，δ=0.556，</p><p>　　B’／B=0，ω=1 600 r／min時，中心轉子B1和D1點的軌跡是不穩(wěn)定的如圖6所示。當ψ=0.003，μ=0.0287

36、 Pa·s，δ=0.556，B’／B=0，,</p><p>　　ω=1 800 r（轉子的角速度ω=6Oπ rad／s），B1和D1點的混亂運動軌跡如圖7所示。圖8表示連續(xù)的時間yD1和同步能量譜處于混亂狀態(tài)。能量譜通常用于低能量級的測量—混亂的一個重要特征。</p><p><b>　?。╝）B1點的軌跡</b></p><p>

37、　　B1點在Poincare截面的軌跡</p><p>　　圖5 當ω=1 580 r／min，δ=0.556，ψ=0.003，B’／B=0.4時，中心轉子在B1點的半周期運動軌跡以及其在Poincare截面投影于x-y 軸的軌跡圖</p><p><b>　　（a）B1點的軌跡</b></p><p><b>　?。╞）D1點的軌

38、跡</b></p><p>　　圖6 當ω=1 600 r／min，δ=0.556，ψ=0.003，B’／B=0時，中心轉子在B1，D1點的軌跡</p><p><b>　?。╝）B1點的軌跡</b></p><p><b>　?。╞）D1點的軌跡</b></p><p>　　圖7

39、當ω=1 800 r／min，δ=0.556，ψ=0.003，B’／B=0時，中心轉子在B1，D1點的混亂運動軌跡</p><p>　　D1在y點的連續(xù)時間</p><p>　?。╞）D1點在y軸的能量譜</p><p>　　圖8當ω=1 800 r／min，δ=0.556，ψ=0.003，B’／B=0時，D1在y連續(xù)的時間圖和同步能量譜</p>&l

40、t;p><b>　　結論</b></p><p>　　軸被描述為多自由度有限元素使用了2個結點帶有8個自由度的Timoshenko 梁軸有限元模型。流體擾性轉子系統(tǒng)應考慮陀螺儀和剪切力的作用。由于軸承轉子系統(tǒng)的局部非線性特征。一個修改過的帶自由接口的模態(tài)綜合技術描述了減少擾性轉子系統(tǒng)模型的自由度。為了改變坐標引起的誤差，非線性影響仍然存在于物理空間。確保了系統(tǒng)的非線性分析的準確性和節(jié)約

41、了計算工作量。在轉子軸承支撐實驗中，油膜的氣蝕領域改變是由于中心轉子的移動和速度擾動引起的。根據(jù)油膜的物理特性，改變約束的方法是通過引入持續(xù)改變的每一階段的動態(tài)的集成與重復的雷諾茲方程式。雷諾茲邊值的流體問題得以解決是由于使用了八交點方法的有限元素分析而沒有增加計算量。非線性油膜力及其Jacobian矩陣可同時計算和得到同樣的精確度。如果系統(tǒng)元素u（例如轉子角速度ω，軸承直徑寬度比d/B，周長寬度比B’／B和距離半徑比ψ，橢圓度δ）當作

42、控制參數(shù)分配，綜合軌跡預測追蹤的延續(xù)算法和PNF方法被稱為連續(xù)周期方法問題導致計算周期運動和分岔現(xiàn)象。軸承轉子系統(tǒng)的這種混亂狀態(tài)是從能量譜中調(diào)查出來的。許多的例子顯示這項研究能節(jié)省計算量而且具有很好的精密度。</p><p><b>　　Reference</b></p><p>　　1 Brancati R，Rocca E，Rosso M，et a1．Jouma1

43、oribt and their stability for</p><p>　　rigid unbalance rotors ASME Journa1 of Tribology,1995，117：709～716</p><p>　　2 Della P L，De R E，Rossi C．Static and dynamic behavior of a rigid rotor</p>

44、;<p>　　on iourna1 bearings．Meceanica，1991，26：229-245</p><p>　　3 Kim YB，Noah S Bifurcation analysis for a modified Jeffcott rotor with</p><p>　　bearing clearance．Nonlinear Dynamics．1990,1：

45、221～241</p><p>　　4 Ch0i S K．Noach S T. Mode-locking and chaos in a Jeffcott rotor with</p><p>　　bearing clearance．ASME Jouma1 of Applied Mechanics． 1994,61：I3I～I38</p><p>　　5 Child

46、D W. Fractiona1．frequency rotor motion due to nonsymmetric clear．</p><p>　　ance effects．ASME Jouraa1 of Energy and Power,1982,104：533～541</p><p>　　6 Ehrich F F．High ordel subharmonic response of

47、 high speed rotors in bearing</p><p>　　clearance．ASME Journa1 of Vibration．Acoustics．Stress, and Reliability in Design,1988．1l0：9～16</p><p>　　7 Saito S．Calculation of nonlinear un balance respon

48、se of horizonta1 Jeffcott</p><p>　　rotors supported by ball bearings with radial clearance．ASME Paper No．85-DET-33,1985</p><p>　　8 Shiau T N，Jean A N．Prediction of steady state response of flexi

49、ble rotor</p><p>　　system with nonlinear supports：a new technique．ASM E Jouma1 0f Vibration and Acoustics．1990, 1l2：5Ol～507</p><p>　　9 Huang J L，Shiau，T N．An application of the generalized polyn

50、omial expansion method to nonlinear rotor beating system．ASME Journa1 of Vibration an d Acoustics．1994,1 13：299～3O8</p><p>　　10 Tiesheng z，Norio H．An efficient analysis of high order dynamic system</p>

51、<p>　　with 1pea1 nonlinearity．ASME Journa1 Vibration and Acoustics．1999,121:</p><p><b>　　408～416</b></p><p>　　l1 Zhang J z．Xu Q Y Zheng T S．A method for determining the period

52、ic</p><p>　　solution and its stability of a dynamic system with 1pea1 nonlinearities</p><p>　　Aeta Mechanica Sinica，l998，30(5)：572～579(In Chinese)</p><p>　　12 Nelson H D．A finite ro

53、tating shaft element using Timoshenko beam</p><p>　　theory．ASME Journa1 ofMechanica1 Design．1980，102(10)：793～803</p><p>　　13 Craig J R R．A review of time-domain and frequency-domain component<

54、;/p><p>　　modes synthesis methods．Combined Expcrimental/Analytical Modeling</p><p>　　of Dynamic Structural Systems Using Substructure Synthesis．1985：1～3l</p><p>　　14 Kinderlenhrer D．St

55、anpacchia G．An introduction to variational inequalities and their applications．New York：Academic Press,1980</p><p>　　l5 Sundarajan P,Noah S Dynamics of forced nonlinear systems using shooting／arc length cont

56、inuation method— app1ic on to rotor system．Jouma1 of Vibration and Acoustics，1997，l19(1)：10--20</p><p>　　I6 Iooss G，Joseph D D．Elementary Stability and Bifurcation Theory．New York：Springer-Veflag lnc.,1980&l

57、t;/p><p>　　17 Parker T S．Chua L O．Practica1 NUlrlerica1 Algorithms for Caotic System．</p><p>　　New York：Springer-Verlag,1989</p><p>　　18 Seydel R．Practica1 Bifurcation and Stability An