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

圓柱-分離盤結構的流致旋擺響應數值研究

2023-08-03 13:53:14朱紅鈞周新宇陳泉宇
空氣動力學學報 2023年6期
關鍵詞:結構

唐 濤,朱紅鈞,周新宇,陳泉宇

(西南石油大學 石油與天然氣工程學院,成都 610500)

0 引言

流致旋擺現象普遍存在于生活和工程中,如海洋風機、水車等。Maxwell[1]最早研究了流致旋轉現象,指出質量中心和氣動力中心不重合導致的扭矩作用是引起結構物旋轉的原因。早期學者主要研究的是平板流致旋轉[2-5],20 世紀90 年代開始才出現使用自由旋轉分離盤控制圓柱旋渦脫落及水動力的研究,其中最廣泛、最核心的是對分岔(bifurcation)現象的探究和解釋。物體旋轉會產生不對稱尾流,進而造成橫向作用力,稱為馬格努斯效應[6]。分岔現象就是馬格努斯效應的典型例子:反對稱破壞機制。即從反對稱尾渦到不對稱尾渦的轉變[7],造成近尾流區(qū)作用在圓柱-分離盤結構上的壓力關于來流不對稱,從而產生凈升力作用,進一步導致凈流體力矩,使分離盤偏移到一側;其旋轉平衡位置與來流存在非零夾角θmean,該夾角與無量綱分離盤長度L*(L*=L/D,其中L為分離盤長度,D為圓柱直徑)、雷諾數Re密切相關。眾多研究[8-11]表明:在亞臨界雷諾數范圍內,θmean與雷諾數無關,主要受分離盤長度影響,可分為三個階段:L*≤ 1 時,隨L*增大θmean驟降;1 <L*≤ 4 時,θmean緩慢減小至零;L*> 4 時,θmean基本為零,此時分離盤旋轉平衡位置與來流平行。但在低雷諾數Re=50 時,Xu 等[12]觀察到θmean在L*≈ 1.7D時就減小到零。

除分岔現象外,分離盤在新平衡位置的擺動大小也是一個研究重點。Xu 等[13]研究了低雷諾數下L*=1 時分離盤的無量綱擺幅大小(定義為盤尖端擺動弧度與圓柱直徑之比),發(fā)現分離盤的無量綱擺幅很小,僅有1 ×10?3量級;Re< 48 時,擺幅為零,分離盤處于靜止狀態(tài);48 <Re< 70 時,擺幅快速增大;Re> 70 時,擺幅緩慢增大,而后基本趨于平穩(wěn)。Shukla等[14]較為系統(tǒng)地分析了L*和Re對分離盤擺動的影響,發(fā)現在Re< 4 000 時,分離盤擺幅隨雷諾數增大而增大;Re超過4 000 后,擺幅基本不變,表明雷諾數的影響減小。相比于低雷諾數,分離盤在高雷諾數范圍的擺幅顯著增大,擺動響應可以劃分為兩個部分,且二者之間存在跳躍現象:第一個部分對應L*≤ 3,擺幅較大,且以單一頻率擺動;第二個部分對應L*≥4,擺幅較小,多頻響應明顯。

至今為止,關于圓柱-分離盤結構的流致旋擺響應研究多見于高雷諾數的實驗分析,對于低雷諾數范圍內的旋擺響應研究相對較少:Xu 等[12-13]探究了Re=40~100 范圍內的圓柱-分離盤結構的擺動特性,但未見流場細節(jié)和水動力的變化規(guī)律;Lu 等[7]雖然分析了流場的變化,但雷諾數僅限于Re=100。因此,針對圓柱-分離盤結構的流致旋擺響應,拓寬雷諾數范圍并進一步闡明其中的流動現象和機理具有顯著意義?;诖耍疚拈_展了低雷諾數Re=40~160,L*=0.5、1.0、1.5、2.0 時,圓柱-分離盤結構的流致旋擺響應數值研究,重點分析了結構旋擺特性、流場細節(jié)及水動力系數。

