無人帆船翼帆桅桿根部應力影響因素及其影響規律研究

牛溫超 孫朝陽 盛忠起 俞建成 趙寶德 趙文濤

1.中國科學院沈陽自動化研究所機器人學國家重點實驗室,沈陽,110016 2.中國科學院機器人與智能制造創新研究院,沈陽,110169 3.中國科學院大學機器人與智能制造學院,北京,100049 4.東北大學機械工程與自動化學院,沈陽,110819

0 引言

無人帆船是一種以海洋環境能源為驅動的新型長續航力海洋機器人[1-2],它通過風帆將海洋風能直接轉化為航行驅動力[3-4],通過太陽能電池板獲取電能給控制系統和載荷使用,突破了傳統海洋機器人能源供給的限制。與燃油無人艇或電動無人艇相比,無人帆船可以實現長航時、大范圍、低成本、高時空分辨率的海洋觀測[5],是國際海洋機器人研究的前沿熱點。

剛性翼帆是一種受飛機機翼啟發而研制的適用于無人帆船的新型風帆,它憑借空氣動力學性能穩定、便于自動化控制等突出優點,已成為無人帆船的主流風帆[6]。剛性翼帆是無人帆船的唯一動力來源,其結構安全對無人帆船的航行能力具有重要影響。剛性翼帆的失效形式主要包括帆面破損、帆面變形、桅桿折斷等,其中以桅桿斷裂最為嚴重。無人帆船在海上航行時,船體的橫搖和俯仰等運動會在翼帆桅桿根部形成應力,在應力的長期作用下,桅桿根部會發生疲勞損傷甚至疲勞斷裂,因此,分析桅桿根部的受力對開展翼帆的可靠性設計、保障無人帆船的航行安全具有重要意義。

目前在無人帆船風帆研究方面,學者們大多聚焦于不同帆裝形式[7]、翼帆的翼型[8-9]、翼帆三維構型[10-12]和翼帆控制[13-14]等方面,極大促進了無人帆船風帆技術與理論的發展,如MEILHAN[15]研究了惡劣天氣下帆船碳纖維桅桿的失效問題;PELLICIOLI等[16]研究了無人帆船航行過程中風帆的受力情況。但翼帆結構安全方面的研究報道較少。本文研究了翼帆桅桿根部應力的影響因素及其影響規律,為無人帆船高性能翼帆的設計和長期航行目標的實現提供理論支撐。

1 翼帆等效模型概述

剛性翼帆結構與機翼結構類似,由蒙皮、肋板和桁架等構成[17-18],但兩者也存在明顯區別：機翼與機身固定連接,而翼帆通過桅桿與船身轉動連接,即翼帆與船身之間是一個僅存在轉動自由度的剛性連接,如圖1中四種無人帆船。通常將桅桿置于翼帆的空氣動力學中心位置(即距離翼帆前緣1/4弦長位置處),以降低轉帆耗能,但會造成翼帆的質心與桅桿軸線的偏離,如圖2所示。

(a)Salidrone(美國) (b)Sailbuoy(挪威)

圖2 翼帆結構Fig.2 Structure of wing sail

無人帆船在海上航行時在風浪的作用下會發生橫搖俯仰等運動,前期海上試驗發現橫搖運動最為劇烈,因此,本文主要探究橫搖運動過程中翼帆桅桿根部應力的影響因素及其影響規律。橫搖運動中,桅桿受慣性載荷的影響最大,而慣性載荷主要與翼帆質量、質心位置、橫搖加速度等相關,暫不考慮橫搖過程中的空氣阻尼和風載荷。對翼帆桅桿根部進行受力分析時對翼帆進行簡化,并抽象設計等效構件,如圖3所示。等效構件中桅桿結構保持不變,桅桿高度為等效后質心高度,在桅桿頂端設置配重,使等效構件的總質量、質心位置與真實翼帆參數一致。等效構件頂端安裝法蘭盤,可以調節配重相對于桅桿的安裝角度以模擬翼帆轉動后的重心位置。

圖3 翼帆等效構件Fig.3 Equivalent component of wing sail

2 桅桿根部應力數值計算模型的建立

