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

浮力修正湍流模型在航空發動機火災模擬中的應用

2017-11-10 09:23:50李松陽邵興晨
航空發動機 2017年1期
關鍵詞:模型

李松陽,邵興晨

(中國航發商用航空發動機有限責任公司,上海200241)

浮力修正湍流模型在航空發動機火災模擬中的應用

李松陽,邵興晨

(中國航發商用航空發動機有限責任公司,上海200241)

由于發動機艙的火災是典型的熱驅動的浮力羽流,從探索浮力羽流的模擬方法出發,針對熱羽流的基準試驗,比較驗證了3種基于浮力修正的2個方程湍流模型;利用prePD F燃燒模型,模擬驗證了Purdue甲烷火燃燒試驗;以RR公司Trent 800發動機的1/2縮比短艙著火試驗器為原型,采用RA N S方法對由燃油泄漏引起的油池火進行了模擬計算,重現了短艙火災的主要物理過程,并與試驗測量的速度及溫度結果進行了對比,驗證了計算方法的準確性,并進一步分析了影響模擬結果的主要原因:湍流模型與燃燒模型能否準確計算近火源區域的火焰鋒面狀態,直接影響空氣卷吸及下游火羽流的溫度與速度。應用CFD技術,可以從防火設計的角度優化通風系統及短艙附件的布局。

防火設計;火災數值模擬;浮力羽流;湍流燃燒;航空發動機

0 引言

商用航空發動機的核心機艙位于核心機機匣與外涵道內壁之間的環腔區域,而風扇艙位于外涵道外壁與短艙外罩之間的環腔區域,這些艙室內布置了大量的控制設備與燃、滑油管路。為保持艙室內的溫度以及防止可燃氣體的聚集,艙內設計了通風冷卻系統。充足的氧氣、局部的高溫、潛在的燃滑油泄漏風險,使得這2個艙室成為發動機的主要火區[1-2]。這是發動機防火設計重點關注的區域。

發動機艙室內傳統的防火設計多依賴于設計經驗與部件的防火試驗,但隨著現代發動機結構日益復雜,采用了大量新材料和新技術,以及對構型的修改,用傳統的經驗和標準試驗很難預測具有復雜結構和氣動環境的短艙內的火災特性。為此,國外各大發動機公司近年大力發展基于CFD技術的短艙防火設計方法[3]。與一般的湍流燃燒模擬不同,短艙中的火災在近場是浮力驅動的熱羽流,在遠場又受到復雜通風冷卻氣流的影響,且由于空間受限,燃燒不充分,更增加了數值模擬的難度。因此,需要結合發動機艙的火災試驗,發展適用于其場景的火災數值模擬方法。

短艙火災屬于熱驅動的浮力羽流,具有很強的流動不穩定性,并且浮力在很大程度上影響著湍動能產生和耗散[4]。因此,模擬短艙火災需要準確預測浮力對湍流的影響。直接數值模擬(DNS)及大渦模擬(LES)雖然能較好地模擬湍流的內部結構,但需要大量的計算資源,目前工程應用還不常見。因此,雷諾平均(RANS)湍流模擬依然是目前常用的方法。其中,k-ε模型是RANS方法中應用最廣泛的湍流模型。不過,標準的k-ε模型需要進行修正來模擬浮力對湍動能的影響。Standard Gradient Diffusion Hypothesis(SGDH)浮力修正模型被耦合到了多種k-ε的輸運方程中[5-6],如 Standard k-ε、RNG k-ε及 Realizable k-ε模型。研究表明,經過浮力修正的k-ε模型對浮力羽流的預測明顯強于未修正的模型[7]。然而,對應這幾種浮力修正的k-ε模型之間的預測能力,還未有研究。

短艙火災模擬需要解決的另一個問題是受限空間內的湍流燃燒問題。湍流、燃燒、輻射等模型的耦合、非充分燃燒過程的描述等問題一直是國際上研究的前沿[8]。在工程計算中往往需要對燃燒反應進行一定的簡化。目前,火災模擬常用的燃燒模型包括體積熱源、渦破碎及基于混合分數的概率密度(PrePDF)模型。H.Xue研究比較了這幾個模型在受限空間火災中的預測能力,PDF模型比前二者更好[9]。然而,在發動機艙火災中的預測能力還有待進一步驗證。

