雷曉珊 ，馬 勇 *，藺世杰

LEI Xiaoshan1，2，MAYong1，2*，LIN Shijie1，2

Neil Pryde RS：X（NP）級別帆板在2008年時正式成為夏季奧林匹克運(yùn)動會比賽項(xiàng)目，其主要動力來源于風(fēng)作用在帆翼上的空氣動力，帆翼空氣動力的大小直接影響比賽成績（羅曉川，2010）。帆翼空氣動力學(xué)研究的方法有理論分析、實(shí)驗(yàn)研究和數(shù)值計(jì)算。目前為止，傳統(tǒng)的計(jì)算流體動力學(xué)（computational fluid dynamics，CFD）仿真方法通常將帆翼視為剛體，只考慮氣流對帆翼氣動性能的影響，很少考慮帆翼的變形以及相應(yīng)的流固耦合問題，忽略了帆翼結(jié)構(gòu)變形和振動與氣流之間的耦合作用，進(jìn)而導(dǎo)致在氣動特性和結(jié)構(gòu)方面的計(jì)算結(jié)果與帆翼實(shí)際狀況存在一定的偏差（Augier，2012；Viola，2011；Persson et al.，2016）。在奧運(yùn)會帆船帆板比賽時，常用航線為奧林匹克航線（藺世杰等，2017），比賽中隨著風(fēng)速的變化，帆翼的結(jié)構(gòu)形狀變化較大，甚至產(chǎn)生振動，流固耦合效應(yīng)更加顯著，直接影響到帆板高速航行時帆翼空氣動力性能（Bak et al.，2013；Durand et al.，2014）。

許多學(xué)者在運(yùn)動帆翼的空氣動力性能方面做了大量研究。王樹杰等（2009）、何海峰（2012）利用數(shù)值模擬方法，研究了奧運(yùn)級別NP帆板帆翼分別在穩(wěn)態(tài)和瞬態(tài)下空氣動力性能，分析討論了不同狀態(tài)下轉(zhuǎn)角和頻率的變化對帆翼空氣動力性能的影響。Durand等（2014）基于非線性有限元方法進(jìn)行耦合，對IMOCA60帆翼順風(fēng)航行時的空氣動力性能進(jìn)行了研究，考慮了帆翼褶皺、振動等問題，為帆翼的推進(jìn)性和穩(wěn)定性優(yōu)化設(shè)計(jì)打下了基礎(chǔ)。馬勇等（2012，2013a，2013b，2013c，2016）基于風(fēng)洞試驗(yàn)與CFD方法，對奧運(yùn)會帆船帆翼的空氣動力學(xué)特性進(jìn)行研究，討論了不同參數(shù)對帆翼空氣動力性能的影響，得到了不同航向角下帆翼的升阻力、側(cè)向力和橫傾力矩等特性與帆翼攻角變化的關(guān)系。Sacher等（2017）利用ISIS-CFD與ARA耦合，對IMOCA主帆進(jìn)行流固耦合計(jì)算，發(fā)現(xiàn)帆翼的氣動力特性與試驗(yàn)結(jié)果較吻合。Deparday等（2018）對風(fēng)角為50°～140°下的三角帆帆翼進(jìn)行流固耦合研究，分析了帆翼的壓力，發(fā)現(xiàn)來流帆緣氣動彈性相對較弱且非穩(wěn)定性振動發(fā)生在前帆。目前，基于流固耦合方法對帆翼空氣動力性能的研究較多關(guān)注大三角帆，而小帆船帆翼研究主要集中在通過傳統(tǒng)CFD方法進(jìn)行空氣動力性能研究，奧運(yùn)會帆船帆板等帆翼流固耦合問題有待進(jìn)一步研究。

本文以奧運(yùn)女子級別的NP帆板帆翼為研究對象，基于ANSYS Workbench平臺，采用單/雙向流固耦合方法，對不同風(fēng)速下帆翼進(jìn)行三維數(shù)值模擬，對比雙向流固耦合與單向流固耦合仿真計(jì)算結(jié)果在氣動力、流場和結(jié)構(gòu)等方面的差異，探討單/雙向流固耦合仿真方法的差異性，并揭示流固耦合效應(yīng)對帆翼氣動性能的影響。

1 研究方法

1.1 流體動力學(xué)數(shù)值求解

流固耦合計(jì)算中，由流體動力學(xué)求解提供了帆翼表面的氣動載荷，在歐拉坐標(biāo)系中，根據(jù)三大守恒定律，可得到流體控制方程：