1 物理模型

如圖1 所示,分離盤固結于圓柱后方,二者在流體力作用下同步旋擺。圓柱直徑為D,分離盤長度L*=0.5、1.0、1.5、2.0,分離盤厚度w=0.2D。初始時刻,分離盤長邊與來流平行;發(fā)生旋擺后,分離盤與初始位置的夾角定義為θ(t),t為流動時間。本文采用矩形計算域,其中上游速度入口邊界、兩側對稱邊界和下游壓力出口邊界與圓心的距離分別為10D、10D和30D。本文采用重疊網格方法將計算域切分為背景網格部分和前置網格部分,其中前置網格部分是直徑為15D的同心圓,圓的外側邊界條件設置為交界線(overset),用于插值和數據傳遞。

圖1 物理模型與計算域Fig.1 Skematic of the model and computational domain

2 數值方法

2.1 控制方程

低雷諾數繞流流場由非定常不可壓縮Navier-Stokes(N-S)方程求解,包括連續(xù)性方程(1)和動量方程(2),其無量綱形式為[15]:

式中:u*為笛卡爾坐標系下的無量綱流動速度,包括流向速度u*和橫向速度v*,u*=u/U,v*=v/U;t*為無量綱流動時間,t*=Ut/D;p*為無量綱壓力,p*=p/ρU2,其中p為實際壓力,ρ為流體密度;Re=UD/ υ,其中U為來流速度,υ為流體運動黏度。采用有限體積法求解上述控制方程,對流項采用二階迎風格式離散,壓力和速度耦合采用Coupled 算法。

基于牛頓第二定律,不考慮剛度和阻尼作用,結構流致旋擺響應的控制方程可寫為[13]:

式中:I*為結構的無量綱質量慣性矩,I*=I/(ρD4),其中I為結構的單位長度質量慣性矩,本文取I*=10;M*為作用于結構的無量綱扭矩,M*=M/(ρU2D2),其中M為扭矩。

迭代計算時,結構的運動控制通過用戶自定義函數(user defined function,UDF)實現:在一個時間步內,首先利用FLUENT 軟件求解N-S 方程獲得流場信息,而后通過積分得到作用于圓柱-分離盤結構上的力矩,根據UDF 執(zhí)行當前時間步內發(fā)生的旋轉響應,網格自適應調整,準備下一次更新迭代。

2.2 計算網格

如前所述,本文采用重疊網格技術,以L*=2 為例,圖2 展示了本文的網格劃分策略。整個計算域均以四邊形單元填充,交界線附近的兩套網格尺寸接近,保障了子網格間的信息傳遞準確率。在圓柱-分離盤結構表面,進行網格加密處理,而在遠場則采用較為稀疏的網格單元。

圖2 數值計算網格Fig.2 Computational mesh

在計算前,首先進行無關性驗證。選取L*=2、Re=160 時的圓柱-分離盤結構,觀測其旋擺平衡角θmean和升力系數均方根值CL,rms的變化。如表1 所示,隨著網格總數增加(增加圓柱表面節(jié)點數),θmean和CL,rms逐步收斂于M3;繼續(xù)增加網格,M3 和M4 之間的結果相對誤差已小于1%??紤]到計算成本,選擇M3 網格。在M3 基礎上,繼續(xù)開展第一層網格高度、計算域大小和時間步長的無關性驗證,結果分別如表2、表3 和表4 所示。以兩次結果間的相對誤差小于1%為判據,同時兼顧計算成本,得到本文使用的圓柱表面第一層網格高度為0.02D(y+=0.5);相應地,分離盤表面的第一層網格沿側邊逐漸增大,但最大不超過0.1D。經過驗證,本文的計算域選擇為40D× 20D,滿足阻塞率小于6%的要求[16-19];計算時間步長為0.002 s,滿足最大庫朗數小于0.2 的要求。最后,如表5 所示,將無關性驗證得到的結果與前人研究結果進行對比,發(fā)現百分比差異不超過5%,證明本文無關性驗證結果是可靠的。