本文比較了浮力修正的Standard k-ε、RNG k-ε及Realizable k-ε3個模型在浮力羽流模擬中的適用性,并結合prePDF燃燒模型和DO輻射模型,模擬了甲烷火焰和發動機機艙火災,重現了火災中的主要物理過程,驗證了計算方法的準確性,并進一步分析了影響模擬結果的主要原因。

1 物理模型

1.1 流動方程

對于熱驅動的浮力羽流,其流動守恒方程包括連續、動量和能量方程,考慮燃燒反應的情況下還包括組分濃度方程。針對熱羽流的特殊性,該守恒方程將引入低馬赫數假設,即不考慮壓縮性,密度變化源于溫度而不是壓力,平均密度通過理想氣體的狀態方程來計算。此外,本文的模擬僅考慮穩態過程,因此守恒方程的笛卡爾張量表示形式為[10]

連續方程

動量方程

式中:ρ為流體密度;P為靜壓減去流體靜壓,流體靜壓通過參考密度來體現;ui為速度矢量;xi為空間矢量;τij為粘性應力張量;gi為重力加速度矢量;φ為焓與組分濃度;Sφ為熱釋放速率與化學反應中的組分生成和消耗。

在動量方程中的雷諾應力項需要通過模型來進行封閉。在經典的Boussinesq假設中,引入了湍流黏性系數,而k-ε模型正是通過求解湍流動能k與湍流耗散率ε的輸運方程來計算湍流黏性系數[11]

能量及組分方程

式中:Cμ為系數。對于考慮浮力修正的Standard k-ε模型

式中:σk和σε分別為k和ε的湍流Prandtl數,為常系數;SK和Sε為湍流動能與湍流耗散率輸運方程的源項,表示為

式中:Gk和Gb分別為黏性應力和浮力對湍流動量能的產生項

式中:Sij為應變率張量,Prt為湍流 Prandtl數,C1ε和C2ε為影響湍流耗散率的常系數,C3ε在浮力對湍耗散的影響中發揮著更大的作用

式中:v為與重力矢量相平行的速度分量;u為與重力矢量相垂直的速度分量。

因此,在浮力剪切層中當平均流方向與浮力方向平行時,C3ε=1;當平均流方向與與浮力方向垂直時,C3ε=0。

對于考慮浮力修正的RNG k-ε模型,其輸運方程的形式與Standard k-ε模型類似,除Cμ、σk和σε的系數取值與Standard k-ε模型有差異外,最主要的不同在系數C2ε上

式中:η=Sk/ε,β 為常系數,η0為常數。與 Standard k-ε模型相比,RNG k-ε模型更適合于低Re及剪切流的計算。

對于考慮浮力修正的Realizable k-ε模型,其主要差異在于渦耗散方程中的源項

式中:C1為η的函數,C2為常系數,v為動力粘性系數。此外,計算湍流黏性系數的Cμ不再為常數,而是通過求解有關應變率張量及渦量張量的函數進行計算。因此,Realizable k-ε模型更適合計算含有大的應變率的自由剪切流。

1.2 湍流燃燒模型

對于火災模擬,最簡單的燃燒模型是體積熱源模型,即不考慮詳細的燃燒反應過程,而將反應產生的熱量以1個體積熱源的方式設置到計算區域內,來模擬熱羽流的輸運及傳熱過程。體積熱源模型不包含組分方程,因而無法預測流場中的組分分布。本文純浮力羽流算例采用該模型。

渦破碎模型需要求解反應物及產物的組分方程,并且需要顯性地輸入單步或多步的化學反應機理。在該模型中,作為反應物的燃料與氧氣一旦接觸便進行快速燃燒反應。模型中燃燒反應速率正比于湍耗散與湍動能的比值ε/k,這樣將反應物與產物的反應速率與湍流的影響關聯起來。燃燒反應對流動作用主要通過在能量方程及組分方程中加入源項的方式體現[12]。

