• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    典型金屬低溫固體界面接觸熱阻數(shù)值模擬及關(guān)聯(lián)式

    2022-11-04 09:17:48王干緣黃永華
    低溫工程 2022年5期
    關(guān)鍵詞:黃銅蒙特卡羅熱導(dǎo)率

    王干緣 沈 逸 黃永華

    (上海交通大學(xué)制冷與低溫工程研究所 上海 200240)

    1 引言

    兩個(gè)相互接觸的固體表面,相互之間的導(dǎo)熱實(shí)際上僅通過在總面積遠(yuǎn)小于名義接觸面積的離散接觸點(diǎn)完成,從而形成了接觸熱阻。接觸熱阻的存在嚴(yán)重影響接觸界面間熱量的傳遞能力,在實(shí)際工程應(yīng)用中往往產(chǎn)生消極影響。因此,準(zhǔn)確計(jì)算界面間的接觸熱阻并提供經(jīng)驗(yàn)關(guān)聯(lián)式,對(duì)涉及接觸熱阻問題的固體導(dǎo)熱相關(guān)應(yīng)用設(shè)計(jì)十分必要。

    以往對(duì)接觸熱阻的研究可從實(shí)驗(yàn)測(cè)量和數(shù)值模擬兩方面歸類。Alcock 和Mech[1]于1943 年首先提出界面熱阻的點(diǎn)接觸理論。隨后,Carslaw 和Jaeger[2]建立了圓柱(圓盤)點(diǎn)接觸模型,包括單點(diǎn)接觸模型和多點(diǎn)接觸模型,并認(rèn)為多點(diǎn)接觸模型可視作多個(gè)單點(diǎn)接觸模型的并聯(lián)。Fenech 和Rohsenow[3]進(jìn)一步細(xì)化圓柱(圓盤)點(diǎn)接觸模型,并給出了圓柱接觸點(diǎn)的換熱系數(shù)和接觸熱阻計(jì)算方法。Cooper 等[4]根據(jù)接觸面的溫度分布和邊界條件給出了更加詳細(xì)的圓柱模型接觸熱阻計(jì)算方法。考慮到實(shí)際接觸表面上各接觸點(diǎn)可能呈現(xiàn)一定的弧度,Clausing 和Chao[5]引入了半球模型,并給出了半球模型接觸熱阻表達(dá)式。Bahrami 等[6]假設(shè)接觸面的觸點(diǎn)發(fā)生彈性變形,建立了適合預(yù)測(cè)粗糙表面接觸壓力分布的半球模型。Williams 和Major[7]則將接觸點(diǎn)假設(shè)為錐型,提出了圓錐計(jì)算模型。Yovanovich[8]以及Shai 和Santo[9]對(duì)比了上述接觸模型,指出圓錐模型更加符合實(shí)際情況。常用的商用仿真軟件也被很多學(xué)者用作計(jì)算接觸熱阻的工具。如Rostami 等[10]利用ANSYS 軟件建立了接觸面的三維模型,分析了接觸點(diǎn)幾何形狀、尺寸等因素對(duì)接觸熱阻的影響。沈軍等[11]利用Delaunay 三角剖分算法對(duì)溫度場(chǎng)進(jìn)行了非結(jié)構(gòu)預(yù)處理,并計(jì)算了常溫下的接觸熱阻,發(fā)現(xiàn)采用30°錐角的圓錐單點(diǎn)模型時(shí)模擬結(jié)果精度最高。吳登倍[12]采用Msc.Patran/Nastran 軟件對(duì)0—300 ℃下接觸表面形貌和接觸熱阻大小之間的關(guān)系進(jìn)行了數(shù)值模擬分析,結(jié)果表明接觸熱阻大小和粗糙表面的高度成正比,和溫度成反比。陸敏恂等[13]基于ANSYS 軟件選用矩形微凸體模型對(duì)25 ℃下MOSFET 功率管和散熱片之間的接觸熱阻進(jìn)行了模擬計(jì)算,也證實(shí)接觸熱阻隨著粗糙度的增大而增加,隨著間隙的導(dǎo)熱介質(zhì)的熱導(dǎo)率增大而減小。

    文獻(xiàn)調(diào)研表明,以往研究主要是針對(duì)室溫以上的接觸熱阻數(shù)值模擬,鮮有對(duì)低溫溫區(qū)接觸熱阻進(jìn)行計(jì)算的報(bào)導(dǎo),更缺乏低溫下典型材料的計(jì)算關(guān)聯(lián)式。本研究分別使用蒙特卡羅法、數(shù)值積分法和有限元法,對(duì)10—30 K 溫區(qū)內(nèi)低溫工程上常用的黃銅、304 不銹鋼和紫銅3 種材料的接觸熱阻值進(jìn)行計(jì)算,結(jié)合與試驗(yàn)數(shù)據(jù)的對(duì)比分析,嘗試給出相關(guān)經(jīng)驗(yàn)關(guān)聯(lián)式。

    2 接觸熱阻實(shí)驗(yàn)測(cè)量數(shù)據(jù)

    為了評(píng)估后文所開展的多種數(shù)值計(jì)算方法的準(zhǔn)確性,本研究搭建了低溫接觸熱阻實(shí)驗(yàn)系統(tǒng)并進(jìn)行了數(shù)據(jù)測(cè)量。下面簡(jiǎn)要介紹該實(shí)驗(yàn)系統(tǒng)的構(gòu)成,各參量測(cè)量精度、實(shí)驗(yàn)樣品以及相關(guān)實(shí)驗(yàn)數(shù)據(jù)。

    2.1 實(shí)驗(yàn)系統(tǒng)

    所搭建的實(shí)驗(yàn)系統(tǒng)如圖1 所示。其中測(cè)試段基于穩(wěn)態(tài)疊片法測(cè)量接觸熱阻,在其兩端分別設(shè)置熱源和冷源。熱源為25Ω 電熱絲,冷源則選用G-M 制冷機(jī)。各儀器設(shè)備如表1 所示。

    表1 主要儀器設(shè)備情況Table 1 Primary equipment and devices

    圖1 測(cè)量低溫接觸熱阻的實(shí)驗(yàn)系統(tǒng)Fig.1 Experimental system for measuring TCR

    2.2 實(shí)驗(yàn)樣品

    實(shí)驗(yàn)被測(cè)樣品材質(zhì)為黃銅、304 不銹鋼以及紫銅。其中,黃銅和紫銅組包含3 種粗糙度,分別為Ra1.6、Ra3.2 以及Ra6.3,304 不銹鋼組有Ra3.2 一種粗糙度;Ra3.2 黃銅組包含0.8—1.4 N·m 4 個(gè)不同螺栓扭矩工況,其他兩個(gè)粗糙度以及紫銅、304 不銹鋼組的螺栓扭矩均固定為1.0 N·m。接觸界面溫度范圍為10—30 K。共計(jì)30 組實(shí)驗(yàn)工況。

    2.3 測(cè)量結(jié)果

    上述各工況下的接觸熱阻測(cè)量值如表2 所列。需要說明的是,紫銅材料由于其性質(zhì)特殊,從測(cè)量結(jié)果發(fā)現(xiàn)穩(wěn)態(tài)疊片法不適用。

    表2 接觸熱阻實(shí)驗(yàn)測(cè)量值Table 2 Experimentally measured value of TCR

    3 低溫接觸熱阻的數(shù)值計(jì)算

    3.1 蒙特卡羅法

    蒙特卡羅法是把概率作為研究對(duì)象的數(shù)值模擬方法,通過構(gòu)建一個(gè)符合系統(tǒng)概率特性的概率模型,對(duì)離散系統(tǒng)進(jìn)行隨機(jī)性試驗(yàn)并獲得其特性。將這一思想應(yīng)用于接觸熱阻計(jì)算,即:將具有一定粗糙度的固體表面視為若干隨機(jī)高度粗糙峰的集合,粗糙峰的高度符合一定的分布規(guī)律。當(dāng)兩個(gè)平面接觸時(shí),統(tǒng)計(jì)上下表面的粗糙峰接觸情況。若給定一個(gè)界面間距,則高度大于這一間距的粗糙峰發(fā)生接觸產(chǎn)生形變,而高度小于這一間距的粗糙峰則不發(fā)生接觸。對(duì)發(fā)生接觸的兩個(gè)粗糙峰之間的接觸熱阻采用單點(diǎn)模型進(jìn)行計(jì)算,并將所有單點(diǎn)接觸熱阻視為并聯(lián),從而可得到整個(gè)界面的總接觸熱阻。應(yīng)當(dāng)指出,上述方法需遵循以下基本假設(shè):

    (1) 熱流垂直通過粗糙接觸表面且恒定不變,名義接觸面積為Am,不考慮徑向漏熱;

    (2) 兩粗糙表面接觸時(shí),將其等效為一粗糙表面和一剛性平面接觸,且接觸時(shí)發(fā)生純塑性變形,如圖3a 所示;

    (3) 粗糙峰分布均勻,且形變和傳熱過程互相獨(dú)立,即視為一個(gè)離散系統(tǒng);

    (4) 上下表面粗糙峰斜率分別為a1、a2,材料硬度分別為H1、H2,兩粗糙面粗糙度分別為σ1、σ2,粗糙峰峰點(diǎn)密度分別為η1、η2,峰高z的概率密度函數(shù)φ(z)服從正態(tài)分布。每個(gè)粗糙峰的單點(diǎn)模型如圖2 所示。

    圖2 蒙特卡羅法接觸熱阻模型Fig.2 Monte Carlo thermal contact resistance model

    蒙特卡羅法計(jì)算接觸熱阻的步驟如下:

    由于將兩粗糙表面接觸等效為一粗糙表面和一剛性平面接觸,因此需要計(jì)算等效后的接觸面的等效參數(shù)。由文獻(xiàn)[14]可知,各等效參數(shù)的計(jì)算方法如下:

    等效粗糙度:

    式中:σ、σ1、σ2分別表示接觸面等效粗糙度和上、下表面粗糙度,μm。

    等效斜率:

    式中:a、a1、a2分別表示接觸面粗糙峰等效斜率和上、下表面粗糙峰斜率。

    等效硬度:

    式中:H、H1、H2分別表示接觸面等效布氏硬度和上、下表面布氏硬度,×9.8 MPa。

    根據(jù)假設(shè),粗糙峰高度z的概率密度函數(shù)φ(z)服從正態(tài)分布:

    式中:位置參數(shù)μ為0,即圖2b 中的x軸,表示粗糙面輪廓高度均值。尺度參數(shù)σ即前述等效粗糙度,為采樣范圍內(nèi)粗糙峰高度z的均方差:

    式中:z為輪廓高度,m;m為輪廓高度均值,m。

    在實(shí)際應(yīng)用中,一般采用輪廓算數(shù)平均偏差Ra表示粗糙度:

    由文獻(xiàn)[15]的研究結(jié)果可知,在進(jìn)行數(shù)值模擬估算時(shí),σ可以近似為Ra的倍,如圖3 所示。

    圖3 σ 和Ra 隨磨損時(shí)間變化曲線[16]Fig.3 σ and Ra vs.wear time[16]

    用隨機(jī)方法模擬N個(gè)服從φ(z)分布的隨機(jī)數(shù)zi作為等效界面粗糙峰高度。其中,N可以用面積乘以粗糙峰密度得到:

    式中:Am為名義接觸面積,m2;η為粗糙峰密度,個(gè)/m2。

    由文獻(xiàn)[16]可知,粗糙峰密度可按照式(8)計(jì)算:

    式中:a為接觸面粗糙峰等效斜率;σ為等效粗糙度,μm。

    粗糙峰峰頂與剛性平面的距離δi為:

    式中:d為兩個(gè)接觸界面的間隙距離,m;zi為粗糙峰高度,m。

    當(dāng)δi>0 時(shí),此粗糙峰與剛性平面發(fā)生接觸,接觸區(qū)圓半徑ri和接觸面積Aci以及接觸載荷Fi分別為:

    式中:H為接觸面等效布氏硬度,×9.8 MPa。

    在獲得每個(gè)接觸點(diǎn)的幾何參數(shù)后,便可使用單點(diǎn)接觸模型計(jì)算每個(gè)接觸點(diǎn)的接觸熱阻。需要指出的是,圖2b 中的a就是上文中的接觸區(qū)圓半徑r。

    根據(jù)上述計(jì)算方法,采用MATLAB 編寫計(jì)算程序來獲取與實(shí)驗(yàn)工況對(duì)應(yīng)的接觸熱阻值。圖4a 為同一工況下采用圓柱模型、半球模型和圓錐模型3 種單點(diǎn)模型的計(jì)算結(jié)果比較,圖4b 為采用圓柱模型時(shí)3種不同壓力和粗糙度工況下接觸熱阻隨溫度的變化曲線。處于篇幅考慮,此處省略用圓柱模型、半球模型和圓錐模型3 種幾何模型的數(shù)學(xué)表征,詳細(xì)可參見文獻(xiàn)[4]、[6-7]。由圖4 可見:(1)在所涉及的10—30 K 溫區(qū)實(shí)驗(yàn)工況范圍內(nèi),使用蒙特卡羅法計(jì)算得到的結(jié)果在10-4—10-3量級(jí)間;(2)接觸熱阻始終隨溫度升高而減小,且溫度越高,減小幅度越小,最終趨于平緩;(3)同一單點(diǎn)模型下,壓力變化對(duì)接觸熱阻的影響不明顯,而粗糙度的影響比較顯著。更大的粗糙度意味著更小的實(shí)際接觸面積。壓力影響不明顯的原因可能是實(shí)驗(yàn)壓力取值范圍較小,對(duì)應(yīng)的界面形變不充分;(4)同一工況下,圓柱模型和半球模型的計(jì)算結(jié)果十分接近,而圓錐模型計(jì)算值則整體偏小1倍左右。

    圖4 蒙特卡羅法計(jì)算黃銅接觸熱阻變化曲線Fig.4 Brass thermal contact resistance change curves obtained by using Monte Carlo method

    3.2 積分法

    積分法和蒙特卡羅法類似,需要構(gòu)建粗糙峰高度正態(tài)分布概率密度函數(shù)。假設(shè)兩個(gè)接觸界面之間的間隙距離為d,當(dāng)粗糙峰高度z>d時(shí),該峰將和剛性平面發(fā)生接觸。根據(jù)上述概率密度函數(shù)可以積分計(jì)算粗糙峰發(fā)生接觸的概率,與名義接觸面積的乘積即為產(chǎn)生接觸的粗糙峰總數(shù),進(jìn)而計(jì)算總接觸面積和總接觸載荷。與蒙特卡羅法不同的是,積分法將所有接觸峰視為一致,通過總接觸面積除以接觸峰數(shù),獲得平均的接觸半徑和熱流通道半徑等幾何參數(shù)。通過單點(diǎn)接觸熱阻模型計(jì)算平均點(diǎn)接觸熱阻,并聯(lián)各點(diǎn)接觸熱阻即可得到總接觸熱阻。由于將所有接觸點(diǎn)視為同樣大小,該方法不具備一定的隨機(jī)性和普適性。圖5給出兩種數(shù)值計(jì)算方法在同一工況下隨溫度變化的接觸熱阻曲線??梢钥闯?積分法計(jì)算得到的接觸熱阻變化規(guī)律和蒙特卡羅法基本一致,但計(jì)算數(shù)值更大。由于蒙特卡羅法考慮到了每個(gè)粗糙峰之間的隨機(jī)性影響,因此認(rèn)為相對(duì)于積分法更加符合實(shí)際情況。但不管是何種數(shù)值計(jì)算方法它們都更適合用作初步預(yù)測(cè),無法準(zhǔn)確模擬低溫工況下各物理量的變化情況,因?yàn)樗鼈兾纯紤]溫度降低所帶來的影響。

    圖5 同一工況下蒙特卡羅法和積分法計(jì)算黃銅接觸熱阻結(jié)果比較Fig.5 Comparison between brass thermal contact resistance of Monte Carlo and integral method under same condition

    4 低溫接觸熱阻的有限元仿真

    選用COMSOL 軟件的固體傳熱物理場(chǎng)模塊進(jìn)行有限元仿真計(jì)算,所建立的模型和網(wǎng)格劃分如圖6 所示。該形狀結(jié)構(gòu)和實(shí)驗(yàn)器件保持一致,其中部分細(xì)節(jié)處如螺帽、螺紋等忽略。上下兩段工件和四周的絲桿選用材料庫內(nèi)置的304 不銹鋼材料,中間兩塊樣品為自行設(shè)置的黃銅材料,接觸面設(shè)置熱接觸物理場(chǎng)用以模擬接觸熱阻。由于材料的熱導(dǎo)率和溫度有關(guān),因此采用美國NIST 材料低溫物性庫給出的式(13)[17]計(jì)算,對(duì)于黃銅材料式中各項(xiàng)系數(shù)見表3。模型同時(shí)考慮了材料隨溫度變化的熱膨脹/收縮效應(yīng)。網(wǎng)格劃分采用物理場(chǎng)控制的自由四面體網(wǎng)格,并對(duì)絲桿和連接處進(jìn)行了細(xì)化處理。經(jīng)網(wǎng)格無關(guān)性驗(yàn)證后,選取的完整網(wǎng)格包含150 611 個(gè)域單元、39 710 個(gè)邊界元和4 084個(gè)邊單元。

    表3 黃銅熱導(dǎo)率計(jì)算公式系數(shù)Table 3 Brass thermal conductivity calculation formula coefficient

    圖6 有限元模型和網(wǎng)格劃分Fig.6 Finite element simulation model and model meshing

    該接觸熱阻模擬主要涉及以下多個(gè)影響因素:溫度、氣隙熱導(dǎo)率、輻射傳導(dǎo)率、表面粗糙度、表面粗糙峰斜率、接觸壓力、布氏硬度、表面輻射率以及熱摩擦。需要說明的是,氣隙熱導(dǎo)率一般在常溫常壓實(shí)驗(yàn)或者各類特殊氣體環(huán)境中需要考慮,但本研究計(jì)算模型所引用的實(shí)驗(yàn)不僅是低至5 K 的低溫環(huán)境,而且處于10-4量級(jí)的高真空狀態(tài),因此氣隙熱傳導(dǎo)可以忽略不計(jì)。輻射傳導(dǎo)率設(shè)定為灰體漫射平行面,表面輻射率來自材料物性。熱摩擦采用默認(rèn)的Charron 關(guān)系。選用的樣品材料分別為黃銅、304 不銹鋼以及紫銅,各邊界條件參數(shù)如表4 所示,其中粗糙峰斜率和布氏硬度均為查閱資料選取的較為合適的值。使用參數(shù)化掃描對(duì)所有組合情況進(jìn)行計(jì)算。

    表4 邊界條件參數(shù)設(shè)置Table 4 Simulation parameter setting after pressure correction

    圖7 為放大顯示的黃銅樣品部位溫度分布云圖。可以發(fā)現(xiàn),當(dāng)存在接觸熱阻時(shí),樣品的接觸表面上下會(huì)形成明顯的溫度階躍。當(dāng)熱流方向?yàn)樽陨隙聲r(shí),頂部的夾持工件(上熱流計(jì))溫度整體明顯高于底部工件(下熱流計(jì)),引起溫度變化的原因主要發(fā)生在樣品接觸面上。

    圖7 樣品部位局部溫度分布Fig.7 Local temperature distribution at sample site

    使用COMSOL 的派生值功能計(jì)算接觸表面的平均接觸熱導(dǎo)(W/(m2·K)),即接觸熱阻的倒數(shù)。在研究溫度、粗糙度以及壓力等影響因素時(shí),選定接觸面硬度為140 ×9.8 MPa 和粗糙峰斜率為。據(jù)此,所研究工況范圍內(nèi)的黃銅接觸熱阻計(jì)算值在1.07 ×10-3—1.31 ×10-2m2·K/W 之間。圖8 給出了同一工況下有限元仿真結(jié)果和前述兩種數(shù)值計(jì)算方法結(jié)果的對(duì)比。不難看出,不同方法計(jì)算的接觸熱阻隨溫度的變化規(guī)律基本一致,但相較于數(shù)值計(jì)算方法,有限元仿真得出的結(jié)果要大得多,與試驗(yàn)數(shù)據(jù)的比較情況將在后文闡述。

    圖8 黃銅接觸熱阻數(shù)值計(jì)算和有限元仿真結(jié)果比較Fig.8 Comparison of numerical calculation and simulation results of brass thermal contact resistance

    通過仿真模型具體分析各個(gè)因素對(duì)接觸熱阻的影響。圖9a 為保持k=,P=0.288 MPa,Hb=140 ×9.8 MPa 不變,不同粗糙度下接觸熱阻隨溫度的變化情況;圖9b 為在此基礎(chǔ)之上Thot=40 K 時(shí)接觸熱阻隨粗糙度的直接變化曲線??梢钥闯?(1) 接觸熱阻隨著溫度的升高而減小,這是因?yàn)樵诘蜏卣婵窄h(huán)境下,固體材料間的接觸熱阻主要和材料自身熱導(dǎo)率、硬度、彈性模量等物理性質(zhì)相關(guān)。隨著溫度升高,黃銅的熱導(dǎo)率增大,同時(shí)硬度和彈性模量減小,樣品自身也會(huì)產(chǎn)生一定熱膨脹,使得實(shí)際接觸面積增大,在兩者的共同作用下,接觸熱阻顯著減小。(2) 接觸熱阻隨粗糙度增大而增大,在1.6—6.3 μm 范圍內(nèi)隨粗糙度呈現(xiàn)近似線性變化趨勢(shì)。這是因?yàn)殡S著粗糙度的增大,接觸界面越發(fā)不平整,實(shí)際接觸面積越小,熱流通道面積越小,從而導(dǎo)致傳熱能力降低,接觸熱阻相應(yīng)增大。Ra6.3 粗糙度下的接觸熱阻比Ra1.6粗糙度高出近2 倍,可見粗糙度對(duì)接觸熱阻影響程度顯著,但仍不及溫度的影響。

    圖9 溫度和粗糙度對(duì)黃銅接觸熱阻的影響Fig.9 Influence of T and σ on brass thermal contact resistance in simulation

    圖10a 為保持k=,σ=3.2 μm,Hb=140 ×9.8 MPa 不變,不同壓力下接觸熱阻隨溫度的變化情況,圖10b 為在此基礎(chǔ)之上Thot=40 K 時(shí)接觸熱阻隨壓力的直接變化曲線??梢园l(fā)現(xiàn):(1) 隨著壓力改變,接觸熱阻隨溫度的變化趨勢(shì)也始終保持一致;(2) 接觸熱阻隨壓力的增大而減小,這是因?yàn)榻佑|界面受到壓力越大,擠壓產(chǎn)生的形變?cè)矫黠@,實(shí)際接觸面積越大,但其變化幅度相較于粗糙度要小很多。同時(shí),其減小速率隨著壓力的增大逐漸減緩,這是因?yàn)殡S著壓力的增大,實(shí)際接觸面積越來越大,可供接觸面形變的空間也越來越少。在本研究引述的實(shí)驗(yàn)工況范圍內(nèi),壓力引起的變化幅度也小于溫度。

    圖10 溫度和壓力對(duì)黃銅接觸熱阻的影響Fig.10 Influence of T and P on brass thermal contact resistance in simulation

    以上研究的樣品均為黃銅材料,為了探究比較不同材料的接觸熱阻特性,還對(duì)304 不銹鋼和紫銅材料進(jìn)行了模擬計(jì)算。保持其他因素一致,即控制k==3.2 μm,P=0.288 MPa 不變。表5 給出了3 種材料的接觸熱阻計(jì)算值??梢园l(fā)現(xiàn),在同一工況下,304 不銹鋼的接觸熱阻值是黃銅的數(shù)倍之多。此外,紫銅的接觸熱阻值趨近于0,遠(yuǎn)小于前兩種材料,且接觸熱阻值幾乎不隨溫度產(chǎn)生變化。綜合比對(duì)3 種材料,認(rèn)為材料自身的熱導(dǎo)率在其中起到了關(guān)鍵作用。在10—30 K 溫度范圍內(nèi),黃銅的熱導(dǎo)率是304 不銹鋼的5—6 倍,而相應(yīng)地,仿真計(jì)算得到的黃銅的接觸熱阻正好是304 不銹鋼的1/5—1/6。而紫銅的熱導(dǎo)率要遠(yuǎn)高于不銹鋼和黃銅,且在20 K 左右其數(shù)值存在一個(gè)峰值,接近1 500 W/(m·K),因此在該溫區(qū)內(nèi)測(cè)量紫銅材料的接觸熱阻遠(yuǎn)比其它材料要困難。

    表5 3 種典型金屬材料接觸熱阻有限元仿真結(jié)果對(duì)比Table 5 Simulation comparison of thermal contact resistances of three materials

    5 低溫接觸熱阻的經(jīng)驗(yàn)關(guān)聯(lián)式建立

    有限元仿真建模計(jì)算能夠給出較為準(zhǔn)確的數(shù)據(jù)并且使得研究不同因素對(duì)接觸熱阻的影響成為可能,但其缺點(diǎn)是建模和計(jì)算費(fèi)時(shí)費(fèi)力,不利于工程實(shí)際應(yīng)用。為了更加方便地獲得接觸熱阻數(shù)據(jù),經(jīng)驗(yàn)關(guān)聯(lián)式是一種受歡迎的表達(dá)形式。下文將給出黃銅材料在10—30 K 溫區(qū)內(nèi)的接觸熱阻計(jì)算關(guān)聯(lián)式。

    從仿真結(jié)果發(fā)現(xiàn),接觸熱阻隨溫度的變化較好地遵循Rc=y0+A1e-T/t1規(guī)律,且其3 個(gè)參數(shù)y0,A1及t1中,無論工況如何變化,t1基本保持不變,t1的取值越大,曲線越平緩。其它2 個(gè)參數(shù)y0和A1均隨著各個(gè)影響因素的變化而改變。y0為T趨向于無窮大時(shí)Rc最終趨近于的數(shù)值;A1為曲線下降的速率,A1越大,接觸熱阻減小得越快。

    結(jié)合仿真數(shù)據(jù)和實(shí)驗(yàn)結(jié)果的規(guī)律分析,參數(shù)A1和t1均為常數(shù),只需確定y0的表達(dá)式。利用仿真數(shù)據(jù)分別繪制y0關(guān)于4 個(gè)影響因素的圖像,并分別擬合得到式(14)—(17)。除了粗糙度σ為線性表達(dá)式外,其余3 個(gè)變量均使用y0=y+A1e-x/t函數(shù),擬合效果如圖11 所示。

    圖11 y0 和各變量擬合關(guān)系Fig.11 Fitting relationship between y0 and each variable

    代入其它數(shù)據(jù)確定并修正未知參數(shù),最終得到的黃銅接觸熱阻計(jì)算關(guān)聯(lián)式如下:

    式中:T為接觸面溫度,K;σ為接觸面粗糙度,μm;P為載荷壓力,MPa;Hb為接觸面布氏硬度,×9.8 MPa;a為接觸面粗糙峰斜率。

    使用該關(guān)聯(lián)式可以計(jì)算各實(shí)驗(yàn)工況下黃銅材料的接觸熱阻值,實(shí)驗(yàn)和關(guān)聯(lián)式計(jì)算結(jié)果比對(duì)如圖12 所示。兩者誤差最大偏差24.15%,平均偏差11.41%,基本滿足該溫區(qū)內(nèi)黃銅低溫接觸熱阻的工程設(shè)計(jì)需求。

    圖12 黃銅接觸熱阻實(shí)驗(yàn)和關(guān)聯(lián)式計(jì)算結(jié)果對(duì)比Fig.12 Comparison of experimental and correlation calculation results of brass thermal contact resistance

    式(18)適用于黃銅材料在10—30 K 溫區(qū)內(nèi)的接觸熱阻預(yù)測(cè)。根據(jù)上文所述,在相同的溫度和壓力工況下,粗糙度和硬度接近的材料,其接觸熱阻和熱導(dǎo)率基本成反比關(guān)系。常溫下黃銅的布氏硬度為100—180 ×9.8 MPa,而304 不銹鋼的布氏硬度一般略小于187 ×9.8 MPa,兩者接近。因此,在黃銅接觸熱阻關(guān)聯(lián)式的基礎(chǔ)上,對(duì)于粗糙度相同、硬度接近的304 不銹鋼,其接觸熱阻可以通過下式進(jìn)行間接計(jì)算:

    式中:R,k,Rb,kb分別為目標(biāo)材料和黃銅的接觸熱阻和熱導(dǎo)率。根據(jù)美國NIST 數(shù)據(jù)庫給出的兩種材料的熱導(dǎo)率計(jì)算公式,可以進(jìn)一步得到熱導(dǎo)率比值和溫度之間的關(guān)系,如式(20)所示。

    因此,用于計(jì)算10—30 K 溫區(qū)內(nèi)304 不銹鋼接觸熱阻的經(jīng)驗(yàn)關(guān)聯(lián)式如下:

    圖13 304 不銹鋼接觸熱阻實(shí)驗(yàn)測(cè)量和關(guān)聯(lián)式計(jì)算結(jié)果對(duì)比Fig.13 Comparison of experimental measurement and correlation calculation results of 304 stainless steel thermal contact resistance

    6 結(jié)論

    本研究使用蒙特卡羅法、積分?jǐn)?shù)值方法和有限元仿真方法分別計(jì)算了不同材料在10—30 K 溫區(qū)內(nèi)的接觸熱阻數(shù)值,比較和分析了它們的適用條件和相互差異?;谟邢拊ǐ@得的仿真數(shù)據(jù)和分離變量法推導(dǎo),給出了適用于該溫區(qū)黃銅和304 不銹鋼的接觸熱阻經(jīng)驗(yàn)關(guān)聯(lián)式,并經(jīng)試驗(yàn)驗(yàn)證。得到如下結(jié)論:

    (1)蒙特卡羅法和積分?jǐn)?shù)值方法之間的接觸熱阻計(jì)算值偏差大,僅可用于初步的估算,并不能準(zhǔn)確地計(jì)算低溫接觸熱阻。

    (2)相較于常規(guī)數(shù)值計(jì)算,有限元仿真可以考慮接觸壓力、材料物性、硬度和熱輻射等作用因素,得出的接觸熱阻值更接近實(shí)際測(cè)量結(jié)果。

    (3)所給出的黃銅接觸熱阻計(jì)算經(jīng)驗(yàn)關(guān)聯(lián)式在多種工況下的最大誤差和平均誤差分別為24.15%和11.41%,基本滿足該溫區(qū)內(nèi)黃銅接觸熱阻的預(yù)測(cè)需求。

    (4)真空低溫環(huán)境下,硬度、粗糙度等機(jī)械性質(zhì)接近的材料,其接觸熱阻與熱導(dǎo)率近似成反比關(guān)系。基于此得出的304 不銹鋼在同溫區(qū)內(nèi)的接觸熱阻經(jīng)驗(yàn)關(guān)聯(lián)式與實(shí)測(cè)數(shù)據(jù)最大偏差小于9.95%。

    猜你喜歡
    黃銅蒙特卡羅熱導(dǎo)率
    空位缺陷對(duì)單層石墨烯導(dǎo)熱特性影響的分子動(dòng)力學(xué)
    連續(xù)碳纖維鋁基復(fù)合材料橫向等效熱導(dǎo)率的模擬分析
    Si3N4/BN復(fù)合陶瓷熱導(dǎo)率及其有限元分析
    利用蒙特卡羅方法求解二重積分
    黃銅容器
    黃銅筆頭的一場(chǎng)革命
    中國制筆(2017年2期)2017-07-18 10:53:08
    用于鑄造電氣元件的鋁黃銅新材料研究
    熱鍛黃銅閥門的缺陷和失效研究
    金屬熱導(dǎo)率的第一性原理計(jì)算方法在鋁中的應(yīng)用
    探討蒙特卡羅方法在解微分方程邊值問題中的應(yīng)用
    成人毛片a级毛片在线播放| 欧美极品一区二区三区四区| 国产av一区在线观看免费| 成人美女网站在线观看视频| 在线观看66精品国产| 欧美+亚洲+日韩+国产| 国产精品国产三级国产av玫瑰| 午夜亚洲福利在线播放| 亚洲人成伊人成综合网2020| 看片在线看免费视频| 久久香蕉精品热| 美女免费视频网站| 国产伦人伦偷精品视频| 干丝袜人妻中文字幕| 成人综合一区亚洲| 精品不卡国产一区二区三区| 久久精品夜夜夜夜夜久久蜜豆| 91久久精品国产一区二区三区| 禁无遮挡网站| 国产色爽女视频免费观看| 久久久久国产精品人妻aⅴ院| 亚洲精品亚洲一区二区| 色综合站精品国产| 乱系列少妇在线播放| 国产色爽女视频免费观看| 亚洲色图av天堂| 国产精品福利在线免费观看| 乱人视频在线观看| 国产成人av教育| 国产成人a区在线观看| 亚洲精品一卡2卡三卡4卡5卡| 国产不卡一卡二| 欧美一区二区亚洲| 色视频www国产| 国产伦精品一区二区三区四那| 国产av麻豆久久久久久久| 国产精品伦人一区二区| 99热网站在线观看| 哪里可以看免费的av片| 色综合站精品国产| 欧美一区二区亚洲| 一个人免费在线观看电影| 亚洲欧美清纯卡通| 干丝袜人妻中文字幕| 99热精品在线国产| 亚洲国产精品合色在线| 国产乱人视频| 成人三级黄色视频| 免费不卡的大黄色大毛片视频在线观看 | 亚洲精品亚洲一区二区| 国内少妇人妻偷人精品xxx网站| 又紧又爽又黄一区二区| 亚洲最大成人中文| 日本色播在线视频| 99久久无色码亚洲精品果冻| 99久久成人亚洲精品观看| 色综合色国产| 日韩欧美一区二区三区在线观看| 18禁裸乳无遮挡免费网站照片| 欧美国产日韩亚洲一区| 国产熟女欧美一区二区| 亚州av有码| 国产成人影院久久av| 久久精品国产鲁丝片午夜精品 | 亚洲五月天丁香| 成熟少妇高潮喷水视频| 一级黄色大片毛片| 深夜精品福利| 最近在线观看免费完整版| 国产综合懂色| 欧美日韩综合久久久久久 | 精品无人区乱码1区二区| 午夜视频国产福利| 色在线成人网| 欧美一区二区亚洲| 久久精品国产99精品国产亚洲性色| 国产亚洲精品综合一区在线观看| 深夜a级毛片| 深爱激情五月婷婷| 国产精品久久久久久久久免| 亚洲性久久影院| 自拍偷自拍亚洲精品老妇| 亚洲精品国产成人久久av| 国产美女午夜福利| 伦精品一区二区三区| 午夜免费成人在线视频| 一个人免费在线观看电影| 国产精品伦人一区二区| 在线观看66精品国产| 久久久成人免费电影| 国产精品久久久久久av不卡| 亚洲在线自拍视频| 日韩欧美国产一区二区入口| 亚洲av一区综合| 免费人成在线观看视频色| 少妇高潮的动态图| 九九久久精品国产亚洲av麻豆| 国产精品久久久久久av不卡| 国产在线精品亚洲第一网站| 22中文网久久字幕| 琪琪午夜伦伦电影理论片6080| 国产淫片久久久久久久久| 成人一区二区视频在线观看| 自拍偷自拍亚洲精品老妇| 超碰av人人做人人爽久久| 免费看av在线观看网站| .国产精品久久| 草草在线视频免费看| 国产精品久久久久久av不卡| 黄色丝袜av网址大全| 少妇高潮的动态图| 在线天堂最新版资源| 天天躁日日操中文字幕| 亚洲国产精品成人综合色| 亚洲av熟女| 亚洲久久久久久中文字幕| 国产精品98久久久久久宅男小说| 天堂√8在线中文| 亚洲成人久久爱视频| 亚洲中文日韩欧美视频| 日韩中字成人| 91久久精品国产一区二区成人| 欧美极品一区二区三区四区| 国产一级毛片七仙女欲春2| 欧美精品国产亚洲| 99精品久久久久人妻精品| 波多野结衣巨乳人妻| 午夜免费成人在线视频| 欧美xxxx黑人xx丫x性爽| 久久国产精品人妻蜜桃| 国产精品久久久久久亚洲av鲁大| 久久久久久久精品吃奶| 国产精品永久免费网站| av在线蜜桃| 国产一区二区激情短视频| 在线看三级毛片| 久久久久久久久久成人| 国产精品av视频在线免费观看| 中国美女看黄片| 九色国产91popny在线| 三级国产精品欧美在线观看| 制服丝袜大香蕉在线| 欧美成人免费av一区二区三区| 婷婷精品国产亚洲av| 丰满人妻一区二区三区视频av| 久久午夜亚洲精品久久| 国产伦在线观看视频一区| 日本色播在线视频| 成人无遮挡网站| 亚洲人成网站高清观看| 自拍偷自拍亚洲精品老妇| 日日干狠狠操夜夜爽| 久久精品国产鲁丝片午夜精品 | 嫩草影院入口| 成人一区二区视频在线观看| 欧美bdsm另类| 日本一二三区视频观看| 国产亚洲精品综合一区在线观看| 久久久久久伊人网av| 亚洲欧美日韩高清专用| 亚洲国产精品合色在线| 亚洲精品粉嫩美女一区| 婷婷丁香在线五月| 1000部很黄的大片| 日日干狠狠操夜夜爽| 天堂影院成人在线观看| 午夜福利在线在线| 一区二区三区四区激情视频 | 欧美3d第一页| 国产乱人伦免费视频| 久久久久久久久久黄片| 久久精品国产99精品国产亚洲性色| 黄色日韩在线| 十八禁网站免费在线| 在线看三级毛片| 伦理电影大哥的女人| 久久久久久久久大av| 麻豆国产av国片精品| 中国美白少妇内射xxxbb| 黄片wwwwww| 国模一区二区三区四区视频| 国产在线男女| 日韩欧美精品v在线| 久久6这里有精品| 别揉我奶头 嗯啊视频| 欧美+日韩+精品| 成人国产综合亚洲| 高清毛片免费观看视频网站| 久久亚洲真实| 国产69精品久久久久777片| 午夜日韩欧美国产| 亚洲成人久久性| 天堂动漫精品| 日韩高清综合在线| 中出人妻视频一区二区| 波多野结衣巨乳人妻| 亚洲电影在线观看av| 黄片wwwwww| 日本免费a在线| 俺也久久电影网| a在线观看视频网站| 亚洲成a人片在线一区二区| 亚洲在线自拍视频| 免费无遮挡裸体视频| 成人国产一区最新在线观看| 国产精品av视频在线免费观看| 嫩草影院精品99| 男人和女人高潮做爰伦理| 欧美高清性xxxxhd video| 国产午夜精品论理片| 久久精品人妻少妇| 国产精品,欧美在线| 少妇猛男粗大的猛烈进出视频 | 亚洲美女搞黄在线观看 | 高清日韩中文字幕在线| 人妻少妇偷人精品九色| 精品不卡国产一区二区三区| 淫秽高清视频在线观看| 99热这里只有是精品50| 无遮挡黄片免费观看| 午夜福利视频1000在线观看| 中文资源天堂在线| 亚洲熟妇熟女久久| 亚洲第一电影网av| 精品无人区乱码1区二区| 日韩国内少妇激情av| 天堂动漫精品| 别揉我奶头 嗯啊视频| 精品人妻偷拍中文字幕| 日韩高清综合在线| 黄色一级大片看看| 久久精品综合一区二区三区| 国产高清三级在线| 午夜福利欧美成人| 在线看三级毛片| 国产熟女欧美一区二区| 欧美zozozo另类| 国产精品野战在线观看| 黄色女人牲交| 国产v大片淫在线免费观看| 精品一区二区三区人妻视频| 久久午夜福利片| 在线播放无遮挡| 免费电影在线观看免费观看| 男人的好看免费观看在线视频| 国产亚洲欧美98| 国国产精品蜜臀av免费| videossex国产| 亚洲美女搞黄在线观看 | 女人被狂操c到高潮| 十八禁国产超污无遮挡网站| 国产白丝娇喘喷水9色精品| 国产蜜桃级精品一区二区三区| 亚洲av熟女| 简卡轻食公司| 久久精品国产清高在天天线| 69av精品久久久久久| 男人狂女人下面高潮的视频| 午夜日韩欧美国产| 婷婷精品国产亚洲av在线| 亚洲精品粉嫩美女一区| 嫁个100分男人电影在线观看| 桃色一区二区三区在线观看| 国产激情偷乱视频一区二区| 午夜福利高清视频| 成人特级av手机在线观看| 在线免费观看不下载黄p国产 | 91午夜精品亚洲一区二区三区 | 精品久久久久久久久久免费视频| 国产熟女欧美一区二区| 日日撸夜夜添| 俺也久久电影网| 男女那种视频在线观看| 女同久久另类99精品国产91| 亚洲av成人av| 春色校园在线视频观看| 婷婷精品国产亚洲av| 国产一区二区亚洲精品在线观看| 亚洲不卡免费看| 最好的美女福利视频网| 午夜激情福利司机影院| 一进一出抽搐gif免费好疼| 女人十人毛片免费观看3o分钟| 亚洲av中文字字幕乱码综合| 九九久久精品国产亚洲av麻豆| 免费在线观看成人毛片| 在线观看av片永久免费下载| 99久久精品国产国产毛片| 久久精品国产亚洲av涩爱 | 欧美bdsm另类| 51国产日韩欧美| 一个人看的www免费观看视频| 久久久精品欧美日韩精品| 亚洲人成网站高清观看| 免费电影在线观看免费观看| 久久精品国产鲁丝片午夜精品 | 午夜久久久久精精品| 亚洲成人精品中文字幕电影| 国产视频内射| 神马国产精品三级电影在线观看| 在线免费十八禁| 变态另类成人亚洲欧美熟女| 丝袜美腿在线中文| 看免费成人av毛片| 我要搜黄色片| 91狼人影院| 国产一区二区亚洲精品在线观看| 欧美最黄视频在线播放免费| 五月玫瑰六月丁香| 日本一本二区三区精品| 国产精品免费一区二区三区在线| 日本在线视频免费播放| 草草在线视频免费看| 久久精品国产99精品国产亚洲性色| 日本免费一区二区三区高清不卡| 日本欧美国产在线视频| 亚洲第一区二区三区不卡| 性插视频无遮挡在线免费观看| 国产精品99久久久久久久久| 欧美一区二区亚洲| 国产伦一二天堂av在线观看| 听说在线观看完整版免费高清| 国内精品久久久久精免费| 男女做爰动态图高潮gif福利片| 成年免费大片在线观看| 日本与韩国留学比较| 国产淫片久久久久久久久| 久久香蕉精品热| 毛片一级片免费看久久久久 | 床上黄色一级片| 男女那种视频在线观看| 亚洲av免费在线观看| 久久久久久久久久久丰满 | 中文字幕精品亚洲无线码一区| 免费观看在线日韩| 国产精华一区二区三区| 成年女人看的毛片在线观看| 高清毛片免费观看视频网站| 中文亚洲av片在线观看爽| 亚洲黑人精品在线| 国产精品国产高清国产av| 一个人免费在线观看电影| 亚洲精华国产精华液的使用体验 | 久久亚洲精品不卡| 内地一区二区视频在线| 久久亚洲精品不卡| 欧美国产日韩亚洲一区| 免费av毛片视频| 亚洲男人的天堂狠狠| 五月伊人婷婷丁香| 精品99又大又爽又粗少妇毛片 | 老师上课跳d突然被开到最大视频| 午夜福利欧美成人| 国产精品久久电影中文字幕| 男女做爰动态图高潮gif福利片| 精品日产1卡2卡| 2021天堂中文幕一二区在线观| 亚洲成人久久爱视频| 两个人的视频大全免费| 能在线免费观看的黄片| 又粗又爽又猛毛片免费看| 1024手机看黄色片| 午夜福利高清视频| 99热6这里只有精品| 在线观看av片永久免费下载| 在线免费观看不下载黄p国产 | 女的被弄到高潮叫床怎么办 | 中文在线观看免费www的网站| 婷婷六月久久综合丁香| 国内久久婷婷六月综合欲色啪| 在线观看一区二区三区| 国产精品永久免费网站| 色吧在线观看| 天堂网av新在线| 日本精品一区二区三区蜜桃| 天堂网av新在线| 日本精品一区二区三区蜜桃| 成人一区二区视频在线观看| 久久99热这里只有精品18| 精品日产1卡2卡| 欧美日韩精品成人综合77777| 男人狂女人下面高潮的视频| 午夜福利欧美成人| 精品99又大又爽又粗少妇毛片 | 亚洲欧美日韩高清在线视频| 久久精品影院6| 啦啦啦观看免费观看视频高清| 国产高清视频在线观看网站| 啦啦啦观看免费观看视频高清| 十八禁网站免费在线| 天堂网av新在线| 波多野结衣高清作品| 精品久久久噜噜| 欧美精品国产亚洲| 大又大粗又爽又黄少妇毛片口| 久久中文看片网| 亚洲天堂国产精品一区在线| 日韩欧美三级三区| 国产又黄又爽又无遮挡在线| 美女高潮喷水抽搐中文字幕| 国产精品电影一区二区三区| 女生性感内裤真人,穿戴方法视频| 日韩一区二区视频免费看| 久久久久久久久久黄片| 美女被艹到高潮喷水动态| 男女边吃奶边做爰视频| 国产精品野战在线观看| 国产精品精品国产色婷婷| 久久久久免费精品人妻一区二区| 欧美日韩乱码在线| 黄色一级大片看看| 欧美zozozo另类| 黄色一级大片看看| 日韩高清综合在线| 免费在线观看影片大全网站| 午夜激情欧美在线| 午夜影院日韩av| 国产v大片淫在线免费观看| 久久精品国产自在天天线| 欧美高清性xxxxhd video| 日韩欧美一区二区三区在线观看| 99热这里只有是精品50| 欧美日韩中文字幕国产精品一区二区三区| 色5月婷婷丁香| 日本a在线网址| 白带黄色成豆腐渣| 美女cb高潮喷水在线观看| 身体一侧抽搐| 三级毛片av免费| 久久久午夜欧美精品| 啦啦啦啦在线视频资源| 成人国产麻豆网| 一本精品99久久精品77| 国产精品免费一区二区三区在线| 欧美潮喷喷水| 中文字幕熟女人妻在线| 国产精品亚洲美女久久久| 欧美成人性av电影在线观看| 在线观看66精品国产| 国产亚洲av嫩草精品影院| 91久久精品电影网| 色视频www国产| 亚洲精品成人久久久久久| 校园春色视频在线观看| 久久久久久伊人网av| 欧美一级a爱片免费观看看| 国产精品日韩av在线免费观看| 色综合婷婷激情| 俄罗斯特黄特色一大片| 国产白丝娇喘喷水9色精品| 午夜福利在线观看吧| 午夜福利在线观看免费完整高清在 | 亚洲真实伦在线观看| 最新在线观看一区二区三区| 少妇丰满av| 乱码一卡2卡4卡精品| 长腿黑丝高跟| 色5月婷婷丁香| 免费av毛片视频| 亚洲第一电影网av| 亚洲内射少妇av| 国产精品不卡视频一区二区| 欧美在线一区亚洲| 国产精品久久久久久久电影| 国产成人aa在线观看| 国产三级在线视频| 蜜桃久久精品国产亚洲av| 亚洲一级一片aⅴ在线观看| 成人综合一区亚洲| 欧美一区二区亚洲| 欧美最新免费一区二区三区| 国产淫片久久久久久久久| 亚洲精品在线观看二区| 欧美极品一区二区三区四区| 免费不卡的大黄色大毛片视频在线观看 | 亚洲av电影不卡..在线观看| 国产精品久久久久久精品电影| 人妻丰满熟妇av一区二区三区| 色综合站精品国产| 不卡视频在线观看欧美| 色精品久久人妻99蜜桃| 悠悠久久av| 亚洲人与动物交配视频| 亚洲精品日韩av片在线观看| 在线播放国产精品三级| 国产一区二区亚洲精品在线观看| 天天一区二区日本电影三级| 一区二区三区免费毛片| 国产高清视频在线播放一区| 国产亚洲欧美98| 最新中文字幕久久久久| 啦啦啦观看免费观看视频高清| 亚洲性夜色夜夜综合| 又黄又爽又免费观看的视频| 毛片女人毛片| 99国产极品粉嫩在线观看| 精品99又大又爽又粗少妇毛片 | 精品久久久久久久人妻蜜臀av| 国产伦人伦偷精品视频| 夜夜看夜夜爽夜夜摸| 99热6这里只有精品| 少妇丰满av| 天堂网av新在线| 欧美性感艳星| 女同久久另类99精品国产91| 热99re8久久精品国产| 国产不卡一卡二| 午夜精品一区二区三区免费看| 99九九线精品视频在线观看视频| 又紧又爽又黄一区二区| 美女 人体艺术 gogo| 美女免费视频网站| 国产精品一区二区免费欧美| 亚洲狠狠婷婷综合久久图片| 欧美最新免费一区二区三区| 一进一出抽搐动态| 男人舔女人下体高潮全视频| 啦啦啦韩国在线观看视频| 校园人妻丝袜中文字幕| 国产极品精品免费视频能看的| 99热这里只有是精品在线观看| 久久99热这里只有精品18| 免费一级毛片在线播放高清视频| 69av精品久久久久久| 两个人的视频大全免费| 如何舔出高潮| 国产欧美日韩精品亚洲av| 熟妇人妻久久中文字幕3abv| 春色校园在线视频观看| ponron亚洲| 精品无人区乱码1区二区| 国产精品久久久久久亚洲av鲁大| 在线观看美女被高潮喷水网站| 国产亚洲av嫩草精品影院| 国产精品98久久久久久宅男小说| 一区二区三区四区激情视频 | 在线观看美女被高潮喷水网站| 国产午夜精品论理片| 久久精品影院6| 国产久久久一区二区三区| 最新中文字幕久久久久| 国产在视频线在精品| 欧美精品国产亚洲| 天堂网av新在线| 丰满人妻一区二区三区视频av| 九色国产91popny在线| 草草在线视频免费看| 成人亚洲精品av一区二区| eeuss影院久久| 人妻久久中文字幕网| 亚州av有码| 国产免费男女视频| 真人一进一出gif抽搐免费| 亚洲国产精品sss在线观看| 久久久久久久亚洲中文字幕| 午夜亚洲福利在线播放| 国产精品一区二区性色av| 中出人妻视频一区二区| 亚洲熟妇熟女久久| 免费无遮挡裸体视频| 狂野欧美白嫩少妇大欣赏| 99久久精品热视频| 国产精品不卡视频一区二区| 午夜免费激情av| 国产真实伦视频高清在线观看 | 国产精品精品国产色婷婷| 熟女电影av网| 欧美3d第一页| 十八禁网站免费在线| 亚洲国产精品合色在线| 亚洲成人久久爱视频| 搡女人真爽免费视频火全软件 | 亚洲18禁久久av| 熟妇人妻久久中文字幕3abv| 久久久色成人| 五月伊人婷婷丁香| 免费一级毛片在线播放高清视频| 亚洲五月天丁香| 婷婷精品国产亚洲av在线| 国产 一区 欧美 日韩| 欧美日韩乱码在线| 国产女主播在线喷水免费视频网站 | 熟女人妻精品中文字幕| 国产伦精品一区二区三区视频9| 亚洲av中文字字幕乱码综合| 欧美成人一区二区免费高清观看| 久久久午夜欧美精品| 久久人妻av系列| 69人妻影院| 在线观看av片永久免费下载| 尤物成人国产欧美一区二区三区| 欧美性感艳星| 精品一区二区三区人妻视频| 亚洲久久久久久中文字幕| or卡值多少钱| 亚洲第一区二区三区不卡| 夜夜夜夜夜久久久久| 国产男人的电影天堂91| 亚洲色图av天堂| 在线看三级毛片| 亚洲av.av天堂| 国产精品99久久久久久久久| 精品一区二区三区视频在线| 国产主播在线观看一区二区| 亚洲午夜理论影院| 一本一本综合久久| eeuss影院久久| 午夜日韩欧美国产| 国内精品宾馆在线| 色噜噜av男人的天堂激情| 国产精品人妻久久久久久| 久久国内精品自在自线图片| 一级a爱片免费观看的视频| 麻豆一二三区av精品| 桃红色精品国产亚洲av|