彈性楔形體各狀態(tài)參數(shù)對(duì)入水運(yùn)動(dòng)性能的影響

王文華，黃一，王言英，翟鋼軍，黃亞南，2（1大連理工大學(xué)a.工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室；b.船舶工程學(xué)院；c.深海工程研究中心，遼寧大連116024；2大連海洋大學(xué)航海與船舶工程學(xué)院，遼寧大連116024）

彈性楔形體各狀態(tài)參數(shù)對(duì)入水運(yùn)動(dòng)性能的影響

王文華1a，b，黃一1b，c，王言英1b，翟鋼軍1c，黃亞南1b，2
（1大連理工大學(xué)a.工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室；b.船舶工程學(xué)院；c.深海工程研究中心，遼寧大連116024；2大連海洋大學(xué)航海與船舶工程學(xué)院，遼寧大連116024）

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

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

1 引言

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

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

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

2 CFD流場(chǎng)計(jì)算理論

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

式中：ρ為流體密度，vx和vy分別為x、y方向的流體速度，p為流體壓強(qiáng)，Bx和By分別為x、y方向的質(zhì)量力加速度，在重力場(chǎng)中，Bx=0，By=-9.8m/s2。

計(jì)算域邊界條件設(shè)置如下：在外邊界處，為了允許流體可以無(wú)反射地自由出入外邊界，流體速度和密度滿足法向梯度為零的條件，壓強(qiáng)設(shè)為定值。在物體邊界處,流體速度滿足不可穿透條件，密度滿足法向梯度為零的條件。對(duì)于流體壓強(qiáng)來(lái)說，將動(dòng)量方程（2）和（3）投影到物體邊界的法向上可得：

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

此外，采用直角切割網(wǎng)格作為空間離散網(wǎng)格，處理物體入水的動(dòng)邊界問題。此網(wǎng)格系統(tǒng)作為一種新興的網(wǎng)格技術(shù)，具有生成簡(jiǎn)單、自適應(yīng)能力強(qiáng)和易于實(shí)現(xiàn)自動(dòng)化等優(yōu)點(diǎn)，并且系統(tǒng)大多數(shù)為形狀比較規(guī)則的矩形網(wǎng)格，能夠?yàn)楦呔群透叻直媛矢袷教峁┓奖愕膽?yīng)用平臺(tái)。另外，在處理動(dòng)邊界問題時(shí)，對(duì)于直角切割網(wǎng)格系統(tǒng)來(lái)說，只需局部更新一些切割網(wǎng)格的信息即可，從而避免了傳統(tǒng)的非結(jié)構(gòu)網(wǎng)格重構(gòu)所帶來(lái)的誤差。網(wǎng)格生成算法和數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)的具體描述可以參照文獻(xiàn)[12]和[13]。

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

3 結(jié)構(gòu)有限元和板條梁?jiǎn)卧碚?/h2>

在流場(chǎng)CFD算法的基礎(chǔ)上，本文進(jìn)一步采用基于板條梁?jiǎn)卧Ｐ偷挠邢拊椒〝?shù)值求解和分析自由入水楔形體彈性邊界的結(jié)構(gòu)變形和響應(yīng)。

首先，建立固定在地面的總體坐標(biāo)系和隨著楔形體一起運(yùn)動(dòng)的局部坐標(biāo)系。將彈性楔形體的邊界離散成許多等間距板條梁?jiǎn)卧藛卧梢酝瑫r(shí)承受彎曲和拉壓作用。此外，在局部坐標(biāo)系中，將單元k的節(jié)點(diǎn)平動(dòng)位移u1，w1，u2，w2和節(jié)點(diǎn)轉(zhuǎn)動(dòng)位移θ1，θ2描述如下：

然后，在局部坐標(biāo)系中，根據(jù)板條梁?jiǎn)卧匦越卧|(zhì)量矩陣和剛度矩陣，并將流場(chǎng)計(jì)算所得到的壓強(qiáng)分布轉(zhuǎn)化到單元節(jié)點(diǎn)上。再利用坐標(biāo)變換的方法將單元局部坐標(biāo)下的質(zhì)量矩陣、剛度矩陣和節(jié)點(diǎn)載荷向量變到總體坐標(biāo)系上，并裝配在一起形成總體坐標(biāo)系下的整體動(dòng)力平衡方程如下：

此處，忽略結(jié)構(gòu)阻尼；F｛｝是總體坐標(biāo)系下的節(jié)點(diǎn)載荷向量，具體計(jì)算方法可以參考本文第4.2.1節(jié)；在總體質(zhì)量陣M［、和剛度陣K［、中，結(jié)構(gòu)單元k所對(duì)應(yīng)單元質(zhì)量陣M［、k和剛度陣K［、k表示如下：

