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

基于離散模塊梁單元水彈性理論的復雜連接處建模方法

2022-03-19 08:40:32陳永強張宇張顯濤
中國艦船研究 2022年1期
關鍵詞:結構方法

陳永強,張宇,張顯濤*,2

1 上海交通大學 海洋工程國家重點實驗室,上海 200240

2 上海交通大學 三亞崖州灣深海科技研究院,海南 三亞 572024

0 引 言

布置于海上的超大型浮體具有多種用途,可作為浮式機場、浮橋和海洋可再生能源的浮式基礎等[1]。由于超大型浮體的水平尺寸與波長相當,且遠大于垂向尺寸,在波浪激勵下浮體的變形影響不可忽略。因此,計及流體與結構變形耦合影響的水彈性方法被用于分析在波浪作用下超大型浮體動力響應[2]。傳統的水彈性分析方法通常分為直接法和模態疊加法。其中,直接法是直接針對浮體濕表面建立流體壓力和浮體變形的耦合方程[3],而模態疊加法通常將水彈性問題分為3個步驟處理:1) 確定浮式結構的固有彈性振動模態;2)針對每個模態建立輻射和繞射問題;3)將所有模態疊加并最終求解問題。

很多研究人員已應用模態疊加法來處理超大型浮體在不同場景下的水彈性響應問題[4-6]。Lu等[7]提出基于離散模塊?梁單元(discrete-module-beam,DMB)的水彈性分析方法,在頻域范圍內處理波浪與連續彈性浮體的相互作用。在DMB水彈性理論(以下也稱DMB理論或方法)框架下,將連續彈性浮體離散成若干個剛性子模塊來模擬等效水動力性能,并假設這些剛性子模塊由等效梁單元連接,以此達到近似模擬結構變形的目的。與傳統的三維水彈性方法相比,DMB方法不需要預先確定結構的彈性模態,更適于處理一些連接形式復雜的超大型浮體動力響應問題,其原因在于,這種連接形式復雜的結構整體變形與連接處的強耦合使得彈性模態的確定變得比較困難。研究表明,DMB方法計算得到的水彈性響應結果與實際測量和模態疊加法計算得到的結果吻合較好。鑒此,很多研究人員將DMB方法推廣應用于兩模塊鉸接超大型浮體[8]和幾何形狀復雜的連續超大型浮體[9]結構的彈性分析。例如,Zhang等[10]和Wei等[11]發展了時域DMB方法,分別應用于研究超大型浮體在非穩態的外部載荷和非均勻海洋環境下的動力響應,Jin等[12]進一步考慮系泊系統并應用DMB方法分析了浮體?連接件?系泊系統耦合動力響應特性。

對于超大型浮體,更傾向于采用多模塊連接的設計方案。在采用模態疊加法對其進行分析時,一般首先基于有限元方法來確定包含連接形式的浮體彈性振動模態,然后建立水彈性運動方程進行求解[13],而在采用DMB方法進行分析時,其處理連接處的方式與模態疊加法有所不同。首先,DMB方法無需顯式求解連接多個浮體結構的振動模態;其次,DMB方法是在離散子模塊重心處建立運動方程,故對于連接處的建模,該方法的核心是將連接處引起的位移限制和力整合到浮體運動方程中。文獻[8]假設連接處及其相鄰的兩個離散子模塊重心處的結構滿足剛性條件,先通過限制矩陣法將鉸接處的力等效到相鄰的兩個子模塊重心上,再建立運動方程求解;文獻[10]則先在模塊連接處兩端引入兩個節點,建立節點力與位移的關系式,得到對應于連接處的剛度矩陣,并將該剛度矩陣整合到結構總剛度陣中,再建立運動方程進行求解。結果表明,在DMB理論框架下,上述兩種連接處建模方法得到的結果與傳統的三維水彈性計算結果吻合較好[13]。