表1 關于圓柱表面節(jié)點數的網格無關性驗證Table 1 Independence verification of the mesh size

表2 關于第一層網格高度的網格無關性驗證Table 2 Independence test of the height of the first layer

表3 關于計算域大小的網格無關性驗證Table 3 Independence test of the computational domain

表4 關于時間步長的網格無關性驗證Table 4 Independence test of time step

表5 本文網格無關性驗證結果與前人結果對比Table 5 Comparison of the equilibrium angle with reported results

2.3 模型驗證

首先以Re=160 時的圓柱-分離盤結構為對象,驗證了旋擺平衡角θmean隨分離盤長度L*的變化規(guī)律。Cimbala 和Chen[10]、Xu 等[13]的研究結果表明,當80 <Re< 2 × 105、L/D≤ 2 時,雷諾數對旋擺平衡位置幾乎沒有影響,所以本文驗證時選擇Re=160 具有合理性。如圖3 所示,驗證結果與前人的研究結果吻合較好。隨后,在低雷諾數范圍內驗證了旋擺平衡角θmean、擺角均方根θrms與雷諾數的關系,結果如圖4所示。值得注意的是,本文與Xu 等[12-13]在計算時采用了同樣的分離盤(L*=1),結構的流致旋擺響應控制方程也一致,因此二者結果具有相似性,而其中的差異主要來自于Xu 等[12-13]并未給出確切的質量慣性矩數值。綜上,以上結果證實了本文數值模型的可靠性。

圖3 旋擺平衡角與盤長的變化關系Fig.3 Variation of equilibrium angle with plate length

圖4 旋擺平衡角與擺角均方根隨雷諾數的變化(L*=1.0)Fig.4 Variations of the equilibrium angle and root-meansquared rotary angle against Re (L*=1.0)

3 結果分析

3.1 流致旋擺響應

初始時刻,分離盤長邊與來流平行,隨著流場不斷演變,在不穩(wěn)定渦流作用下,圓柱-分離盤結構發(fā)生偏轉[8-9],調穩(wěn)后將圍繞新的平衡位置往復旋擺,圖5 展示了旋擺平衡角隨雷諾數、盤長的變化。結果表明,當L*=2.0(Re=40、50)和L*=1.5(Re=40)時,結構的旋擺平衡角度為零,此時未發(fā)生分岔現象;而在其他工況下則發(fā)生了分岔,即結構的旋擺平衡位置與來流存在夾角。隨著雷諾數增大,旋擺平衡角度的變化經歷兩個階段:先增大,而后基本保持不變;這表明雷諾數的影響在逐步減弱。Xu 等[12]的實驗也證實了這種旋擺平衡位置與雷諾數的變化關系。本文中,兩個階段的拐點雷諾數Re*隨盤長增加而增大,分別為Re*=60(L*=0.5)、70(L*=1.0)、80(L*=1.5)、90(L*=2.0)。本文還注意到,分岔現象的臨界雷諾數Rec也與盤長有關。本文結果表明,分離盤越長,該臨界雷諾數越大。當L*=2.0 時,Rec=50;當L*=1.5 時,Rec=40;而對于更短的分離盤,在本文雷諾數范圍內暫未出現臨界雷諾數。對比盤長對旋擺平衡位置的影響發(fā)現,分離盤越短,旋擺平衡角度越大,且隨著盤長增加,θmean的縮減速率在逐漸減弱。這是因為當分離盤較短時,結構兩側的剪切層依然能夠在近尾流區(qū)相互作用,交替脫落的旋渦提供了較大的壓差,導致結構偏轉;而較長的分離盤則可以有效阻斷兩側流體的交互作用,將旋渦脫落向下游推移,從而減弱旋擺[8]。

圖5 旋擺平衡角隨雷諾數、盤長的變化Fig.5 Variations of the equilibrium angle versus Re and L*