根據物體動平衡原理,船舶的平衡條件為∑M=0[19],船舶橫搖運動方程可寫為

(1)

式中,M為船舶所承受的合力矩;φ為翼帆等效構件橫搖角;φ0為橫搖幅值(極限橫傾角);v為等效構件桅桿頂部的速度;ω為桅桿頂部的角速度;t為船舶橫搖過程的時刻。

船舶在規則波中的橫搖僅由波浪的強迫橫搖所決定[19],具體形式如下：

φ=φ0sin(ω0t+a)+A

(2)

(3)

式中,ω0為橫搖固有角速度;TRoll為帆船的橫搖周期;A、a為常數。

建立等效構件的坐標系如圖4所示。等效構件繞x軸的旋轉角度φ為帆船的橫搖角,繞z軸旋轉角度β為翼帆的偏轉角。

圖4 翼帆等效構件的坐標系Fig.4 Coordinate system of the equivalent model of wing sail

運用理論力學整體隔離法[20]分別進行受力分析,將配重和桅桿視為兩個整體,配重在橫搖過程中的受力如圖5a所示。此時,翼帆等效構件的偏轉角為β,如圖5b所示。翼帆隨帆船橫搖運動的過程可以抽象為翼帆等效構件繞橫搖軸(x軸)的橫搖運動。

(a)等效構件的受力分析 (b)偏轉示意圖圖5 翼帆等效構件橫搖過程中的受力分析Fig.5 Force analysis during the roll of the equivalent model of wing sail

定義繞原點O順時針方向的切向力為正,指向圓心O的法向力為正。將桅桿視為剛體,對配重進行受力分析,得

F2sinφ=F3cosφ

(4)

-F2cosφ-F3sinφ=G

(5)

式中,F2為桅桿對配重塊的支撐力;F3為桅桿與配重塊之間的摩擦力;G為配重的重力。

聯立式(4)和式(5),得

(6)

(7)

由式(2)、式(3)可知,翼帆等效構件的橫搖角度φ是關于時間t、帆船橫搖幅值φ0、橫搖周期TRoll的正弦函數。配重隨翼帆等效構件的橫搖運動可以看作是繞原點O做周期性的變加速圓周運動。桅桿頂部的角速度ω和角加速度α分別為

(8)

(9)

由變加速圓周運動的計算公式可求得配重質心B處的切向加速度和法向加速度：

at=αlOC

(10)

an=ω2lOC

(11)

式中,lOC為桅桿底面中心O與桅桿重心C之間的距離。

則配重所受切向力Ft和法向力Fn分別為

Ft=mat=mαlOC

(12)

Fn=man=mω2lOC

(13)

式中,m為配重質量。

如圖6所示,配重質心位置B處的切向合力FT和法向合力FN分別為

圖6 載荷作用示意圖Fig.6 Schematic diagram of load action

FT=F3+Ft

(14)

FN=F2+Fn

(15)

聯立式(6)～式(15),得

(16)

(17)

由圖6可知,切向合力FT和法向合力FN的作用點位于等效模型配重的重心位置,并且相對于桅桿處于偏心狀態。在翼帆橫搖和偏轉的耦合運動過程中,切向合力FT在桅桿頂部產生繞x軸的彎矩Mx和繞z軸的扭矩N分別為

Mx=FTl1cosβ

(18)

N=FTlBCcosβ

(19)

采用微元法[21]分析,任取桅桿上一微分單元。法向合力FN在桅桿頂部產生繞y軸的彎矩為

My=FTl1sinβ+FNlBC

(20)

式中,l1為該單元與桅桿頂部中心位置C的距離;lBC為配重質心位置B至C點的距離。

將彎矩Mx和My合成可得合彎矩M：

(21)

合彎矩M在桅桿根部產生的彎曲應力為

(22)

式中,y為該單元至桅桿軸心線的距離,ymax=D/2;D為桅桿外圓直徑;d為桅桿內圓直徑,λ=d/D。

由式(22)可知,彎曲應力σa與合彎矩M、y均成正比。因此,桅桿根部表面的單元處的正應力最大,此時l1=lOC,y=ymax。

同時,法向合力FN擠壓桅桿產生壓應力