然而,關于DMB理論框架下浮體連接處的建模有一些問題仍有待進一步分析,原因是:1)目前針對不同方式連接處的建模方法缺乏系統的比較;2)現有的連接處建模方法均會增加運動方程的變量和系數矩陣維數。因此,尋求不改變變量和系數矩陣維數的連接處建模方法也有一定的意義,尤其是針對多模塊超大型浮體連接處的動力響應分析。

綜上,本文將在DMB理論框架下討論多模塊超大型浮體連接處的建模方法及其適用性,新提出虛擬剛度法和子結構法這兩種新的建模方法,詳述兩種新方法及已有的限制矩陣法和全節點法, 并對這4種方法進行系統的對比分析。

1 基于DMB水彈性分析方法簡介

1.1 子模塊重心處的位移與受力求解

本節將簡述文獻[7]和文獻[10]提出的DMB水彈性分析方法。如圖1所示,以某個連續的彈性浮體結構為例,采用DMB方法將此浮體離散成多個剛性子模塊(旨在用于近似考慮水動力),而這些剛性子模塊則通過梁單元(旨在用于近似考慮結構變形)連接。關于DMB方法的控制方程和邊界條件可詳見文獻[10]。

圖1 在DMB理論框架下的連續型彈性浮體建模原理圖Fig.1 Modelling of a continuous elastic floating structure in the framework of DMB hydroelasticity theory

具體而言,在DMB理論框架下,一個連續的彈性浮體可被離散成N個子模塊,其中N的值取決于收斂性分析。首先,假設這些子模塊滿足剛性條件且可自由漂浮于水面,此時可運用多剛體水動力學理論求解所有子模塊的水動力參數,包括:波浪激勵力矩陣FE、附加質量力(和力矩)FA=ω2A(ω)ξ 、輻射阻尼力FRd=iωB(ω)ξ(其中:ω為波浪頻率;ξ為所有子模塊重心處的位移,其維度為 6N×1;A(ω)為所有子模塊的附加質量矩陣,其維度為 6N×6N;B(ω)為所有子模塊的輻射阻尼矩陣,其維度為 6N×6N)、所有子模塊受到的靜水回復力FHS=?Cξ (其中C為靜水回復力系數矩陣,維度為 6N×6N)、所有子模塊的慣性力FIn=ω2Mξ ( 其中M為 質量矩陣,其維度為6N×6N)。

上述所有力都假設作用在子模塊重心處,故所有子模塊被簡化成集中質量,其受力表示為

然后,假設相鄰的集中質量通過梁單元連接,基于歐拉?伯努利梁和圣維南扭曲理論來考慮彈性浮體的變形影響。因彈性浮體的質量都集中于集中質量上,故假設梁單元的質量為0。另外,梁單元的結構屬性與真實的彈性浮體的屬性保持一致。對于某個梁單元,其兩端的位移和結構力通過梁單元的剛度矩陣KE建 立聯系,KE的表達式可詳見文獻[10]。

基于有限元理論的標準方法組裝梁單元剛度陣,可得到總的結構剛度矩陣KSt(其維度為6N×6N),進而可建立作用于集中質量上的總的結構力FSt與集中質量(或子模塊重心處)位移ξ的關系,即FSt=?KStξ。

根據集中質量處的合力為0的原則,可得到:

式(1)可進一步表示如下運動方程:

定 義 初 始 的 剛 度 矩 陣Kˉo=?ω2(M+A(ω))?iωB(ω)+C+KSt,式(2)可重寫為

通過求解式(2)或式(3),所有子模塊重心處(或集中質量)的位移都可求解得到。

1.2 VLFS任意位置的位移和結構力求解

對于DMB方法,考慮鉸接的影響,求解式(2)或式(3),可得到集中質量(或每個子模塊重心處)的位移,即圖2中ξp和 ξq。而對于因結構變形所引起的梁單元兩端的結構力,即和,則可通過式(4)計算。

圖2 梁單元e與子單元f示意圖Fig.2 Schematic of a beam element e and a sub-element f

式中,KE,e為梁單元的剛度矩陣。

