999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于多物理場耦合分析的機載天線仿真技術研究

2015-02-24 07:07:01賈云峰楊柳胡修魏鴻浩邱琳
電波科學學報 2015年6期

賈云峰 楊柳 胡修 魏鴻浩 邱琳

(北京航空航天大學電子信息工程學院,北京 100191)

?

基于多物理場耦合分析的機載天線仿真技術研究

賈云峰楊柳胡修魏鴻浩邱琳

(北京航空航天大學電子信息工程學院,北京 100191)

摘要為了預測、評估對物理場影響較為敏感的微帶陣列天線在復雜的溫度、應力等多物理場效應影響下的性能變化規律,本文從多場耦合的機理出發,立足于機載陣列天線裝機狀態下的復雜物理環境,深入分析陣列天線多物理場效應.通過將求解電大尺寸具有優勢的多層快速多極子算法和計算溫度場及結構場問題的有限元數值算法的聯合運用,成功突破模型關聯技術、網格共享技術、多物理場協同仿真控制等關鍵技術,形成了一套較為成熟的解決天線多物理場仿真中典型電熱結構耦合問題的仿真流程.

關鍵詞多物理場耦合;微帶陣列;機載天線;天線間隔離度

資助項目: 航空基金(20112051017)

聯系人: 楊柳 E-mail:yangliu2116@126.com

引言

多物理場耦合(Multi-Physics Coupling, MPC)問題在自然界和工程應用中廣泛存在,其表現形式及種類繁多[1].由于多場耦合問題內在的復雜性,研究人員在遇到多場耦合問題時一般都是對問題做出較大的簡化,只考慮一個主要物理場的效應,而忽略其它多場耦合的影響.這種分析方法忽略了多物理場之間的相互耦合,所得到的結果與實際情況偏差較遠,無法用于天線性能的精確預測[2-3].尤其是在對許多高頻天線進行輻射特性分析時,發現其電性能對結構、應力等物理場極其敏感.其中比較具有代表性的就是微帶陣列天線[4].微帶陣列天線由于具有高增益、高功率、低旁瓣、波束可電控等特性,通常作為雷達的發射天線廣泛應用于航空、航天領域的氣象探測、火控瞄準、武器引導、目標搜索、空域預警、導航、數據傳輸、通信等方面[5-6].

在影響微帶陣列天線性能指標的諸多多物理場效應中,熱應力所帶來的形變是其中具有代表性的問題[7-12],因此,其多物理場耦合問題可以用電熱耦合來表征.電熱耦合關系是指陣列天線工作時熱效應及其帶來的結構形變與電磁場之間的相互影響、相互制約關系,具有典型的多學科交叉的特點[13].

目前考慮多物理場效應往往還是集中在高壓高功率、大電流相關領域,而針對天線的多物理場效應的分析相對較少:葉菁等人從空饋體制相控陣雷達天線的熱設計方案著手,提出了用模型試驗與計算機仿真相結合的方法對天線熱設計方案進行驗證[12];陳杰等人根據星載合成孔徑雷達(Synthetic Aperture Radar, SAR)相控陣天線結構特點,建立了天線熱變形誤差的數學模型,提出了定量分析星載SAR相控陣天線熱變形誤差對模糊性能影響的有效方法[13];陳國強提出了一種對相控陣天線進行流-熱-結構耦合分析的流程,得出了影響天線陣電性能的初步結論[14];段寶巖院士團隊分析了矩形陣面相控陣微帶天線輻射單元位置誤差對天線電場和方向性的影響[15].

由于算法的限制,以上方法均只適用于對天線自由空間的輻射特性進行討論,無法應用于天線裝機性能的分析.該文依托于航空基金(20112051017)天線陣列多物理場仿真技術研究,針對天線裝機性能分析的需求,運用模型修復等手段實現在不同算法間的模型轉換,實現了機載天線在多物理場影響下電性能的評估與分析.

1多物理場耦合仿真分析

1.1 天線建模