PrePDF模型首先求解混合分數和混合分數的方差的輸運方程,然后通過混合分數來計算組分的濃度、反應熱等。由于化學反應的時間尺度遠小于湍流混合的時間尺度,該模型也假設反應速率無限快[13-14]。prePDF模型也會考慮湍流對燃燒反應的影響。對于1個簡單的燃料與氧氣的反應系統,混合分數可以表示為

式中:σt為湍流 Prandtl數;Cg、Cd為常系數。混合分數方差主要用于描述湍流對化學反應的影響,并用于封閉方程。一旦該守恒變量求解后,其他熱力學標量均可通過混合分數及平衡方程計算出來。但是,對于湍流燃燒反應,不僅需要計算瞬態值,還需要計算時均值。因此,該模型還將引入概率密度函數PDF,來考慮湍流與燃燒反應之間的相互作用。

該PDF函數描述的是混合分數在f與f+Δf之間脈動的時間概率

式中:τi為混合分數f的值出現在Δf區間內的時間分數。該概率密度函數的分布源于混合分數的脈動特性。一旦空間位置上概率密度函數確定后,可作為1個權重函數來計算組分濃度、密度及溫度的時間平均量,如以下積分方程

一般而言,概率密度函數滿足β函數分布。prePDF模型能夠考慮中間產物及湍流對燃燒反應的影響,并且由于其避免了求解每一組分輸運方程,能夠極大提高計算效率。因而,prePDF模型在火災數值模擬中的應用較為廣泛。在本文含有燃燒反應的模擬計算中采用該模型。

1.3 輻射模型

在發動機短艙火災中,輻射傳熱是重要的熱傳導方式。在考慮燃燒反應的情況下,也會考慮輻射對湍流燃燒的影響。這里選擇的輻射模型是Discrete ordinates(DO)模型[15]。該模型在流體網格的基礎上,針對每個網格劃分出若干個體積角,基于體積角求解關于輻射強度的輸運方程。采用Weighted-sum-ofgray-gases model(WSGGM)加權平均的灰氣體假設模型來計算混合物的吸收率。同時,由于火災中的煙顆粒對于散射效果而言較小,因此在模擬中不考慮內散射與外散射項。

2 數值方法

本文的數值計算采用ANSYS Fluent 14.5,該軟件基于非結構化網格采用有限體積方法求解守恒方程。計算采用壓力求解器,壓力速度解耦采用Coupled的數值方法。并采用Pseudo瞬態松弛技術,該技術在穩態計算中對主要標量加入1個瞬態的時間步來控制松弛因子,進而提高數值穩定性[10]。壓力的空間離散采用Body-force weighted格式,其他標量采用2階迎風格式。對于收斂標準,在浮力羽流計算中要求所有變量的殘差小于10-5,在含有燃燒反應的計算中要求所有變量的殘差小于10-3。

3 模擬結果與試驗驗證

3.1 浮力羽流

本節選取George等[16]的浮力羽流試驗作為浮力修正湍流模型的校驗數據。試驗的環境溫度為302 K,環境壓力為101325 Pa,熱源的直徑為6.35 cm,熱源的溫度穩定在573 K,熱源口的流速為67 cm/s,對應的熱源口的Re=870,Fr=1.4。試驗確定在熱源口2倍熱源直徑以上的位置,浮力羽流由層流轉捩為湍流。

在Fluent模擬中,采用2D軸對稱網格,邊界條件的選取如圖1所示。其中,計算區域為3 m×1 m;熱源采用速度入口,湍流強度為0.5%,湍流長度尺度為D0/15(D0為熱源直徑);網格總量為4000,軸向為100個網格,徑向為40個網格。同時,采用200×80、400×160的網格檢驗了網格獨立性,平均軸向速度與溫度分布的差異在三者之間小于1.0%,因此選擇了100×40的網格。

Shabbir與George針對各自試驗中的無量綱平均軸向速度與無量綱浮力,擬合了如下公式