對于梁單元 上的任意一點,其與梁單元左端 之間構成了子單元,且剛度矩陣為 ,則點處的位移 和結構力 滿足如下關系式:er pfKE,f rξr

式中,,K,K和為剛度矩陣KE,f的子矩陣。

進而,可求得點r處的位移ξr和結構力:

在線性假設的框架下,式(7)計算得到的結構力和力矩僅在每個子模塊的端面處是“準確”的,這是因為基于DMB水彈性方法將子模塊受到的分布力等效集中在了子模塊重心處,而這種等效操作僅對子模塊端面處的力沒有影響。在得到各個子模塊端面處的“準確”的結構力和力矩之后,可利用插值法得到彈性浮體任意位置的結構力和力矩,具體可詳見文獻[8]。

2 多模塊彈性浮體連接處建模方法

本節將詳述基于DMB水彈性分析方法的4種針對多模塊彈性連接的浮體連接處的建模方法。浮體的連接方式既可是剛性連接、也可是彈性連接或鉸接等。本文將以鉸接連接方式作為研究對象,期望研究結果可很方便地推廣到其他多模塊超大型浮體及其他連接方式的浮體研究之中。

如圖3所示,整個浮體被離散成N個子模塊,連接處位于第j和j+1個子模塊之間。定義大地坐標系的原點在浮體的最左側,xoy平面在靜水面上,x軸沿浮體的縱向方向,z軸垂直向上為正。在討論的4種不同建模方法中,限制矩陣法和全節點法由文獻[8]和文獻[10]分別提出,虛擬剛度法和子結構法則為本文新提出的方法。

圖3 在DMB理論框架下的兩模塊超大型浮體連接處建模原理圖Fig.3 Schematic of modelling interface of a two-modulesinterconnected VLFS in the framework of DMB hydroelasticity theory

2.1 虛擬剛度法

圖4所示為基于DMB水彈性分析方法的兩模塊超大型浮體鉸接處及其附近的集中質量示意圖。對于虛擬剛度法,浮體的鉸接處由如下剛度矩陣Kcv表示。

圖4 在DMB理論框架下的兩模塊彈性連接處及其附近的集中質量[10]Fig.4 Schematic of the interface of elastically-interconnected modules and its two adjacent lumped masses in the framework of DMB hydroelasticity theory[10]

式中,kncv(n=1,2,···,6)分別為6個自由度方向的剛度。對于鉸接連接方式而言,應賦予kncv很大的值,以使鉸接近似于對應自由度的剛性條件。在鉸接處,因浮體可圍繞平行于y軸的軸自由轉動,k5cv=0 ,故可將其從矩陣Kcv中去除。

關于第j和j+1個集中質量的位移,可分別定義如下:uj,vj,wjxj,yj,zjαj,βj,θjxj,yj,zj

式中: 分別為沿圖中 軸的線位移;分別為繞 軸的角位移。鉸接處在其相鄰的兩個局部坐標系中的坐標分別表示為:

若忽略連接處與其相鄰的集中質量之間的結構變形影響,則鉸接處兩側的位移可用其相鄰的集中質量的位移來表示,即

式中,L為轉換矩陣(或限制矩陣[8]),其表達式如下:

通過補加0元素,轉換矩陣L可 從維度5×12拓展為維度 5×6N,即得到新的矩陣:

鉸接處的受力表示為,其中 分別是沿 軸的受力;分別是繞 軸的力矩。其中, 可通過如下等式計算。Fh=[Fh1,Fh2,Fh3,Fh4,Fh6]TFhk(k=1,2,3)x,y,z Fhk(k=4,6)x,zFh

鉸接處的受力 還可轉換為相鄰的兩個集中質量上的受力 ,即FhFhl

將由鉸接引起的集中質量的受力 加入到浮體運動方程(1)中,可得到兩模塊鉸接彈性浮體的運動方程:Fhl

值得一提的是:虛擬剛度法忽略了連接處及其附近的集中質量之間的結構變形,即采用剛性假設;另外,相比于連續彈性浮體的運動方程,該方法在對連接處建模時并不增加運動方程的維數。