該文以項目提供的某型機多普勒雷達天線參數進行建模和多物理場耦合仿真方法說明.該天線為一個含360單元的微帶陣列天線,其陣列排布為12×30,中心頻率fr=13 GHz,帶寬約為500 MHz.單元結構以及陣列模型如圖1所示,紅色部分即為饋源.

圖1 12×30陣列天線及單元模型

其陣列單元可以分為三個部分:表面金屬貼片、基板以及饋電部分.其中,金屬貼片材質為銅,尺寸為7.5 mm×9.1 mm,基板材質為聚四氟乙烯(FR-4),尺寸為15.1 mm×18.21 mm,饋點位置位于距離貼片邊緣2.3 mm處,基板厚度為0.6 mm.天線單元中心間距為21.14 mm.

為了方便分析,假設該天線各單元在工作時饋點供電、陣面溫度環境等均相同.由于該微帶陣列天線各單元采用同軸饋電,天線各單元間物理結構固定相對獨立,因此可以將天線的溫度、位移場分析轉換為兩個部分:1) 天線單元貼片金屬和饋電部分的溫度場、位移場分析;2) 基片的溫度場、位移場分析.最后將兩個部分得到的形變模型結合,進行電磁場分析.

1.2 溫度場分析

熱分析遵循熱力學第一定律,即能量守恒定律:對于一個封閉的系統(沒有質量的流入或者流出)有關系式為

Q-W=ΔU+ΔEk+ΔEp.

(1)

式中: Q為熱量; W為功;ΔU為系統能量;ΔEk為系統動能;ΔEp為勢能.

對于大多數的工程傳熱問題,有

ΔEk=ΔEp=0,

(2)

通常考慮沒有做功W=0,因此Q=ΔU.

由此可以得出以下結論:在機載天線熱分析的過程中,在不考慮環境對天線影響(做功、熱傳遞)的情況下,天線的熱能完全由天線工作時金屬表面的電流所產生,這樣的簡化為天線本身的溫度場分析提供了便利.下面分別對其最高和最低溫度在溫度場內進行分析.

機載陣列天線所能達到的最低溫度是由其非工作狀態下的自然環境溫度所決定的,該文針對低空飛行的直升機所處的自然環境,其最低溫度通過對國內主要城市極端低溫情況進行統計調研發現約為-40℃.

機載陣列天線的最高溫度是由其滿負荷工作狀態下功率器件發熱所造成的.將貼片與饋源相接觸的面熱流密度設置為最大值20 W/cm2[14],陣面與空氣接觸,考慮通風設施為正常風冷條件下,設空氣溫度為20℃.根據銅的材料特性通過有限元算法仿真計算得到單元貼片在工作狀態下所能達到的最高溫度為130℃.

根據熱傳導的相關特性,在該文的穩態分析中認為基板達到最高溫度與貼片相同為130℃.

因此,通過溫度場的分析得到機載微帶陣列天線所處的溫度范圍是[-40℃,130℃].

1.3 位移場分析

根據溫度場中得到的天線所處的[-40℃,130℃]溫度區間,該文以10℃為步進,考慮18個溫度觀測點的天線結構形變情況.該部分針對材料不同的兩個部分(即饋源與金屬貼片部分以及基板部分)分別進行討論[12-13].

陣列天線陣面由于熱源的存在,陣面溫度升高,將產生熱變形,變形后的陣面又將影響溫度的傳導以及熱交換,熱與結構的這種耦合應該是雙向的.然而,實際應用中,由于天線的陣面熱變形很小,變形后的陣面對溫度場的影響非常小,幾乎可以忽略,因此認為天線只存在由溫度場到結構場的單向耦合.

以金屬貼片的位移場仿真為例進行結構形變的分析.圖2和圖3分別為100℃下的金屬貼片結構形變以及應力分布情況.

圖2 100℃金屬貼片及饋點結構形變圖

圖3 100℃金屬貼片及饋點應力分布圖