(23)

式中,R為桅桿外圓半徑;r為桅桿內圓半徑。

則桅桿根部處的正應力σ為

(24)

扭矩N在桅桿根部產生的剪切應力為

(25)

由式(24)可知,正應力σ主要由彎曲應力σa和擠壓應力σc組成。其中,σa的方向取決于表面單元的位置。桅桿根部在彎矩M的作用下產生彎曲,故桅桿根部一側受拉一側受壓。受拉側表面單元的彎曲應力σa方向與擠壓應力σc相反,而受壓側表面單元的σa方向與σc相同。因此,受壓側表面單元的正應力σ最大。

根據式(22)～式(25),桅桿在擠壓、扭轉和彎曲的組合作用下,最大應力發生于桅桿根部表層。由于從構件表層取出的微分單元體就屬于二向應力狀態[21],故取桅桿根部表層的微分單元體,此時桅桿在彎扭組合作用下處于二向應力狀態。根據解析法[21],所取出微分單元體研究點的主平面不是此時的橫截面,而應是傾斜一定角度的斜截面。此時,有

(26)

(27)

式中,σ1、σ2、σ3分別為該微分單元的第一、二、三主應力。

二向應力狀態下,根據第四強度理論[21],可得等效應力σr4：

(28)

3 模型驗證

為驗證所建數值計算模型計算結果的正確性,利用有限元計算方法計算了桅桿根部應力值,并將該值與采用數值計算模型獲得的應力值進行對比。

分別設定7個不同參數組合進行10組有限元計算,并與所建模型計算得到的數據進行對比。7種不同參數組合見表1。根據式(2),橫搖角度是隨時間實時變化的,在驗證過程中選擇橫搖周期中第一次達到橫搖幅值的時刻。

表1 模型驗證參數組合

根據表1中各組參數組合,采用數值計算模型和有限元計算得到的桅桿根部應力數值及偏差如圖7所示。

(a)應力隨橫搖周期的變化曲線 (b)應力隨橫搖幅值的變化曲線 (c)應力隨橫搖角度的變化曲線

由圖7可以看出,數值計算模型和有限元計算得到的數據之間的誤差均小于3.6%。圖7c中,翼帆橫搖角度為2°時,兩種計算方法得到的數據偏差達20.9%,計算結果間的誤差相較于其他組明顯過大,認為它屬于隨機誤差,忽略不計。由圖7可知,模型計算和有限元計算結果的變化趨勢一致。

為了進一步驗證模型,分別計算了四次驗證工作的均方根誤差(root mean square error,RMSE),用以表征模型計算結果和有限元計算結果的誤差。計算出的RMSE數值見表2。

表2 模型計算結果與有限元計算結果的均方根誤差

根據統計學觀點：均方根誤差一般不超過2是合理的[22]。表2中各均方根誤差均小于2,證明本文模型與有限元計算結果相符,因此,利用本文模型可以正確獲得真實橫搖運動中桅桿根部的應力值。

4 影響因素和影響規律分析

根據式(28),等效應力σr4主要隨正應力σ和切應力τ的變化而變化。根據式(18)～式(25),影響正應力σ大小的參數主要有合彎矩M及法向力FN,而合彎矩M的計算主要依靠FT、FN、lOC、lBC和β,以及桅桿的內外徑尺寸d和D;影響切應力τ大小的參數主要有扭矩N和β,以及桅桿的內外徑尺寸d和D。根據式(16)、式(17),影響切向力FT和法向力FN的參數有G、φ、φ0、T和lOC。

綜上所述,等效應力σr4的值取決于T、φ0、φ、lOC、lBC、m和β,以及桅桿的內外徑尺寸d和D。本文主要探究橫搖周期T、橫搖幅值φ0、橫搖角度φ、翼帆重心高度lOC、翼帆重心B到桅桿重心C的偏心距離lBC、翼帆質量m和偏轉角度β這7個參數對桅桿應力的影響規律和影響機理。利用所建數值計算模型分別計算表3中各參數組合條件下桅桿根部的最大應力值。

表3 參數組合條件

4.1 橫搖周期的影響規律