2.2 限制矩陣法

文獻[8]提出采用限制矩陣法分析兩模塊鉸接超大型浮體的動力響應,此方法忽略了兩模塊連接處與其相鄰的兩個集中質量之間的結構變形。鉸接處的位移ξ滿足如下限制條件:

作用在兩個模塊鉸接處的受力Fh可等效為作用在其相鄰的兩個集中質量上的受力Fhl,即

聯立等式(3),式(15)和式(16),可得到如下所示兩模塊鉸接超大型浮體的運動方程:

通過求解等式(17),可得到所有集中質量的位移和鉸接處的受力。

2.3 子結構法

圖5為在DMB水彈性分析方法框架內用于多模塊鉸接浮體連接處建模的子結構法示意圖。對于子結構法,連接處是通過其兩側的節點(即C1和C2)來表示的,這也與文獻[10]提出的全節點法類似。在介紹該方法之前,這里先給出一些基本定義。如圖5所示,上標 1 代表節點處的力來源于其左側的梁單元(或連接處)的變形,上標2 代表節點處的力來源于其右側的梁單元(或連接處)的變形,F為 作用于節點C1上因左側梁單元變形引起的力,F為 作用于節點C1上因右側連接處引起的力,F和 F的 含義與F和 F類似。例如,第 j 個節點處的力 FSt,j由兩個分量組成,即

圖 5 在DMB理論框架下多模塊鉸接處建模的子結構法原理圖Fig.5 Schematic of modelling the hinged interconnection by substructure approach in the framework of DMB hydroelasticity theory

子結構法的核心思想是將連接處及其附近的梁單元(即eC1和 eC2)視為一個整體結構,以下介紹其數學推導過程。

首先,結構變形引起的力和連接處的力與對應的節點位移有如下關系:

式 中, KeC1和 KeC2分 別 為 梁 單 元eC1和 eC2的 剛 度 矩陣,其維度是1 2×12(具體表達式可詳見文獻[10]中的附錄A),其中上標(1,1),(1,2),(2,1)和(2,2)代表了一個矩陣的子矩陣;KC為對角矩陣。

式(18)可調整為,

式(22)在推導過程中采用了 F?2=0。這是合理的,這是因為在鉸接處兩側的節點C1和C2上沒有受到外力作用的影響。圖5描述了式(22)建立的“子結構”兩端的位移及與其受力的關系,而式(22)中K??( 其維度為1 2×12)則代表了該“子結構”的等效單元剛度。按照有限元理論的剛度矩陣的標準組裝方法,將K??置 入總的結構剛度陣KSt(式(2))中,可求得所有集中質量的位移,而鉸接處兩側的節點位移也可通過式(19)進行求解。

2.4 全節點法

文獻[10]采用全節點法考慮連接處的建模。全節點法有別于限制矩陣法,其將鉸接處的力等效到相鄰的集中質量上,在連接處兩側引入兩個額外的節點來直接考慮連接處的力和位移。如圖5 所示,連接處通過兩個節點 C1和 C2來表示,兩個節點處的位移分別定義為 ξC1和 ξC2。因連接處存在的限制條件,將節點 C1和 C2產生的力定義為 FC1和 FC2。節點處的力與位移有如下關系,

式中,對角矩陣 KC的 維度為 6 ×6,并定義為

對于 KC對角矩陣中的不同元素的取值(即knc(n=1,2,···,6))取決于連接處的具體形式。例如,對于鉸接連接的方式,只允許圍繞平行于 y軸的某一軸自由轉動,此時 k5c=0 且 knc應該賦予很大的值來近似剛性條件。

采用全節點法對連接處進行建模,兩模塊連接的彈性浮體與連續彈性浮體的唯一區別在于新增加了鉸接處兩端的兩個節點。因此,連續彈性浮體承受的結構變形引起的力 FSt(見式(1)),其維度是 6N×1)可推廣到鉸接的兩模塊彈性浮體,如下所示。

