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

    誘餌球鏡面反射內(nèi)壁的輻射傳遞系數(shù)計(jì)算*

    2017-06-27 08:14:35厲夫兵冷俊敏李紅蓮
    現(xiàn)代防御技術(shù) 2017年3期
    關(guān)鍵詞:鏡面反射面元球心

    厲夫兵,冷俊敏,李紅蓮

    (北京信息科技大學(xué) 信息與通信工程學(xué)院,北京 100101)

    誘餌球鏡面反射內(nèi)壁的輻射傳遞系數(shù)計(jì)算*

    厲夫兵,冷俊敏,李紅蓮

    (北京信息科技大學(xué) 信息與通信工程學(xué)院,北京 100101)

    針對(duì)內(nèi)表面具備鏡面反射特性的空間誘餌球,基于蒙特卡羅方法建立漫發(fā)射模型,推導(dǎo)光線(xiàn)的最大反射次數(shù),利用光線(xiàn)在球腔的傳播特性求解光線(xiàn)與球壁的多次反射交點(diǎn),計(jì)算球面微元吸收的光束能量,獲取不同反射率情況下單一球面微元對(duì)全部球面微元的輻射傳遞系數(shù)。計(jì)算結(jié)果表明,單一微元對(duì)球面各微元的輻射傳遞系數(shù)中存在2個(gè)峰值,分別為微元對(duì)其自身的輻射傳遞系數(shù),以及微元對(duì)其球?qū)ΨQ(chēng)微元的輻射傳遞系數(shù)。隨著反射率的增大,微元對(duì)其自身的輻射傳遞系數(shù)逐漸接近于微元對(duì)其球?qū)ΨQ(chēng)微元的輻射傳遞系數(shù);當(dāng)反射率接近于1時(shí),兩者的值近似相等。

    輻射傳遞系數(shù);輻射換熱;鏡面反射;蒙特卡羅;角系數(shù);球

    0 引言

    在彈道目標(biāo)飛行中段釋放氣球誘餌的方法可以降低目標(biāo)的探測(cè)概率[1-2],增加天基紅外監(jiān)視系統(tǒng)對(duì)目標(biāo)的識(shí)別與跟蹤難度[3-4]。公開(kāi)文獻(xiàn)表明,空間氣球多采用鋁膜-聚酯薄膜-鋁膜的超薄層結(jié)構(gòu)[5]。在日光照射下,氣球向陽(yáng)面與背陰面的溫差會(huì)加劇氣球內(nèi)部的輻射換熱,而輻射換熱又反過(guò)來(lái)影響著目標(biāo)的溫度分布。因此,提高輻射換熱的計(jì)算精度對(duì)提高目標(biāo)溫度的計(jì)算準(zhǔn)確度具有重要意義。

    近年來(lái),國(guó)內(nèi)外許多學(xué)者對(duì)影響彈頭、誘餌等空間目標(biāo)的溫度與紅外輻射特性的主要因素進(jìn)行了研究[6-8],獲得了目標(biāo)在日照區(qū)及地球陰影區(qū)的整體溫度或溫度分布隨時(shí)間的變化[9-10]。傳統(tǒng)方法計(jì)算目標(biāo)內(nèi)部輻射換熱時(shí),通常將目標(biāo)內(nèi)壁視為灰體面元[11-12],利用輻射角系數(shù)建立傳熱方程,通過(guò)多次迭代求解灰體面元間的輻射換熱,整個(gè)求解過(guò)程較費(fèi)時(shí)間。張偉清等先利用輻射角系數(shù)求解灰體面元間的輻射傳遞系數(shù)[13],再直接利用輻射傳遞系數(shù)求解灰體面元間的輻射換熱,避免了直接利用角系數(shù)求解輻射換熱時(shí)的迭代計(jì)算過(guò)程,加速了計(jì)算進(jìn)程。值得注意的是,當(dāng)氣球內(nèi)壁鋁膜足夠光滑時(shí),氣球內(nèi)表面將具備良好的鏡面反射特性而非漫反射特性,此時(shí)利用基于面元灰體假設(shè)的計(jì)算方法將會(huì)給輻射換熱計(jì)算引入新的誤差,因此需要計(jì)算鏡反射面元間的輻射傳遞系數(shù)。

    由于鏡面反射的存在,面元的輻射熱量不可避免地在球腔內(nèi)多次反射,因此優(yōu)先采用蒙特卡羅方法計(jì)算面元間的輻射傳遞系數(shù)。經(jīng)典的蒙特卡羅方法利用面元的發(fā)射率隨機(jī)吸射光線(xiàn),通過(guò)統(tǒng)計(jì)面元吸收光線(xiàn)數(shù)目的方法計(jì)算輻射傳遞系數(shù)[14-15],這種隨機(jī)吸收光線(xiàn)的方法在面元發(fā)射率ε較大時(shí)會(huì)導(dǎo)致反射光線(xiàn)數(shù)量的急劇減少。例如,當(dāng)系統(tǒng)發(fā)射率ε=0.9、面元發(fā)射光線(xiàn)數(shù)量N=10 000條時(shí),經(jīng)過(guò)系統(tǒng)的1,2,3次反射后,剩余光線(xiàn)數(shù)量分別約為1 000,100,10。多次反射發(fā)生后,由于剩余光線(xiàn)的樣本數(shù)量太少,難以對(duì)某面元的反射光線(xiàn)在半空間的分布情況進(jìn)行分析。通過(guò)大量增加發(fā)射光線(xiàn)數(shù)量的方法雖然能夠緩解上述問(wèn)題,但同時(shí)也引入了大量的碰撞與求交計(jì)算,限制了蒙特卡羅方法的計(jì)算速度。

    本文通過(guò)建立鏡反射內(nèi)壁面元的蒙特卡羅漫發(fā)射模型,根據(jù)光線(xiàn)的最大反射次數(shù)計(jì)算每次反射事件發(fā)生時(shí)光線(xiàn)損失的能量份額,利用光線(xiàn)在球腔內(nèi)的傳播特性直接求解光線(xiàn)與內(nèi)壁的多次反射交點(diǎn),通過(guò)統(tǒng)計(jì)各次反射事件中面元吸收的光線(xiàn)能量份額計(jì)算面元間的輻射傳遞系數(shù),在此基礎(chǔ)上分析鏡反射內(nèi)壁的輻射傳遞系數(shù)分布特性。

    1 蒙特卡羅漫發(fā)射模型

    1.1 目標(biāo)建模與微元剖分

    氣球目標(biāo)被建模為半徑為R0的球面,并在天頂角和方位角方向進(jìn)行等間隔剖分,獲得一系列球面微元,并利用微元在天頂角、方位角方向的位序確定微元的位置。由球的性質(zhì)可知,內(nèi)壁微元表面任意位置點(diǎn)的法矢量為該點(diǎn)指向球心的單位矢量。

    1.2 光線(xiàn)的隨機(jī)發(fā)射位置、方向與能量

    在蒙特卡羅方法中,面元向半空間的輻射能量被視為從面元隨機(jī)位置、以隨機(jī)方向發(fā)射的光線(xiàn)的集合,每束光線(xiàn)都具備一定的能量[14-15]。假定面元向半空間發(fā)射的光線(xiàn)總數(shù)為N,則每束光線(xiàn)攜帶的初始能量為

    (1)

    式中:T為面元的瞬態(tài)溫度;ρ為面元的反射率;σ為Steven-Boltzmann常數(shù)。

    在球坐標(biāo)系中,光線(xiàn)在球面微元內(nèi)表面的隨機(jī)發(fā)射位置可表示為

    (2)

    式中:θa,θb分別是球面微元天頂角的上、下界;φa,φb分別是微元方位角的上、下界;Rθ,Rφ是在[0,1]區(qū)間滿(mǎn)足均勻分布的隨機(jī)變量。

    光線(xiàn)在隨機(jī)發(fā)射位置點(diǎn)上的隨機(jī)發(fā)射方向的求解過(guò)程可參考文獻(xiàn)[14]。通過(guò)分析面元向半空間范圍內(nèi)投射光束的二維概率密度,可推導(dǎo)光束的二維分布函數(shù)并得到光束的發(fā)射方向與隨機(jī)數(shù)的關(guān)系為

    (3)

    式中:θf(wàn),φf(shuō)分別為發(fā)射光束的天頂角、方位角;Rθ,Rφ為[0,1]區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)。

    2 輻射傳遞系數(shù)的計(jì)算

    2.1 光線(xiàn)在球腔中的傳播

    圖1所示為參考坐標(biāo)系OXYZ中球面微元A上光束的某次隨機(jī)發(fā)射位置p(θp,φp)。為計(jì)算方便,建立與參考坐標(biāo)系OXYZ重合的本體坐標(biāo)系Oxyz,并將本體坐標(biāo)系Oxyz先沿z軸逆時(shí)針旋轉(zhuǎn)角度φp,再沿y軸順時(shí)針旋轉(zhuǎn)角度π-θp,使該坐標(biāo)系的z軸負(fù)半軸恰好通過(guò)光線(xiàn)的發(fā)射位置p。因此,從本體坐標(biāo)系Oxyz向參考坐標(biāo)系OXYZ的坐標(biāo)轉(zhuǎn)換矩陣為

    (4)

    式中:Ly(β),Lz(γ)分別為坐標(biāo)系沿y軸、z軸逆時(shí)針旋轉(zhuǎn)β,γ角度對(duì)應(yīng)的旋轉(zhuǎn)矩陣,即

    圖1 面元隨機(jī)發(fā)射位置p在坐標(biāo)系OXYZ與Oxyz中的位置Fig.1 Random position ‘p’ on a spherical facet A in OXYZ coordinates and Oxyz coordinates

    如圖2所示,在坐標(biāo)系Oxyz中,光線(xiàn)從點(diǎn) p[0,0,-R0] 以角度 (θf(wàn),φf(shuō)) 發(fā)射。光線(xiàn)與球壁第1次相交于點(diǎn)d (1次反射點(diǎn)),并遵循鏡面反射定律繼續(xù)向前傳播。由于球內(nèi)表面上任意點(diǎn)的法矢量都指向球心,故過(guò)點(diǎn)d的反射光線(xiàn)必定在由球心、點(diǎn)p和點(diǎn)d所確定的大圓切平面內(nèi)。容易推得,發(fā)射光線(xiàn)的多次反射光線(xiàn)也全部落在該圓內(nèi)。從圖2中可以看出,發(fā)射光線(xiàn)與球壁的1次反射點(diǎn)d在大圓中的圓心角α1=π-2θf(wàn),光線(xiàn)與球壁的k次反射點(diǎn)對(duì)應(yīng)的圓心角為

    (5)

    式中:%表示求余操作。

    容易求出,k次反射點(diǎn)在坐標(biāo)系Oxyz中的天頂角θk和方位角φk為

    (6)

    將坐標(biāo)系Oxyz中的k次反射點(diǎn)坐標(biāo) (R0,θk,φk) 先轉(zhuǎn)換為直角坐標(biāo),再利用式(4)的坐標(biāo)轉(zhuǎn)換矩陣L求反射點(diǎn)在坐標(biāo)系OXYZ中的直角坐標(biāo),并將該直角坐標(biāo)轉(zhuǎn)換為球坐標(biāo)(R0,θr,φr),最后利用反射點(diǎn)在坐標(biāo)系OXYZ中的天頂角θr和方位角φr求出反射點(diǎn)所在的球面微元。

    圖2 光線(xiàn)在球腔中的多次反射 (坐標(biāo)系Oxyz)Fig.2 Multi-reflections of a light ray inside the sphere in Oxyz coordinates

    2.2 輻射傳遞系數(shù)的計(jì)算

    輻射傳遞系數(shù)是離開(kāi)一個(gè)表面到達(dá)另一個(gè)表面的能量份額,其中包含了從任意面元反射后到達(dá)目標(biāo)表面的能量。因此,球面微元的輻射傳遞系數(shù)需要考慮光線(xiàn)的多次反射貢獻(xiàn)。

    本文中,單束光線(xiàn)攜帶的初始能量e0是被面元依吸收率逐漸吸收而衰減的。當(dāng)光線(xiàn)發(fā)生第1次反射時(shí),反射光束的能量變?yōu)棣裡0,被交點(diǎn)所屬微元吸收的能量為(1-ρ)e0。當(dāng)光線(xiàn)發(fā)生第2次反射時(shí),反射光束的能量變?yōu)棣?e0,被交點(diǎn)所屬微元吸收的能量為(1-ρ)ρe0。以此類(lèi)推,第k次反射時(shí),反射光束的能量變?yōu)棣裬e0,被交點(diǎn)所屬微元吸收的能量為

    (7)

    給定閾值η,當(dāng)?shù)趉次反射發(fā)生時(shí),若反射光束的能量ρke0與初始能量e0之比低于閾值η時(shí)視為光束能量被完全吸收,即ρk≤η,則最小反射次數(shù)nr的取值可表示為

    (8)

    式中:「?表示上取整。

    綜上所述,對(duì)于球面微元的任意一束初始發(fā)射光線(xiàn),利用式(5)~(7)可一次求得其nr個(gè)反射點(diǎn)的所屬微元及微元吸收能量。對(duì)微元的全部N束發(fā)射光線(xiàn)進(jìn)行遍歷,統(tǒng)計(jì)每束光線(xiàn)的全部nr個(gè)反射點(diǎn)的所屬微元,并計(jì)算微元吸收能量,可求得該微元對(duì)其它微元的輻射傳遞系數(shù)為

    (9)

    式中:Ej為微元i的輻射能量中最終被微元j吸收的部分,包括微元i對(duì)微元j的直接入射能量和被其它面元反射后到達(dá)微元j的能量。

    3 數(shù)值計(jì)算結(jié)果與分析

    本例的參數(shù)設(shè)置如下。氣球半徑R0=1 m,內(nèi)壁鏡面鋁膜的反射率為ρ,球面沿天頂角剖分18份 (dθ=10°),沿方位角方向剖分36份 (dφ=10°)。對(duì)特定的球面微元,隨機(jī)發(fā)射光線(xiàn)的數(shù)目N=100萬(wàn)條。光束被多次吸收、反射后,當(dāng)反射光束的能量與初始能量之比低于閾值η=0.01時(shí),輻射能量被視為全部吸收,光線(xiàn)的多次反射計(jì)算過(guò)程結(jié)束。

    3.1 輻射角系數(shù)

    當(dāng)球內(nèi)壁為黑體時(shí) (ε=1,ρ=0),輻射傳遞系數(shù)等于輻射角系數(shù)。球內(nèi)壁任意兩球面微元m,n之間的角系數(shù)理論值為[15]

    (10)

    式中:Sn為微元n的面積;θa,θb分別為微元n的天頂角上、下界;φa,φb分別為微元的方位角上、下界。

    由式(10)可知,球面微元間的角系數(shù)僅與微元n的面積有關(guān),與發(fā)射微元m的位置和面積無(wú)關(guān)。

    圖3a)為球面發(fā)射微元A(天頂角、方位角中心為 (135°,115°)) 對(duì)球內(nèi)壁全部微元的角系數(shù)計(jì)算值二維展開(kāi)圖,垂直方向?yàn)榍騼?nèi)壁微元的天頂角,水平方向?yàn)榉轿唤?。圖3b)為3個(gè)特定的方位角下 (對(duì)應(yīng)于圖3a)方框) 的輻射角系數(shù)計(jì)算值與理論值對(duì)比,圖3c)為計(jì)算值與理論值的相對(duì)誤差。

    圖3 球面微元的角系數(shù)理論值與計(jì)算值Fig.3 Theoretical and calculated shape factors for sphere facets

    從圖3a)可以看出,該微元對(duì)球面各微元的輻射角系數(shù)在方位角方向保持相對(duì)穩(wěn)定,僅在天頂角方向變化。從圖3b)~c)可以看出,輻射角系數(shù)計(jì)算值與理論值很接近。在本例的參數(shù)設(shè)置下,相對(duì)誤差大部分處于±5%的范圍內(nèi),最大相對(duì)誤差不超過(guò)10%。由理論式(10)可知,當(dāng)球面沿方位角方向均勻剖分時(shí),φb,φa為常數(shù),角系數(shù)與方位角無(wú)關(guān),由此證明了本文方法計(jì)算輻射角系數(shù)的正確性。

    3.2 輻射傳遞系數(shù)

    圖4a)~c)分別為鏡面反射率ρ=0.40,0.70,0.99時(shí) (反射次數(shù)nr=6,13,459),發(fā)射微元A(135°,115°) 對(duì)全體微元的輻射傳遞系數(shù)矩陣。從圖中可以看出,與圖3a)的黑體面元不同,鏡面反射時(shí),輻射傳遞系數(shù)沿方位角方向也發(fā)生變化,并呈現(xiàn)出幾個(gè)特性。①隨著反射率的增大,輻射傳遞系數(shù)存在兩個(gè)較為明顯的峰值區(qū)域,一個(gè)恰為微元A(135°,115°)自身,另一個(gè)是其球心對(duì)稱(chēng)微元B(45°,295°) (注:球面兩點(diǎn)a,b關(guān)于球心對(duì)稱(chēng)時(shí),其天頂角滿(mǎn)足θa+θb=180°,方位角滿(mǎn)足 |φa-φb|=180°)。②當(dāng)反射率較低時(shí),微元A的輻射傳遞系數(shù)小于微元B,隨著反射率的增加,其輻射傳遞系數(shù)逐漸接近微元B。當(dāng)反射率接近于1時(shí),兩者的值近似相等。

    圖4 發(fā)射微元A對(duì)球面微元的輻射傳遞系數(shù)Fig.4 Radiative heat transfer coefficient between a manua-lly selected facet A and the whole spherical facet

    圖5a)~b)為反射率ρ=0.99時(shí) (nr=459),2個(gè)不同位置的發(fā)射微元A(5°,55°),(85°,85°) 對(duì)全體球面微元的輻射傳遞系數(shù)矩陣。從圖5中可以看出,與圖4c)類(lèi)似,輻射傳遞系數(shù)矩陣中的兩個(gè)峰值所在的微元在空間上關(guān)于球心對(duì)稱(chēng)。由此可以推知,當(dāng)反射率較大時(shí),任意位置的發(fā)射微元對(duì)全體球面微元的輻射傳遞系數(shù)均存在2個(gè)峰值,且分別位于微元自身及其球心對(duì)稱(chēng)微元。

    接下來(lái)首先分析圖5中輻射傳遞系數(shù)存在2個(gè)關(guān)于球心對(duì)稱(chēng)的峰值的原因。在圖2所示的坐標(biāo)系Oxyz中,假設(shè)微元A,B分別為發(fā)射微元及其球心對(duì)稱(chēng)微元。當(dāng)反射率ρ增大時(shí),單束光線(xiàn)的反射次數(shù)將增加;當(dāng)ρ接近于1時(shí),反射次數(shù)nr將非常大,各次反射點(diǎn)將散布于大圓圓周上。雖然單一光線(xiàn)的多次反射點(diǎn)在大圓圓周上的散布位置是確定的;但對(duì)于發(fā)射天頂角θf(wàn)∈[0,π/2]的大量隨機(jī)發(fā)射光線(xiàn)而言,其多次反射點(diǎn)落在大圓圓周各處的概率將近似相等。由于大圓切割球心對(duì)稱(chēng)微元A,B形成的弧長(zhǎng)相等,因此,在方位角φ=φf(shuō)上,反射點(diǎn)落在微元A上弧線(xiàn)的概率等于其落在微元B上弧線(xiàn)的概率。特別地,由圖2可知,由于任意方位角φ∈[0,2π)的大圓都經(jīng)過(guò)微元A,B,因此反射點(diǎn)落在微元A,B上的概率遠(yuǎn)高于其他微元。即反射率接近于1時(shí),微元A,B的輻射傳遞系數(shù)近似相等且遠(yuǎn)高于其他球面微元。

    圖5 ρ=0.99時(shí)不同位置的發(fā)射微元A對(duì)應(yīng)的輻射傳遞系數(shù)矩陣Fig.5 Radiative heat transfer coefficient matrix corresponding to two different manually selected facets A when reflectivity ρ is 0.99

    對(duì)圖4a)~c)中微元A的輻射傳遞系數(shù)隨反射率ρ的增大而逐漸逼近微元B的原因進(jìn)行分析時(shí),需要對(duì)微元A,B吸收的光線(xiàn)能量成份進(jìn)行研究。通過(guò)對(duì)微元A的每條發(fā)射光線(xiàn)進(jìn)行追蹤,記錄微元上的各個(gè)反射點(diǎn)都源于光線(xiàn)的第幾次反射行為,然后統(tǒng)計(jì)微元上各次反射行為發(fā)生的次數(shù)。圖6所示為在圖4參數(shù)設(shè)置下,反射率ρ=0.40,0.70,0.83時(shí),微元A,B上光線(xiàn)的第k次反射行為 (k=1,2,…,nr) 發(fā)生的次數(shù)統(tǒng)計(jì)。

    圖6 不同反射率時(shí)面元A,B上各次反射發(fā)生的次數(shù)Fig.6 Statistics of multi-reflection of rays on facet A and its spherically symmetric facet B

    當(dāng)η=0.01,ρ=0.4時(shí),由式(8)可知光線(xiàn)發(fā)生6次反射即被完全吸收。如圖6a)所示,當(dāng)k=1時(shí)微元A,B上光線(xiàn)的第1次反射行為發(fā)生的次數(shù)近似相等,即微元A發(fā)出的光線(xiàn)中,與微元A,B直接發(fā)生碰撞的光線(xiàn)數(shù)量近似相等。由于微元A,B是關(guān)于球心對(duì)稱(chēng)的微元,面積相等,因此光線(xiàn)與兩面元發(fā)生碰撞的概率相等,這與節(jié)3.1中輻射角系數(shù)的計(jì)算結(jié)論是一致的。

    k=2表示光線(xiàn)的第2次反射行為,即發(fā)射微元A發(fā)出的光線(xiàn)經(jīng)過(guò)球面微元的1次鏡面反射后再與微元A或B相交。由圖2可知,光線(xiàn)由微元A發(fā)出后,只有那些能夠與微元B(或其附近微元) 發(fā)生直接碰撞的光線(xiàn)才有可能經(jīng)1次反射后與微元A相交,此部分光線(xiàn)數(shù)量較少。但是,光線(xiàn)由微元A發(fā)出后,所有能夠與“赤道”微元 (與Oxy平面相交的球面微元) 發(fā)生直接碰撞的光線(xiàn)經(jīng)過(guò)1次鏡面反射后,都有可能與微元B相交。因此,2次反射點(diǎn)落在微元B上的光線(xiàn)數(shù)量要遠(yuǎn)大于落在微元A上的光線(xiàn)數(shù)量,即微元B吸收的光線(xiàn)的2次碰撞能量遠(yuǎn)大于微元A,直接導(dǎo)致此時(shí)微元B的輻射傳遞系數(shù)較微元A大。

    由圖6a)~c)可以看出,當(dāng)反射率ρ增大時(shí),反射次數(shù)增多,且光線(xiàn)的第k次反射行為 (k≥3) 發(fā)生在微元A,B上的數(shù)目趨于相等。此時(shí),雖然微元B吸收的光線(xiàn)的2次碰撞能量仍較微元A高,但相對(duì)于總能量而言,其所占的能量份額隨反射率的增大而逐漸減小。因此,如圖4所示,隨著反射率的增大,微元A的輻射傳遞系數(shù)逐漸接近于微元B。

    在本例中,光線(xiàn)被多次反射后,其攜帶的能量會(huì)逐漸降低,若反射光射的能量與初始能量之比低于閾值η=0.01,則光線(xiàn)能量被視為全部吸收,即光束被吸收的能量份額為0.99。當(dāng)該閾值減小時(shí),計(jì)算精度將會(huì)提高。圖7所示為η=0.000 1、其他設(shè)置與圖5相同時(shí),得到的輻射傳遞系數(shù)矩陣。

    圖7 η=0.000 1,ρ=0.99時(shí)不同位置的微元A對(duì)應(yīng)的輻射傳遞系數(shù)矩陣 Fig.7 Radiative heat transfer coefficient matrix corresponding to two different manually selected facets when η=0.000 1,ρ is 0.99

    從圖7與圖5可以看出計(jì)算結(jié)果非常相近。這是因?yàn)楫?dāng)閾值η=0.000 1時(shí),光束被吸收的能量份額為0.999 9;而當(dāng)η=0.01時(shí),光束被吸收的能量份額為0.99。由此可見(jiàn),兩者的能量變化不足1%,因此圖7得到的結(jié)果與圖5類(lèi)似,但數(shù)據(jù)精度更高。

    4 結(jié)束語(yǔ)

    本文利用蒙特卡羅方法計(jì)算了球內(nèi)壁的輻射角系數(shù)及不同鏡面反射率情況下的輻射傳遞系數(shù)。輻射角系數(shù)的計(jì)算值與理論值一致。輻射傳遞系數(shù)的結(jié)果表明,首先,微元對(duì)全部球面微元的輻射傳遞系數(shù)矩陣中存在2個(gè)峰值,即微元對(duì)其自身、以及微元對(duì)其球心對(duì)稱(chēng)微元的輻射傳遞系數(shù)。其次,隨著反射率的增大,微元對(duì)其自身的輻射傳遞系數(shù)逐漸接近于微元對(duì)其球心對(duì)稱(chēng)微元的輻射傳遞系數(shù)。當(dāng)反射率接近于1時(shí),兩者的輻射傳遞系數(shù)值近似相等。由此可知,對(duì)于誘餌球鏡反射內(nèi)壁微元而言,當(dāng)某微元的輻射出射能量一定時(shí),反射率越大,輻射能量越傾向于向該微元及其球?qū)ΨQ(chēng)微元集中。本文的結(jié)論有助于分析空間誘餌球向陽(yáng)面高溫微元與背陰面低溫微元之間的輻射換熱與溫度變化過(guò)程。

    [1] 酈蘇丹,任萱.大氣層外誘餌釋放研究[J].宇航學(xué)報(bào),2001,22(2):100-104. LI Su-dan,REN Xuan.Research of Releasing Decoy Outside Atmosphere[J].Journal of Astronautics,2001,22(2):100-104.

    [2] 林兩魁,謝愷,徐暉,等.中段彈道目標(biāo)群的紅外成像仿真研究[J].紅外與毫米波學(xué)報(bào),2009,28(3):218-223. LIN Liang-kui,XIE Kai,XU Hui,et al.Research on Infrared Imaging Simulation of Midcourse Ballistic Object Target Complex[J].Journal of Infrared and Millimeter Waves,2009,28(3):218-223.

    [3] 王碩,張奕群,孫冰巖.紅外點(diǎn)目標(biāo)跟蹤方法綜述[J].現(xiàn)代防御技術(shù),2016,44(2):124-134. WANG Shuo,ZHANG Yi-qun,SUN Bing-yan.Review of Point Target Tracking Methods Based on Infrared Sensor[J].Modern Defence Technology,2016,44(2):124-134.

    [4] 浦甲倫,崔乃剛,郭繼峰.天基紅外預(yù)警衛(wèi)星系統(tǒng)及其探測(cè)能力分析[J].現(xiàn)代防御技術(shù),2008,36(4):68-72. PU Jia-lun,CUI Nai-gang,GUO Ji-feng.Space-Based Infrared System and the Analysis of Its Detecting Capability[J].Modern Defence Technology,2008,36(4):68-72.

    [5] SESSLER A M,CORNWALL J M,DIETZ B,et al.Countermeasures:a Technical Evaluation of the Operational Effectiveness of the Planned US National Missile Defense System[R].Cambridge,MA:Union of Concerned Scientists and MIT Security Studies Program,April 2000.

    [6] 姚連興,侯秋萍,羅繼強(qiáng).彈道導(dǎo)彈中段目標(biāo)表面溫度與紅外突防研究[J].航天電子對(duì)抗,2005,21(2):5-16. YAO Lian-xing,HOU Qiu-ping,LUO Ji-qiang.Research on Surface Temperature and Infrared Signature of the Ballistic Warhead in Midcourse[J].Aerospace Electronic Warefare,2005,21(2):5-16.

    [7] 汪洪源,陳赟.天基空間目標(biāo)紅外動(dòng)態(tài)輻射特性建模與仿真[J].紅外與激光工程,2016,45(5):0504002. WANG Hong-yuan,CHEN Yun.Modeling and Simulation of Infrared Dynamic Characteristics of Space-Based Targets[J].Infrared and Laser Engineering,2016,45(5):0504002.

    [8] 敬韓博,王英瑞.臨近空間高超聲速目標(biāo)天基紅外探測(cè)技術(shù)研究[J].現(xiàn)代防御技術(shù),2016,44(6):80-84. JING Han-bo,WANG Ying-rui.Space-Based Infrared Detection for Near Space Hyptrsonic Targets[J].Modern Defence Technology,2016,44(6):80-84.

    [9] ZHU D Q,SHEN W T,CAI G B,et al.Numerical Simulation and Experimental Study of Factors Influencing the Optical Characteristics of a Spatial Target[J].Applied Thermal Engineering,2013(50):749-762.

    [10] 張駿,楊華,凌永順,等.彈道導(dǎo)彈中段彈頭表面溫度場(chǎng)分布理論分析[J].紅外與激光工程,2005,34(5):583-586. ZHANG Jun,YANG Hua,LING Yong-shun,et al.Theoretical Analysis of Temperature Field on the Surface of Ballistic Missile Warhead in Midcourse[J].Infrared and Laser Engineering,2005,34(5):583-586.

    [11] LI F B,XU X J.Modeling Time-Evolving Infrared Characteristics for Space Object with Mircromotions[J].IEEE Transactions on Aerospace and Electronic Systems,2012,48(4):3567-3577.

    [12] 厲夫兵,許小劍.計(jì)算超薄多層介質(zhì)誘餌溫度的快速算法[J].紅外與激光工程,2011,40(6):986-991.

    LI Fu-bing,XU Xiao-jian.Fast Algorithm for Temperature Calculation of Multi-Layered Ultra-Thin Decoy[J].Infrared and Laser Engineering,2011,40(6):986-991.

    [13] 張偉清,宣益民,韓玉閣.單元表面間輻射傳遞系數(shù)的新型計(jì)算方法[J].宇航學(xué)報(bào),2005,26(1):77-81. ZHANG Wei-qing,XUAN Yi-min,HAN Yu-ge.A New Method of Radiative Transfer Coefficient between Unit Surfaces[J].Journal of Astronautics,2005,26(1):77-81.

    [14] HOWELL J R,PERLMUTTER M.Monte Carlo Solution of Thermal Transfer Through Radiant Media between Gray Walls[J].Journal of Heat Transfer,1964,86(1):116-122.

    [15] 卞伯繪.輻射換熱的分析與計(jì)算[M].北京:清華大學(xué)出版社,1988:81-83,144-145. BIAN Be-hui.Analysis and Calculation of Radiative Heat Transfer[M].Beijing:Tsinghua Universtiy Press,1988:81-83,144-145.

    Radiative Heat Transfer Coefficient Calculation of Inflatable Sphere with Specular Inner Surface

    LI Fu-bing,LENG Jun-min,LI Hong-lian

    (Beijing Information Science & Technology University,School of Information & Communication Engineering,Beijing 100101,China)

    Monte Carlo technique is used to calculate the radiative heat transfer coefficient for an inflatable sphere with specular reflection inner surface. Light rays are emitted diffusely from a spherical facet and being traced. The maximum number of times a ray can be reflected is derived, and the multi-intersections of a ray with the inner side of spherical surface are obtained. Radiative heat transfer coefficients between a manually selected facet and the whole spherical facet are calculated at different reflectivities. The results indicate that there are two peaks in the calculated coefficients, the smaller one is the coefficient between the selected facet and itself (i.e., the self-to-self coefficient), and the larger one is between the selected facet and its spherically symmetric facet (i.e., the self-to-spherical-symmetry coefficient). The self-to-self radiative heat transfer coefficient gradually approaches and finally equals the self-to-spherical-symmetry coefficient as the reflectivity increases to 1.

    radiation transfer coefficient;radiation heat transfer;specular reflection;Monte Carlo;shape factor;sphere

    2016-12-01;

    2017-02-11

    北京市自然科學(xué)基金資助項(xiàng)目(3164043)

    厲夫兵(1982-),男,山東日照人。講師,博士,主要從事目標(biāo)紅外輻射特性建模研究。

    10.3969/j.issn.1009-086x.2017.03.030

    TK124;TP391.9

    A

    1009-086X(2017)-03-0193-07

    通信地址:100101 北京市朝陽(yáng)區(qū)北四環(huán)中路35號(hào)北京信息科技大學(xué)信息與通信工程學(xué)院

    E-mail:lifubing@bistu.edu.cn

    猜你喜歡
    鏡面反射面元球心
    隨機(jī)粗糙面散射中遮蔽效應(yīng)算法的改進(jìn)
    直擊多面體的外接球的球心及半徑
    光滑物體表面反射光偏振特征分析及反射光分離技術(shù)*
    基于最短路徑的GNSS-R鏡面反射點(diǎn)算法
    ?如何我解決幾何體的外接球問(wèn)題
    例析確定球心位置的策略
    基于改進(jìn)Gordon方程的RCS快速算法
    畫(huà)好草圖,尋找球心
    面元細(xì)分觀測(cè)系統(tǒng)應(yīng)用分析
    化工管理(2014年14期)2014-08-15 00:51:32
    一種基于Kd-tree 射線(xiàn)追蹤法的衛(wèi)星RCS 預(yù)估方法
    一个人观看的视频www高清免费观看| 欧美性猛交╳xxx乱大交人| 日韩有码中文字幕| 久久中文看片网| 动漫黄色视频在线观看| 在线播放无遮挡| 最近最新中文字幕大全电影3| 97热精品久久久久久| 一个人免费在线观看电影| 亚州av有码| 99久久九九国产精品国产免费| 麻豆一二三区av精品| 非洲黑人性xxxx精品又粗又长| 制服丝袜大香蕉在线| 国产真实伦视频高清在线观看 | 好男人在线观看高清免费视频| 日韩欧美精品免费久久 | 99热精品在线国产| 中出人妻视频一区二区| 有码 亚洲区| 精品人妻视频免费看| 久久久久久久久大av| 亚洲真实伦在线观看| 99视频精品全部免费 在线| 可以在线观看毛片的网站| 精品一区二区三区av网在线观看| 亚洲成av人片在线播放无| 精品不卡国产一区二区三区| 日本与韩国留学比较| 每晚都被弄得嗷嗷叫到高潮| 波多野结衣高清无吗| 中文字幕av成人在线电影| 国产一区二区三区视频了| 亚洲va日本ⅴa欧美va伊人久久| 亚洲欧美日韩东京热| 成人特级黄色片久久久久久久| 亚洲成人久久性| 此物有八面人人有两片| 人妻制服诱惑在线中文字幕| 亚洲午夜理论影院| 日日夜夜操网爽| 日本在线视频免费播放| 国产精品电影一区二区三区| 性欧美人与动物交配| 免费观看的影片在线观看| 毛片女人毛片| 亚洲美女黄片视频| 欧美精品啪啪一区二区三区| 亚洲专区国产一区二区| 国产高潮美女av| 日韩免费av在线播放| 精品人妻视频免费看| 国产一级毛片七仙女欲春2| 成年版毛片免费区| 三级男女做爰猛烈吃奶摸视频| 小蜜桃在线观看免费完整版高清| 色av中文字幕| 亚洲成a人片在线一区二区| 在现免费观看毛片| 免费搜索国产男女视频| 99国产综合亚洲精品| 亚洲av.av天堂| 欧美bdsm另类| 亚洲成人精品中文字幕电影| 婷婷精品国产亚洲av| 一个人看视频在线观看www免费| 人人妻人人看人人澡| 中文字幕人妻熟人妻熟丝袜美| 尤物成人国产欧美一区二区三区| 香蕉av资源在线| 国产一区二区激情短视频| 舔av片在线| 最近中文字幕高清免费大全6 | 亚洲一区二区三区不卡视频| 久久久久免费精品人妻一区二区| 精品日产1卡2卡| 国产单亲对白刺激| 人妻久久中文字幕网| 欧美成人免费av一区二区三区| 欧美乱妇无乱码| 久久久久性生活片| 成年女人毛片免费观看观看9| 精品99又大又爽又粗少妇毛片 | 亚洲avbb在线观看| 日韩欧美在线二视频| 日本一本二区三区精品| 亚洲精华国产精华精| 久久久久久久精品吃奶| 欧美激情在线99| 日本 欧美在线| 可以在线观看的亚洲视频| av女优亚洲男人天堂| 好男人在线观看高清免费视频| 亚洲五月婷婷丁香| 一级作爱视频免费观看| a在线观看视频网站| 一个人观看的视频www高清免费观看| 久久性视频一级片| 日本撒尿小便嘘嘘汇集6| 黄色丝袜av网址大全| 免费观看精品视频网站| 久久国产精品影院| 亚洲av二区三区四区| 亚洲av二区三区四区| 欧美最新免费一区二区三区 | 国产一区二区三区视频了| 搞女人的毛片| 免费无遮挡裸体视频| 日本熟妇午夜| 国产白丝娇喘喷水9色精品| 国产午夜福利久久久久久| 色哟哟哟哟哟哟| 国产伦一二天堂av在线观看| 午夜精品一区二区三区免费看| 免费看日本二区| 午夜精品在线福利| 岛国在线免费视频观看| 欧美一区二区亚洲| 国产精品久久久久久人妻精品电影| 亚洲av一区综合| 真人一进一出gif抽搐免费| 午夜福利在线观看吧| 亚洲精品影视一区二区三区av| 欧美zozozo另类| 亚洲av第一区精品v没综合| 在线观看午夜福利视频| 国产亚洲精品久久久com| 高清毛片免费观看视频网站| 欧美色欧美亚洲另类二区| 欧美中文日本在线观看视频| 精品久久久久久久久久久久久| 亚洲久久久久久中文字幕| 午夜久久久久精精品| 亚洲久久久久久中文字幕| 亚洲精品456在线播放app | 国产亚洲av嫩草精品影院| 国产精品免费一区二区三区在线| 日韩欧美一区二区三区在线观看| 国产成人福利小说| 免费大片18禁| 最后的刺客免费高清国语| bbb黄色大片| 中文字幕精品亚洲无线码一区| 在现免费观看毛片| 1024手机看黄色片| 欧美成人免费av一区二区三区| 男人狂女人下面高潮的视频| 国产69精品久久久久777片| 欧美丝袜亚洲另类 | 婷婷色综合大香蕉| 国产真实伦视频高清在线观看 | 夜夜看夜夜爽夜夜摸| 成人毛片a级毛片在线播放| 亚洲精品一区av在线观看| 日本五十路高清| 午夜a级毛片| 国产精品,欧美在线| 18禁黄网站禁片免费观看直播| 一区福利在线观看| 免费看a级黄色片| 成人永久免费在线观看视频| 久久久久国产精品人妻aⅴ院| 亚洲精品一卡2卡三卡4卡5卡| 亚洲电影在线观看av| 日本黄大片高清| 国产一区二区三区视频了| 好男人电影高清在线观看| 白带黄色成豆腐渣| 欧美最新免费一区二区三区 | 午夜精品一区二区三区免费看| 国产精华一区二区三区| av视频在线观看入口| 欧美bdsm另类| 男人舔女人下体高潮全视频| 五月伊人婷婷丁香| 九九久久精品国产亚洲av麻豆| 久久久久久久久久黄片| 亚洲精品色激情综合| 婷婷六月久久综合丁香| 三级国产精品欧美在线观看| 亚洲第一区二区三区不卡| 成人毛片a级毛片在线播放| 亚洲av成人不卡在线观看播放网| 国产精品久久久久久精品电影| 婷婷六月久久综合丁香| 国产精品三级大全| 51午夜福利影视在线观看| 亚洲熟妇中文字幕五十中出| 精品久久国产蜜桃| 亚洲性夜色夜夜综合| 啪啪无遮挡十八禁网站| 男女之事视频高清在线观看| 国产精品精品国产色婷婷| 国产精品一区二区三区四区久久| 久久精品国产亚洲av香蕉五月| 午夜福利欧美成人| 男女下面进入的视频免费午夜| 男女做爰动态图高潮gif福利片| 亚洲av二区三区四区| 好看av亚洲va欧美ⅴa在| 老司机深夜福利视频在线观看| 在线a可以看的网站| 中文字幕av成人在线电影| 精品不卡国产一区二区三区| 亚洲中文字幕一区二区三区有码在线看| 大型黄色视频在线免费观看| av在线观看视频网站免费| 久久亚洲精品不卡| 色综合站精品国产| 成年女人永久免费观看视频| 长腿黑丝高跟| 观看美女的网站| www.www免费av| 美女 人体艺术 gogo| 又爽又黄a免费视频| 在线播放无遮挡| 午夜两性在线视频| 日本黄大片高清| 亚洲专区中文字幕在线| 自拍偷自拍亚洲精品老妇| 午夜福利欧美成人| 精品午夜福利在线看| 国产欧美日韩精品亚洲av| 嫩草影院精品99| 哪里可以看免费的av片| 精品乱码久久久久久99久播| 亚洲国产欧美人成| 女人十人毛片免费观看3o分钟| 国产真实乱freesex| 男插女下体视频免费在线播放| 男人的好看免费观看在线视频| 午夜老司机福利剧场| 国产人妻一区二区三区在| 日韩欧美精品免费久久 | 天堂网av新在线| 亚洲无线观看免费| av福利片在线观看| 中文字幕av成人在线电影| 男女下面进入的视频免费午夜| 99国产综合亚洲精品| 毛片女人毛片| 极品教师在线视频| av在线老鸭窝| 日日摸夜夜添夜夜添av毛片 | netflix在线观看网站| 国产伦一二天堂av在线观看| 伦理电影大哥的女人| 久久国产精品影院| h日本视频在线播放| av在线蜜桃| 他把我摸到了高潮在线观看| 亚洲av.av天堂| 亚洲黑人精品在线| 嫩草影院精品99| 男女视频在线观看网站免费| 欧美乱妇无乱码| 国产69精品久久久久777片| 国产精华一区二区三区| 国产麻豆成人av免费视频| 精品日产1卡2卡| 国产探花在线观看一区二区| 欧美激情国产日韩精品一区| 一个人看的www免费观看视频| 国产午夜福利久久久久久| 麻豆成人av在线观看| 久久欧美精品欧美久久欧美| 国产伦精品一区二区三区视频9| 国产精品一区二区三区四区久久| 观看美女的网站| 国产视频内射| 国产av一区在线观看免费| 特级一级黄色大片| 国产一区二区三区视频了| 国产蜜桃级精品一区二区三区| 极品教师在线免费播放| 免费搜索国产男女视频| 国产精品久久视频播放| 国内精品久久久久久久电影| 九色国产91popny在线| 日本一本二区三区精品| 在线国产一区二区在线| 日韩欧美在线乱码| 又紧又爽又黄一区二区| 91狼人影院| 网址你懂的国产日韩在线| av中文乱码字幕在线| 国产精品免费一区二区三区在线| 欧美性猛交╳xxx乱大交人| 99riav亚洲国产免费| x7x7x7水蜜桃| 中文在线观看免费www的网站| 赤兔流量卡办理| 又紧又爽又黄一区二区| 人人妻人人澡欧美一区二区| 九九热线精品视视频播放| 国产日本99.免费观看| 国产精品人妻久久久久久| 国产免费一级a男人的天堂| 在线免费观看的www视频| 欧美成人性av电影在线观看| 麻豆成人午夜福利视频| 日本三级黄在线观看| 亚洲av五月六月丁香网| 国产男靠女视频免费网站| 国产成年人精品一区二区| 最后的刺客免费高清国语| 国产久久久一区二区三区| h日本视频在线播放| 97人妻精品一区二区三区麻豆| 日本免费a在线| 狂野欧美白嫩少妇大欣赏| 成人欧美大片| 国产爱豆传媒在线观看| 久久人人爽人人爽人人片va | 亚洲欧美激情综合另类| 精品久久久久久久人妻蜜臀av| 日本黄大片高清| 国产精品三级大全| 午夜福利成人在线免费观看| 午夜日韩欧美国产| 免费在线观看影片大全网站| 他把我摸到了高潮在线观看| а√天堂www在线а√下载| 在线看三级毛片| 可以在线观看毛片的网站| 亚洲av成人不卡在线观看播放网| 日本免费a在线| 两人在一起打扑克的视频| 成人精品一区二区免费| 我的老师免费观看完整版| 国产一区二区在线av高清观看| 成人午夜高清在线视频| 欧美丝袜亚洲另类 | 国产在线精品亚洲第一网站| 国产精品国产高清国产av| 午夜福利18| 男女做爰动态图高潮gif福利片| 久99久视频精品免费| 亚洲欧美日韩东京热| 自拍偷自拍亚洲精品老妇| 国产精品,欧美在线| 国产伦人伦偷精品视频| 高清毛片免费观看视频网站| 极品教师在线视频| 中文字幕高清在线视频| 中文亚洲av片在线观看爽| 老熟妇仑乱视频hdxx| 午夜两性在线视频| 久久亚洲精品不卡| 搡老妇女老女人老熟妇| 精品久久久久久久末码| 中亚洲国语对白在线视频| 最后的刺客免费高清国语| 久9热在线精品视频| 精品熟女少妇八av免费久了| 亚洲成a人片在线一区二区| 午夜福利成人在线免费观看| 亚洲一区高清亚洲精品| 国产久久久一区二区三区| 一本精品99久久精品77| 国产又黄又爽又无遮挡在线| 国产激情偷乱视频一区二区| av专区在线播放| 韩国av一区二区三区四区| 亚洲av五月六月丁香网| 国产一区二区三区视频了| 三级男女做爰猛烈吃奶摸视频| 国产高清视频在线观看网站| 男人狂女人下面高潮的视频| 一区二区三区激情视频| 久久久精品欧美日韩精品| 性欧美人与动物交配| 日韩欧美在线乱码| 3wmmmm亚洲av在线观看| 中国美女看黄片| netflix在线观看网站| 久久精品国产亚洲av香蕉五月| 日韩有码中文字幕| 国产精品国产高清国产av| 日韩欧美三级三区| 国产真实乱freesex| 久久久久久久久久成人| 精华霜和精华液先用哪个| 深爱激情五月婷婷| 国产成人啪精品午夜网站| 老熟妇乱子伦视频在线观看| 五月伊人婷婷丁香| 男女视频在线观看网站免费| 色尼玛亚洲综合影院| 亚洲成人久久性| 成人特级av手机在线观看| 夜夜夜夜夜久久久久| 成人欧美大片| 午夜福利在线在线| 性色av乱码一区二区三区2| 午夜精品久久久久久毛片777| 国产欧美日韩精品一区二区| 亚洲乱码一区二区免费版| 免费人成视频x8x8入口观看| 免费观看的影片在线观看| 丁香六月欧美| 少妇丰满av| 午夜精品在线福利| 欧美高清成人免费视频www| 热99在线观看视频| 18禁黄网站禁片午夜丰满| 99热6这里只有精品| 天堂网av新在线| 免费在线观看影片大全网站| 国产黄色小视频在线观看| 久久久精品大字幕| 此物有八面人人有两片| 国产高清有码在线观看视频| 精品久久久久久久久久久久久| 一区二区三区激情视频| 国产精品av视频在线免费观看| 1024手机看黄色片| 国语自产精品视频在线第100页| 欧美性猛交黑人性爽| 成人一区二区视频在线观看| 男人舔女人下体高潮全视频| 亚洲一区高清亚洲精品| 午夜两性在线视频| 欧美激情在线99| 午夜视频国产福利| 99热只有精品国产| 男人狂女人下面高潮的视频| 成人特级黄色片久久久久久久| av在线天堂中文字幕| 欧美高清性xxxxhd video| 91久久精品电影网| 俄罗斯特黄特色一大片| 国产在视频线在精品| 两个人的视频大全免费| 亚洲最大成人av| 一卡2卡三卡四卡精品乱码亚洲| 麻豆成人av在线观看| 国产一区二区激情短视频| 久久精品国产99精品国产亚洲性色| 欧美日本视频| av在线蜜桃| 亚洲激情在线av| 51午夜福利影视在线观看| 亚洲人成伊人成综合网2020| 亚洲中文日韩欧美视频| 三级国产精品欧美在线观看| 国产精品精品国产色婷婷| 欧美日韩亚洲国产一区二区在线观看| 国产午夜精品论理片| 亚洲美女视频黄频| 麻豆av噜噜一区二区三区| 亚洲专区国产一区二区| 日本免费a在线| 国产av麻豆久久久久久久| 国产精品久久久久久精品电影| 在线免费观看不下载黄p国产 | 51国产日韩欧美| 成年女人看的毛片在线观看| 色播亚洲综合网| 国产伦一二天堂av在线观看| 欧美一区二区国产精品久久精品| 老司机福利观看| 国产视频一区二区在线看| 久久久久久久精品吃奶| 真实男女啪啪啪动态图| 午夜影院日韩av| 性色avwww在线观看| 搞女人的毛片| 深爱激情五月婷婷| 成人精品一区二区免费| 国产私拍福利视频在线观看| 亚洲精品乱码久久久v下载方式| 久久国产乱子伦精品免费另类| 两人在一起打扑克的视频| 成年女人看的毛片在线观看| 97热精品久久久久久| 亚洲经典国产精华液单 | 国产精品乱码一区二三区的特点| a级毛片a级免费在线| 观看免费一级毛片| 黄片小视频在线播放| 一本精品99久久精品77| 色噜噜av男人的天堂激情| 老司机午夜十八禁免费视频| 欧美日韩黄片免| 天天躁日日操中文字幕| 亚洲国产日韩欧美精品在线观看| 精品乱码久久久久久99久播| 国产在视频线在精品| 久久99热这里只有精品18| 国产精品永久免费网站| 一区福利在线观看| 狠狠狠狠99中文字幕| 窝窝影院91人妻| 九九热线精品视视频播放| 真人做人爱边吃奶动态| 久久人人爽人人爽人人片va | 久久精品91蜜桃| 中国美女看黄片| 小说图片视频综合网站| 午夜免费激情av| 亚洲精品在线美女| 国产亚洲精品久久久com| 国产极品精品免费视频能看的| 性色avwww在线观看| 夜夜看夜夜爽夜夜摸| 男插女下体视频免费在线播放| av福利片在线观看| 久久天躁狠狠躁夜夜2o2o| 变态另类成人亚洲欧美熟女| 午夜福利在线观看吧| 久久精品91蜜桃| 啪啪无遮挡十八禁网站| 亚洲专区国产一区二区| 69av精品久久久久久| 国产精品1区2区在线观看.| 一本久久中文字幕| 不卡一级毛片| 亚洲人与动物交配视频| 人人妻人人澡欧美一区二区| 91麻豆av在线| 国产免费av片在线观看野外av| av黄色大香蕉| 国产精品一及| 国产探花在线观看一区二区| 亚洲国产欧美人成| 色5月婷婷丁香| 欧美成人性av电影在线观看| 国产精品1区2区在线观看.| 给我免费播放毛片高清在线观看| 午夜福利高清视频| 欧美色欧美亚洲另类二区| 久久国产乱子免费精品| 91麻豆av在线| 我的女老师完整版在线观看| 真实男女啪啪啪动态图| 亚洲国产欧美人成| 欧美激情国产日韩精品一区| 美女xxoo啪啪120秒动态图 | 日日夜夜操网爽| 久久久久久国产a免费观看| 在线国产一区二区在线| 欧美一区二区亚洲| 亚洲,欧美,日韩| 欧美成人a在线观看| 久久香蕉精品热| 嫁个100分男人电影在线观看| 国产免费一级a男人的天堂| 热99re8久久精品国产| 最好的美女福利视频网| 校园春色视频在线观看| 久久99热这里只有精品18| 国产又黄又爽又无遮挡在线| 波多野结衣高清作品| 婷婷丁香在线五月| 婷婷六月久久综合丁香| 国产又黄又爽又无遮挡在线| 97超视频在线观看视频| 99精品在免费线老司机午夜| 一个人免费在线观看的高清视频| 欧美高清成人免费视频www| 亚洲内射少妇av| 精品久久久久久久久久免费视频| 一二三四社区在线视频社区8| 亚洲在线自拍视频| 搡女人真爽免费视频火全软件 | a级一级毛片免费在线观看| 18禁黄网站禁片免费观看直播| 久久精品国产99精品国产亚洲性色| 亚洲自拍偷在线| 国产单亲对白刺激| 久久精品人妻少妇| aaaaa片日本免费| 在线观看一区二区三区| 亚洲精品456在线播放app | 人人妻,人人澡人人爽秒播| 日韩有码中文字幕| 成人av一区二区三区在线看| 一边摸一边抽搐一进一小说| 伦理电影大哥的女人| 一区二区三区四区激情视频 | 色av中文字幕| 成人三级黄色视频| 成年女人毛片免费观看观看9| www日本黄色视频网| 欧美3d第一页| 午夜免费成人在线视频| 日本撒尿小便嘘嘘汇集6| 波多野结衣高清无吗| 欧美日韩福利视频一区二区| 两个人的视频大全免费| 欧美色视频一区免费| 亚洲成人免费电影在线观看| 一区二区三区免费毛片| 国产亚洲精品久久久久久毛片| 给我免费播放毛片高清在线观看| 日本免费一区二区三区高清不卡| 热99re8久久精品国产| 免费看a级黄色片| 国产亚洲欧美98| av国产免费在线观看| 欧美另类亚洲清纯唯美| 久久精品夜夜夜夜夜久久蜜豆| 国产综合懂色| 精品乱码久久久久久99久播| 亚洲在线自拍视频| 丰满乱子伦码专区| 免费看光身美女| 免费人成视频x8x8入口观看| .国产精品久久| 99riav亚洲国产免费| 亚洲国产高清在线一区二区三| 久久久久九九精品影院| 日本黄色视频三级网站网址| 99在线视频只有这里精品首页|