在表3中第1組參數組合條件下,計算得到的桅桿根部應力值隨無人帆船橫搖周期T的變化情況如圖8所示。由圖8可知,隨著橫搖周期的延長,桅桿根部的應力值減小;同時,結合圖7a中的變化曲線可知,應力值隨橫搖周期呈現遞減的指數變化趨勢。出現這種變化趨勢的原因是：式(3)中橫搖角速度與橫搖周期為反比例關系,式(16)、式(17)中FT和FN與橫搖角速度近似為負二次冪函數關系。同時,FT和FN的大小直接影響桅桿根部的彎矩M、扭矩N和壓應力σc,進而影響正應力σ和切應力τ,最終導致應力值與橫搖周期呈現遞減的負二次冪函數關系。

圖8 應力值隨橫搖周期的變化Fig.8 The stress value varies with the roll cycle

4.2 φ0、lOC、lAC和m的影響規律

在表3中第2～5組參數組合條件下,分別計算得到的桅桿根部應力值隨橫搖幅值φ0、翼帆重心高度lOC、偏心距離lBC和翼帆質量m的變化情況,如圖9～圖12所示?？芍?應力值分別隨φ0、lOC、lBC和m的增大而增大;參考圖7b、圖7e～圖7g中變化曲線可知,應力值隨φ0、lOC、lBC與m的增大均呈現正比例增大的趨勢。

圖9 應力值隨橫搖幅值的變化Fig.9 The stress value varies with roll amplitude

圖10 應力值隨翼帆重心高度的變化Fig.10 The stress value varies with height of mast’s center of gravity

圖11 應力值隨偏心距離的變化Fig.11 The stress value varies with eccentricity distance

圖12 應力值隨翼帆質量的變化Fig.12 The stress value varies with the mass of wing sail

由式(16)、式(17)可知,切向力FT和法向力FN與橫搖幅值φ0成正比,導致桅桿根部的彎矩M、扭矩N和壓應力σc同樣正比例增大,進而引起正應力σ和切應力τ正比例變化,最終使等效應力σr4隨橫搖幅值φ0呈現正比例增大的趨勢。

應力值隨lOC正比例增大的原因主要有兩方面：①與橫搖幅值φ0對等效應力σr4的影響機理相同,根據式(16)、式(17),切向力FT和法向力FN是隨翼帆重心高度lOC正比例增大的,從而引起等效應力σr4的正比例變化;②根據式(18)和式(20),彎矩Mx和My與lOC均為正比例關系,從而導致正應力σ正比例增大。

應力值隨偏心距離lBC正比例增大的原因是：式(19)中扭矩N與偏心距離lBC、式(20)中彎矩My與偏心距離lBC均為正比例關系,從而引起切應力τ與正應力σ正比例增大。最終,導致應力值與偏心距離lBC成正比例的關系。

分析翼帆質量m對應力值的影響機理,整理式(16)、式(17),得

(29)

(30)

由式(29)、式(30)可知,切向力FT和法向力FN與翼帆質量m成嚴格的正比例關系,引起彎矩M、扭矩N和壓應力σc的正比例變化,進而引起正應力σ和切應力τ的正比例變化,從而使應力值與翼帆質量也呈現正比例關系。

4.3 翼帆轉角的影響規律

在表3中第6組參數組合條件下,計算得到的桅桿根部應力值隨翼帆偏轉角度β的變化如圖13所示。由圖13可知,應力值隨著偏轉角度β的增大也表現出逐漸增大的趨勢。由圖7d可以看出,應力值與偏轉角度之間的變化曲線并非正比例增大的曲線,而是近似于定義域在0°～90°內的正弦曲線。根據式(18)、式(20),彎矩Mx和扭矩N隨偏轉角度β的變化曲線近似于余弦曲線,而彎矩Mx隨偏轉角度β的變化曲線近似于正弦曲線。因此,以表1中第1組參數組合條件所得數值為例,分別計算Mx和My的合彎矩M、扭矩N與桅桿根部應力值,得到的數值變化曲線如圖14所示。

圖13 應力值隨偏轉角度的變化Fig.13 The stress value varies with deflection angle

圖14 數值變化曲線Fig.14 The curve of values