式中：t為時間；ρf為流體的密度→為速度矢量為流體介質(zhì)的體積力矢量；τf為剪切應(yīng)力張量；p為流體的壓力；μ為動力粘度；e為速度應(yīng)力張量。

1.2 固體模型的有限元分析

帆翼結(jié)構(gòu)的動力響應(yīng)可以用下面的運(yùn)動方程描述：

式中：M為質(zhì)量矩陣；C為阻尼矩陣；K為剛度矩陣；F為由流體導(dǎo)致的作用在結(jié)構(gòu)上的外力；u?，u?分別為速度和加速度。

1.3 流固耦合交界面的動力學(xué)條件

根據(jù)已知的邊界條件，通過逆解法和半逆解法進(jìn)行求解，得到應(yīng)變、應(yīng)力和位移的結(jié)果數(shù)據(jù)。應(yīng)力邊界條件：

式中：n為單位法向量分量；τf為流體應(yīng)力；τs為結(jié)構(gòu)應(yīng)力。

位移邊界條件：

式中：df為流體位移；ds為固體位移。

2 數(shù)值模擬的建立

2.1 計(jì)算模型與網(wǎng)格

本文以奧運(yùn)會女子級別的NP帆板帆翼為研究對象，實(shí)際航行如圖1（A），計(jì)算速度為2～12 m·s-1，涵蓋比賽過程中最常用的風(fēng)速。帆翼主要規(guī)格參數(shù)如表1，通過對NP帆翼外形的測繪得到（羅曉川，2010）?；贏NSYS Workbench平臺，采用單、雙流固耦合方法對奧運(yùn)女子級別的NP帆板帆翼進(jìn)行迭代計(jì)算，在建模軟件Design Modeler進(jìn)行帆翼模型和流動區(qū)域的構(gòu)建，得到三維模型?？紤]到計(jì)算量、流場充分發(fā)展和模擬的真實(shí)性（馬勇等，2012；Ma et al.，2016），確定采用以帆翼底邊弦長為基準(zhǔn)向前、左右和上各延伸6倍，向后延伸10倍的長方體為計(jì)算區(qū)域，帆翼和計(jì)算域位置設(shè)置如圖1（B）。

帆板運(yùn)動處于近海平面且空氣流動速度不高，該流場可被視為不可壓縮流體。通過對比不同湍流模型的優(yōu)缺點(diǎn)，采用Realizable k-ε模型來計(jì)算湍流度（馬勇等，2013a），在近壁面區(qū)使用壁面函數(shù)法進(jìn)行處理。入口邊界條件采用速度入口，入口速度按照工況要求而設(shè)定。出口條件采用自由流出，出口處速度滿足流動充分發(fā)展。

圖1 （A）實(shí)帆圖（B）NP帆翼模型圖Figure 1. （A）ASail Wing（B）Model ofASail Wing for NPClass

表1 女子NP帆板帆翼主要規(guī)格Table 1 Main Parameters of the Sail Wing for Neil Pryde RS：X Class

考慮流固耦合計(jì)算耗時較多，在保證足夠計(jì)算精度的基礎(chǔ)上，應(yīng)充分發(fā)揮有限的硬件條件。本文進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證，通過改變帆翼表面網(wǎng)格尺寸和帆翼表面第1層網(wǎng)格高度來改變網(wǎng)格數(shù)量，在不同網(wǎng)格數(shù)下對帆翼模型采用相同的算法進(jìn)行模擬。由圖2不同網(wǎng)格方案下升力、阻力系數(shù)對比可看出，當(dāng)網(wǎng)格單元數(shù)大于209萬時，網(wǎng)格單元數(shù)對升阻力系數(shù)影響極小。考慮到計(jì)算能力，同時兼顧計(jì)算精度，本研究選擇第2種方案進(jìn)行劃分，整個流體域的網(wǎng)格單元總數(shù)為209萬個。

圖2 不同網(wǎng)格方案下升力系數(shù)、阻力系數(shù)對比Figure 2.Comparison of Lift Coefficient and Drag Coefficient for Different Mesh schemes

2.2 計(jì)算方法與條件設(shè)置

在實(shí)際航行過程中，在不同風(fēng)速下伴隨著帆翼表面發(fā)生變形、搖擺等復(fù)雜的運(yùn)動，需要采用動網(wǎng)格模型來解決這一問題（郭澤宇等，2018）。而帆翼由桅桿和帆面成拱形結(jié)構(gòu)，且受風(fēng)影響沿一個方向形變，屬于小變形疊加狀態(tài)，可通過彈簧光順（霍世慧等，2011）和局部網(wǎng)格重構(gòu)（龐宇飛等，2017）2種動網(wǎng)格方法混合，能夠高精度和高效率地實(shí)現(xiàn)帆翼流固耦合過程。