在圖2中,形變最大的部位(紅色部位)分別位于金屬貼片偏離饋源距離較遠的頂角部位.在不考慮貼片與基板間應力影響的情況下,金屬板受熱膨脹時在水平方向沒有阻力存在,其形變值由饋點中心向四周呈現遞增的趨勢;由于金屬片下表面與基板共面,金屬片無法進行豎直方向的形變.因此,其膨脹形變最大值出現于金屬片四周頂角上端符合實際情況.而由于饋點與其下端饋線是相對固定的,形變最小值出現于饋點底部,這與實際情況也是相符合的.

由圖3可以看出,由于金屬貼片水平方向的形變相對無阻力,因此整個貼片部位未出現應力(均為藍色).而由于饋點與其下端饋線的相對固定,饋點的膨脹與下端固定部位間產生了較大的熱應力,其應力最大部位為饋點下端附近,約190 Mpa,這對饋點部分的材料屬性提出了更高的要求.該文中僅考慮天線的電磁特性,在此處認為該饋點處可以承受該應力并保持正常工作.

以20℃為理想環境溫度,基于熱脹冷縮原理,當溫度升高時,模型受熱膨脹,將膨脹后形變值定義為正;當溫度降低時,模型緊縮,將緊縮后模型的形變值定義為負.統計18個溫度觀測點的形變仿真結果,如表1所示.

表1 金屬貼片結構形變仿真結果

表1顯示:當溫度降低時,隨著溫度相對于理想環境溫度越低,陣列單元模型形變絕對值越大,最大收縮值達到7.28 μm;當溫度升高時,隨著溫度相對于理想環境溫度越高,陣列單元模型形變絕對值越大,最大膨脹值達到13.3 μm.

采用同樣的方法對基板進行結構形變仿真,材料為FR-4,仿真結果如表2所示.

表2 基板結構形變仿真結果

對比表1與表2中對應溫度的形變值數據,基板的結構形變值遠大于金屬貼片,因此由基板形變導致的陣列分布的變化將對天線的電性能產生更大的影響.

1.4 電磁場分析

該文中所重點關注的是機載天線裝機狀態下的電磁特性.在溫度場和位移場中,前文運用了有限元數值算法對問題進行仿真求解.如果繼續使用該算法對中心頻率達到10 GHz的模型進行電磁仿真計算,在保證計算精度的前提下需要將網格單元尺寸設置為5 mm左右,則單位體積內的網格數量將達到106~107.而對于動輒數米甚至數十米的飛機機身而言,其龐大的網格數量已遠遠超出了大部分服務器工作站的計算能力.即使在犧牲仿真精度的前提下削減部分網格數量,但仿真計算所需長達數周甚至數月的計算時間也是大部分工程項目所不允許的.

該文中采用求解電大尺寸問題具有較大優勢的多層快速多極子算法(Multilevel Fast Multipolc Method,MLFMM)對形變天線的裝機特性進行分析.因此,如何將有限元算法得到的體網格模型轉換為可用于MLFMM算法計算的面網格問題成為了研究中的一個難點.

該文采用網格模型重構的方法解決了這一難題.利用ICEM軟件與ANSYS軟件間的接口程序,將形變后的四面體網格模型進行模型重構,形成實體模型,對該實體模型進行修模后重新劃分為MLFMM中可識別的三角形面網格模型.

在整個模型的提取、重構和剖分過程中,尤其是重構時,經常會發生實體模型缺少線、點等模型信息,因此需要靈活采用面切割等操作對模型進行修復,此過程較為繁瑣,但有利于提高仿真精度.

將金屬貼片以及基板的形變模型導入電磁仿真軟件中,首先進行單個單元的電磁仿真分析.網格材料分別設置為理想金屬導體以及介質材料;并將饋電形式設置為同軸饋電;仿真頻率設置為13 GHz;由于該文重點觀察微帶陣列天線的方向圖,此處饋源幅值為1 V.

隨著溫度的變化,對天線輻射特性有重要影響的介電常數也變化,根據文獻[17]中研究的A類介質基板介電常數隨溫度升高逐漸減小的變化趨勢,該微帶陣列天線介質基板FR-4的介電常數溫變化如圖4所示.計算天線增益時,對應不同溫度點設置該基板對應的介電常數.