式 中: FSt,j和 ξj(j=1,2,3,···,N)分 別 為 第 j個 節 點(或集中質量)的受力和位移; F?St, ξ? 和 K?St分別為考慮了鉸接的新的結構力、位移和剛度矩陣,其維度分別是 (6N+12)×1, (6N+12)×1和 (6N+12)×(6N+12); K?St為 按 照有 限元 理 論 的 總 結 構 剛 度矩陣標準組裝程序將梁單元剛度矩陣 KE和連接處剛度矩陣 KˉC組裝得到的矩陣。

同樣的,其他矩陣如 M , A(ω), B(ω), C 和 FE也可通過補充0 元素擴展到考慮鉸接的情況,與此對應的新的矩陣分別為 M? , A?(ω), B?(ω) , C? 和 F?E(具體表達式可詳見文獻[13]的附錄C)。因此,采用全節點法可建立如下兩模塊連接的超大型浮體運動方程:

通過求解式(26),所有集中質量和連接處兩側節點的位移都可求得。可以看出,對于全節點法,兩模塊連接的超大型浮體的運動方程(即式(26))與連續超大型浮體的運動方程(即式(2))形式上基本一樣,唯一的區別是前者的矩陣維度更大。另外,全節點法還考慮了連接處及其附近的集中質量之間的結構變形。

2.5 4 種不同的連接處建模方法簡要總結

表1 給出了在DMB 水彈性分析方法框架下的4 種不同連接處的建模方法。首先,相對于連續彈性浮體的運動方程而言,限制矩陣法和全節點法都增加了運動方程的系數矩陣的維度,但本文新提出的虛擬剛度法和子結構法并未增加運動方程系數矩陣的維度。在這4 種方法中,全節點法運動方程系數矩陣維度最高,限制矩陣法和虛擬剛度法針對連接處及其附近的集中質量之間的結構都采用了剛性假設,子結構法和全節點法考慮的是結構變形的影響。總體上,子結構法和全節點法相比限制矩陣法和虛擬剛度法更“準確”。

表 1 4 種不同的連接處建模方法比較Table 1 Comparison of modelling the interconnections by different approaches

此外,限制矩陣法可以被認為是虛擬剛度法在其元素取較大值的一個特例,而虛擬剛度法還可用于分析除鉸接形式以外的其他連接方式。本文新提出的子結構法則是在全節點法的基礎上利用靜力凝聚方法推導得到。而且,因其不需要如全節點法一樣在系數矩陣中補0 元素,所以子結構法的運動方程系數矩陣是滿秩的,可方便用于基于Cummins 方程的時域問題求解。

3 結果和討論

3.1 連接處的不同建模方法驗證分析

圖6所示為文獻[13]建立的兩模塊鉸接超大型浮體模型。圖中,紅線表示鉸接,其位置在x=L/2=150 m處,在此處超大型浮體被分為左、右兩個模塊。本文將對建模方法進行驗證分析。表2給出了該模型的系數和波浪參數。對于基于DMB水彈性分析方法,圖6所示超大型浮體的左、右模塊都被離散成4個子模塊,其結果已收斂。

表2 兩模塊鉸接超大型浮體的尺度及波浪參數 [13]Table 2 Parameters of waves and two-modules-hinged VLFS [13]

本文將4種方法得到的兩模塊鉸接的超大型浮體位移和彎矩與文獻[13]基于模態疊加法的計算結果進行了對比,如圖7所示。由圖7可見,不同方法計算得到的結果總體上吻合得較好,這證實了4種方法在基于DMB水彈性分析方法下處理連接處建模的有效性。圖中,w為超大型浮體垂向位移(其值利用入射波幅A進行了無量綱化處理),λ為波長,My為彎矩。

3.2 連接處的不同建模方法的對比分析