式中:W為平均軸向速度;η=r/z,為相似變量;z為軸向高度;r為距羽流中心軸線的距離;ΔT為羽流中溫度與環境溫度的差;β為熱膨脹系數;F0為熱源處的特征浮力

該擬合公式適用于浮力羽流的自相似區,區域內羽流特征參數不受軸向位置影響。因此,本節選擇該區域內的z=1.75 m位置作為比較截面,比較浮力修正的 Standard k-ε、RNG k-ε 及 Realizable k-ε 模型的模擬結果與擬合曲線的差異。

無量綱平均軸向速度與無量綱浮力在z=1.75 m位置上的徑向分布分別如圖2、3所示。對軸向速度的比較,Standard k-ε模型與RNG k-ε模型的預測結果相近,徑向擴展上的結果好于Realizable k-ε模型,但Realizable k-ε模型對軸線上峰值的預測結果最好;對反映溫度差的無量綱浮力的比較,RNG k-ε模型與Realizable k-ε模型在徑向擴展上的預測結果較好,但對軸線峰值的預測,Realizable k-ε模型模擬值比試驗擬合值高,Standard k-ε模型與RNG k-ε模型模擬值比擬合值低。

3.2 燃燒器火焰

本節采用Purdue大學甲烷燃燒器的試驗數據[17-18],比較驗證prePDF模型在浮力驅動的火焰模擬上的適用性。該擴散燃燒器的直徑D=7.1 cm,燃料的組分為92.2%的甲烷、3.3%的乙烯、3.9%的氮氣、0.6%的二氧化碳,燃料的質量流量為84.3 mg/s,燃燒器出口流速為3.14 cm/s,因而相應的Fr=0.109。試驗中的環境溫度為288 K,可見火焰高度為36.4 cm,在燃料充分燃燒的情況下火源的熱釋放速率為4.2 kW。

在Fluent數值模擬中,采用3D結構化網格,火源為質量入口,其他邊界為壓力出口邊界;雖然在浮力羽流模擬中Realizable k-ε模型的預測結果與試驗值更貼近,但在考慮燃燒過程的算例中Realizable k-ε模型沒有獲得收斂結果,因此,本算例中湍流模型采用浮力修正的RNG k-ε模型;燃燒模型為prePDF模型、輻射模型為DO模型;網格數量為130萬,最小網格尺度為1 mm,計算區域為直徑30 cm、高40 cm的圓柱體。

在Z/D=0.07、1.41時軸向截面上的周向平均溫度與周向平均混合分數的徑向分布分別圖4、5所示。通過對比可知,在近火源截面上(Z/D=0.07),軸線中心的溫度較低,火焰面在0.030~0.035 m區域,模擬結果中的中心溫度及火焰面位置與試驗值符合,但最高溫度低于試驗值;在遠火源截面上(Z/D=1.41),火焰中心溫度較高,模擬值與試驗值符合較好。混合分數的分布決定了火焰的熱釋放速率以及溫度的分布,在遠火源截面上的模擬結果與試驗值的符合性好于近火源截面上的結果,這也與溫度的比較結果相對應。

周向平均軸向速度分量的徑向分布比較(不同軸向截面)如圖6所示。火焰及熱羽流區域的溫度較高,由于密度差誘導的浮力推動羽流向上運動,一方面促使軸向速度加速,另一方面從周圍環境卷吸空氣,不同截面上軸向速度有一定的對稱性。圖6比較了Z/D=0.42、0.85時軸向截面上的軸向速度,中軸線上的速度從1.0 m/s加速到1.5 m/s,模擬值與試驗值基本符合,但是模擬結果低估了中軸線上的速度,主要由于中軸線上溫度的模擬值較低。

4.3 Trent 800發動機風扇艙火災