圖4 FR-4介電常數隨溫度變化關系曲線

為了對比該微帶陣列天線在不同溫度條件下的輻射特性,選擇4個具有代表性的溫度點下陣列單元在φ=0°平面的方向圖曲線進行對比,結果如圖5所示.

圖5 不同溫度下陣列單元φ=0°平面方向圖曲線

隨著溫度的變化,微帶天線結構和介電常數發生變化,導致其輻射特性有所變化,主瓣輻射增益變化約3 dB(100℃時增益最大為10.33 dB和50℃時增益最小為7.35 dB),這是因為基板和金屬片的結構模型以及介電常數隨溫度變化后已不滿足最初設計時理想常溫達到最大輻射增益,而且隨著介電常數隨溫度的變化,在100℃時的天線結構和介電常數恰好匹配良好,主瓣輻射電平值達到最大.

采用等效源的仿真方法對360單元微帶陣列天線的自由空間輻射特性進行仿真計算,常溫條件(20℃)天線輻射方向圖與測試數據比較如圖6所示.

圖6 20℃天線輻射方向圖與測試數據比較

由圖6對比可得仿真的結果可靠,將主瓣放大,可以看到測試數據比仿真數據偏小,有可能是此時測量天線進入穩定工作狀態產生一定熱量,從而導致主瓣增益略小于理想20℃仿真情況.

將不同溫度下微帶陣列天線在φ=0°平面的方向圖曲線進行對比,結果如圖7所示.并將θ=0°和θ=90°方位隨溫度變化的電平值列于表3和表4.

圖7 12×30陣列天線E面方面圖曲線

表3 θ=0°方位電平隨溫度變化

表4 θ=90°方位電平隨溫度變化

由圖7可知:不同溫度下,天線在z軸方向的增益變化較小(表3顯示變化小于1 dB:130 ℃時增益最大為33.690 dB,-40 ℃時增益最小為33.614 dB);較為明顯的改變在于其旁瓣電平,尤其是在θ=90°附近變化較為劇烈,部分溫度點對應的方向圖曲線在該方向上形成了新的副瓣.

對微帶天線的主瓣與副瓣的幅值變化進行進一步分析發現,該天線主瓣增益的變化主要源于陣列單元的方向圖最大增益隨溫度的變化和介質基板介電常數隨溫度變化,而θ=90°處副瓣增益的變化主要源于基板伸縮造成的陣列單元間距的變化,導致陣列單元不滿足其最佳耦合條件[16],致使副瓣電平抬高,以至于產生新的柵瓣.

表4數據顯示:θ=90°方位輻射增益最大變化值達到67 dB(最大變化為:20 ℃時增益最小為-57.5 dB,130 ℃時增益最大為9.50 dB).該方位電平幅值的變化使得該天線與同一飛機平臺上的其他天線間相互耦合干擾情況急劇惡化.

2多物理場耦合結果分析

2.1 隔離度分析

通過對該天線的裝機位置進行研究,發現其θ=90°平面附近存在數部重要天線,其中,高度表發射信號(4.3 GHz)的三次諧波正好落在該微帶陣列天線的工作頻帶范圍內,因此必須對兩部天線間的隔離度進行分析.

高度表天線模型及其三次諧波輻射方向圖如圖8所示.

多普勒雷達天線位于坐標(0,0,-22)處,高度表天線位于坐標(-460,0,-22)處,兩部天線間距離0.46 m,兩部天線的相對位置如圖9所示.

圖9 天線裝機位置與局部示意圖

由于該文重點關注微帶陣列天線在溫度影響下的電熱耦合特性,此處為了簡化分析,假設高度表天線處于理想狀態.

通過仿真計算微帶陣列天線接收到的功率來分析兩部天線間的耦合情況.天線間隔離度的計算公式為

ISO=-LC

(3)

式中:ISO為天線間隔離度,dB;LC為天線耦合系數,dB;Pr為接收功率,W;Pt為發射天線功率,W.將高度表的干擾作為等效源時,Pt為1 W,通過公式(3)可分別將對應溫度點下的接收功率轉換為耦合值,不同溫度觀測點下兩部天線的隔離度數據如表5所示.