子結構法和全節點法這兩種方法考慮了鉸接處及其附近的集中質量之間的結構變形,虛擬剛度法和限制矩陣法忽略了鉸接處附近的結構變形。但是,從圖7可看出,4種方法計算得到的鉸接的浮體位移和結構力誤差很小,這可能是由于超大型浮體的剛度相對較大,導致鉸接處及其附近的集中質量之間的結構彈性變形受到的影響相對較小所致。

圖7 不同方法計算得到的結果對比Fig.7 Comparison of calculated results by different methods

為了驗證上述解釋,本文選取一組不同的彎曲剛度k·EI(k=0.01, 0.06, 0.1, 0.6, 1.0, 1.6),且保持其他參數不變,重新計算了兩模塊超大型浮體的動力響應,結果如圖8所示。圖中,黑色虛線表示8個子模塊重心位置,選取的超大型浮體的尺度參數中除了彎曲剛度0.01改為4.77×109N·m2外,其他不變(具體見表2)。

由于限制矩陣法和全節點法的結果與虛擬剛度法和子結構法的結果非常接近,為保持圖線清晰,本節僅展示了虛擬剛度法和子結構法的結果。從圖8可看出,虛擬剛度法和子結構法這兩種方法計算得到的位移和彎矩在鉸接處附近有明顯的差異,且該差異隨著彎曲剛度的增加而減小。

圖8 不同彎曲剛度下超大型浮體垂向位移和彎矩分布Fig.8 Distribution of vertical displacement and bending moment along the hinged VLFS under different bending stiffnesses

此外,從1.2節還可知,計算超大型浮體任意一個位置處動力響應的源頭在于每個子模塊重心處的位移與受力,所以無論采取剛性還是彈性的假設,本質上,這4種方法都是對鉸接位置所在的桿件剛度矩陣和受力矩陣進行修正,以及尋找新的力與位移的關系。因此,虛擬剛度法與子結構法計算結果間的差異其根源在于鉸接處相鄰的兩個子模塊重心位置的動力響應,而各子模塊重心位置的動力響應是經過求解矩陣方程式(2)得到的,在采用不同方法計算得到的方程存在差異的情況下,很難使用理論推導給出公式中存在的差異。盡管如此,仍可探究在彎曲剛度發生變化時虛擬剛度法和子結構法兩者的計算結果在子模塊重心位置的差值。

圖9所示為彎曲剛度在0.01~2.0變化時,兩種方法計算得到的鉸接處相鄰子模塊重心的動力響應之間的差值。從圖9可明顯看出,各差值隨著彎曲剛度的增加逐漸趨于0。可見,當超大型浮體的彎曲剛度較大時,鉸接處可近似作為剛體處理而不至于產生大的誤差,但若彎曲剛度較小,鉸接處附近的變形便不可忽略。

圖9 鉸接處相鄰子模塊重心處的垂向位移、剪力和彎矩之間的差值隨彎曲剛度的變化曲線Fig.9 Variation of the D-values between two adjacent sub-modules gravity center's vertical displacement, shear force and bending moment

如圖8(c)和圖8(d)所示,對于相對較小的彎曲剛度,虛擬剛度法得到的在鉸接處附近的彎矩相比子結構法的結果偏大,若離鉸接處有一定的的距離,兩者的差別基本可忽略。對上述現象的解釋如下:

如圖10所示,鉸接處左端的彎矩由兩個因素決定,即鉸接處的力Fh3和 分布載荷q(x1),其中x1是 與鉸接處的距離;如式(2)所示,q(x1)可能包括慣性力、靜水力和水動力等。在鉸接處附近,彈性浮體的彎矩主要受到鉸接處的力的影響。

圖10 鉸接連接處的結構受力示意圖Fig.10 Schematic of the loads near the structure at the hinged position

