• <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ā)動機軸承試驗器支承剛度研究
    欧美97在线视频| 国产成人精品无人区| 亚洲精品,欧美精品| 国产亚洲欧美精品永久| 亚洲av免费高清在线观看| 丝瓜视频免费看黄片| 亚洲,欧美精品.| 国产福利在线免费观看视频| 大片电影免费在线观看免费| 午夜激情av网站| 亚洲欧美色中文字幕在线| 成人免费观看视频高清| 日韩伦理黄色片| 看免费av毛片| 黄网站色视频无遮挡免费观看| 美女中出高潮动态图| 精品国产一区二区三区久久久樱花| 26uuu在线亚洲综合色| 国产精品女同一区二区软件| 亚洲三区欧美一区| 亚洲成av片中文字幕在线观看 | 日日撸夜夜添| 高清不卡的av网站| 飞空精品影院首页| 成人免费观看视频高清| 一级毛片 在线播放| 美女福利国产在线| 九九爱精品视频在线观看| 亚洲久久久国产精品| 亚洲精品久久久久久婷婷小说| 亚洲内射少妇av| 亚洲av免费高清在线观看| tube8黄色片| 黑丝袜美女国产一区| 黄片无遮挡物在线观看| 人人妻人人添人人爽欧美一区卜| 青春草国产在线视频| 亚洲人成网站在线观看播放| 免费高清在线观看视频在线观看| www.熟女人妻精品国产| 国产精品国产三级国产专区5o| 久久精品aⅴ一区二区三区四区 | 午夜91福利影院| 亚洲美女搞黄在线观看| 国产熟女欧美一区二区| 中文精品一卡2卡3卡4更新| 婷婷色综合www| 日日啪夜夜爽| 久久国产精品大桥未久av| 最新中文字幕久久久久| 人妻人人澡人人爽人人| 欧美97在线视频| 男女下面插进去视频免费观看| 一本色道久久久久久精品综合| 国产免费一区二区三区四区乱码| 国产一区有黄有色的免费视频| 叶爱在线成人免费视频播放| 99热网站在线观看| 老司机影院成人| 国产又爽黄色视频| 亚洲精品自拍成人| 欧美另类一区| 一区二区三区四区激情视频| 黄片无遮挡物在线观看| 午夜免费男女啪啪视频观看| 亚洲成av片中文字幕在线观看 | 男的添女的下面高潮视频| 91午夜精品亚洲一区二区三区| 26uuu在线亚洲综合色| 免费久久久久久久精品成人欧美视频| 91久久精品国产一区二区三区| 免费av中文字幕在线| 国产精品人妻久久久影院| 成年女人毛片免费观看观看9 | 国产高清国产精品国产三级| 亚洲婷婷狠狠爱综合网| av线在线观看网站| 精品国产露脸久久av麻豆| 99热网站在线观看| 午夜老司机福利剧场| 欧美av亚洲av综合av国产av | 大片电影免费在线观看免费| 性色av一级| 欧美+日韩+精品| av福利片在线| 嫩草影院入口| 免费黄频网站在线观看国产| 久久免费观看电影| 久久99热这里只频精品6学生| 飞空精品影院首页| 90打野战视频偷拍视频| 啦啦啦在线观看免费高清www| a级毛片黄视频| 日韩一区二区视频免费看| 老鸭窝网址在线观看| 国产精品一区二区在线不卡| 天堂俺去俺来也www色官网| 亚洲,欧美精品.| 国产熟女午夜一区二区三区| av电影中文网址| 亚洲国产毛片av蜜桃av| freevideosex欧美| 精品卡一卡二卡四卡免费| 精品亚洲成a人片在线观看| 永久网站在线| 色哟哟·www| 侵犯人妻中文字幕一二三四区| 18禁国产床啪视频网站| 精品国产一区二区久久| av女优亚洲男人天堂| 人妻 亚洲 视频| 婷婷色麻豆天堂久久| 黑人欧美特级aaaaaa片| 女的被弄到高潮叫床怎么办| 午夜福利在线观看免费完整高清在| 中文精品一卡2卡3卡4更新| 久久精品夜色国产| 国产毛片在线视频| 视频区图区小说| 国产免费一区二区三区四区乱码| 欧美变态另类bdsm刘玥| 日韩精品有码人妻一区| 亚洲一区中文字幕在线| 黄色一级大片看看| 男女边吃奶边做爰视频| 在线亚洲精品国产二区图片欧美| 制服丝袜香蕉在线| 午夜免费鲁丝| 飞空精品影院首页| 亚洲av日韩在线播放| 又粗又硬又长又爽又黄的视频| 少妇人妻精品综合一区二区| 午夜影院在线不卡| 永久网站在线| 欧美日韩亚洲国产一区二区在线观看 | 国产色婷婷99| 亚洲精品日韩在线中文字幕| 国产成人免费无遮挡视频| 丝袜喷水一区| 久久久亚洲精品成人影院| 国产日韩一区二区三区精品不卡| 欧美成人午夜免费资源| 久久精品国产亚洲av高清一级| 亚洲欧洲日产国产| 亚洲精品国产av蜜桃| 国产精品av久久久久免费| 国产成人精品婷婷| 国产成人欧美| 狠狠婷婷综合久久久久久88av| www.精华液| 另类精品久久| 黑丝袜美女国产一区| 日韩av免费高清视频| 青春草国产在线视频| 亚洲经典国产精华液单| 中文精品一卡2卡3卡4更新| 麻豆精品久久久久久蜜桃| 久久久久国产网址| 午夜激情久久久久久久| 18禁动态无遮挡网站| 狠狠婷婷综合久久久久久88av| www.熟女人妻精品国产| 丁香六月天网| 男女啪啪激烈高潮av片| xxxhd国产人妻xxx| 国产xxxxx性猛交| 18在线观看网站| 亚洲 欧美一区二区三区| 日日啪夜夜爽| 视频在线观看一区二区三区| 亚洲婷婷狠狠爱综合网| 欧美老熟妇乱子伦牲交| 日韩大片免费观看网站| 晚上一个人看的免费电影| 国产av码专区亚洲av| 亚洲欧洲日产国产| 涩涩av久久男人的天堂| 日本色播在线视频| 下体分泌物呈黄色| 成人免费观看视频高清| 国产日韩一区二区三区精品不卡| 欧美精品高潮呻吟av久久| 亚洲av国产av综合av卡| 性色avwww在线观看| 久久久精品免费免费高清| 日本色播在线视频| 最近的中文字幕免费完整| 成人免费观看视频高清| 在线观看国产h片| www.自偷自拍.com| 街头女战士在线观看网站| 精品国产国语对白av| 亚洲人成网站在线观看播放| 各种免费的搞黄视频| 欧美成人午夜精品| 亚洲av日韩在线播放| 日本av手机在线免费观看| 国产野战对白在线观看| 国产精品一国产av| 夜夜骑夜夜射夜夜干| 国产毛片在线视频| 国产精品三级大全| 亚洲国产精品国产精品| 两性夫妻黄色片| 美国免费a级毛片| 久久毛片免费看一区二区三区| 免费在线观看黄色视频的| 日本欧美国产在线视频| 性高湖久久久久久久久免费观看| 亚洲精品国产色婷婷电影| 丝瓜视频免费看黄片| 亚洲av国产av综合av卡| 免费日韩欧美在线观看| 色94色欧美一区二区| 欧美国产精品一级二级三级| 精品亚洲成国产av| 中文字幕人妻熟女乱码| 国产精品免费视频内射| 夫妻性生交免费视频一级片| 亚洲 欧美一区二区三区| 美女中出高潮动态图| 久久精品久久久久久噜噜老黄| 曰老女人黄片| www.av在线官网国产| av网站在线播放免费| 欧美亚洲日本最大视频资源| 日韩人妻精品一区2区三区| 卡戴珊不雅视频在线播放| 久久久国产欧美日韩av| 国产极品粉嫩免费观看在线| 三级国产精品片| 精品国产露脸久久av麻豆| 人人澡人人妻人| 少妇精品久久久久久久| 国产麻豆69| 亚洲五月色婷婷综合| 丝袜美腿诱惑在线| 中文字幕av电影在线播放| 考比视频在线观看| 日日撸夜夜添| 国产成人欧美| 亚洲一区二区三区欧美精品| 下体分泌物呈黄色| 制服丝袜香蕉在线| av电影中文网址| 精品人妻熟女毛片av久久网站| 国产欧美日韩综合在线一区二区| 久久午夜综合久久蜜桃| 建设人人有责人人尽责人人享有的| 免费少妇av软件| 免费女性裸体啪啪无遮挡网站| 免费不卡的大黄色大毛片视频在线观看| 亚洲欧美一区二区三区久久| 欧美成人午夜精品| 欧美激情 高清一区二区三区| 女人精品久久久久毛片| 欧美成人午夜免费资源| 中文字幕人妻熟女乱码| 久久久久久免费高清国产稀缺| 男女下面插进去视频免费观看| 精品少妇内射三级| 亚洲精品美女久久av网站| 97在线视频观看| 极品少妇高潮喷水抽搐| 欧美bdsm另类| 欧美成人午夜精品| 国产成人精品久久二区二区91 | 国产精品久久久久成人av| 久久久久精品人妻al黑| 一级毛片电影观看| 国产日韩欧美亚洲二区| 视频区图区小说| 久久精品久久久久久久性| 亚洲精品aⅴ在线观看| 少妇猛男粗大的猛烈进出视频| 久久毛片免费看一区二区三区| 久久久久视频综合| 十八禁高潮呻吟视频| 91aial.com中文字幕在线观看| 在线 av 中文字幕| 狠狠婷婷综合久久久久久88av| 久久韩国三级中文字幕| 国产一区二区激情短视频 | 国产高清不卡午夜福利| 最近中文字幕高清免费大全6| 久久这里有精品视频免费| 性少妇av在线| 欧美日韩亚洲高清精品| 91国产中文字幕| 九色亚洲精品在线播放| 纵有疾风起免费观看全集完整版| 超碰成人久久| 欧美日韩综合久久久久久| 黄色怎么调成土黄色| 秋霞伦理黄片| 国产成人精品婷婷| 精品一区二区三卡| 国产精品人妻久久久影院| 亚洲伊人色综图| 婷婷色综合大香蕉| 亚洲精品中文字幕在线视频| 只有这里有精品99| 一区二区三区激情视频| av在线观看视频网站免费| 美女国产高潮福利片在线看| 国产成人精品无人区| 日韩免费高清中文字幕av| 女性生殖器流出的白浆| 国产1区2区3区精品| 国产成人一区二区在线| 亚洲国产精品一区三区| 午夜福利,免费看| 男女国产视频网站| 两个人免费观看高清视频| 大码成人一级视频| 欧美av亚洲av综合av国产av | 亚洲欧美清纯卡通| 国产精品av久久久久免费| 国产1区2区3区精品| 精品人妻熟女毛片av久久网站| 久久影院123| 国产av一区二区精品久久| 成年女人毛片免费观看观看9 | 欧美国产精品va在线观看不卡| 国产精品一区二区在线观看99| 午夜免费鲁丝| 亚洲成av片中文字幕在线观看 | 精品视频人人做人人爽| www.精华液| 男人添女人高潮全过程视频| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲伊人色综图| 国产国语露脸激情在线看| 婷婷色av中文字幕| 观看av在线不卡| 建设人人有责人人尽责人人享有的| 麻豆av在线久日| 国产成人aa在线观看| 国产女主播在线喷水免费视频网站| 人妻人人澡人人爽人人| 大陆偷拍与自拍| 好男人视频免费观看在线| videosex国产| 卡戴珊不雅视频在线播放| www.av在线官网国产| 精品少妇久久久久久888优播| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 伊人久久大香线蕉亚洲五| 欧美日韩一级在线毛片| 欧美成人午夜免费资源| 亚洲av男天堂| 亚洲av在线观看美女高潮| 欧美bdsm另类| 久久久久久人妻| 观看av在线不卡| 欧美成人精品欧美一级黄| 免费观看a级毛片全部| 激情视频va一区二区三区| 只有这里有精品99| 国产精品久久久久久久久免| 美女主播在线视频| 国产深夜福利视频在线观看| 久久精品国产自在天天线| 国产一级毛片在线| 视频区图区小说| 天天躁夜夜躁狠狠躁躁| av片东京热男人的天堂| 一个人免费看片子| 熟女少妇亚洲综合色aaa.| 国产又爽黄色视频| 97在线视频观看| 午夜福利一区二区在线看| 亚洲内射少妇av| 久久精品aⅴ一区二区三区四区 | 成年动漫av网址| 中文精品一卡2卡3卡4更新| 欧美97在线视频| 久久久国产一区二区| 国产xxxxx性猛交| 2018国产大陆天天弄谢| 午夜免费观看性视频| 午夜久久久在线观看| 国产精品香港三级国产av潘金莲 | 一本—道久久a久久精品蜜桃钙片| 十八禁高潮呻吟视频| 色吧在线观看| 一级爰片在线观看| 中文字幕色久视频| 久久久精品免费免费高清| 国产片内射在线| 国产精品不卡视频一区二区| 久久精品久久精品一区二区三区| 国产av一区二区精品久久| 午夜激情久久久久久久| 国产精品国产av在线观看| 精品一区在线观看国产| 99国产精品免费福利视频| 色视频在线一区二区三区| 日韩,欧美,国产一区二区三区| 免费观看无遮挡的男女| 精品国产一区二区三区久久久樱花| 亚洲精品国产av成人精品| 欧美xxⅹ黑人| 最近2019中文字幕mv第一页| 国产极品天堂在线| 黄色 视频免费看| 最新中文字幕久久久久| 日本91视频免费播放| 在线观看三级黄色| 免费看不卡的av| 久久久久久久久免费视频了| 性高湖久久久久久久久免费观看| 亚洲综合精品二区| 久久久久久久久久久久大奶| 国产精品国产三级国产专区5o| 香蕉国产在线看| 久久久久久久国产电影| 色婷婷av一区二区三区视频| 狠狠精品人妻久久久久久综合| 又大又黄又爽视频免费| av在线app专区| 一本—道久久a久久精品蜜桃钙片| 美女午夜性视频免费| 黑人猛操日本美女一级片| 亚洲,欧美精品.| 久久99热这里只频精品6学生| 久久久国产精品麻豆| 国语对白做爰xxxⅹ性视频网站| 久久精品国产亚洲av高清一级| 午夜激情av网站| 久久99蜜桃精品久久| 亚洲国产欧美日韩在线播放| 亚洲av男天堂| 亚洲精品久久成人aⅴ小说| 亚洲欧洲精品一区二区精品久久久 | 久久久国产欧美日韩av| 男女边吃奶边做爰视频| 亚洲在久久综合| 日日摸夜夜添夜夜爱| 午夜福利在线免费观看网站| 亚洲一区中文字幕在线| 黄片无遮挡物在线观看| 久久精品国产鲁丝片午夜精品| 久久久久国产一级毛片高清牌| 久久精品熟女亚洲av麻豆精品| 在线观看国产h片| 国产欧美亚洲国产| 欧美人与性动交α欧美软件| 亚洲精品日韩在线中文字幕| 午夜精品国产一区二区电影| 99久久中文字幕三级久久日本| 久久婷婷青草| 人人妻人人添人人爽欧美一区卜| 欧美日韩综合久久久久久| 美女国产高潮福利片在线看| 香蕉国产在线看| 亚洲伊人色综图| 国产亚洲欧美精品永久| 欧美人与善性xxx| 99久久综合免费| 在线看a的网站| 在线免费观看不下载黄p国产| 91aial.com中文字幕在线观看| 久久精品国产亚洲av天美| 日韩成人av中文字幕在线观看| 久久久精品免费免费高清| 日韩成人av中文字幕在线观看| 亚洲人成77777在线视频| 亚洲四区av| 亚洲精品视频女| 亚洲精品中文字幕在线视频| 最近最新中文字幕免费大全7| 一区福利在线观看| 久久久久久久精品精品| 欧美日韩综合久久久久久| 9热在线视频观看99| 只有这里有精品99| 亚洲少妇的诱惑av| 色哟哟·www| 国产在线视频一区二区| 2021少妇久久久久久久久久久| 亚洲精品日韩在线中文字幕| 亚洲美女搞黄在线观看| 日本猛色少妇xxxxx猛交久久| 熟女av电影| 精品少妇内射三级| 欧美精品人与动牲交sv欧美| 18禁动态无遮挡网站| 少妇的丰满在线观看| 国产不卡av网站在线观看| 少妇猛男粗大的猛烈进出视频| 久久精品国产亚洲av高清一级| 一区二区三区激情视频| 日日撸夜夜添| 亚洲精华国产精华液的使用体验| 91精品国产国语对白视频| 国产福利在线免费观看视频| 黄色 视频免费看| 国产精品亚洲av一区麻豆 | 免费高清在线观看视频在线观看| 久久精品国产鲁丝片午夜精品| 国产精品三级大全| 女人久久www免费人成看片| 久久久久网色| 亚洲色图 男人天堂 中文字幕| 免费黄频网站在线观看国产| 丰满迷人的少妇在线观看| 国产高清不卡午夜福利| 少妇熟女欧美另类| 精品国产一区二区久久| 欧美变态另类bdsm刘玥| 午夜日韩欧美国产| 男女边摸边吃奶| 伊人久久大香线蕉亚洲五| 欧美中文综合在线视频| 性高湖久久久久久久久免费观看| 欧美日韩精品网址| 欧美日韩一级在线毛片| 一区福利在线观看| 美女福利国产在线| 麻豆av在线久日| 黑人巨大精品欧美一区二区蜜桃| 国产精品免费大片| 性高湖久久久久久久久免费观看| 777久久人妻少妇嫩草av网站| 日韩伦理黄色片| 精品一品国产午夜福利视频| 少妇猛男粗大的猛烈进出视频| 亚洲av欧美aⅴ国产| 丝袜在线中文字幕| 日韩中文字幕视频在线看片| 亚洲男人天堂网一区| 中文字幕精品免费在线观看视频| 日本av免费视频播放| 国产一区有黄有色的免费视频| 欧美日韩视频高清一区二区三区二| 免费观看av网站的网址| 欧美亚洲日本最大视频资源| av网站在线播放免费| 岛国毛片在线播放| 777久久人妻少妇嫩草av网站| 亚洲欧洲国产日韩| 精品人妻在线不人妻| 丝袜美足系列| av在线观看视频网站免费| 在线观看人妻少妇| 欧美日韩视频高清一区二区三区二| 亚洲第一区二区三区不卡| 久久久精品区二区三区| 亚洲国产欧美在线一区| 捣出白浆h1v1| 中文乱码字字幕精品一区二区三区| 最黄视频免费看| 日韩一区二区视频免费看| 精品少妇内射三级| 亚洲人成电影观看| videos熟女内射| 成人18禁高潮啪啪吃奶动态图| 中文字幕另类日韩欧美亚洲嫩草| 欧美亚洲日本最大视频资源| 免费久久久久久久精品成人欧美视频| 久久久国产精品麻豆| 久久人人爽人人片av| 七月丁香在线播放| 国产毛片在线视频| 国产老妇伦熟女老妇高清| 久久精品久久久久久噜噜老黄| 赤兔流量卡办理| 一级毛片我不卡| 欧美日韩成人在线一区二区| 在线观看一区二区三区激情| 亚洲精品久久午夜乱码| 老鸭窝网址在线观看| 日日啪夜夜爽| 国产成人av激情在线播放| 亚洲一码二码三码区别大吗| 天天躁夜夜躁狠狠躁躁| 久久久久人妻精品一区果冻| 天堂中文最新版在线下载| 久久久久久久大尺度免费视频| 欧美另类一区| 丝袜美腿诱惑在线| 如何舔出高潮| 男女高潮啪啪啪动态图| 中文字幕人妻丝袜制服| 国产成人欧美| 欧美少妇被猛烈插入视频| 国产高清不卡午夜福利| 超碰97精品在线观看| 国产男人的电影天堂91| 亚洲经典国产精华液单| 精品国产超薄肉色丝袜足j| 亚洲国产av新网站| 国产精品国产av在线观看| 99久久精品国产国产毛片| 亚洲天堂av无毛| 日韩精品免费视频一区二区三区| 精品卡一卡二卡四卡免费| 亚洲美女黄色视频免费看| 国产xxxxx性猛交| 精品少妇内射三级| 欧美精品一区二区免费开放| 一区福利在线观看| 国语对白做爰xxxⅹ性视频网站| 国产熟女午夜一区二区三区| 少妇的逼水好多| 午夜福利影视在线免费观看| 日韩中字成人| 宅男免费午夜| a级毛片黄视频| 亚洲av.av天堂| 久久精品国产亚洲av涩爱| 久久精品夜色国产| 久久久久国产一级毛片高清牌| 国产免费现黄频在线看|