表5 多普勒雷達與高度表天線間隔離度

隨著微帶陣列天線在θ=90°處副瓣幅值的升高,兩部天線間的隔離度也隨之相對應的趨勢減小.隔離度下降幅度最大達到62.11 dB(最大變化為:20℃時隔離度最大為90.16 dB,120℃時隔離度最小為28.15 dB).

2.2 安全余量分析

當隔離度低于某一限值時,將帶來嚴重的干擾問題,接收機的安全余量A/dB的計算公式如(4)所示:

A=Si-Pt+ISO-ηr-ηt+Lcab+Lob.

(4)

式中:Si為接收機靈敏度,dB/dBm;Pt為發射機功率,dB/dBm;ηr和ηt分別為接收天線和發射天線的效率,dB;Lcab為線纜衰減,dB;Lob為接收機帶外衰減,dB.

根據GJB151A-97中規定,該多普勒雷達接收機安全余量需大于6 dB,其靈敏度為-90.1 dBm;接收天線效率為-8.792 dB;線纜衰減為3 dB;發射天線效率為-1.2 dB;發射機功率為27 dBm;帶外衰減取值按GJB151A-97中的CE106規定:“二次和三次諧波抑制50+10lgP(P為基波峰值輸出功率,W)或80 dB,取抑制要求較小者.”此處雷達的基波峰值輸出功率為5 W,因此發射機帶外衰減為57 dB.帶入公式(4)計算得到當留有6 dB的安全余量時,天線隔離度不得低于53.108 dB.

將隔離度隨溫度變化趨勢以曲線的形式畫出,并與允許的最低隔離度比較,結果如圖10所示.在溫度為80℃、90℃、100℃、110℃、120℃、130℃時微帶陣列天線接收機安全余量將小于6 dB,低于本項目中該系統的安全余量最低限值,該多普勒雷達有可能受擾.

圖10 12×30陣列天線隔離度分布曲線圖

3結論

該文通過對多物理場耦合機理的研究,以機載多普勒雷達為例,運用ANSYS、ICEM以及FEKO等軟件,實現了天線的電-熱-結構三場耦合仿真技術.通過討論該陣列天線在-40℃~130℃溫度區間范圍內熱應力場變化導致輻射特性變化的情況,得到如下結論:

1) 其主要輻射方向上,輻射增益變化較小,變化小于1 dB,其影響因素為陣列單元的主要輻射方向輻射增益隨溫度的變化和介質基板介電常數隨溫度的變化.

2) 其θ=90°方位電平幅值隨溫度變化明顯,電平變化最大值達到67 dB,主要源于基板伸縮造成的陣列單元間距的變化,導致天線陣列排布情況發生改變,進而影響其柵瓣數量.

3) 天線裝機后由于其輻射特性在多物理場影響下發生改變,影響其自身工作性能.該天線與高度表天線之間的隔離度隨溫度變化明顯,隔離度最大下降62.11 dB,性能降級嚴重,在溫度80℃到130℃時,多普勒雷達安全余量將小于6 dB,有可能受擾.

參考文獻

[1]葉鳴, 賀永寧, 崔萬照. 基于電熱耦合效應的微帶線無源互調機理研究[J]. 電波科學學報, 2013, 28(3): 220-225.

YE Ming, HE Yongning, CUI Wanzhao. Passive intermodulation mechanism of microstrip lines based on the electro-thermal coupling effect[J]. Journal of Radio Science, 2013, 28(3): 220-225. (in Chinese)

[2] HECKLER M V T, DREHER A. Performance of micro-strip antenna arrays installed on aircraft[J]. Aerospace Science and Technology, 2013, 26(1): 235-243.

[3] RIFAI S M, FERENCZ R M, WANG W, et al. The Role of Multiphysics Simulation in Multidisciplinary Analysis[R]. New York: National Aeronautics and Space Administration, Lewis Research Center, 1998.