由圖14可知,扭矩N隨β的變化曲線近似于余弦曲線,彎矩M隨β的變化曲線近似于正弦曲線,并且與彎矩M相比,扭矩N較小。而應力值隨β的變化曲線也近似于正弦曲線,同時表明應力值受彎矩的影響較為顯著。

4.4 無人帆船橫搖角度的影響規律

無人帆船的橫搖角度φ是隨時間變化的,因此圖中時間坐標軸又可表示橫搖角度坐標軸,0～1 s內桅桿橫搖至左側幅值后返回中位,1～2 s內桅桿橫搖至右側幅值后返回中位。

分析圖8～圖13中時間坐標軸與應力坐標軸可以看出,每半個周期又可看作一個小周期;橫搖角度值越大,應力值越大;應力值在一個橫搖周期內隨時間(橫搖角度)的變化曲線近似于正弦曲線。出現這種變化趨勢的原因是：式(16)、式(17)中,一個橫搖周期中μ和ρ隨時間的變化曲線如圖15a和圖15b所示。

(a)μ值隨時間變化曲線

圖15a中μ隨時間的變化曲線為一正弦形波,圖15b中ρ隨時間的變化曲線為一余弦形波。圖15c所示為圖15a和圖15b中正弦波與余弦波疊加后波形,為正弦波形。根據式(2)和圖15c,t在0～0.5 s時橫搖角度與μ和ρ疊加后的數值均逐漸增大到最大值。因橫搖角度并未規定正方向,因此,在計算中橫搖角度φ均為正值。

圖15d所示為圖12中翼帆質量m=12 kg時的數值案例所得到的應力值隨時間的變化曲線。根據圖15c和圖15d,應力值對應于μ和ρ疊加后的正弦波形,并且根據圖15d中的變化曲線與式(2),應力值和橫搖角度與時間均成正弦變化關系,結合圖7c的變化曲線可知,應力值與橫搖角度成正比。

5 桅桿根部結構優化與實驗

5.1 針對性結構優化

由前文計算和仿真可知,桅桿最大應力位于桅桿根部位置。為減小桅桿根部應力,進行了針對性的結構優化：在桅桿根部添加了四個均布的加強筋,同時為了使加強筋與船體的連接牢固,在加強筋下方添加了一個法蘭盤,如圖16a所示。分別采用表1中第2、4、6、7組的參數組合進行了無加強筋與添加加強筋之后桅桿應力的有限元仿真計算,以驗證該結構優化的效果。此外,為了減少變量,在桅桿根部同樣添加了相同的法蘭盤,如圖16b所示。

(a)有加強筋的桅桿根部 (b)無加強筋的桅桿根部圖16 結構優化前后的桅桿根部Fig.16 Mast roots before and after structural optimization

通過有限元仿真計算發現,桅桿僅添加法蘭盤之后產生最大應力的位置仍為桅桿根部位置,如圖17a所示;而在桅桿根部添加法蘭盤和加強筋之后產生最大應力的位置變為加強筋頂端與桅桿表面接觸的位置,如圖17b所示。

(a)無加強筋時的仿真結果 (b)有加強筋時的仿真結果圖17 結構優化前后的仿真結果對比Fig.17 Comparison of simulation results before and after structural optimization

在表4中各參數組合條件下,分別對添加加強筋前后的翼帆桅桿進行有限元計算,得到桅桿根部的最大等效應力值數據對比,如圖18所示。

表4 優化對比仿真參數組合條件

(a)橫搖幅值不同時

由圖18可知,相同條件下桅桿根部添加加強筋后仿真計算得到的最大應力值均小于無加強筋的最大應力值,證明在桅桿根部添加加強筋的方式可以有效減小在不同受力條件下桅桿根部的最大應力值,從而保證翼帆的正常工作。同時可以看出,在相同的參數條件變化下,添加加強筋與無加強筋時桅桿的最大應力數值的變化趨勢基本一致。

5.2 優化設計結果實驗驗證