Cranfield大學針對RR公司的Trent 800型航空發動機的風扇艙,搭建了1/2縮比尺寸的試驗臺[19-20]。該試驗臺內外2個圓桶構成1個環腔,內壁與外壁的直徑分別為1220、1520 mm,軸向長度為850 mm。試驗臺配備了水冷系統,保證了試驗臺在火災試驗中的耐久性。對風扇艙內部附件與管路進行了簡化,僅包含了齒輪箱、滑油換熱器與ECU。通風口在風扇艙上部,左右各1個成45°入射角,排風口在底部偏20°方位,火源在底部正中央,面積為128 cm2,燃料為丙烷,出口流速約為20 cm/s,火源功率相當于100 kW,用于模擬齒輪箱滑油泄漏后形成的油池火。試驗中在距風扇艙底部順時針方向 60°、120°、150°3 個位置上,布置了熱電偶陣列,用于測量這3個周向截面上的溫度分布。試驗臺外觀與結構如圖7所示。

在Fluent數值模擬中,采用3D非結構化網格,網格數量為120萬,最小網格尺度為2 mm,模型的選擇與Purdue大學甲烷燃燒火焰模擬一致。在試驗中由于艙內附件布局的不對稱性,火焰及熱羽流沿順時針方向傾斜,高溫羽流最遠撞擊到滑油換熱器一側,模擬結果與試驗觀察相符。

平均溫度在軸向上的分布比較(不同測量站)如圖 8所示。圖中比較了 60°、120°、150°3個測量站分別距風扇艙內壁徑向距離為37.5、75 mm位置上平均溫度的試驗值與模擬值。在靠近火源的60°測量站上,模擬結果在軸向上的分布與試驗值相符,但是高估了火焰中心位置上的溫度,并且模擬的火焰寬度低于試驗值。在120°、150°2個測量站上,模擬結果與試驗結果符合較好,尤其在軸向位置大于0.3 m的區域,但在小于0.3 m的高溫區域,模擬結果比試驗值低。其主要原因在于火源區域火焰的擴展上模擬結果低于試驗值,這與羽流對空氣的卷吸相關。同時,浮力修正的湍流模型的預測能力不足,以及在模擬中未考慮煙氣模型,可能是導致模擬偏差的主要原因。

總之,模擬結果能夠基本預測風扇艙發生火災后的溫度、速度等主要參數的空間分布,并描述主要的物理特征,該模擬方法將為實際工程設計提供一定的模擬數據支撐。

4 結論

利用Fluent模擬了熱羽流、甲烷燃燒器火焰及RR公司Trent 800發動機風扇艙火災,并與試驗值進行了比較驗證,得到了如下結論:

(1)對于熱羽流的模擬,經過浮力修正的Standard k-ε、RNG k-ε 及 Realizable k-ε 模型均能較好地預測羽流速度、溫度的分布,相比而言,Realizable k-ε模型預測結果更貼近試驗值;

(2)對甲烷燃燒器火焰的模擬,prePDF燃燒模型能夠較好地預測混合分數、溫度及速度的分布,不過在近火源區預測的溫度峰值與試驗值差距較大,進而也低估了中軸線上的速度;

(3)在Trent 800發動機風扇艙火災的模擬中,能夠基本模擬火焰及煙氣的溫度分布,在遠離火源的區域模擬結果較好,但是在近火源區域模擬值與試驗值仍有較大偏差;湍流、燃燒、輻射模型及數值方法的組合在一定程度上能夠描述發動機艙火災的主要物理過程,為發動機艙防火通風系統的設計及附件的布局提供技術支持;但是,由于發動機艙環境復雜,要想準確模擬其火災發生發展過程,仍具有較大的技術挑戰;

下一步將著重研究浮力修正湍流模型和燃燒模型,更好地描述該受限空間中的非充分燃燒過程,提高對發動機艙火災的預測能力。同時,在火災模擬中考慮煙氣對燃燒及熱輻射的影響,進一步提高模擬預測的精度。

[1]中國民用航空局.CCAR-33-R2航空發動機適航規定 [S].北京:中國民用航空局,2011:8-9.Civil Aviation Administration of China.CCAR-33-R2 Aircraft engine airworthiness regulations[S].Beijing:Civil Aviation Administration,2011:8-9.(in Chinese)