式中：對(duì)于板條梁結(jié)構(gòu)單元來(lái)說，ρs是結(jié)構(gòu)材料密度，h是結(jié)構(gòu)單元厚度，l是結(jié)構(gòu)單元長(zhǎng)度，A=1.0*h是結(jié)構(gòu)單元的橫剖面面積，I=1.0*h3/12是結(jié)構(gòu)單元的慣性矩，E1=E/1-μ2（）是板條梁的等效彈性模量，E是結(jié)構(gòu)材料的彈性模數(shù)，μ是結(jié)構(gòu)材料的泊松比。

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

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

對(duì)于彈性楔形體自由入水模型來(lái)說，因?yàn)榱鞴恬詈蠁栴}涉及到彈性楔形體的整體運(yùn)動(dòng)、局部變形和入水流場(chǎng)的相互作用，所以如何將流場(chǎng)和結(jié)構(gòu)的求解聯(lián)系起來(lái)是非常重要的。因此，這里基于第2章的流場(chǎng)計(jì)算和第3章的結(jié)構(gòu)有限元分析理論討論彈性楔形體自由入水的流固耦合策略，以及流場(chǎng)壓強(qiáng)和結(jié)構(gòu)單元節(jié)點(diǎn)位移、速度在彈性邊界處的數(shù)據(jù)傳遞方法。

4.1 流固耦合策略

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

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

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

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

另一方面，在彈性楔形體的中低速入水過程中，彈性楔形體的局部變形很小，從而局部結(jié)構(gòu)變形對(duì)同一物理時(shí)刻流場(chǎng)的反作用可以忽略不計(jì)。因此，為了提高計(jì)算效率，這里采用單向耦合策略求解流場(chǎng)和局部結(jié)構(gòu)變形的相互作用。

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

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

在上述流固耦合算法中，涉及到流場(chǎng)和結(jié)構(gòu)計(jì)算結(jié)果的數(shù)據(jù)傳遞問題。根據(jù)流場(chǎng)的直角切割網(wǎng)格和彈性結(jié)構(gòu)的板條梁?jiǎn)卧奶攸c(diǎn)，兩種網(wǎng)格單元在耦合界面上并不直接匹配，從而不能直接進(jìn)行數(shù)據(jù)交換，因此這里需要推導(dǎo)出適合本文算法的邊界數(shù)據(jù)傳遞方法。因?yàn)樵诹鞴恬詈辖缑嫔蟼鬟f的數(shù)據(jù)包括流場(chǎng)壓強(qiáng)和結(jié)構(gòu)單元節(jié)點(diǎn)位移、速度，所以這里將數(shù)據(jù)傳遞分兩個(gè)方面加以討論，一方面是將分布式流場(chǎng)壓強(qiáng)轉(zhuǎn)化為節(jié)點(diǎn)集中力，另一方面是將結(jié)構(gòu)響應(yīng)反饋為流場(chǎng)計(jì)算的邊界條件。

4.2.1 分布式流場(chǎng)壓強(qiáng)轉(zhuǎn)化為節(jié)點(diǎn)集中力

圖2 流場(chǎng)切割和結(jié)構(gòu)單元不同幾何關(guān)系的示意圖Fig.2 Different position relations of cut cell and structural elements

根據(jù)力和力矩的平衡原則，將流場(chǎng)計(jì)算獲得的物面分布式壓強(qiáng)向結(jié)構(gòu)單元節(jié)點(diǎn)傳遞，得到節(jié)點(diǎn)等效載荷向量。根據(jù)流場(chǎng)直角切割網(wǎng)格和彈性結(jié)構(gòu)板條梁?jiǎn)卧奈恢藐P(guān)系，大致分為三類情況進(jìn)行討論：（1）結(jié)構(gòu)單元L1L2完整包含切割單元的物體邊界AB，如圖2（a）所示；（2）物體邊界AB完整包含結(jié)構(gòu)單元L1L2，如圖2（b）所示；（3）物體邊界AB與結(jié)構(gòu)單元L1L2部分相交。這里，僅以圖2（a）的結(jié)構(gòu)單元L1L2為例說明如何根據(jù)分布在物體邊界AB上的流場(chǎng)壓強(qiáng)得到節(jié)點(diǎn)L1、L2的等效載荷F1和F2。

首先，以節(jié)點(diǎn)L1為坐標(biāo)原點(diǎn)在結(jié)構(gòu)單元L1L2上建立局部坐標(biāo)系x’o’y’，其中坐標(biāo)軸o’x’沿物面逆時(shí)針方向，由L1指向L2；坐標(biāo)軸o’y’垂直物面，由流體指向結(jié)構(gòu)。因此o’x’和o’y’的單位向量為：

其中：（x1，y1）和（x2，y2）分別為節(jié)點(diǎn)L1和L2的位置坐標(biāo)。再根據(jù)流場(chǎng)切割單元中心O點(diǎn)的坐標(biāo)（xo，yo），位置向量OL1可以寫成：