這里,以彎曲剛度0.01EI的情況為例,分別采用虛擬剛度法和子結構法計算鉸接處的力,結果如圖11所示。由圖11可見,在鉸接處位置x=150 m,采用虛擬剛度法計算得到剪力結果(800 kN)大于子結構法計算的結果(約265 kN),這就造成了虛擬剛度法計算得到的鉸接處附近的彎矩大于子結構法計算的結果。但是,若在遠離鉸接處的位置,兩種方法計算得到的子模塊端面處的結構力(用于插值求得任意位置的彎矩)基本相同,這也就解釋了在遠離鉸接處不同方法計算得到的彎矩差別很小的原因。

4 結 語

在DMB水彈性理論框架下,本文提出了模擬連接形式復雜的超大型浮體連接處的虛擬剛度法和子結構法,并與已有的限制矩陣法和全節點法進行了對比分析。結果表明,本文新提出的方法并沒有增加運動方程的系數矩陣的維度,計算也更方便。當彈性浮體的彎曲剛度相對較小時,因為子結構法和全節點法考慮了鉸接處附近的結構變形,所以更適合該情況下的動力響應計算。

對于多鉸接或具有斷裂位置等的非連續超大型浮體,研究者只需要定義連接處的剛度陣,就可按照本文所述步驟構建新的力與位移的關系,然后進行水彈性分析。

最后,值得一提的是本文結果并沒有采用合適的實驗數據與其進行對比,因此,在未來的研究中,將針對不同的建模方法開展實驗驗證,并考慮在更復雜的應用場景(例如非均勻海洋環境、非均勻海底以及有越浪)下的適用性。

猜你喜歡
結構方法
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
學習方法
論《日出》的結構
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 97超碰精品成人国产| 国产精品视频免费网站| 国产真实乱人视频| 91视频免费观看网站| 最新精品国偷自产在线| 91欧美亚洲国产五月天| 国产欧美日韩综合在线第一| 亚洲视频欧美不卡| 欧美日本二区| 国产在线自在拍91精品黑人| 一区二区三区成人| 国产高清色视频免费看的网址| 在线精品亚洲国产| 四虎永久免费地址在线网站| 伊人久久久大香线蕉综合直播| 国产99在线| 制服丝袜一区二区三区在线| 国产精品深爱在线| 极品私人尤物在线精品首页 | 国产精品成| 国产一级二级三级毛片| 99在线观看精品视频| 欧美性猛交一区二区三区| 久久精品国产在热久久2019| 中文字幕 91| 99久视频| 中文国产成人久久精品小说| 综合成人国产| 国产一二三区视频| 青青草综合网| 国产成人高清精品免费5388| 国内老司机精品视频在线播出| 亚洲狠狠婷婷综合久久久久| 另类重口100页在线播放| 欧美综合成人| 国产肉感大码AV无码| 欧美在线黄| 白浆免费视频国产精品视频| 五月天福利视频| 午夜性刺激在线观看免费| 欧美中文字幕无线码视频| 亚洲中文精品久久久久久不卡| 91色爱欧美精品www| 99精品国产自在现线观看| 欧美在线视频不卡第一页| 久久这里只有精品免费| 五月丁香在线视频| 亚洲天堂网2014| 99热这里只有精品国产99| 国产成人AV大片大片在线播放 | 99国产精品免费观看视频| 亚洲美女一区| 久青草网站| 国产粉嫩粉嫩的18在线播放91| 亚洲精品动漫| 国产在线精品人成导航| 欧美午夜网| 91视频99| 在线亚洲小视频| 久草中文网| 国产在线观看成人91| 人人看人人鲁狠狠高清| 蜜桃视频一区二区三区| 一区二区三区四区精品视频 | 九九热视频精品在线| 特级做a爰片毛片免费69| 青青草原国产免费av观看| 91 九色视频丝袜| 国产91蝌蚪窝| 亚洲国内精品自在自线官| 99re热精品视频国产免费| 亚洲日韩国产精品无码专区| 999精品免费视频| 久久人人妻人人爽人人卡片av| 91视频区| 亚洲国产成人综合精品2020 | 亚洲欧美另类色图| 精品三级在线| 午夜不卡福利| 91尤物国产尤物福利在线| 永久免费av网站可以直接看的| 久久综合婷婷|