為驗證有限元分析結果的準確性,同時也為了驗證添加加強筋后桅桿根部的應力安全性,進行了翼帆桅桿等效構件的橫搖實驗。橫搖實驗借助北京航空精密機械研究所研制的T1601三軸轉動平臺(以下簡稱平臺),如圖19a所示,操作界面如圖19b所示。平臺最大負載40 kg,角度控制可精確到0.001°。選擇正弦搖擺運動方式,使翼帆桅桿等效構件以標準正弦波的運動方式擺動,并在桅桿根部監視點處粘貼應變片,以監視各點的應力變化。

(a)三軸轉動平臺 (b)平臺的操作界面圖19 橫搖試驗平臺Fig.19 Roll test platform

翼帆桅桿橫搖實驗過程如圖20所示。主要步驟如下：①根據翼帆幾何參數模型設計翼帆等效構件;②將加工好的翼帆等效構件進行安裝;③在桅桿根部設置監視點粘貼應變片;④連接應變儀并進行應變儀的校準和調試;⑤將翼帆等效構件安裝到T1601三軸轉動平臺上;⑥啟動三軸轉動平臺,調整翼帆等效構件的橫搖角度和橫搖幅值;⑦記錄桅桿根部的應力變化曲線;⑧檢查桅桿根部的變形情況以及裂紋位置。

圖20 翼帆桅桿橫搖實驗步驟Fig.20 The experimental steps of wing sail mast roll

在T=2 s、不同橫搖幅值和偏轉角度的條件下,針對桅桿根部應力值分別進行了有限元仿真計算。翼帆桅桿在不同橫搖幅值、翼帆轉角下的試驗測試結果和有限元仿真結果對比如圖21和圖22所示。由圖21、圖22可知,實驗測試和有限元計算得到的桅桿根部監視點的應力變化曲線的趨勢一致。

(a)橫搖幅值為5°

(a)翼帆偏轉角度為0°

設定T=2 s,β=30°,在不同橫搖幅值條件下,實驗測試(test)和有限元仿真分析(finite element analysis, FEA)桅桿根部監視點所得應力值結果之間的誤差統計如圖23所示,最大應力偏差范圍為±4.18%,最小應力偏差范圍為±4.57%。設定T=2 s,φ0=15°,同理,在不同偏轉角度條件下,兩種方式所得應力值結果之間的誤差統計如圖24所示,最大應力偏差范圍為±2.60%,最小應力偏差范圍為±2.61%。

(a)最大等效應力變化曲線

(a)最大等效應力變化曲線

在進行桅桿的應力應變測量的過程中,由于實驗環境、條件、人的操作、實驗工藝等方面的因素,實驗結果與真實值之間存在一些誤差[23]。此外,應變儀在測試桅桿應力過程中容易出現零點漂移的現象,造成測量結果的穩定性誤差。因此,該測試結果包含了桅桿應力測量的實驗誤差。測試結果說明,通過有限元仿真計算可以較為準確地評估桅桿根部的應力變化。

6 結論

(1)本文基于對無人帆船翼帆等效構件的力學分析,建立了翼帆桅桿根部應力數值計算模型,并通過有限元仿真和實驗的方式驗證了該數值計算模型的準確性,實現了桅桿根部應力計算的參數化。在此基礎上,通過該模型提取了在帆船橫搖運動過程中影響桅桿根部應力的7個參數。

(2)利用該模型采用實例計算了桅桿根部應力值,分析了各參數的影響規律和機理。得到應力值與橫搖周期成遞減的負二次冪函數關系,與橫搖幅值、翼帆重心高度、偏心距離和翼帆質量成正比例關系,應力值隨偏轉角度的變化曲線與正弦曲線相吻合,并且在一個橫搖周期內應力值曲線呈現正弦波形的結論。

(3)針對桅桿根部進行了初步的結構優化,并采用實驗驗證了有限元分析用以反映真實工況下桅桿應力變化的準確性,證明利用在桅桿根部添加加強筋的方法可以有效減小桅桿根部應力。

(4)為進一步減小桅桿根部應力,在實際工程應用中可以從兩方面開展深入研究：一是在船舶設計方面要降低無人帆船的橫搖幅值,增大橫搖周期;二是在翼帆本體結構設計方面要開展翼帆輕量化研究,減小翼帆質量,降低翼帆整體重心。