帆面和桅桿使用高分子和碳纖維復(fù)合材料制造，材料屬性在ANSYS自帶的材料屬性數(shù)據(jù)庫中不存在，因此對實(shí)際比賽中的NP帆翼進(jìn)行了相關(guān)的應(yīng)力松弛實(shí)驗(yàn)，確定了各個材料的參數(shù)。為了簡化計(jì)算的結(jié)構(gòu)模型，假定NP帆翼模型帆面為各項(xiàng)同性材料，帆翼參數(shù)是：帆面彈性模量為 14.4 GPa、泊松比為 0.33、密度為 100.86 kg·m-3；桅桿彈性模量為161.7 GPa、泊松比為0.23。

3 結(jié)果與討論

3.1 升力和阻力以及壓力場分析

單/雙向流固耦合效應(yīng)的2種數(shù)值仿真方法，在求解得到的升力系數(shù)、阻力系數(shù)以及單/雙向流固耦合時差異如圖3。單/雙向流固耦合數(shù)值計(jì)算得到的NP帆翼升力系數(shù)和阻力系數(shù)隨著風(fēng)速的改變而變化的趨勢一致，都隨著風(fēng)速的增加而降低，且在風(fēng)速為2 m·s-1時升、阻力系數(shù)都是最大。而雙向流固耦合計(jì)算下，升、阻力系數(shù)變化過渡更圓滑，更接近實(shí)際情況。從氣動力數(shù)據(jù)的結(jié)果對比發(fā)現(xiàn)，2種仿真方法的差異率都在9%以內(nèi)，但是在風(fēng)速為4 m·s-1時升力系數(shù)和阻力系數(shù)差異比其他風(fēng)速時都要大，說明單/雙向耦合對帆翼空氣動力性能有不同程度的影響，而帆翼空氣動力的大小直接影響到帆板的推進(jìn)性能，關(guān)系到帆翼行駛的操縱穩(wěn)定性。而考慮流固耦合效應(yīng)的仿真方法更加貼切運(yùn)動實(shí)際，有利于帆板操縱性的準(zhǔn)確控制。

圖3 單/雙向流固耦合時（A）升力系數(shù)（B）阻力系數(shù)差異對比Figure 3.Comparison of（A）Lift Coefficient and（B）Drag Coefficient in One-way/Two-way Fluid-Structure Coupling

帆翼空氣動力的產(chǎn)生，主要是由于帆翼前后表面存在壓力差，本文進(jìn)一步對比了6種風(fēng)況下壓力場的差異。通過對比2種仿真方法帆翼表面的壓力數(shù)據(jù)發(fā)現(xiàn)，帆翼迎風(fēng)面和背風(fēng)面的壓力都有一定差異。圖4對比了流速從2～12 m·s-1仿真計(jì)算收斂后的2種耦合方式下Z=1.6 m橫截面處流場壓力的分布情況。對比顯示，由于雙向流固耦合效應(yīng)作用下帆翼面發(fā)生變形，帆翼后方負(fù)壓區(qū)增大。

圖4 不同風(fēng)速作用下Z=1.6 m截面處流場壓力分布對比Figure 4. Comparison of Pressure Distribution of Flow Field at Z=1.6 m Section at Different Speeds

隨著氣流沿著帆翼表面流動，在背風(fēng)面后部會逐漸脫落，當(dāng)風(fēng)速較高時會在背風(fēng)面產(chǎn)生明顯的渦脫現(xiàn)象，帆翼表面持續(xù)受到周期性的載荷作用而發(fā)生變形和振動，產(chǎn)生風(fēng)激振現(xiàn)象，也會使帆翼周圍流場發(fā)生變化。因此，考慮帆翼與流場的耦合效應(yīng)十分有必要。

3.2 單/雙向耦合時應(yīng)變分析

6種風(fēng)況下，帆翼的等效應(yīng)變云圖如圖5，帆翼最大應(yīng)變處都在桅桿上部附近和尾緣中部到下部，且其在氣動力的作用下極易發(fā)生變形。隨著風(fēng)速的增大，帆面最大應(yīng)變量變大，最大應(yīng)變范圍也相應(yīng)地?cái)U(kuò)大。以風(fēng)速12 m·s-1為例，采用單/雙向流固耦合計(jì)算方法，帆翼最大變形量分別為4.65和7.95 mm，皆小于弦長的1%，即23.3 mm，在其允許的最大變形量內(nèi)，符合帆翼設(shè)計(jì)的變形要求（圖6）。