[2]中國民用航空局.CCAR-25-R4中國民用航空規章運輸類飛機適航標準[S].北京:中國民用航空局,2011:123-127.Civil Aviation Administration of China.CCAR-25-R4 Civil aviation regulations of China,transport aircraft airworthiness standards[S].Beijing:Civil Aviation Administration 2011:123-127.(in Chinese)

[3]Mullender A J,Coney M H.The numerical simulation and experimental validation of ventilation flow and fire events in a Trent nacelle fire zone[C]//ICAS Congress,International Council of the Aeronautical Sciences,2000.

[4]Dewan A.Tackling turbulent flows in engineering[M].Berlin:Springer-Verlag,2011:81-85.

[5]Markatos N C,Malin M R,Cox G.Mathematical modeling of buoyancy-induced smoke flow in enclosures[J].International Journal of Heat and Mass Transfer,1982(25):63-75.

[6]Nam S,Bill R G.Numerical simulation of thermal plumes[J].Fire Safety Journal,1993(21):231-256.

[7]Kumar R,Dewan A.Assessment of buoyancy-corrected turbulence models for thermal plumes[J].Engineering Applications of Computational Fluid Mechanics,2013(7):239-249.

[8]Alexandre N.CFD fire simulation in enclosures[D].Lisboa:Universidade Tecnica de Lisboa,2007.

[9]Xue H,Ho J C,Cheng Y M.Comparison of different combustion models in enclosure fire simulation[J].Fire Safety Journal,2011(36):37-54.

[10]ANSYS Fluent 14.5 theory guide documentation [R/OL].http://support.ansys.com/docdownloads.ANSYS,Inc.,2012.

[11]Rodi W.Turbulence models and their applications in hydraulics:a state-of-the-art review[M].Netherlands:International Association for Hydraulic Research,1984:23-25.

[12]Magnussen B F,Hjertager B H.Mathematical models of turbulent combustion with special emphasis on soot formation and combustion[C]//16th Symposium(International)on Combustion,Cambridge,MA,1976:719-729.

[13]Sivathann Y R,Faeth G M.Generalized state relationship for scalar properties in non-premixed hydrocarbon/air flame[J].Combustion and Flame,1990(82):211-230.

[14]Jones W P,Whitelaw J H.Calculation methods for reacting turbulent flows:a review[J].Combustion and Flame,1982(48):1-26.

[15]Chui E H and Raithby G D.Computation of radiant heat transfer on a non-orthogonal mesh using the finite-volume method[J].Numerical Heat Transfer,Part B.1993(23):269-288.

[16]A.Shabbir,W.George.Experiments on a round turbulent buoyant plume[J].Journal of Fluid Mechanics.l994,275:1-32.

[17]Xin Y,Gore J P,McGrattan K B,et al.Fire dynamics simulation of a turbulent buoyant flame using a mixture-fraction-based combustion model[J].Combustion and Flame 2005(141):329-335.

[18]Xin Y,Gore J P,McGrattan K B,et al.Large eddy simulation of buoyant turbulent pool fires[J].Proceedings of the Combustion Institute,2002(29):259-266.

[19]Fasquelle A,Rubini P A.Experiments and transient numerical simulations of fire events in a gas turbine nacelle [C]//2nd International Conference,Fire Bridge.Belfast,2005.

[20]Moss J B,Rubini P A,Sanderson V,et al.The simulation of engine casing fires:computational modelling and experiment [C]//NATO RTO Fire and Survivability,Specialists Meeting.Aalborg,Denmark:Applied Vehicle Technology Panel,2002:23-26.

Application on Buoyancy Corrected Turbulent Model in Nacelle Fire Modeling

LI Song-yang,SHAO Xing-chen
(AECC Commercial Aircraft Engine Co.,Ltd,Shanghai 200241,China)