[4] WANG H S C. Performance of phased-array antennas with mechanical errors[J]. IEEE Transactions on Aerospace and Electronic Systems, 1992, 28(2): 535-545.

[5] 林昌錄. 天線工程手冊[M]. 北京:電子工業出版社, 2002: 463-530.

[6] CHUNG D C, CHOI S Y, KO Y H, et al. Circularly polarized HTS microstrip antenna array[J]. IEEE Transactions on Applied Superconductivity, 2003, 13(2): 301-304.

[7] MASAYUKI N, EIHISA M, YOSHISDA K, et al. Development of thermal control for phased array antenna[C]//21st International Communications Satellite Systems Conference and Exhibit, AIAA-2003-2226: 1-10.

[8] HOFF N J. Thermal buckling of supersonic wing panels[J]. Journal of the Aeronautical Sciences, 2012, 23(11): 1019-1028.

[9] 李維特, 黃保海, 畢仲波.熱應力理論分析及應用[M]. 北京:中國電力出版社, 2004.

[10]BOLEY B A, WEINER J H. Theory of Thermal Stresses[M]. New York: Courier Dover Publications, 2012.

[11]GOMEZ D, MATEO P D, ALTET S J. Electro-thermal coupling analysis methodology for RF circuits[J]. Microelectronics Journal, 2012, 43(9): 633-641.

[12]葉菁. 相控陣雷達天線的熱設計[J]. 電子機械工程, 2001, 89(1): 42-46.

YE Jing. Thermal design of phased array antenna[J]. Electro-mechanical Engineering, 2001, 89(1): 42-46. (in Chinese)

[13]陳杰, 周蔭清. 星載 SAR 相控陣天線熱變形誤差分析[J]. 北京航空航天大學學報, 2004, 30(9): 839-843.

CHEN Jie, ZHOU Yinqing. Ambiguity performance analysis of spaceborne SAR with thermal mechanism errors in phased array antenna[J]. Journal of Beijing University of Aeronautics and Astronautics, 2004, 30(9): 839-843. (in Chinese)

[14]陳國強. 有源相控陣天線流-熱-結構耦合研究[D]. 西安: 西安電子科技大學, 2008.

CHEN Guoqiang. Research of Fluid-Thermal-Structural Coupling for Active Phased Array Antenna[D]. Xi’an: Xidian University, 2008. (in Chinese)

[15]WANG C S, DUAN B Y, ZHANG F S, et al. Coupled structural-electromagnetic-thermal modeling and analysis of active phased array antennas[J]. IET Microwaves, Antennas & Propagation, 2010, 4(2): 247-257.

[16]RAO Jiaren, PENG Zong, BECERRA D R. Mutual coupling effect of microstrip antenna array[J]. Procedia Engineering, 2012, 29: 1984-1988.

[17]KABACIK P, BIALKOWSKI M E. The temperture dependence of substrate parameters and their effect on microstrip antenna performance[J]. IEEE Transactions on Antennas and Propagation, 1999, 47(6): 1042-1049.

賈云峰 (1975-),男,湖北人,副教授,研究方向為電磁場、電磁兼容.

楊柳(1992-),女,河北人,碩士研究生,研究方向為電磁兼容.

胡修(1987-),男,江蘇人,碩士研究生,研究方向為電磁兼容.

郭小路, 陶海紅, 黨博. 基于矩陣分解的多輸入多輸出雷達解相關算法[J]. 電波科學學報,2015,30(6):1033-1038. 10.13443/j.cjors. 2015011601

GUO Xiaolu, TAO Haihong, DANG Bo. Matrix decomposition based de-correlation algorithm in MIMO radars[J]. Chinese Journal of Radio Science,2015,30(6):1033-1038. (in Chinese). doi: 10.13443/j.cjors. 2015011601

賈云峰, 楊柳, 胡修, 等. 基于多物理場耦合分析的機載天線仿真技術研究[J]. 電波科學學報,2015,30(6):1025-1032. doi: 10.13443/j.cjors. 2014120101