最后，根據(jù)力和力矩的平衡原則將物體表面AB上的分布式壓強(qiáng)向結(jié)構(gòu)單元節(jié)點(diǎn)L1和L2轉(zhuǎn)化，從而獲得結(jié)構(gòu)板條梁?jiǎn)卧狶1L2的節(jié)點(diǎn)等效載荷F1和F2：

式中：Po和Δp=pxi+pxj是通過流場(chǎng)CFD算得的切割單元中心O點(diǎn)處的流場(chǎng)壓強(qiáng)和壓強(qiáng)梯度。

同理，采用與圖2（a）類似的方法和思路進(jìn)行推導(dǎo)可以獲得其他兩種情況（圖2（b）和圖2（c））下結(jié)構(gòu)板條梁?jiǎn)卧狶1L2的節(jié)點(diǎn)等效載荷。

4.2.2 結(jié)構(gòu)響應(yīng)反饋為CFD邊界條件

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

首先，根據(jù)結(jié)構(gòu)單元L1L2的形函數(shù)可以獲得L1L2上物體邊界端點(diǎn)A在局部坐標(biāo)系x’o’y’沿o’x’和o’y’軸的的速度和w˙A：

式中：N1~N4和N5、N6分別是結(jié)構(gòu)單元彎曲和拉壓的形函數(shù)。對(duì)于梁?jiǎn)卧獊?lái)說，L1A/L1L2和=2（L1A-L1L2/2）/L1L2是結(jié)構(gòu)單元L1L2的自然坐標(biāo)。同理，可以獲得物體邊界另一端點(diǎn)B的速度和B。

然后，物體邊界AB的平均法向速度vsn和切向速度vst可以計(jì)算如下：

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

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

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

5.1 數(shù)值算法驗(yàn)證

為了驗(yàn)證本文數(shù)值方法的準(zhǔn)確性和合理性，這里基于孫輝（2003）[5]的彈性楔形體入水實(shí)驗(yàn)進(jìn)行建模并加以數(shù)值計(jì)算，然后將三種彈性楔形體入水模型的數(shù)值結(jié)果和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較。基于實(shí)驗(yàn)水池參數(shù)建立0.8m×1.45m方形計(jì)算域，水深為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 計(jì)算域內(nèi)部的整體網(wǎng)格和物面邊界附近的局部網(wǎng)格Fig.3 Globalmesh(a)and local grid(b)near solid boundary

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

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

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

5.2 彈性楔形體各狀態(tài)參數(shù)對(duì)入水運(yùn)動(dòng)性能的影響

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

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

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

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

6 結(jié)論

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

本文為彈性楔形結(jié)構(gòu)的自由入水問題的求解提供了一種新的CFD算法，并且所得結(jié)論能夠?yàn)檫M(jìn)一步研究船舶砰擊問題、指導(dǎo)船舶航行操縱和結(jié)構(gòu)設(shè)計(jì)提供一些理論基礎(chǔ)。此外，在以后的研究中可以進(jìn)一步考慮將流固相互作用的緊耦合策略和結(jié)構(gòu)變形的幾何非線性理論引入算法中，從而擴(kuò)展本文數(shù)值方法的應(yīng)用領(lǐng)域和范圍。

[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]李國(guó)鈞,黃震球.平底物體對(duì)水面的斜向沖擊[J].華中理工大學(xué)學(xué)報(bào),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]孫輝,盧熾華,何友聲.二維楔形體沖擊入水時(shí)的流固耦合響應(yīng)的實(shí)驗(yàn)研究[J].水動(dòng)力學(xué)研究與進(jìn)展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]顧懋祥,程貫一,張效慈.平頭旋轉(zhuǎn)殼撞水水彈性效應(yīng)的研究[J].水動(dòng)力學(xué)研究與發(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]盧熾華,何友聲.二維彈性結(jié)構(gòu)入水沖擊過程中的流固耦合效應(yīng)[J].力學(xué)學(xué)報(bào),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]陳占暉,盧永錦.高速運(yùn)動(dòng)物體砰擊水面后的動(dòng)力學(xué)特性研究[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]王文華.物體入水問題的分析研究及其在船舶與海洋工程中的應(yīng)用[D].大連:大連理工大學(xué),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

國(guó)家自然科學(xué)創(chuàng)新研究群體基金資助項(xiàng)目“海洋環(huán)境災(zāi)害作用與結(jié)構(gòu)安全防護(hù)”（50921001）；國(guó)家青年科學(xué)基金項(xiàng)目“運(yùn)動(dòng)物體穿過深海密度分界面的水動(dòng)力特性研究”（11202047）

王文華（1981-），男，博士后，E-mail：wangwenhua0411@yahoo.cn；黃一（1964-），男，大連理工大學(xué)教授，博士生導(dǎo)師。