圖5 不同風(fēng)速下帆翼的等效應(yīng)變云圖Figure 5. Equivalent Strain Cloud Diagram of Sail Wing at Different Wind Speeds

由于氣動力的作用，帆翼在經(jīng)過一段時間的大幅度波動后，在2.0 s后逐漸平穩(wěn)。表2為基于單/雙向流固耦合計(jì)算的帆翼等效應(yīng)變量的對比表，Do、Dt分別為單、雙向流固耦合條件下的帆翼等效應(yīng)變量。

3.3 單/雙向耦合時應(yīng)力分析

基于不同流固耦合方法計(jì)算時，帆翼應(yīng)變分布圖略有差異，但在不同風(fēng)況下，帆翼最大應(yīng)力都位于桅桿附近上部和尾緣下方區(qū)域。以風(fēng)速12 m·s-1為例，采用單/雙向流固耦合計(jì)算方法，帆翼最大應(yīng)力值分別為4.109和7.237 MPa（圖7）。表3為基于單/雙向流固耦合計(jì)算的帆翼等效應(yīng)變量和等效應(yīng)力大小的對比表，σo、σt分別為單、雙向流固耦合條件下的帆翼等效應(yīng)力值。

圖6 風(fēng)速12 m·s-1時（A）單向耦合（B）雙向耦合時帆翼的變形圖Figure 6. Deformation Diagram of the Sail Wing under（A）Oneway Fluid-Structure and（B）Two-way Fluid-Structure Interaction at 12 m·s-1

表2 單/雙向流固耦合計(jì)算的帆翼應(yīng)變量比較Table 2 Comparison of the Strains of the Sail Wing under Oneway and Two-way Fluid-Structure Interaction

表3 單/雙向流固耦合計(jì)算的帆翼應(yīng)力大小比較Table 3 Comparison of the Stress of the Sail Wing under Oneway and Two-way Fluid-Structure Interaction

通過對比基于單/雙向流固耦合方法的結(jié)構(gòu)特性計(jì)算結(jié)果可以發(fā)現(xiàn)，相同工況下，雙向流固耦合計(jì)算得到的帆翼的應(yīng)力值略大，隨著風(fēng)速增加，兩者的差值逐漸變大。結(jié)合圖6帆翼變形云圖可發(fā)現(xiàn)，雙向流固耦合數(shù)值計(jì)算得到的應(yīng)變更符合帆翼實(shí)際的變形情況?？紤]到雙向流固耦合計(jì)算時固體結(jié)構(gòu)與流體之間的相互作用，及基于雙向流固耦合方法下帆翼的變形量與應(yīng)力值比單向流固耦合計(jì)算時過渡的更圓滑，從結(jié)構(gòu)強(qiáng)度校核角度考慮，使用雙向流固耦合的計(jì)算結(jié)果更保守、更接近實(shí)際情況。

圖7 風(fēng)速12 m·s-1時帆（A）單向耦合（B）雙向耦合翼的應(yīng)力云圖Figure 7. Stress Diagram of the Sail Wing under（A）One-way Fluid-Structure and（B）Two-way Fuid-Structure Interaction at 12 m·s-1

4 結(jié)論

流固耦合效應(yīng)對帆翼的升力和阻力產(chǎn)生不同程度的影響，考慮流固耦合效應(yīng)的仿真方法對帆翼空氣動力預(yù)測更加準(zhǔn)確，有利于帆板操縱性的準(zhǔn)確控制；單/雙向耦合下帆翼最大應(yīng)變處都在桅桿上部附近和尾緣中部到下部，在氣動力的作用下極易發(fā)生變形；雙向流固耦合計(jì)算得到的帆翼的應(yīng)力值略大，隨著風(fēng)速增加，兩者的差值逐漸變大。結(jié)合帆翼變形云圖分析，雙向流固耦合數(shù)值計(jì)算得到的應(yīng)變更符合帆翼實(shí)際的變形情況。因此，研究帆翼空氣動力性能時，流固耦合效應(yīng)不可忽略，研究結(jié)果為運(yùn)動員在訓(xùn)練和實(shí)際比賽中科學(xué)調(diào)整帆翼提供依據(jù)。