JIA Yunfeng, YANG Liu, HU Xiu, et al. The technology of airborne antenna simulation based on analysis of coupled multi-field[J]. Chinese Journal of Radio Science,2015,30(6):1025-1032. (in Chinese). doi:10.13443/j.cjors. 2014120101

The technology of airborne antenna simulation based on

analysis of coupled multi-field

JIA YunfengYANG LiuHU XiuWEI HonghaoQIU Lin

(SchoolofElectronicandInformationEngineering,BeihangUniversity,Beijing100191,China)

AbstractIn order to assess the performance pattern of the environmentally sensitive micro-strip antenna array within multi-field such as complex temperature and stress field, based on the multi-field coupling mechanism, this paper analyzes the multi-field effects on the micro-strip antenna array under the circumstances of the airborne antenna application within complicated physical environment. By combining multilevel fast multi-pole method, which has the advantage of calculating electrically large size antenna, and finite element method, which is a classic method for calculating temperature field and structure field, this paper is able to achieve a breakthrough in terms of some of the key techniques such as the parallel modeling technique, the grid sharing technique, as well as the multi-field collaborating simulation technique, and to establish a mature simulation process in terms of solving the electric-thermal coupling problems during the antenna multi-physic field simulation.

Key wordsMPC; micro-strip antenna array; airborne antenna; antenna isolation

作者簡介

收稿日期:2014-12-01

中圖分類號TN957

文獻標志碼A

文章編號1005-0388(2015)06-1025-08

主站蜘蛛池模板: 手机精品福利在线观看| 欧美日韩成人| 国产在线精品香蕉麻豆| 久久美女精品国产精品亚洲| 青青国产在线| 性色一区| 成年A级毛片| 亚洲精品国产首次亮相| 中文字幕无线码一区| 在线欧美国产| 天堂网亚洲综合在线| 亚洲欧美综合在线观看| 国产丝袜无码精品| 精品福利视频网| 成人免费视频一区| 99热这里只有精品免费| 日韩专区欧美| 一本一本大道香蕉久在线播放| 色悠久久久久久久综合网伊人| 重口调教一区二区视频| 成人小视频网| 国产成人在线无码免费视频| av天堂最新版在线| 国产男女XX00免费观看| 免费在线色| 国产在线一区视频| 国产青榴视频在线观看网站| 日韩国产精品无码一区二区三区| 91小视频在线| 久久久久亚洲av成人网人人软件| 伊人久久精品无码麻豆精品| 麻豆精品国产自产在线| 国产精品久久精品| 精品国产网| 亚洲AV电影不卡在线观看| 国产三级国产精品国产普男人| 国产性生大片免费观看性欧美| 免费人成视网站在线不卡| 思思热精品在线8| 国产96在线 | 97在线免费视频| 色天天综合| 国产高清不卡| 青青草原国产一区二区| 综1合AV在线播放| 国产主播一区二区三区| 精品久久综合1区2区3区激情| 毛片最新网址| 毛片基地视频| 激情无码视频在线看| 国产剧情国内精品原创| 在线免费不卡视频| 国产精品毛片一区视频播| 成人在线观看一区| 国产精品无码一二三视频| 永久免费av网站可以直接看的| 中国丰满人妻无码束缚啪啪| 国产精品浪潮Av| 国产一区二区人大臿蕉香蕉| 色亚洲成人| 全色黄大色大片免费久久老太| 国产免费久久精品99re不卡| 国产精品美女免费视频大全 | 小说区 亚洲 自拍 另类| 欧美啪啪精品| 色婷婷综合激情视频免费看| 日韩成人免费网站| 免费在线观看av| 97在线视频免费观看| 亚洲精品国产综合99久久夜夜嗨| 乱人伦视频中文字幕在线| 中美日韩在线网免费毛片视频| 国产真实自在自线免费精品| 成人福利在线免费观看| h视频在线观看网站| 亚洲成人动漫在线观看| 美美女高清毛片视频免费观看| 精品久久久久久中文字幕女| 久久亚洲综合伊人| 99视频在线看| 中国丰满人妻无码束缚啪啪| 亚洲有无码中文网|