除分岔現象外,旋擺幅度也是表征圓柱-分離盤結構流致旋擺響應的一個重要參數,本文選用的是結構達到旋擺平衡后,旋擺角度均方根值θrms。如圖6所示,隨著雷諾數增大,θrms首先保持為零(結構為靜止或基本靜止狀態(tài)),而后快速增大,最后緩慢增大或基本保持不變。與高雷諾數結果[14]相比,低雷諾數的擺幅更小,但變化趨勢相近。盤長不僅對θrms≠ 0的起始雷諾數產生影響,還對擺幅大小產生影響??傮w而言,分離盤越長,該起始雷諾數越大,結構的擺幅也越大。此外我們還注意到,L*=2.0、Re> 100時,θrms近乎呈線性增大;而盤長減小,擺幅增大的速率也減小。

圖6 擺角均方根隨雷諾數、盤長的變化Fig.6 Variations of the root-mean-squared rotary angle versus Re and L*

探究結構的擺動頻率有助于理解結構的穩(wěn)定性能。圖7 展示了結構的無量綱旋擺頻率f*(f*=fD/U,其中f為實際的旋擺頻率)。結果表明,隨著雷諾數、盤長的變化,當L*=0.5 和1.0 時,無量綱旋擺頻率隨雷諾數增加呈近似線性增長,這表明相較于分離盤長度,雷諾數對旋擺頻率的影響更明顯。繼續(xù)增加盤長,發(fā)現旋擺頻率在Re=40~130 范圍內依舊呈近似線性增長,但之后增長速度變小。此外,縱向對比盤長對旋擺頻率的影響發(fā)現,分離盤越長,無量綱旋擺頻率越小,這是因為較長的分離盤與流體作用面積更大,在橫向上受到的阻力越大。

圖7 無量綱旋擺頻率隨雷諾數、盤長的變化Fig.7 Variations of the non-dimensional rotation frequencies versus Re and L*

3.2 流場特性

根據本文結果,可以將圓柱-分離盤結構的流致旋擺響應分為3 類:分岔后大幅度往復擺動、分岔后基本靜止、不分岔。圖8~圖11 分別展示了3 類響應模式的流場特性。圖8 以L*=1.0、Re=100 為例,展示了結構旋擺發(fā)生分岔后一個周期內的流場演變。在Ⅰ時刻,結構旋擺至最大角度,一個小尺寸順時針旋轉的旋渦A 停留在分離盤上側;同時在分離盤下側觀察到另一個較大尺寸逆時針旋轉的旋渦B。從Ⅰ時刻到Ⅳ時刻,旋渦B 從分離盤下側泄放并逐漸向下游遷移,而旋渦A 始終位于分離盤上側;在這一過程中,分離盤上側的壓力幾乎不改變,而下側的壓力變大,由此旋擺速度逐漸減小,并在Ⅳ時刻發(fā)生轉向。此外,在Ⅳ時刻,一個新的順時針旋轉的旋渦C 在分離盤尾端出現;旋渦C 是由于邊界層從圓柱表面分離后再附著于分離盤上側,最后在板尖端發(fā)生二次分離引起的[7]。在后半個旋擺周期內,旋渦C 尺寸逐漸增大并從板尖端脫落,并伴隨有在分離盤下側形成的、新的逆時針旋轉旋渦D。在一個旋擺周期內,旋渦A 始終存在,這實質上是圓柱上側分離點與分離盤上再附著點之間的再循環(huán)區(qū)域。圖8 的結果還顯示了在一個周期內發(fā)生了兩次旋渦脫落;但與裸圓柱不一樣的是,在圓柱-分離盤結構中兩個旋渦分別從盤尖端和盤下側脫落,并非從圓柱兩側交替脫落。這種旋渦脫落模式與文獻[7]的結果一致,因此在普遍意義上也可以稱為“2S”模式。

圖8 分岔后大幅度往復擺動時流場演變(L*=1.0、Re=100)Fig.8 Post-bifurcation evolution of flow fields around the model with large oscillations (L*=1.0 and Re=100)

