亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        電磁軸承-柔性轉子系統(tǒng)多目標加權的主動振動控制

        2016-12-05 10:30:18蔣科堅祝長生
        浙江大學學報(工學版) 2016年10期
        關鍵詞:振動

        蔣科堅, 祝長生

        (1. 浙江理工大學 信息學院,浙江 杭州 310018; 2. 浙江大學 電氣工程學院,浙江 杭州 310027)

        ?

        電磁軸承-柔性轉子系統(tǒng)多目標加權的主動振動控制

        蔣科堅1,2, 祝長生2

        (1. 浙江理工大學 信息學院,浙江 杭州 310018; 2. 浙江大學 電氣工程學院,浙江 杭州 310027)

        為了避免影響系數法求解出一個過大的校正力計算結果,解決由此造成的電磁力飽和以及轉子懸浮失穩(wěn)問題,拓展基于最小二乘法的影響系數法,提出多目標價值函數混合加權的柔性轉子振動主動控制方法.利用多個目標價值函數的加權設計,不僅實現了多個節(jié)點振動之間按權值實現振動抑制,而且能夠在不同量綱的振動幅值和所需電磁力之間實現權衡.算例仿真表明,提出的控制方法可以在實現轉子振動控制的同時,顧及到電磁軸承(AMB)電磁力的能力極限.在轉子振動能夠控制在允許的范圍內,適當減小對電磁力的需求,降低對電磁軸承的功率要求.

        柔性轉子;電磁軸承(AMB);影響系數法;主動振動控制

        主動電磁軸承(active magnetic bearing, AMB,以下簡稱電磁軸承) 是一種新型的轉子磁懸浮支承技術.電磁軸承作為轉子振動主動控制器件,可以通過施加電磁力對轉子振動進行主動控制.電磁軸承的懸浮支承為轉子系統(tǒng)動力學和轉子振動抑制提出了新的研究內容.

        從原理上區(qū)分,柔性轉子的振動控制方法可以分為兩大類:模態(tài)法(modal method)和影響系數法(influence coefficient method).模態(tài)法的基本思想是根據轉子不平衡量可以按振動模態(tài)分解成多個正交不平衡分量,每一分量只能激勵起相應的一個振動模態(tài)的動力學理論.影響系數法直接以激勵力和轉子節(jié)點振動之間的關系實現控制.與模態(tài)法不同,影響系數法的最大優(yōu)點是與轉子的振動模態(tài)信息無關,不依賴于控制對象的數學模型.實踐證明,在轉子振動的線性范圍內,采用影響系數法可以獲得較好的控制效果.

        柔性轉子動力學及其振動控制研究已有較長的歷史,但有關采用電磁軸承抑制柔性轉子振動的研究相對較少.Stephen[1]指出對柔性轉子的振動控制,無法通過轉子上有限數量的位置(離散的,有限個節(jié)點)施加電磁力,消除整個連續(xù)轉子的振動,并提出多目標節(jié)點的柔性轉子振動抑制方法.Arias等[2]為雙盤柔性轉子建立有限元模型,每個軸端視為歐拉梁模型,開展不平衡振動的建模和控制.Tseng等[3]提出根據影響系數法計算不平衡量,以卡爾曼濾波器實現測量噪聲消除的轉子實時動平衡方案.Kang等[4]采用基于影響系數法的自動平衡方法,其中不平衡質量修正用神經PD控制.謝振宇等[5]通過實驗發(fā)現,調節(jié)電磁軸承的支承剛度和阻尼能夠明顯改變柔性轉子的平動和錐動兩個剛體臨界轉速,但對一階彎曲固有頻率的影響很小.Xu等[6]開展基于電磁軸承的柔性轉子過彎曲臨界的實驗,其團隊在“十五”和“十一五”863重點項目支持下,完成了10 MW高溫氣冷堆發(fā)電機磁軸承實驗臺架的設計工作.李紅偉等[7]采用影響系數法和振型平衡法,分別對轉子的剛性模態(tài)和前兩階撓性模態(tài)進行該機動平衡.

        目前,在基于電磁軸承的柔性轉子振動控制中,影響系數法是普遍采用的方法.影響系數法往往會得出過大的、超出電磁軸承實際能力的計算結果.由于電磁軸承的懸浮控制和振動控制是疊加在一起的,電磁力一旦飽和,不只是轉子振動控制效果變差,更嚴重的是破壞了轉子懸浮控制的反饋機制,造成轉子失穩(wěn),甚至跌落事故.其次,在轉子振動控制中,目標節(jié)點振動最小化不是評判最優(yōu)的唯一標準.控制設計必須考慮電磁軸承的功率極限和成本.在振動可接受的范圍內,解決好所需電磁力要求和振動控制效果的靈活權衡.

        已有研究在影響系數法中引入最小二乘法,目的是解決經典影響系數法中目標節(jié)點數量不能大于激勵節(jié)點數量的控制問題.本文拓展了最小二乘法的思路,在價值函數中不僅包含節(jié)點振動,還加入了電磁力的因素,使得控制方法不僅可以抑制轉子各節(jié)點的振動,還可以在電磁力和轉子振動之間實現靈活調控.以較小的電磁力,將轉子振動抑制在允許范圍內.

        本文的仿真算例以KEOGH[8-12]實驗室的電磁軸承-柔性轉子實驗平臺為對象.轉子有限元建模軟件和模態(tài)分析軟件來源于該實驗室.以有限元轉子模型為核心,本文構建提出轉子振動抑制方法的控制模型,驗證了控制效果.

        1 影響系數法及最小二乘法簡介

        轉子上需要進行振動抑制的節(jié)點稱為目標節(jié)點;施加矯正力的節(jié)點稱為激勵節(jié)點.在線性范圍內,目標節(jié)點的運動與外激勵力之間可以表示為如下關系:

        (1)

        式中:n為目標節(jié)點的數量,A為目標節(jié)點的當前振動,B為目標節(jié)點的原始振動,m為激勵節(jié)點的數量,W為各個激勵節(jié)點上施加的矯正力;C為影響系數矩陣,是維復數矩陣,描述了在m個激勵節(jié)點施加單位激勵后,在n個目標節(jié)點上引起響應的矩陣關系.

        假設在激勵W作用下,目標節(jié)點的振動響應為A1;激勵改變?yōu)閃2,目標節(jié)點響應變?yōu)锳2,則有

        A1=CW1+B ,

        A2=CW2+B .

        (2)

        為了使激勵W2能夠完全抑制目標節(jié)點的振動,即A2=0,可得

        W2=W1-C-1A1.

        (3)

        可見,在系統(tǒng)線性且C已知的前提下,可以通過式(3)準確地求得激勵力W2,使得激勵后的轉子振動為零.以上所述是影響系數法的基本思想.

        式(3)中存在C-1項,因此必須n=m,C為方陣時式(3)才能求解;否則,C-1不存在.當n≠m時,式(3)的求解方法無效.

        結合最小二乘法的影響系數法被提出,可以使得目標節(jié)點的數量不受限制.

        最小二乘法的控制目標是追求在某個轉速下所有目標節(jié)點振動的平方和為最小.首先,根據每個目標節(jié)點的振動要求,定義每個目標節(jié)點振動的評價權值,是n×n對角矩陣Q,

        (4)

        式中:元素qi為第i目標節(jié)點振動的權值.

        最小二乘法以各節(jié)點振動幅值的平方和加權最小化為價值函數,即

        (5)

        式中:A為目標節(jié)點的振動;*表示矩陣的Hermitian變換,即轉置加共軛.式中的1/2是為了推導公式的格式整齊,無特別意義.下標k和k+1為離散控制中每個離散時刻的順序編號.

        式(5)中,C、Wk及Ak都為已知量.為了在k+1時刻的激勵力Wk+1能夠使目標函數J為最小極值,可以通過對J求Wk+1的偏導為零得到,即

        (6)

        通過式(6)可以解得,使價值函數J為最小極值,用Wopt表示:

        (7)

        式(7)為采用最小二乘法的影響系數法公式.可知,倘若C為方陣,且Q為單位矩陣,式(7)與(3)是完全等價的.

        最小二乘法解決了經典影響系數法中目標節(jié)點的數量限制問題,但采用最小二乘法控制是有殘余振動的.最小二乘法追求目標節(jié)點振動的加權平方和為最小,即價值函數最小.

        2 多價值函數混合加權的控制方法

        對于電磁軸承的工業(yè)應用,有限的電磁力是主要的技術瓶頸之一.要產生更大電磁力,必然需要更大體積的電磁軸承和更高驅動能力的控制系統(tǒng).對于轉子振動的最優(yōu)控制,目標節(jié)點振動最小化不是唯一的最優(yōu)評判標準,所需電磁力最小化是一個重要的考慮因數.為此,提出多價值函數混合加權的控制方法,在價值函數中加入電磁力的因素.

        與節(jié)點振動權值矩陣Q類似,電磁力權值矩陣R為m′m對角矩陣,

        (8)

        式中:m為激勵節(jié)點數量,ri為根據各個電磁軸承的驅動能力確定的電磁力之間的權值.當各個電磁軸承相同時,R可以取為單位矩陣.

        建立一個轉子振動最小化和電磁力最小化的綜合最優(yōu)控制價值函數:

        (9)

        如果直接使用式(9)為價值函數,由于轉子振動和電磁力的物理量綱不一致,即A*A和W*W兩個量的數值可能不在一個數量級上,那么數值較大的量在綜合價值函數中占絕對主導地位,而數值較小量則無關緊要.必須增加一個混合權值參數:

        (10)

        式中:p和q分別為節(jié)點振動和電磁力的混合權值參數,作用是將兩個不同物理量調整至相同的數量級,使之在數值上具有可比性.

        價值函數J與轉子振動Ak、激勵力Wk及Wk+1的關系為

        (11)

        對目標函數J求Wk+1的偏導為零時,目標函數J為最小值,即

        (12)

        根據式(12),可以解得目標函數J為最小極值的Wopt:

        Wopt=Wk+1=

        (13)

        式(13)是提出的能綜合權衡轉子振動和所需電磁力的控制計算公式.

        3 主動控制仿真

        3.1 仿真算例

        圖1 柔性轉子算例的結構和尺寸Fig.1 Example of flexible rotor system

        仿真算例的基本結構如圖1所示.它由一根均質等直徑軸和4個剛性圓盤組成.軸的長度為1 930 mm,軸的直徑為50 mm.每個圓盤的質量約為10 kg.轉子采用有限元法建模,電磁軸承的轉子軸頸部分視為剛性圓盤.具體的建模方法請參閱有關轉子動力學教科書,本文限于篇幅,不詳細敘述.根據轉子的結構分為6個節(jié)點,如圖1所示.仿真假設,轉子的不平衡質量只集中在4個剛性圓盤上,忽略其余部分的不平衡質量.4個剛性圓盤的不平衡質量都設置為1 g·mm,角度分別為4個盤從左到右[0,p/2,p,-p/2].仿真中,影響系數采用事先離線測定,控制中實時查表獲取.

        電磁軸承懸浮控制采用4個徑向自由度離散的PID控制.仿真中,電磁軸承工作在4.3 A的偏置電流下,位移剛度系數為2×106N/M,電流剛度系數為487 N/A.

        在電磁軸承等效剛度為106N/M,等效阻尼為4 800 N·s/m的條件下,采用模態(tài)軟件分析算例轉子的前四階模態(tài)振型和對應的臨界轉速,如圖2所示.其中,錐形渦動和平行渦動為剛體振型,一階彎曲和二階彎曲為柔性彎曲振型.以下所有仿真都是在該支承條件下進行的.

        圖2 算例轉子的前四階模態(tài)振型Fig.2 First four vibration modes of example

        3.2 最小二乘法的控制仿真

        圖3 1~300 rad/s轉子各節(jié)點的原始振動Fig.3 Original vibration of rotor in 1-300 rad/s

        檢測并記錄轉子所有6個節(jié)點在1~300 rad/s轉速下的原始振動,如圖3所示.圖中,n為轉速.由圖3可見,由于平動和錐動兩個剛體臨界距離較近,兩個峰值區(qū)別不明顯,模糊為一個峰值.一階彎曲臨界非常清楚,在178.6 rad/s處峰值較高.

        采用經典的影響系數法,以2個電磁軸承作為激勵節(jié)點,實現2個目標節(jié)點的振動抑制.電磁軸承在2號和5號節(jié)點.

        開展2個仿真過程.如圖4所示為當節(jié)點2和5既作為激勵節(jié)點,又作為目標節(jié)點時,轉子全部6個節(jié)點的振動變化情況.如圖5所示為以節(jié)點3、4為目標節(jié)點時,轉子各節(jié)點的振動變化情況.

        從圖4、5可見,采用影響系數法的振動控制對指定目標節(jié)點的振動抑制效果十分有效,能夠很好地抑制目標節(jié)點的振動.需要特別指出的是,圖4、5都顯示在一階彎曲臨界轉速位置有一點點小凸起,

        圖4 以2、5為目標節(jié)點的經典影響系數法控制Fig.4 Vibration control for node 2 and node 5 with classic influence coefficient method

        圖5 以3、4為目標節(jié)點的經典影響系數法控制Fig.5 Vibration control for node3 and node4 with classic influence coefficient method

        表明還有一些殘余振動.經反復驗證和數據核對發(fā)現,由于仿真采用階梯形階躍升速,階躍幅度為1 rad/s,每次階躍穩(wěn)定時間為3 s,然后記錄振動數據.由于在一階臨界轉速附近,轉子振動變化劇烈,穩(wěn)定時間相對不夠長,造成了殘余振動.依據影響系數法理論可知,只要穩(wěn)定時間足夠長,可以完全抑制目標節(jié)點的振動,沒有殘余振動.

        從圖4、5可以看出,由于影響系數法是以消除目標節(jié)點的振動為控制目標,沒有顧及其他節(jié)點的振動;其他非目標節(jié)點上的振動可能大幅增強.更糟的情況是影響系數法的控制結果在未知的轉速下引起了新的振動峰值,如圖4所示大約在300 rad/s,圖5所示在250 rad/s左右都有一個新的振動峰值.分析后發(fā)現,由于電磁軸承在這一轉速對目標節(jié)點的控制能力較弱,即影響系數較小,必須產生很大的電磁力才能抑制目標節(jié)點的振動.此時,這個巨大的電磁力已在其他非目標節(jié)點上產生不期望的強烈振動.

        采用結合最小二乘法的影響系數法進行仿真,檢驗該方法能夠兼顧多個目標節(jié)點,甚至轉子的整體振動抑制的能力.

        以電磁軸承節(jié)點2、5為激勵節(jié)點,以全部節(jié)點1~6為目標節(jié)點,取各目標節(jié)點的振動權值都相等,即Q為單位矩陣.如圖6所示為在1~300 rad/s速度下的振動控制效果.

        從圖6可見,控制效果與圖4、5有很大差異.在最小二乘法策略下,對減小轉子整體振動的效果十分有效.最重要的是因為最小二乘法能夠兼顧每個節(jié)點的振動,不會發(fā)生在非目標節(jié)點引起額外振動的情況.

        圖6 結合最小二乘法的影響系數法仿真Fig.6 Vibration control for all nodes with least squaremethod

        圖7 節(jié)點2的兩種控制效果比較Fig.7 Vibration contrast of node 2 between two controlmethods

        將圖4中節(jié)點2的振動曲線與圖6中節(jié)點2的振動曲線進行比較,按同樣比例顯示在圖7中.可見,經典影響系數法對目標節(jié)點振動的抑制效果十分理想,只有因為穩(wěn)定時間不夠引起的微量殘余振動.追求目標節(jié)點的振動完美抑制,可能需要很大的電磁力,而過大的電磁力會加劇其他節(jié)點的振動.可能只要對目標節(jié)點振動做很小的犧牲,就可以換來對整個轉子其他節(jié)點振動的整體大幅下降.最小二乘法能夠很好地實現這一點.

        3.3 多價值函數混合加權的控制仿真

        在圖6的最小二乘法的仿真中,兩個電磁軸承在y方向(垂直方向)所產生的電磁力變化幅值被記錄,如圖8所示.在超越一階彎曲臨界轉速的時刻,左、右兩個電磁軸承所需產生的電磁矯正力波動幅值達到的最大值約為190和150 N.

        圖8 最小二乘法振動控制的電磁力變化Fig.8 Magnetic force in least square method

        在進行此項仿真前,先討論式(10)中混合權值p和q的物理含義和取值依據.p為節(jié)點振動幅值平方的權值,根據已有仿真數據,轉子振動大致在0.1 mm數量級,即p為10-8m2數量級.q為電磁力平方的權值,電磁力大致在10 N數量級,即q為102N2數量級.當p∶q=1010時,表示轉子0.1 mm的振動和10 N的電磁力是等權的.通俗地講,在控制中允許以轉子振動增大0.1 mm為代價,使得電磁力減小10 N;相反,可以增大10 N電磁力為代價,換來轉子振動減小0.1 mm的改善.

        采用p∶q=1010的混合加權控制進行仿真,Q、R都取相應維數的單位矩陣.電磁力變化的仿真結果如圖9所示.轉子6個節(jié)點的振動的變化如圖10所示.將圖9、10分別與圖8、6比較,采用電磁力與轉子振動混合權值最優(yōu)控制后,左、右兩個電磁軸承所需產生的最大矯正力波動幅值分別減小到125和120 N,而各節(jié)點的振動都普遍小幅增大.

        增加電磁力權重,取p∶q=0.2×1010進行仿真,電磁力變化和各節(jié)點振動情況記錄如圖11、12所示.可見,由于電磁力的權重的增大,左、右兩個電磁軸承所需產生的最大矯正力波動幅值分別減小到約60和60 N,但各節(jié)點的振動進一步增大.

        仿真算例表明,采用電磁力與轉子振動混合權值最優(yōu)控制方法,可以根據產生電磁力的能力和轉子振動的實際允許范圍,靈活調配電磁力需求與轉子振動的權衡比例.

        圖9 p∶q=1010時的電磁力變化Fig.9 Magnetic force with p∶q=1010

        圖10 p∶q=1010時的各節(jié)點振動情況Fig.10 Vibration control with p∶q=1010

        圖11 p∶q=0.2×1010時的電磁力變化Fig.11 Magnetic force with p∶q=0.2×1010

        圖12 p∶q=0.2×1010時的節(jié)點振動情況Fig.12 Vibration control with p∶q=0.2×1010

        4 結 論

        (1) 結合最小二乘法的影響系數法,解決了目標節(jié)點數量不能大于激勵節(jié)點數量的矛盾,使得目標節(jié)點數量不受限制.它的實質是對某些節(jié)點的振動作少量犧牲,換來對整個轉子其他節(jié)點振動的整體大幅下降.

        (2) 提出多價值函數混合加權的控制方法,不僅可以在多個節(jié)點振動之間按權值比例實現振動調配,而且可以對節(jié)點振動和所需電磁力之間實現權值調控.目的是可以在轉子振動允許范圍內,適當降低對電磁力的要求,用較小的電磁力,得到可以接受的轉子振動抑制效果.

        (3) 提出的轉子振動抑制方法是建立在影響系數法基礎上的,所以應用效果和應用局限與影響系數法一樣,必須滿足轉子振動的激勵和響應在線性范圍或近似線性范圍內.

        [1] STEPHEN W. Adaptive optimal control of active balancing systems for high-speed rotating machinery [D]. Michigan: The University of Michigan, 1999.

        [2] ARIAS M, SILVA G. Finite element modeling and unbalance compensation for a two disks asymmetrical rotor system [C]∥ 5th International Conference on Electrical Engineering, Computing Science and Automatic Control. Mexico City: [s.n.], 2008: 386-391.

        [3] TSENG C, SHIH T, LIN J. A Kalman filter-based automatic rotor dynamic balancing scheme for electric motor mass production [C]∥ The Proceedings of Materials Science Forum(Vols: 505-507). Switzerland: [s.n.], 2006: 997-1002.[4] KANG Y, LIN T, CHU M. Design and simulation of a neural-pd controller for automatic balancing of rotor[C]∥ 3rd International Symposium on Neural Networks. Chengdu: [s.n.], 2006: 1104-1109.

        [5] 謝振宇,徐龍祥,李迎,等.控制參數對磁懸浮軸承轉子系統(tǒng)動態(tài)特性的影響[J].航空動力學報,2004, 19(2): 174-178.

        XIE Zhen-yu, XU Long-xiang, LI Ying, et al. Influence of control parameters on dynamic characteristics ofactive magnetic bearing system [J]. Journal of Aerospace Power, 2004, 19(2): 174-178.

        [6] XU Y, YU H, ZHAO L. Development of a magnetic bearing rotor dynamics analysis software package [C]∥ Proceedings of the 12th International Symposium on Magnetic Bearings. Arlington: [s. n.], 2010: 53-70.

        [7] 李紅偉,徐旸,谷會東.電磁軸承-撓性轉子系統(tǒng)的本機動平衡方法[J].中國機械工程,2008, 19(12): 1419-1428.

        LI Hong-wei, XU Yang, GU Hui-dong. Field dynamic balancing method in AMB—flexible rotor system [J]. China Mechanical Engineering, 2008, 19(12): 1419-1428.

        [8] CADE I, KEOGH P, SAHINKAYA M. Fault identification in rotor / magnetic bearing systems using discrete time wavelet coefficients [J]. IEEE/ASME Transactions on Mechatronics, 2005, 10(6): 648-657.

        [9] SCHLOTTER M, KEOGH P. Synchronous positionrecovery control for flexible rotors in contact with auxiliary bearings [J]. Journal of Vibration and Acoustics, 2007, 129(5): 550-558.

        [10] SAHINKAYA M, ABULRUB A, KEOGH P, et al. Multiple sliding and rolling contact dynamics for a flexible rotor/magnetic bearing system [J]. IEEE/ASME Transactions on Mechatronics, 2007, 12(2): 179-189.

        [11] KEOGH P. Transient rotor/active magnetic bearing control using sampled wavelet coefficients [J]. Journal of Engineering for Gas Turbines and Power, 2006,192(2): 549-555.

        [12] COLE M, KEOGH P, SAHINKAYA M, et al.Towards fault-tolerant active control of rotor-magnetic bearing systems [J]. Control Engineering Practice, 2004, 12(4): 491-501.

        Vibration suppressing with mixed weight for multi-targets in active magnetic bearing-flexible rotor system

        JIANG Ke-jian1,2, ZHU Chang-sheng2

        (1.CollegeofInformaticsandElectronics,ZhejiangSci-TechUniversity,Hangzhou310018,China;2.CollegeofElectricalEngineering,ZhejiangUniversity,Hangzhou310027,China)

        An exceedingly great magnetic force can often be required by the influence coefficient, so that the magnetic force saturating would destroy the stability of rotor levitating. A vibration suppressing method with mixed weight for multi-cost functions was proposed. Not only the rotor vibration can be given attention, but also the required magnetic force is considered in the cost function. The simulation results indicate that the proposed method can mediate not only among a group of the target nodes vibration, but between the rotor vibration and the required magnetic force. The advantage of the proposed method means that the requirement for the magnetic force of active magnetic bearing (AMB) can be appropriately decreased as long as the rotor vibration is within the allowable range. The method provides the feasible way for the low-power control of AMB.

        flexible-rotor; active magnetic bearing (AMB); influence coefficient method; active vibration control

        2015-09-19.

        國家自然科學基金資助項目(51477155,11272288,11172261);浙江省自然科學基金資助項目(LZ13E070001);浙江省公益技術應用研究資助項目(2015C31063);先進航空發(fā)動機協(xié)調創(chuàng)新中心資助項目.

        蔣科堅(1972—),男,教授,從事電磁軸承支承特性、轉子振動主動控制的研究. E-mail: jkjofzju@163.com

        祝長生,男,教授,博導.

        10.3785/j.issn.1008-973X.2016.10.014

        TH 133

        A

        1008-973X(2016)10-1946-06

        浙江大學學報(工學版)網址: www.zjujournals.com/eng

        猜你喜歡
        振動
        振動的思考
        科學大眾(2023年17期)2023-10-26 07:39:14
        某調相機振動異常診斷分析與處理
        大電機技術(2022年5期)2022-11-17 08:12:48
        振動與頻率
        天天愛科學(2020年6期)2020-09-10 07:22:44
        This “Singing Highway”plays music
        具非線性中立項的廣義Emden-Fowler微分方程的振動性
        中立型Emden-Fowler微分方程的振動性
        基于ANSYS的高速艇艉軸架軸系振動響應分析
        船海工程(2015年4期)2016-01-05 15:53:26
        主回路泵致聲振動分析
        UF6振動激發(fā)態(tài)分子的振動-振動馳豫
        計算物理(2014年2期)2014-03-11 17:01:44
        帶有強迫項的高階差分方程解的振動性
        人妻中文字幕乱人伦在线| 男女发生关系视频网站| 日本一区二区视频免费在线观看| 强d乱码中文字幕熟女免费| 午夜精品久久久久久99热 | 亚洲免费不卡| 日韩精品极品视频在线观看蜜桃| 加勒比东京热一区二区| 国内精品伊人久久久久网站| 99re久久精品国产| 制服丝袜人妻中出第一页| va精品人妻一区二区三区| 国产精品福利一区二区| 亚洲精品无码mv在线观看| 99精品国产闺蜜国产在线闺蜜| 日本大片一区二区三区| 亚洲av中文无码乱人伦在线咪咕 | 五月综合高清综合网| 午夜日本理论片最新片| 亚洲国产精品高清一区| 国产成人一区二区三区影院动漫 | 性无码国产一区在线观看| 国产91精品在线观看| 亚洲av无码国产精品草莓在线| 国产一在线精品一区在线观看| 午夜视频免费观看一区二区| 国产三级不卡一区不卡二区在线| 在线天堂www中文| 国产精品入口牛牛影视| 中文字幕一区二区在线| 亚洲最大成人综合网720p| 伊人狠狠色丁香婷婷综合| 国产精品白浆视频免费观看| 亚洲六月丁香色婷婷综合久久| 国产乱对白刺激视频| 婷婷综合久久中文字幕蜜桃三电影 | 亚洲六月丁香色婷婷综合久久| 日本熟妇色xxxxx日本妇| 国产精品露脸视频观看| 日本高清中文字幕二区在线| 色视频网站一区二区三区|