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

    彈性楔形體各狀態(tài)參數(shù)對入水運動性能的影響

    2014-06-12 12:13:20王文華黃一王言英翟鋼軍黃亞南1大連理工大學a工業(yè)裝備結構分析國家重點實驗室船舶工程學院深海工程研究中心遼寧大連1160242大連海洋大學航海與船舶工程學院遼寧大連116024
    船舶力學 2014年11期
    關鍵詞:結構單元楔形流場

    王文華,黃一,王言英,翟鋼軍,黃亞南,2(1大連理工大學a.工業(yè)裝備結構分析國家重點實驗室;b.船舶工程學院;c.深海工程研究中心,遼寧大連116024;2大連海洋大學航海與船舶工程學院,遼寧大連116024)

    彈性楔形體各狀態(tài)參數(shù)對入水運動性能的影響

    王文華1a,b,黃一1b,c,王言英1b,翟鋼軍1c,黃亞南1b,2
    (1大連理工大學a.工業(yè)裝備結構分析國家重點實驗室;b.船舶工程學院;c.深海工程研究中心,遼寧大連116024;2大連海洋大學航海與船舶工程學院,遼寧大連116024)

    為了分析砰擊問題的各種狀態(tài)參數(shù)對船舶入水運動性能的影響,文章采用一種新的CFD(Computational Fluid Dynamic)方法動態(tài)數(shù)值模擬了二維彈性楔形結構的自由入水過程。此方法利用液面捕捉法和直角切割網(wǎng)格系統(tǒng)解決入水過程中瞬時移動的自由液面和動邊界問題,結合彈性板條梁結構的有限元理論將算法擴展應用于彈性結構的入水特性分析。其中,根據(jù)流場的直角切割網(wǎng)格和彈性結構的板條梁單元的特點,提出適合該文算法的流固耦合策略和邊界數(shù)據(jù)傳遞方法。通過與實驗數(shù)據(jù)進行比較,證明了文中算法求解彈性楔形體入水問題的正確性和合理性。最后,建立不同狀態(tài)參數(shù)(結構材料屬性、板厚和質量、底面斜升角和入水高度等)的彈性楔形結構自由入水模型,研究了各參數(shù)對自由入水的彈性楔形結構的整體運動性能和局部變形響應的影響。

    彈性楔形結構;自由入水;液面捕捉法;直角切割網(wǎng)格系統(tǒng);有限元理論;數(shù)據(jù)傳遞方法

    1 引言

    在船舶與海洋工程領域中,會遇到流體和彈性結構相互作用而引起的各種物理現(xiàn)象,其中船舶砰擊是比較常見的一類工程問題。在船舶遭遇的各種砰擊問題中都會包含物體入水過程,并且在物體從空氣進入水中的極短時間內會產(chǎn)生高峰值的沖擊壓強,從而導致局部結構發(fā)生變形損壞[1]或者影響整體運動性能[2]。所以,有必要研究彈性結構入水過程中的流場變化和受力狀況,來預報入水結構的局部變形情況和整體運動特性,用以指導船舶的結構設計和航行操縱。

    關于彈性結構入水的研究主要有實驗方法[3-5],解析半解析方法[6]和完全數(shù)值方法[7-10]。其中,由于實驗本身對設備的硬件要求較高,并且需要比較多的人力、物力和財力,因此需要理論研究來進行輔助。在理論方法中,大規(guī)模簡化的流體模型使得解析和半解析法的應用范圍具有較大的局限性,對各種復雜的瞬態(tài)入水現(xiàn)象不能很好地同時把握,因此數(shù)值計算方法對于需要考慮多種因素共同作用的復雜入水問題顯得非常重要。并且隨著計算機硬件技術飛速發(fā)展,強大的計算能力為數(shù)值模擬提供了良好的運算平臺,能夠用來同時分析飛濺、射流、氣墊和水彈性等多種現(xiàn)象,所以目前有關入水問題的研究以數(shù)值模擬方法作為熱點。但是,在入水過程中自由液面的追蹤和動邊界的處理比較復雜,此外如果再計入流固耦合和結構水彈性的影響,那么準確地預測入水過程中彈性結構的整體運動和局部變形性能是具有挑戰(zhàn)性的。因此,到目前為止同時耦合考慮入水流場、彈性結構整體運動和局部變形相互作用的文獻并不多見。

    本文基于作者早期的研究工作采用一種新的CFD方法對彈性結構自由入水問題進行數(shù)值模擬[11]。結合液面捕捉法和直角切割網(wǎng)格系統(tǒng)的優(yōu)勢,采用液面捕捉方法將自由液面作為計算域內的接觸間斷進行處理,建立非均勻不可壓縮的歐拉方程作為流場控制方程組?;谥苯乔懈罹W(wǎng)格離散計算域,通過局部更新物體邊界附近的網(wǎng)格信息來解決動邊界問題。在此基礎上,利用基于單元中心的有限體積法對控制方程進行數(shù)值離散。將彈性板條梁結構的有限元理論與流場算法相結合,將算法擴展應用于彈性結構的入水特性分析。根據(jù)流場的直角切割網(wǎng)格和彈性結構的板條梁單元的特點,推導出適合本文算法的邊界數(shù)據(jù)傳遞方法。建立彈性楔形體的自由入水模型,將本文數(shù)值計算結果與彈性楔形板入水實驗數(shù)據(jù)[5]進行比較,證明了本文算法求解彈性楔形體入水問題的正確性和合理性。最后,通過改變計算模型的狀態(tài)參數(shù)(結構材料屬性、板厚和質量、底面斜升角和入水高度等)分析各種參數(shù)對自由入水的彈性楔形體的整體運動和局部變形的影響,從而為進一步研究船舶砰擊問題、指導船舶航行操縱和結構設計提供了技術基礎。

    2 CFD流場計算理論

    本文采用液面捕捉法追蹤自由液面,將自由液面看作是密度域內的接觸間斷、使用高精度和高分辨率的離散格式對流場密度進行數(shù)值求解,據(jù)此確定任意時刻的自由液面位置。采用此方法不僅可以忽略復雜的界面重構過程,而且能模擬破碎、重新連接等復雜的自由液面擾動問題。在此方法中,計算域包含水和空氣兩種介質,主場控制方程采用守恒形式的二維非均勻不可壓縮歐拉方程組,包括質量守恒方程:

    式中:ρ為流體密度,vx和vy分別為x、y方向的流體速度,p為流體壓強,Bx和By分別為x、y方向的質量力加速度,在重力場中,Bx=0,By=-9.8m/s2。

    計算域邊界條件設置如下:在外邊界處,為了允許流體可以無反射地自由出入外邊界,流體速度和密度滿足法向梯度為零的條件,壓強設為定值。在物體邊界處,流體速度滿足不可穿透條件,密度滿足法向梯度為零的條件。對于流體壓強來說,將動量方程(2)和(3)投影到物體邊界的法向上可得:

    式中:vf=vxi+vyj是物體邊界處的流體速度,n=nxi+nyj是邊界的法向單位向量,vsn是物體速度沿邊界法向的投影。對于彈性結構入水問題來說,方程(5)中各項的求解將在本文第4.2.1節(jié)的彈性邊界數(shù)據(jù)傳遞方法中進一步說明。

    此外,采用直角切割網(wǎng)格作為空間離散網(wǎng)格,處理物體入水的動邊界問題。此網(wǎng)格系統(tǒng)作為一種新興的網(wǎng)格技術,具有生成簡單、自適應能力強和易于實現(xiàn)自動化等優(yōu)點,并且系統(tǒng)大多數(shù)為形狀比較規(guī)則的矩形網(wǎng)格,能夠為高精度和高分辨率格式提供方便的應用平臺。另外,在處理動邊界問題時,對于直角切割網(wǎng)格系統(tǒng)來說,只需局部更新一些切割網(wǎng)格的信息即可,從而避免了傳統(tǒng)的非結構網(wǎng)格重構所帶來的誤差。網(wǎng)格生成算法和數(shù)據(jù)存儲結構的具體描述可以參照文獻[12]和[13]。

    在此基礎上,采用基于單元中心的有限體積法對流場控制方程組進行數(shù)值離散。在流體網(wǎng)格邊界處,利用Roe的近似Riemann解算子的復合形式求解流函數(shù),邊界兩側的流場變量根據(jù)二階精度Upwind型格式進行重構。此外,在重構表達式中引入Superbee限制器來調節(jié)數(shù)值耗散與色散效應,以保持格式的單調性,達到數(shù)值結果高分辨率和高保真的目的。采用人工壓縮法實現(xiàn)流場速度和壓強的耦合求解,使用一階隱式雙時間推進法進行時間步進,以獲得每一物理時間步的流場變量。

    3 結構有限元和板條梁單元理論

    在流場CFD算法的基礎上,本文進一步采用基于板條梁單元模型的有限元方法數(shù)值求解和分析自由入水楔形體彈性邊界的結構變形和響應。

    首先,建立固定在地面的總體坐標系和隨著楔形體一起運動的局部坐標系。將彈性楔形體的邊界離散成許多等間距板條梁單元,此單元可以同時承受彎曲和拉壓作用。此外,在局部坐標系中,將單元k的節(jié)點平動位移u1,w1,u2,w2和節(jié)點轉動位移θ1,θ2描述如下:

    然后,在局部坐標系中,根據(jù)板條梁單元特性建立單元質量矩陣和剛度矩陣,并將流場計算所得到的壓強分布轉化到單元節(jié)點上。再利用坐標變換的方法將單元局部坐標下的質量矩陣、剛度矩陣和節(jié)點載荷向量變到總體坐標系上,并裝配在一起形成總體坐標系下的整體動力平衡方程如下:

    此處,忽略結構阻尼;F{}是總體坐標系下的節(jié)點載荷向量,具體計算方法可以參考本文第4.2.1節(jié);在總體質量陣M[、和剛度陣K[、中,結構單元k所對應單元質量陣M[、k和剛度陣K[、k表示如下:

    式中:對于板條梁結構單元來說,ρs是結構材料密度,h是結構單元厚度,l是結構單元長度,A=1.0*h是結構單元的橫剖面面積,I=1.0*h3/12是結構單元的慣性矩,E1=E/1-μ2()是板條梁的等效彈性模量,E是結構材料的彈性模數(shù),μ是結構材料的泊松比。

    最后,引入位移邊界條件改寫方程系數(shù)矩陣,利用Newmark方法數(shù)值求解整體動力平衡方程,獲得總體坐標系下的單元節(jié)點位移、速度和加速度,從而求得彈性邊界上各點的應變。然后,將計算結果反饋回CFD求解器,用于下一時刻的流場計算。

    4 流固耦合策略和邊界數(shù)據(jù)傳遞方法

    對于彈性楔形體自由入水模型來說,因為流固耦合問題涉及到彈性楔形體的整體運動、局部變形和入水流場的相互作用,所以如何將流場和結構的求解聯(lián)系起來是非常重要的。因此,這里基于第2章的流場計算和第3章的結構有限元分析理論討論彈性楔形體自由入水的流固耦合策略,以及流場壓強和結構單元節(jié)點位移、速度在彈性邊界處的數(shù)據(jù)傳遞方法。

    4.1 流固耦合策略

    在本文算法中,將彈性楔形體自由入水模型的流固耦合問題分為入水流場和楔形體整體運動以及入水流場和楔形體局部彈性變形兩種相互作用,具體過程如圖1所示。

    一方面,為了保證足夠的數(shù)值穩(wěn)定性,采用雙向耦合方法(Kleefsman等人(2005))[14]分析入水流場和楔形體整體運動的相互作用,楔形體速度計算公式如下:

    式中:Δt是物理時間步長,m是楔形體質量,vns是楔形體在物理時間n時的垂向整體速度。ω是用來控制子迭代過程穩(wěn)定性的松弛因子,數(shù)值從0變化到1。ω越大,子迭代收斂越快。在物理時間n時,迭代初值可以取為。對于子迭代第k步,根據(jù)楔形體垂向速度可以獲得流場壓強,并且通過沿物體邊界積分可以算得流體作用力,代入到方程(10)中進行流固耦合子迭代分析。當d V=-<10-3時,子迭代過程收斂。

    圖1 流固耦合策略的流程圖Fig.1 Flow chartof iteration procedure for the fluid-structure coupled solution strategy

    另一方面,在彈性楔形體的中低速入水過程中,彈性楔形體的局部變形很小,從而局部結構變形對同一物理時刻流場的反作用可以忽略不計。因此,為了提高計算效率,這里采用單向耦合策略求解流場和局部結構變形的相互作用。

    因此,在每一物理時間步上,首先,采用雙向耦合迭代算法計算彈性楔形體自由入水的整體運動,然后使用單向耦合策略實現(xiàn)流場和彈性邊界局部變形的相互耦合,獲得彈性邊界的應變和應力,并且同時將物面邊界信息反饋到下一物理時間步的CFD求解中。

    4.2 彈性邊界的數(shù)據(jù)傳遞方法

    在上述流固耦合算法中,涉及到流場和結構計算結果的數(shù)據(jù)傳遞問題。根據(jù)流場的直角切割網(wǎng)格和彈性結構的板條梁單元的特點,兩種網(wǎng)格單元在耦合界面上并不直接匹配,從而不能直接進行數(shù)據(jù)交換,因此這里需要推導出適合本文算法的邊界數(shù)據(jù)傳遞方法。因為在流固耦合界面上傳遞的數(shù)據(jù)包括流場壓強和結構單元節(jié)點位移、速度,所以這里將數(shù)據(jù)傳遞分兩個方面加以討論,一方面是將分布式流場壓強轉化為節(jié)點集中力,另一方面是將結構響應反饋為流場計算的邊界條件。

    4.2.1 分布式流場壓強轉化為節(jié)點集中力

    圖2 流場切割和結構單元不同幾何關系的示意圖Fig.2 Different position relations of cut cell and structural elements

    根據(jù)力和力矩的平衡原則,將流場計算獲得的物面分布式壓強向結構單元節(jié)點傳遞,得到節(jié)點等效載荷向量。根據(jù)流場直角切割網(wǎng)格和彈性結構板條梁單元的位置關系,大致分為三類情況進行討論:(1)結構單元L1L2完整包含切割單元的物體邊界AB,如圖2(a)所示;(2)物體邊界AB完整包含結構單元L1L2,如圖2(b)所示;(3)物體邊界AB與結構單元L1L2部分相交。這里,僅以圖2(a)的結構單元L1L2為例說明如何根據(jù)分布在物體邊界AB上的流場壓強得到節(jié)點L1、L2的等效載荷F1和F2。

    首先,以節(jié)點L1為坐標原點在結構單元L1L2上建立局部坐標系x’o’y’,其中坐標軸o’x’沿物面逆時針方向,由L1指向L2;坐標軸o’y’垂直物面,由流體指向結構。因此o’x’和o’y’的單位向量為:

    其中:(x1,y1)和(x2,y2)分別為節(jié)點L1和L2的位置坐標。再根據(jù)流場切割單元中心O點的坐標(xo,yo),位置向量OL1可以寫成:

    最后,根據(jù)力和力矩的平衡原則將物體表面AB上的分布式壓強向結構單元節(jié)點L1和L2轉化,從而獲得結構板條梁單元L1L2的節(jié)點等效載荷F1和F2:

    式中:Po和Δp=pxi+pxj是通過流場CFD算得的切割單元中心O點處的流場壓強和壓強梯度。

    同理,采用與圖2(a)類似的方法和思路進行推導可以獲得其他兩種情況(圖2(b)和圖2(c))下結構板條梁單元L1L2的節(jié)點等效載荷。

    4.2.2 結構響應反饋為CFD邊界條件

    同樣以圖2(a)(結構單元L1L2完全包含物體邊界AB)為例,說明如何將通過有限元分析獲得的單元L1L2的結構響應(局部坐標系x’o’y’下結構單元的節(jié)點速度L1)和L2(u,))反饋到物體邊界AB的CFD邊界條件。

    首先,根據(jù)結構單元L1L2的形函數(shù)可以獲得L1L2上物體邊界端點A在局部坐標系x’o’y’沿o’x’和o’y’軸的的速度和w˙A:

    式中:N1~N4和N5、N6分別是結構單元彎曲和拉壓的形函數(shù)。對于梁單元來說,L1A/L1L2和=2(L1A-L1L2/2)/L1L2是結構單元L1L2的自然坐標。同理,可以獲得物體邊界另一端點B的速度和B。

    然后,物體邊界AB的平均法向速度vsn和切向速度vst可以計算如下:

    與此同時,物體邊界AB的法向速度梯度▽vsn=(?vsn/?x )i+(?vsn/?y )j可以獲得,

    最后,通過公式(19)和(20)可以獲得物體邊界AB的平均切向速度、法向速度和梯度,然后將其應用到切割單元新時刻的邊界條件(5)式中,用以進行流場計算。同樣,上述方法對于其他兩種情況(圖2(b)和圖2(c))也同樣適用,區(qū)別在于物體邊界的兩個端點所在的結構單元發(fā)生了變化。

    5 數(shù)值算例分析和討論

    5.1 數(shù)值算法驗證

    為了驗證本文數(shù)值方法的準確性和合理性,這里基于孫輝(2003)[5]的彈性楔形體入水實驗進行建模并加以數(shù)值計算,然后將三種彈性楔形體入水模型的數(shù)值結果和實驗數(shù)據(jù)進行比較。基于實驗水池參數(shù)建立0.8m×1.45m方形計算域,水深為1.1m。這里,彈性楔形體的材料屬性如下:密度1.059×103kg/m3,彈性模量2.67×109kg/ms2(),泊松比0.357。此外,三種不同楔形體模型的幾何參數(shù)如表1所示。

    表1 三種不同彈性楔形體模型的幾何參數(shù)Tab.1 Geometric parameters of three different elastic wedges

    圖3 計算域內部的整體網(wǎng)格和物面邊界附近的局部網(wǎng)格Fig.3 Globalmesh(a)and local grid(b)near solid boundary

    圖4 三種不同彈性楔形體的整體運動加速度和局部彈性應變的時間歷程曲線Fig.4 Time history of global acceleration and local strain of observation pointon sloping edge for three elastic wedges

    在數(shù)值計算中,初始時刻(t=0)選為楔形體底部頂點距離初始自由液面0.1m的位置,然后彈性楔形體開始自由下落。物理計算間隔Δt為0.000 1 s,虛擬時間間隔Δτ為0.01 s,人工壓縮系數(shù)β為500,重力加速度為9.81m/s2,結構單元長度為0.01m。此外,在流場計算域中采用非均勻網(wǎng)格進行劃分,在近壁面處為局部細化網(wǎng)格Δx=Δy=0.01m,然后隨著到楔形體邊界距離逐漸增大,網(wǎng)格單元尺度以一定比例逐步增加最終形成66×90的不等間距網(wǎng)格分布,如圖3所示。

    將表1中編號為1~3的三種不同楔形體自由入水模型的整體加速度和斜邊參考點彈性應變的數(shù)值解和實驗數(shù)據(jù)進行比較展示如圖4。從圖中可以看出,本文數(shù)值計算結果和實驗數(shù)據(jù)的變化趨勢基本一致,能夠較好地描述出入水初期整體運動加速度和局部彈性應變的峰值大小和發(fā)生時間。只不過由實驗裝置自身阻力造成的誤差使得整體加速度的數(shù)值解和實驗數(shù)據(jù)略微有些區(qū)別。此外,如圖4(b)所示,隨著楔形體斜升角變小,斜邊中點處的應變在入水過程出現(xiàn)負值,這說明在入水過程中楔形體斜邊先是向內凹陷而后再向外凸起。與實驗數(shù)據(jù)相比,采用本文算法所得負應變略微偏小,這可能由單向耦合計算過程稍微低估了局部結構變形對流場反作用所引起的。綜上所述,從上述三個彈性楔形體入水模型的數(shù)值結果來看,本文算法基本上可以描述入水過程彈性楔形體的整體運動和局部變形的規(guī)律和變化趨勢。因此,可以采用本文算法求解和分析彈性楔形體自由入水問題。

    5.2 彈性楔形體各狀態(tài)參數(shù)對入水運動性能的影響

    這里,為了討論彈性楔形體的狀態(tài)參數(shù)對入水過程中整體運動和局部變形的影響,建立不同參數(shù)(結構材料泊松比、彈性模量、板厚、入水質量、底面斜升角和入水高度等)的彈性楔形結構自由入水模型,并且采用本文算法分別進行數(shù)值計算和分析。

    計算域和模型的基本參數(shù)如下:在2.0m×2.0 m半水半空氣的方形計算域中,長0.6m,寬0.2m,斜升角45°,泊松比為0.357,密度1 059 kg/m3,彈性模量2.67×109kg/ms2(),斜邊板厚1.5 mm,結構重量7.5 kg,從初始0.1m高度自由下落入水。在數(shù)值計算中,流場計算的物理時間步長Δt為0.000 1 s,虛擬時間步長Δτ為0.01 s,人工壓縮系數(shù)β為500,結構單元長度為0.01m,處理位移約束條件所引用的大數(shù)為1013。計算整體動力方程時,Newmark法中的參數(shù)β和δ分別為0.5和0.25。重力加速度為9.81m/s2,流固耦合算法的迭代松弛因子ω為0.8。分別改變結構材料泊松比、彈性模量、板厚、入水質量、底面斜升角和初始入水高度,建立不同狀態(tài)參數(shù)的彈性楔形體入水模型。采用本文數(shù)值算法進行計算和分析,將不同入水模型所對應的整體運動速度和距離楔形體頂點1/4斜邊長度處的局部彈性應變的時間歷程曲線展示如圖5所示。從圖中可以看出,在楔形體入水初期,彈性邊界受到局部入水沖擊載荷的作用會產(chǎn)生極大的應變峰值,并且楔形結構開始整體做減速運動。然后,隨著物體進一步入水,在忽略結構阻尼的前提下,彈性邊界將根據(jù)其固有頻率作高頻和低頻相結合的“拍”振動。進一步,可以發(fā)現(xiàn)彈性楔形體的各狀態(tài)參數(shù)對其入水整體運動和局部變形的影響各不相同。

    其中,改變底面斜升角、入水高度和質量這三個參數(shù)會引起流場發(fā)生強烈變化,從而改變物面所受的局部沖擊載荷和整體水作用力,最終顯著地影響楔形體局部彈性變形和整體運動性能。如圖5(d)所示,如果楔形體模型初始接觸水面時所受水動力相同,那么質量越小的結構,入水速度減小的越快。進而造成彈性邊界所受到的局部沖擊載荷相應地降低,導致楔形結構入水初期的應變峰值變小。從圖5(e)中可以看出,底面斜升角越小,楔形結構在入水過程中所受的局部沖擊載荷和整體水作用力越大,這將導致楔形結構整體減速越快、物面應變峰值和“拍”振動幅值越大。而在圖5(f)中,隨著初始入水高度增加,楔形結構初始接觸水面的時間和速度將會變大,從而入水初期受到的整體水作用力增大,因此楔形結構整體減速越快,與此同時局部沖擊載荷的增大也會導致應變峰值的加大。另一方面,其余的三個參數(shù)(結構板厚、材料泊松比、彈性模量)雖然對楔形結構的整體運動性能作用很小,但是會影響到彈性邊界的局部變形。根據(jù)圖5(a)~(c)中展示,結構板厚對于彈性應變的作用效果比較顯著。板厚越大,彈性斜邊的剛度也就越大,從而造成入水初期的結構應變峰值越小。此外,材料泊松比和彈性模量對結構局部彈性變形也有一定影響。泊松比越小或是彈性模量越小,板條梁模型的等效彈性模量越小,則結構剛度越小,那么楔形邊界的局部應變峰值也就越小。

    圖5 不同彈性楔形體模型的整體運動速度和距離楔形體頂點1/4斜邊長度處的局部應變Fig.5 Global acceleration and local strain of 1/4 length from vertex on sloping edge versus differentmodels

    6 結論

    本文結合液面捕捉法和直角切割網(wǎng)格系統(tǒng)的優(yōu)勢,采用液面捕捉方法和直角切割網(wǎng)格處理物體入水瞬時移動的自由液面和動邊界問題。在此基礎上引入彈性板條梁結構的有限元理論,并且根據(jù)流場的直角切割網(wǎng)格和彈性結構的板條梁單元的特點,提出適合本文算法的流固耦合策略和邊界數(shù)據(jù)傳遞方法。為了分析船舶砰擊問題,建立了類似船首橫剖面的二維彈性楔形體的自由入水模型,并且通過與實驗數(shù)據(jù)進行比較,證明了本文算法求解彈性楔形體入水問題的正確性和合理性。最后,改變彈性楔形結構的狀態(tài)參數(shù)建立不同入水模型,數(shù)值分析這些參數(shù)對彈性楔形結構的整體運動性能和局部彈性變形的影響,可以發(fā)現(xiàn)對楔形結構整體運動特征有比較明顯作用的參數(shù)主要有底面斜升角、入水高度和質量三個參數(shù),并且這些參數(shù)同時會顯著影響入水初期局部彈性響應峰值,這說明船舶在航行過程中砰擊發(fā)生的位置、砰擊入水前的抬升高度以及砰擊處的質量分布對船舶整體運動性能和局部結構強度有很大的影響。此外,楔形結構的斜邊板厚、泊松比和彈性模量雖然不會使整體運動發(fā)生很大變化,但是會影響到結構的局部變形。其中,結構板厚對彈性應變的作用效果最為明顯,而彈性模量和泊松比相對來說要弱一些。

    本文為彈性楔形結構的自由入水問題的求解提供了一種新的CFD算法,并且所得結論能夠為進一步研究船舶砰擊問題、指導船舶航行操縱和結構設計提供一些理論基礎。此外,在以后的研究中可以進一步考慮將流固相互作用的緊耦合策略和結構變形的幾何非線性理論引入算法中,從而擴展本文數(shù)值方法的應用領域和范圍。

    [1]Yamamoto Y,Iida K,Fukasawa T,et al.Structural damage analysis of a fast ship due to bow flare slamming[J].International Shipbuilding Progress,1985,32(369):124-136.

    [2]Gu X K,Moan T.Long-term fatigue damage of ship structures under non-linear wave loads[J].Marine Technology,2002, 39(2):95-104.

    [3]Chuang SL.Investigation of impact of rigid and elastic bodieswith water[R].NavalWarfare and Marine Eng.,1970.

    [4]李國鈞,黃震球.平底物體對水面的斜向沖擊[J].華中理工大學學報,1995,23(1):145-147. LiG J,Huang ZQ.The slant impact of a flat-bottom body on water[J].Journal of Huazhong University of Science and Technology,1995,23(1):145-147.

    [5]孫輝,盧熾華,何友聲.二維楔形體沖擊入水時的流固耦合響應的實驗研究[J].水動力學研究與進展A輯,2003, 1:104-109. Sun H,Lu Z H,He Y S.Experimental research on the fluid-structure interaction in water entry of 2D elastic wedge[J]. Journal of Hydrodynamics,2003,1:104-109.

    [6]Takagi K.Influence of elasticity on hydrodynamic impact problem[J].JKansai Social Navel Architecture(Japan).1994,222: 97-106.

    [7]顧懋祥,程貫一,張效慈.平頭旋轉殼撞水水彈性效應的研究[J].水動力學研究與發(fā)展,1991,6(1):42-51. Gu M X,Cheng G Y,Zhang X C.The study of hydroelastic effect on the water-entry of a revolutionary shell with a flat nose[J].Journal of Hydrodynamics,1991,6(1):42-51.

    [8]盧熾華,何友聲.二維彈性結構入水沖擊過程中的流固耦合效應[J].力學學報,2000,32(2):129-140. Lu CH,He Y S.Coupled analysis of nonlinear interaction between fluid and structure during impact[J].ACTA MECHANICA SINICA,2000,32(2):129-140.

    [9]陳占暉,盧永錦.高速運動物體砰擊水面后的動力學特性研究[J].船舶工程,2009,3:62-65. Chen Z H,Lu Y J.Study on the dynamic characteristics of high-speed object slamming water[J].Ship Engineering,2009, 3:62-65.

    [10]Sun H.A boundary elementmethod applied to strongly nonlinear wave-body interaction problems[D].Norway:Norwegian University of Science and Technology,2007.

    [11]WangW H,Wang Y Y.An improved free surface capturingmethod based on Cartesian cut cellmesh forwater-entry and -exit problems.Proc.R.Soc.A.,2009,465(2106):1843-1868.

    [12]王文華.物體入水問題的分析研究及其在船舶與海洋工程中的應用[D].大連:大連理工大學,2010. Wang W H.Analysis and research on water entry problem and its application in ship and ocean engineering[D].Dalian: Dalian University of Technology,2010.

    [13]Causon D M,Ingram D M,Mingham CG.A Cartesian cut cellmethod for shallow water flowswith moving boundaries[J]. Advanced Water Resources,2001,24(8):899-911.

    [14]Kleefsman K M T,Fekken G,Veldman A E P,et al.A Volume of Fluid based simulation method for wave impact problems[J].Journal of Computational Physics,2005,206:363-393.

    Effect of status parameters for elastic wedge on dynam ic performance of water-entry

    WANGWen-hua1a,b,HUANG Yi1b,c,WANG Yan-ying1b,ZHAIGang-jun1c,HUANG Ya-nan1b,2
    (a.State Key Laboratory of Structural Analysis for Industrial Equipment;b.School of Naval Architecture and Ocean Engineering;c.Deepwater Engineering Research Center,Dalian University of Technology,Dalian 116024,China; 2 Institute of Navigation and Marine Engineering,Dalian Maritime University,Dalian 116024,China)

    In order to analyze the influence of various status parameters of ship slamming on dynamic performance ofwater-entry,a new CFDmethod was presented to numerically simulate the physical process of 2D elastic wedge entering water.Therein surface capturingmethod and Cartesian cut cellmesh system were used to treatmoving free surface and solid boundary.Furthermore,by combiningwith finite element theory of elastic beam structure,the CFD method was developed to study hydroelastic performance of elastic structure during water-entry process.In thismethod,based on the characteristic of Cartesian cut cellmesh and beam element,a particular data transfermethod on elastic boundary is deduced and described.By comparing with the experimental data,the calculated results show the feasibility and validity of this present method to calculate water entry problem of elastic wedge.Finally,free fallingwater-entrymodels of elastic wedge with various structural parameters(material properties,structural thickness and mass,deadrise angle,initial height,and so on)were created,and then the effects of structural parameters on global rigid mo-tion and local elastic response of free-falling elastic wedgewere studied and discussed.

    elastic wedge structure;free fallingwater entry;surface capturingmethod;cartesian cut cellmesh system;finite element theory;data transfermethod

    U671.5

    A

    10.3969/j.issn.1007-7294.2014.11.007

    1007-7294(2014)11-1320-11

    2014-04-26

    國家自然科學創(chuàng)新研究群體基金資助項目“海洋環(huán)境災害作用與結構安全防護”(50921001);國家青年科學基金項目“運動物體穿過深海密度分界面的水動力特性研究”(11202047)

    王文華(1981-),男,博士后,E-mail:wangwenhua0411@yahoo.cn;黃一(1964-),男,大連理工大學教授,博士生導師。

    猜你喜歡
    結構單元楔形流場
    大型空冷汽輪發(fā)電機轉子三維流場計算
    大電機技術(2021年2期)2021-07-21 07:28:24
    History of the Alphabet
    鋼絲繩楔形接頭連接失效分析與預防
    Eight Surprising Foods You’er Never Tried to Grill Before
    轉杯紡排雜區(qū)流場與排雜性能
    基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計分析
    腹腔鏡下胃楔形切除術治療胃間質瘤30例
    基于瞬態(tài)流場計算的滑動軸承靜平衡位置求解
    一種具有表面活性功能的聚合物及其制備方法和應用
    石油化工(2015年9期)2015-08-15 00:43:05
    基于ANSYS的某型航空發(fā)動機軸承試驗器支承剛度研究
    亚洲av不卡在线观看| 成人午夜高清在线视频| av在线蜜桃| 麻豆成人av在线观看| 日韩精品中文字幕看吧| 欧洲精品卡2卡3卡4卡5卡区| 特大巨黑吊av在线直播| 成人亚洲精品av一区二区| 麻豆成人午夜福利视频| 淫秽高清视频在线观看| 国产主播在线观看一区二区| 亚洲无线观看免费| 午夜福利成人在线免费观看| 精品99又大又爽又粗少妇毛片 | 天堂影院成人在线观看| 国产单亲对白刺激| 97超级碰碰碰精品色视频在线观看| 久久99热6这里只有精品| 夜夜夜夜夜久久久久| 国产麻豆成人av免费视频| 无人区码免费观看不卡| 在线免费观看不下载黄p国产 | 老熟妇仑乱视频hdxx| av在线蜜桃| 日韩欧美国产一区二区入口| 亚洲va日本ⅴa欧美va伊人久久| 老熟妇乱子伦视频在线观看| 国产老妇女一区| 成人av一区二区三区在线看| 男女啪啪激烈高潮av片| 国产精品久久久久久亚洲av鲁大| 精品人妻视频免费看| 亚洲第一区二区三区不卡| 国内精品美女久久久久久| 一区二区三区激情视频| 国产高清视频在线播放一区| 此物有八面人人有两片| 国产av在哪里看| 国产一区二区三区av在线 | 午夜久久久久精精品| 真实男女啪啪啪动态图| 中文字幕人妻熟人妻熟丝袜美| 亚洲四区av| 一本一本综合久久| 非洲黑人性xxxx精品又粗又长| 成人鲁丝片一二三区免费| 91麻豆av在线| 久久精品国产亚洲av天美| 老师上课跳d突然被开到最大视频| 日本三级黄在线观看| 内地一区二区视频在线| 精品久久久久久成人av| 成熟少妇高潮喷水视频| 在线播放国产精品三级| 日韩欧美国产在线观看| 大又大粗又爽又黄少妇毛片口| 成人精品一区二区免费| 日韩国内少妇激情av| 欧美最新免费一区二区三区| 国产精品98久久久久久宅男小说| 欧美日韩综合久久久久久 | 久久精品国产亚洲av香蕉五月| 深夜精品福利| 日韩强制内射视频| 少妇裸体淫交视频免费看高清| 亚洲七黄色美女视频| 一边摸一边抽搐一进一小说| 久9热在线精品视频| 非洲黑人性xxxx精品又粗又长| 人妻久久中文字幕网| 国产色爽女视频免费观看| 国产精品电影一区二区三区| 亚洲精华国产精华液的使用体验 | 国产亚洲av嫩草精品影院| 亚洲欧美激情综合另类| 国产综合懂色| 国产在线男女| 欧美最新免费一区二区三区| 久久精品人妻少妇| 99久久精品热视频| 九九热线精品视视频播放| 一个人观看的视频www高清免费观看| 日韩欧美精品免费久久| 亚洲精华国产精华精| 99久久精品国产国产毛片| 成人性生交大片免费视频hd| 免费看av在线观看网站| 极品教师在线视频| 久久国产精品人妻蜜桃| 精品久久久久久久人妻蜜臀av| 特级一级黄色大片| av.在线天堂| 一个人看的www免费观看视频| 午夜日韩欧美国产| 国产精品99久久久久久久久| 小说图片视频综合网站| 欧美一区二区精品小视频在线| 97人妻精品一区二区三区麻豆| 日韩欧美三级三区| 国产精品爽爽va在线观看网站| 久久午夜福利片| 日韩 亚洲 欧美在线| aaaaa片日本免费| 亚洲第一电影网av| 亚洲三级黄色毛片| 国产黄片美女视频| 国内少妇人妻偷人精品xxx网站| 深夜a级毛片| 深夜精品福利| 天美传媒精品一区二区| 男插女下体视频免费在线播放| 麻豆精品久久久久久蜜桃| 九九爱精品视频在线观看| 亚洲三级黄色毛片| 黄色一级大片看看| 午夜福利在线观看免费完整高清在 | 国产精品1区2区在线观看.| 麻豆久久精品国产亚洲av| 无遮挡黄片免费观看| 久久这里只有精品中国| 国产一区二区激情短视频| 一级毛片久久久久久久久女| 国模一区二区三区四区视频| 午夜精品久久久久久毛片777| 国产精品久久电影中文字幕| 啦啦啦观看免费观看视频高清| 一个人看视频在线观看www免费| 国产精品亚洲一级av第二区| 午夜免费成人在线视频| 国内揄拍国产精品人妻在线| 国产成年人精品一区二区| 日本爱情动作片www.在线观看 | 禁无遮挡网站| 如何舔出高潮| 免费观看精品视频网站| 不卡视频在线观看欧美| 麻豆久久精品国产亚洲av| 国产精品无大码| 波多野结衣高清作品| 一区二区三区四区激情视频 | av在线观看视频网站免费| 久久久成人免费电影| 亚洲狠狠婷婷综合久久图片| 国产精品综合久久久久久久免费| 欧美日韩黄片免| av中文乱码字幕在线| 免费av毛片视频| 校园人妻丝袜中文字幕| 欧美色视频一区免费| 蜜桃亚洲精品一区二区三区| 午夜精品一区二区三区免费看| 搞女人的毛片| 成人av一区二区三区在线看| 亚洲美女视频黄频| 91久久精品国产一区二区三区| 亚洲精品影视一区二区三区av| 午夜激情欧美在线| 欧洲精品卡2卡3卡4卡5卡区| 亚洲国产日韩欧美精品在线观看| 亚洲色图av天堂| 日韩精品青青久久久久久| 99热网站在线观看| 两个人视频免费观看高清| 午夜爱爱视频在线播放| 国产综合懂色| 伊人久久精品亚洲午夜| 一进一出抽搐gif免费好疼| 黄色女人牲交| 一级毛片久久久久久久久女| 麻豆成人午夜福利视频| 久久精品国产亚洲av天美| 免费av观看视频| 亚洲中文字幕一区二区三区有码在线看| 亚洲欧美日韩无卡精品| 长腿黑丝高跟| 亚洲va在线va天堂va国产| 女人被狂操c到高潮| av在线老鸭窝| 白带黄色成豆腐渣| 亚洲精品一区av在线观看| 深夜精品福利| 日本 欧美在线| 亚洲精品456在线播放app | 国产精品一区二区三区四区久久| 观看美女的网站| 男人和女人高潮做爰伦理| 日本在线视频免费播放| 国产精品乱码一区二三区的特点| 欧美一区二区精品小视频在线| 美女大奶头视频| 亚洲av美国av| 一个人观看的视频www高清免费观看| 国产亚洲欧美98| 啦啦啦观看免费观看视频高清| 看十八女毛片水多多多| 少妇人妻一区二区三区视频| 九九热线精品视视频播放| 日韩亚洲欧美综合| 国内久久婷婷六月综合欲色啪| 韩国av一区二区三区四区| 久久久久久九九精品二区国产| 久久中文看片网| 听说在线观看完整版免费高清| 色综合站精品国产| 最后的刺客免费高清国语| 美女高潮的动态| 中文字幕高清在线视频| 国产精品一区二区性色av| 久久婷婷人人爽人人干人人爱| 黄色视频,在线免费观看| 久久国内精品自在自线图片| 亚洲av电影不卡..在线观看| 搡老熟女国产l中国老女人| 嫩草影视91久久| 久久久午夜欧美精品| 床上黄色一级片| 给我免费播放毛片高清在线观看| 国产av一区在线观看免费| 国产免费男女视频| 亚洲成人免费电影在线观看| www日本黄色视频网| 乱人视频在线观看| 夜夜夜夜夜久久久久| 国产欧美日韩精品一区二区| 精品一区二区三区人妻视频| 他把我摸到了高潮在线观看| 国产精品自产拍在线观看55亚洲| 真人一进一出gif抽搐免费| 直男gayav资源| 色综合色国产| av国产免费在线观看| 我要看日韩黄色一级片| 午夜激情福利司机影院| 午夜视频国产福利| 一区二区三区高清视频在线| 十八禁网站免费在线| 禁无遮挡网站| 亚洲av成人精品一区久久| 亚洲成人免费电影在线观看| 久久这里只有精品中国| 女人十人毛片免费观看3o分钟| 国产精品久久久久久久久免| 中文字幕高清在线视频| 国内少妇人妻偷人精品xxx网站| 精品久久久久久久久av| 亚洲欧美清纯卡通| av黄色大香蕉| 精品不卡国产一区二区三区| 日本 av在线| 69av精品久久久久久| 特大巨黑吊av在线直播| 天天一区二区日本电影三级| 色在线成人网| 国产精品三级大全| 人妻少妇偷人精品九色| 亚洲人成网站高清观看| 动漫黄色视频在线观看| 88av欧美| 国产亚洲av嫩草精品影院| 亚洲在线观看片| 精品一区二区三区人妻视频| 欧美精品国产亚洲| 嫁个100分男人电影在线观看| 国内久久婷婷六月综合欲色啪| 国产精品一区二区三区四区免费观看 | 夜夜爽天天搞| 成人鲁丝片一二三区免费| 欧美人与善性xxx| 久久国产乱子免费精品| 国产精品日韩av在线免费观看| 中国美女看黄片| 美女高潮的动态| 久久久久精品国产欧美久久久| 国产成人一区二区在线| 午夜福利在线在线| 十八禁网站免费在线| 中亚洲国语对白在线视频| 黄色丝袜av网址大全| 变态另类成人亚洲欧美熟女| www.色视频.com| 久久久久国内视频| 九色国产91popny在线| av.在线天堂| 久久久成人免费电影| 夜夜爽天天搞| 亚洲天堂国产精品一区在线| 18+在线观看网站| 欧美不卡视频在线免费观看| 精品国产三级普通话版| 日本 欧美在线| 观看免费一级毛片| 国产极品精品免费视频能看的| 香蕉av资源在线| 国产精品久久视频播放| 日韩欧美国产在线观看| 国产极品精品免费视频能看的| 天堂动漫精品| 麻豆国产av国片精品| 亚洲av中文字字幕乱码综合| 波多野结衣高清无吗| 亚洲性夜色夜夜综合| 美女高潮的动态| 一个人看的www免费观看视频| 久久久久久伊人网av| 国产淫片久久久久久久久| 日日摸夜夜添夜夜添小说| 欧美又色又爽又黄视频| 中文字幕熟女人妻在线| 老司机午夜福利在线观看视频| 伦精品一区二区三区| 美女免费视频网站| 亚洲精品影视一区二区三区av| 无人区码免费观看不卡| 国产伦精品一区二区三区四那| 日本与韩国留学比较| 中文字幕熟女人妻在线| 老司机福利观看| 亚洲电影在线观看av| 亚洲内射少妇av| 亚洲成a人片在线一区二区| 黄色一级大片看看| 色综合亚洲欧美另类图片| 成年人黄色毛片网站| 亚洲精华国产精华液的使用体验 | 超碰av人人做人人爽久久| 精品99又大又爽又粗少妇毛片 | 亚洲五月天丁香| 日本在线视频免费播放| av.在线天堂| 亚洲av日韩精品久久久久久密| 99久国产av精品| av.在线天堂| 淫秽高清视频在线观看| 婷婷丁香在线五月| 中文亚洲av片在线观看爽| 亚洲无线在线观看| 精品日产1卡2卡| 国产精品亚洲一级av第二区| 男女边吃奶边做爰视频| 最新在线观看一区二区三区| 国产免费av片在线观看野外av| 国产三级在线视频| 日韩av在线大香蕉| 亚洲不卡免费看| 1000部很黄的大片| 综合色av麻豆| 国产69精品久久久久777片| 成人国产综合亚洲| 我的女老师完整版在线观看| 国产精品嫩草影院av在线观看 | 日韩一区二区视频免费看| 亚洲美女黄片视频| 两个人视频免费观看高清| 午夜福利成人在线免费观看| 婷婷精品国产亚洲av| 亚洲18禁久久av| 春色校园在线视频观看| 亚洲在线自拍视频| 日韩欧美 国产精品| 亚洲国产欧洲综合997久久,| 乱系列少妇在线播放| 日韩欧美在线乱码| 日韩中字成人| 国内精品宾馆在线| 国产激情偷乱视频一区二区| www日本黄色视频网| 高清毛片免费观看视频网站| 精品久久久久久久人妻蜜臀av| 国产不卡一卡二| 一进一出抽搐动态| 在线观看免费视频日本深夜| 99在线人妻在线中文字幕| 最新在线观看一区二区三区| 又紧又爽又黄一区二区| 亚洲无线在线观看| 午夜爱爱视频在线播放| 伦理电影大哥的女人| 变态另类丝袜制服| 男人狂女人下面高潮的视频| 精品不卡国产一区二区三区| 亚洲性夜色夜夜综合| 级片在线观看| 国产三级在线视频| 一本久久中文字幕| 欧美日韩国产亚洲二区| 久久精品国产清高在天天线| 韩国av在线不卡| 国产熟女欧美一区二区| 噜噜噜噜噜久久久久久91| 国产69精品久久久久777片| 国产熟女欧美一区二区| 亚洲av一区综合| 国产欧美日韩精品一区二区| 精品国产三级普通话版| 乱系列少妇在线播放| 成年版毛片免费区| 国产精品三级大全| 国产单亲对白刺激| 日韩中文字幕欧美一区二区| 国产精品自产拍在线观看55亚洲| 蜜桃亚洲精品一区二区三区| 精品一区二区三区人妻视频| 日韩强制内射视频| 欧美国产日韩亚洲一区| 欧美一区二区精品小视频在线| 搡老妇女老女人老熟妇| 日本在线视频免费播放| 国产男靠女视频免费网站| 日韩亚洲欧美综合| 老熟妇乱子伦视频在线观看| 久久中文看片网| 久久久精品欧美日韩精品| 国内久久婷婷六月综合欲色啪| 国产精品福利在线免费观看| 最后的刺客免费高清国语| 日韩国内少妇激情av| 日本精品一区二区三区蜜桃| 校园春色视频在线观看| 一区二区三区高清视频在线| 国产精品久久久久久av不卡| 伦精品一区二区三区| 亚洲内射少妇av| 久久精品国产亚洲av天美| 亚洲乱码一区二区免费版| 亚洲国产色片| 国产乱人伦免费视频| 乱人视频在线观看| 精品久久久久久久久久免费视频| 国产精品福利在线免费观看| 成人午夜高清在线视频| 免费观看精品视频网站| 九九久久精品国产亚洲av麻豆| 干丝袜人妻中文字幕| 中文字幕av成人在线电影| 俺也久久电影网| 搞女人的毛片| 国国产精品蜜臀av免费| 中文在线观看免费www的网站| 中国美女看黄片| 成人综合一区亚洲| 久久精品夜夜夜夜夜久久蜜豆| 成人av在线播放网站| 少妇的逼水好多| 村上凉子中文字幕在线| 久久天躁狠狠躁夜夜2o2o| 一个人看视频在线观看www免费| 久久久久久国产a免费观看| 老师上课跳d突然被开到最大视频| 麻豆成人av在线观看| 日本a在线网址| 赤兔流量卡办理| 99久久中文字幕三级久久日本| 精品午夜福利视频在线观看一区| 神马国产精品三级电影在线观看| 精品欧美国产一区二区三| 国产精品1区2区在线观看.| 日本爱情动作片www.在线观看 | av在线亚洲专区| 色哟哟·www| 啦啦啦观看免费观看视频高清| 高清日韩中文字幕在线| 免费看日本二区| 亚洲精品在线观看二区| 国内精品宾馆在线| 欧美区成人在线视频| 亚洲av熟女| 久久久久免费精品人妻一区二区| 国产一区二区亚洲精品在线观看| 老女人水多毛片| 久久精品国产亚洲av香蕉五月| 久久草成人影院| 亚洲美女黄片视频| 中文字幕精品亚洲无线码一区| 久久热精品热| 国产亚洲精品久久久com| 婷婷精品国产亚洲av| 色播亚洲综合网| 欧美精品啪啪一区二区三区| 国产精品女同一区二区软件 | 国产伦一二天堂av在线观看| 女生性感内裤真人,穿戴方法视频| 天堂网av新在线| 美女高潮喷水抽搐中文字幕| av在线老鸭窝| 香蕉av资源在线| 日韩欧美国产在线观看| 亚洲综合色惰| 欧美国产日韩亚洲一区| 欧美bdsm另类| 成人亚洲精品av一区二区| 亚洲精品日韩av片在线观看| 干丝袜人妻中文字幕| 日韩欧美精品v在线| 桃色一区二区三区在线观看| 一进一出抽搐gif免费好疼| 久久九九热精品免费| 动漫黄色视频在线观看| 高清毛片免费观看视频网站| 女同久久另类99精品国产91| 中亚洲国语对白在线视频| 国产精品98久久久久久宅男小说| 亚洲五月天丁香| 网址你懂的国产日韩在线| 69av精品久久久久久| 国产精品国产高清国产av| 国产伦人伦偷精品视频| 久久久久久久久久成人| 日本爱情动作片www.在线观看 | 精品不卡国产一区二区三区| 夜夜爽天天搞| 国产黄片美女视频| 精品久久久久久久末码| h日本视频在线播放| 天天躁日日操中文字幕| 白带黄色成豆腐渣| 一区二区三区四区激情视频 | 在线免费观看不下载黄p国产 | 一级黄片播放器| 三级国产精品欧美在线观看| 精品一区二区三区视频在线| 成人av一区二区三区在线看| 精品人妻熟女av久视频| 嫩草影视91久久| 狂野欧美白嫩少妇大欣赏| 婷婷精品国产亚洲av在线| 久久精品国产亚洲网站| 老熟妇仑乱视频hdxx| 51国产日韩欧美| 精品人妻1区二区| 非洲黑人性xxxx精品又粗又长| 日日啪夜夜撸| 很黄的视频免费| 亚洲乱码一区二区免费版| 久久精品国产99精品国产亚洲性色| 最好的美女福利视频网| 中出人妻视频一区二区| 婷婷色综合大香蕉| 亚洲不卡免费看| 国产蜜桃级精品一区二区三区| 在线播放国产精品三级| 免费不卡的大黄色大毛片视频在线观看 | 精品人妻视频免费看| 女的被弄到高潮叫床怎么办 | 性欧美人与动物交配| 国产精品福利在线免费观看| 欧美中文日本在线观看视频| 最好的美女福利视频网| 91精品国产九色| www日本黄色视频网| 亚洲avbb在线观看| 色综合亚洲欧美另类图片| 中文在线观看免费www的网站| 亚洲18禁久久av| 内地一区二区视频在线| a在线观看视频网站| 国产熟女欧美一区二区| 禁无遮挡网站| 黄色日韩在线| 午夜老司机福利剧场| 欧美人与善性xxx| 丰满乱子伦码专区| 中文字幕av在线有码专区| 超碰av人人做人人爽久久| 两性午夜刺激爽爽歪歪视频在线观看| 狂野欧美激情性xxxx在线观看| 老师上课跳d突然被开到最大视频| 深夜精品福利| 欧美最黄视频在线播放免费| 成人特级av手机在线观看| 国产精品人妻久久久久久| 18禁在线播放成人免费| 午夜免费男女啪啪视频观看 | 国内毛片毛片毛片毛片毛片| 亚洲国产精品合色在线| 久久6这里有精品| 婷婷精品国产亚洲av| 国产亚洲精品久久久com| 亚洲中文字幕日韩| 国产精品嫩草影院av在线观看 | 国产精品自产拍在线观看55亚洲| 欧美+日韩+精品| а√天堂www在线а√下载| 国产乱人伦免费视频| 美女cb高潮喷水在线观看| 91久久精品电影网| 韩国av在线不卡| 亚洲自拍偷在线| 极品教师在线视频| 亚洲av.av天堂| 大又大粗又爽又黄少妇毛片口| 欧美日韩中文字幕国产精品一区二区三区| 亚洲中文日韩欧美视频| 深爱激情五月婷婷| 精品福利观看| 好男人在线观看高清免费视频| 亚洲精品亚洲一区二区| 久久久国产成人精品二区| 99久久精品热视频| 国内精品久久久久精免费| 亚洲精品国产成人久久av| 日韩,欧美,国产一区二区三区 | 九九爱精品视频在线观看| 免费不卡的大黄色大毛片视频在线观看 | 一进一出好大好爽视频| 免费人成在线观看视频色| 久久久久久大精品| 搞女人的毛片| 免费人成视频x8x8入口观看| 久久香蕉精品热| 国产精品一区www在线观看 | 桃色一区二区三区在线观看| 国内揄拍国产精品人妻在线| 日本一本二区三区精品| 91久久精品国产一区二区三区| 国产探花极品一区二区|