為進一步確認上述旋渦脫落模式,圖9 展示了一個周期內的渦量場分布,選取的代表性時刻Ⅰ~Ⅷ與圖8 中的時刻一一對應,粉紅色線條為u=0 等值線,用以表征回流區(qū)大小及分離點和再附著點位置。圖8中的旋渦A 導致了分離盤上側回流區(qū)始終存在,分離點和再附著點分別位于圓柱上側和分離盤上側。相比之下,下側回流區(qū)尺度增大,再附著點位置不再固定。在前半個周期內(時刻Ⅰ~Ⅳ),下側旋渦從分離盤下側開始脫落并逐漸向下游遷移,此時再附著點位于分離盤尖端。而在后半個周期內,旋渦形成和泄放位置位于分離盤尖端,且分離盤下側的旋渦尚未發(fā)展充分,在二者的夾擊作用下,再附著點位置從分離盤尖端轉移至分離盤下側。在一個旋擺周期內,從分離盤尖端和下側分別脫落了一個旋渦,因此渦脫模式為“2S”。

圖9 一個周期內的渦量場分布(L*=1.0、Re=100)Fig.9 Vorticity fields during one cycle (L*=1.0 and Re=100)

當旋擺發(fā)生分岔但結構基本保持靜止時(擺幅不超過0.01°),流場處于擬穩(wěn)定狀態(tài),這便是第二類響應模式。圖10 展示了該種響應模式下的流場特性,云圖上的粉紅色實線代表流函數等值線,用以表征旋渦的尺寸;藍色點劃線代表u=0 等值線。結果表明,結構周圍始終存在擬靜止的旋渦A、B 和C,分別位于分離盤上側、盤尖端及分離盤下側,尺寸的相對大小關系為C > A > B;在一個旋擺周期內,三個旋渦的位置并未發(fā)生改變,而僅在尺度上有微小變化。因此,由u=0 等值線確定的回流區(qū)大小、分離點和再附著點位置也不變。從力矩角度出發(fā),圓柱上的壓力始終通過轉軸中心,對總力矩無貢獻;分離盤表面的切應力所對應的力臂僅為0.1D,對總力矩貢獻可忽略。因此,本文僅考慮圓柱上的切應力和分離盤上的壓力作用。如圖10(c)所示,分離盤表面的壓力系數Cp(Cp=(p?p)/(0.5ρU2),其中p∞為入口處的參考壓力值)在結構上下兩側的分布存在微小差異,這是引起結構小幅旋擺的主要原因;而圓柱表面的切應力系數 τ*(τ*=τ/(0.5ρU2),τ為結構表面切應力)在兩側的分布則基本相同。

圖10 分叉后基本靜止時壓力場、壓力系數與切應力分布(L*=1.5、Re=60)Fig.10 Post-bifurcation pressure fields,pressure coefficient,and shear stress when the model is quasi static (L*=1.5 and Re=60)

第三類響應模式則是結構不發(fā)生分岔現象,同時擺幅幾乎為零(圖6),這種情況類似于固定的分離盤。如圖11 所示,剪切層在圓柱表面發(fā)生第一次分離后,再附著于分離盤后側。在此過程中,分離盤兩側形成了尺寸相當的回流區(qū),進一步調控了近尾流區(qū)的流線特性,有助于穩(wěn)定流場。壓力系數與切應力分布也表明結構兩側受到的作用力相當,因此結構處于靜止狀態(tài)。這種旋擺響應模式的出現表明,雷諾數越低(流速越低),分離盤越長,越有利于再附著現象在結構兩側同時發(fā)生,進而使流場越穩(wěn)定。

圖11 無分岔時流場分布(L*=2.0、Re=40)Fig.11 Flow field without bifurcation (L*=2.0 and Re=40)