Since nacelle fire is a typical thermal driven buoyancy plume,the study explored the simulation method of buoyancy flow on a benchmark test of thermal plume,compared three kinds of buoyancy corrected two-equation turbulence model.The Purdue methane burning experiment was verified by using one turbulence model together with prePDF combustion model.The fire tests based on 1/2 reduced scaling model of Trent 800 nacelle,where the pool fire was caused by oil leakage using RANS method,and was adopted to verify the accuracy of the modeling method.The simulation reproduced main physical processes of nacelle fire,and a good agreement on hot air velocity and temperature was obtained against testing results.The main factors affecting simulation accuracy were further discussed,which is whether turbulence model together with combustion model can accurately capture the flame front near the fire source plays direct effect on air entrainment as well as the temperature and velocity of the downstream plume.The CFD technology is applicable for optimizing nacelle fire protection design and engine accessories layout.

fire protection;fire modeling;buoyancy plume;turbulent combustion;aeroengine

V 228.6

A

10.13477/j.cnki.aeroengine.2017.01.011

2016-05-21 基金項目:上海市聯合創新研究項目資助

李松陽(1983),男,博士,工程師,從事計算燃燒學、發動機防火研究工作;E-mail:lisongyang@acae.com.cn。

李松陽,邵興晨.浮力修正湍流模型在航空發動機火災模擬中的應用 [J].航空發動機,2017,43(1):58-65.LI Songyang,SHAOXingchen.Applicationonbuoyancycorrectedturbulentmodelinnacellefiremodeling[J].Aeroengine,2017,43(1):58-65.

(編輯:張寶玲)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产美女精品人人做人人爽| 国产成人成人一区二区| 国产一区二区视频在线| 亚洲69视频| 久久久无码人妻精品无码| 5555国产在线观看| 国产簧片免费在线播放| 国产肉感大码AV无码| 久久国产乱子| 无码内射在线| 亚洲制服丝袜第一页| 香蕉综合在线视频91| 久久 午夜福利 张柏芝| 国产香蕉在线视频| 99热这里只有免费国产精品 | 中文字幕无码av专区久久 | 中文字幕久久波多野结衣| 国产精品性| 亚洲综合色婷婷中文字幕| 久久精品亚洲专区| 亚洲视频a| 一区二区午夜| 国产成人精品免费视频大全五级| 欧美一区二区丝袜高跟鞋| 超级碰免费视频91| 欧美日韩成人在线观看| 欧美日韩久久综合| 手机精品视频在线观看免费| 国产美女精品人人做人人爽| 久久天天躁狠狠躁夜夜2020一| 日韩AV无码免费一二三区| 极品尤物av美乳在线观看| 国产精品香蕉在线观看不卡| 乱系列中文字幕在线视频| 日韩av高清无码一区二区三区| 国产高颜值露脸在线观看| 成人中文在线| 日韩成人免费网站| 亚洲开心婷婷中文字幕| 91久久青青草原精品国产| 五月六月伊人狠狠丁香网| 亚洲国产成人超福利久久精品| 亚洲精品无码高潮喷水A| 乱人伦视频中文字幕在线| 亚洲精品制服丝袜二区| 一区二区影院| 熟妇丰满人妻| 91www在线观看| 日本免费一级视频| 免费99精品国产自在现线| 国产青榴视频在线观看网站| 成人亚洲国产| 国产成人乱码一区二区三区在线| 91精品小视频| 国产高清在线观看91精品| 欧美色图久久| 国产精品区视频中文字幕| 国产成人精品一区二区不卡| 中文字幕无线码一区| 乱人伦99久久| 国产午夜福利在线小视频| 日韩精品毛片人妻AV不卡| 欧美综合中文字幕久久| 国产成人亚洲精品色欲AV | 欧美在线精品怡红院| 日韩国产黄色网站| 国产精品视频公开费视频| 91av成人日本不卡三区| 思思热精品在线8| 伊人91视频| 人妻中文久热无码丝袜| 成人亚洲天堂| 亚洲成人福利网站| 国产视频一区二区在线观看| 亚洲激情99| 欧美国产日韩在线播放| 91无码人妻精品一区二区蜜桃| 亚洲人妖在线| 在线观看国产网址你懂的| 2022国产91精品久久久久久| 亚洲天堂网在线视频| 青青草国产免费国产|