上述3 種不同的旋擺模式表明,當近尾流場處于穩(wěn)定狀態(tài)時,分離點和再附著點位置基本不變(如圖10 和圖11);然而存在交替脫落的旋渦時,尾跡被擾亂,分離點和再附著點位置也隨之動態(tài)變化(如圖9)。因此,基于umean=0 等值線的平均流場可以從整體上反映邊界層分離點和再附著位置,但掩蓋了瞬時的流動特征[20-21]。圖12 歸納了本文出現的3 種平均再附著方式。第一種方式的再附著點位于分離盤上下兩側,并且在分離點與再附著點之間形成回流區(qū),這種再附著方式主要發(fā)生在θmean較大(圖5 中的L*=0.5)或者θmean較小但θrms較大時。即,在圖12中,第一種再附著的控制范圍隨著盤長增加逐漸向高雷諾數轉移;也就是旋擺平衡角減小時,需要較大的擺幅來實現第一種再附著現象,此時旋渦僅在下側分離盤的前部形成和脫落,并未觸及整個分離盤。第二種再附著方式與第一種的區(qū)別在于,分離盤下側的再附著點轉移至盤尖端。當分離盤較長(旋擺平衡角較?。┣覕[幅較小時,分離盤兩側的旋渦得以沿盤發(fā)展,同時分離盤的運動不至于擠壓旋渦并使其強制脫落。因此對比第一種方式,第二種再附著模式下,不僅分離盤上側的回流區(qū)更寬,而且整個回流區(qū)的長度也顯著增加。這種模式在小攻角的靜止圓柱-分離盤結構中也有報道[22-23]。第三種再附著方式對應于無分叉時的工況,圓柱-分離盤結構處于靜止狀態(tài),邊界層從圓柱上下表面分離后再附著于分離盤上下兩側,且上下兩側的回流區(qū)尺寸相同。

圖12 再附著模式分區(qū)圖Fig.12 Classification of reattachment behavior

3.3 水動力系數

圖13 顯示了可旋擺圓柱-分離盤結構的水動力系數隨雷諾數、盤長的變化,并與裸圓柱進行了對比。裸圓柱的阻力平均值[24-25]隨雷諾數增加不斷減小并趨于穩(wěn)定,而圓柱-分離盤結構的阻力系數則呈現相反的變化特征,即阻力系數隨雷諾數增大而不斷增大。然而,在本文所研究的Re=40~160 范圍內,旋擺分離盤帶來的減阻效果是顯而易見的,即使在Re=160 時阻力也小于裸圓柱。這可以解釋為分離盤增大了背壓,從而減小了結構前后兩側的壓差,最終減小了阻力[7-8]。此外,在Re=40~100 范圍內盤長對阻力系數幾乎沒有影響,之后隨雷諾數增大,不同盤長結構的阻力系數開始出現差異,但總體呈現的規(guī)律是分離盤越短,阻力越大。

圖13 水動力系數隨雷諾數、盤長的變化Fig.13 Variations of hydrodynamic coefficients versus Reynolds number and plate length

裸圓柱的升力系數均方根值[21,26]隨雷諾數增大幾乎呈線性增長;而圓柱-分離盤結構的升力變化呈現二次函數模式,即在Re=40~80 范圍內增長緩慢,但雷諾數超過80 后增長迅速。與阻力一樣,分離盤有助于減小結構的升力作用,這是因為分離盤在一定程度上減弱了兩側流體的相互作用[27]。同樣的,盤長的影響隨著雷諾數增大得以逐步體現,遵循與阻力變化相同的規(guī)律,即分離盤越短,升力越大。綜合水動力系數的變化來看,雷諾數的影響大于盤長的影響,同時分離盤具有顯著的減阻減升效果。

4 結束語

本文對圓柱-分離盤結構的流致旋擺響應進行了數值模擬研究,分析了雷諾數、分離盤長度對結構旋擺響應、流場特性及水動力系數的影響,主要結論如下:

1)同時觀察到了分岔和不分岔現象,但不分岔現象僅出現在L*=2.0、Re=40 和50 以及L*=1.5、Re=40 時。分岔現象的臨界雷諾數與盤長有關,分離盤越長,此臨界雷諾數越大。隨雷諾數的增大,旋擺平衡角先增大而后基本維持不變,兩個階段的臨界雷諾數也隨盤長的增加而增大。分離盤越短,旋擺平衡角越大。盤長增加時,旋擺平衡角的縮減速率逐漸減小。

2)結構擺幅隨雷諾數增大先保持為零,后快速增大,最后緩慢增大或基本不變。擺幅不為零的起始雷諾數受盤長的影響,盤長越長,該起始雷諾數越大,結構的擺幅越大。雷諾數是影響結構旋擺頻率的主導因素,二者呈正相關關系。

3)發(fā)生分岔并具有大幅度擺動時,在分離盤上側存在一個滯止渦,而在盤尖端和下側各脫落一個旋渦,呈現“2S”模式。發(fā)生分岔但結構基本靜止時,結構上下兩側受到的流體力基本相等,在分離盤上側、盤尖端及盤下側存在穩(wěn)定的回流區(qū),但三者尺度不一。不發(fā)生分岔時,在分離盤兩側形成了相同尺度的回流區(qū),結構與流場更加穩(wěn)定。

4)隨著雷諾數的增大或盤長的減小,圓柱-分離盤結構的水動力系數增大,但小于裸圓柱。

猜你喜歡
結構
DNA結構的發(fā)現
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環(huán)結構謹防“死循環(huán)”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: www.99在线观看| 免费女人18毛片a级毛片视频| 国产在线观看一区精品| 久久精品人人做人人综合试看| 亚洲二区视频| 91久久大香线蕉| 欧美日韩久久综合| 国产麻豆精品在线观看| 久久夜色精品国产嚕嚕亚洲av| 婷婷综合色| 狠狠操夜夜爽| 国产精品大白天新婚身材| 久久99国产乱子伦精品免| 国产成人一区在线播放| 激情亚洲天堂| 亚洲日韩高清无码| 美女被躁出白浆视频播放| 亚洲综合日韩精品| 东京热高清无码精品| 日韩毛片基地| 乱色熟女综合一区二区| 国产乱子伦一区二区=| 精品久久人人爽人人玩人人妻| 欧美高清三区| 国产精品无码作爱| 99偷拍视频精品一区二区| 国产成人精彩在线视频50| 国产va免费精品| 欧美日韩免费在线视频| 国产欧美视频综合二区 | 国产精品男人的天堂| 大香网伊人久久综合网2020| 一本色道久久88| 午夜福利网址| 女人av社区男人的天堂| 免费高清毛片| 91小视频版在线观看www| 亚洲欧美日韩中文字幕一区二区三区| 国产自在自线午夜精品视频| 色一情一乱一伦一区二区三区小说| 99久久性生片| 黄色网在线免费观看| 亚洲欧美日本国产专区一区| 欧日韩在线不卡视频| 国产精品lululu在线观看| 国产一在线观看| 久久久久国产一级毛片高清板| 久久综合亚洲鲁鲁九月天| 亚洲Va中文字幕久久一区| 亚洲色欲色欲www在线观看| 91精品国产一区自在线拍| 永久免费精品视频| 久久综合AV免费观看| 国产成年无码AⅤ片在线| 91网在线| 国产免费久久精品99re丫丫一| 真实国产乱子伦高清| 精品中文字幕一区在线| 亚洲成人一区二区三区| a级毛片免费播放| 亚洲欧美激情小说另类| 日韩一区二区三免费高清| 91精品久久久无码中文字幕vr| 女人18一级毛片免费观看| 欧美中文字幕无线码视频| 亚洲高清在线天堂精品| 久久国产精品嫖妓| 国产精品亚洲αv天堂无码| 91视频区| 456亚洲人成高清在线| 欧美h在线观看| 亚洲天堂福利视频| 国产精品福利导航| 亚洲欧美在线综合一区二区三区| 超薄丝袜足j国产在线视频| 国产在线91在线电影| 免费高清a毛片| 九月婷婷亚洲综合在线| 天天躁日日躁狠狠躁中文字幕| 亚洲国产精品日韩av专区| 亚洲综合专区| 精品国产www|