基于CFD 水動力參數(shù)的水下機器人軌跡跟蹤控制

羅一漢，吳家鳴，周匯鋒

華南理工大學 土木與交通學院，廣東 廣州 510641

0 引 言

隨著海洋資源的開發(fā)和水下領(lǐng)域作業(yè)任務(wù)增加，水下機器人扮演著越來越重要的角色。但機器人在水下環(huán)境工作時，有著非線性強、易受環(huán)境影響的特點。為獲取良好的作業(yè)效果，對水下機器人的運動進行有效控制，使之能較好地跟蹤預先設(shè)計的軌跡變得十分重要，這也使得軌跡跟蹤控制成為水下機器人領(lǐng)域的主要研究熱點之一。同時在對控制器控制效果進行檢驗時，水下機器人水動力參數(shù)的精確與否直接影響控制器控制效果的可信度。

針對水下機器人的水動力參數(shù)獲取和航行軌跡跟蹤控制，國內(nèi)外有著相當豐富的研究成果。在水動力參數(shù)獲取方面，使用CFD 技術(shù)獲取是主要手段。Ramírez-macías 等[1]使用CFD 技術(shù)對某遙控水下機器人（ROV）艏部面罩進行了水動力模型建模，采用黏性流解算器ReFRESCO 求解單相湍流不可壓縮流的穩(wěn)態(tài)和非穩(wěn)態(tài)N-S 方程，對不同來流方向下的ROV 橫搖及其繞X，Y，Z軸旋轉(zhuǎn)時的力和力矩進行了計算。程嘯鵬等[2]基于Fluent 軟件探究某小型水下機器人水動力與速度、角加速度之間的關(guān)系，并通過實體試驗證明了仿真結(jié)果的可信性。趙凱凱等[3]以升力、阻力及縱傾力矩為對比指標進行數(shù)值仿真，對蛇形機器人外形進行優(yōu)化，提高了此水下機器人的水動力性能。吳家鳴等[4]基于動網(wǎng)格和滑移網(wǎng)格技術(shù)，采用Fluent 軟件對潛器做回轉(zhuǎn)運動時在空泡影響下的螺旋槳推力特性進行了研究。陳東軍等[5]基于CFD技術(shù)對帶纜水下機器人系統(tǒng)升沉運動時的螺旋槳推力特性、水動力響應(yīng)進行了分析研究。

在水下機器人軌跡跟蹤方面，薛乃耀等[6]基于神經(jīng)網(wǎng)絡(luò)、滑?？刂茖壽E跟蹤提出一種變增益抗飽和系統(tǒng)并利用Matlab 進行仿真實驗，結(jié)果表明，軌跡跟蹤時推進系統(tǒng)推力飽和持續(xù)時間降低27%，從而驗證了控制策略的可行性。Zhang等[7]提出一種新型自適應(yīng)非線性二階滑?？刂破饔糜谙墩襁\動，并使用自適應(yīng)律實時估計外界干擾的幅值，實模試驗結(jié)果表明，設(shè)計的控制器相較于線性二階滑模面控制器有著更快的收斂速度和抗抖振能力。Liang 等[8]針對參數(shù)不確定性和外部干擾下的三維軌跡跟蹤控制問題，提出了一種基于反步法和滑模的自適應(yīng)魯棒控制策略，使用模糊邏輯逼近非線性函數(shù)。仿真結(jié)果表明，該控制器有著良好的魯棒性。陳偉等[9]針對欠驅(qū)動自主水下航行器基于視線法（LOS）和反步法設(shè)計了一種水平面軌跡跟蹤控制器，并通過數(shù)值仿真驗證了該控制器的有效性。

上述關(guān)于控制策略研究文獻中的水動力參數(shù)來源有：CFD 計算獲取[6]、實模參數(shù)[7-8]以及為驗證控制器的有效性而適當選取[9]。然而，這些水動力參數(shù)存在精度差異，且設(shè)計控制參數(shù)數(shù)量過多也會導致調(diào)參困難的問題。鑒此，本文將使用CFD 軟件STAR-CCM+和ANSYS AQWA 來獲取所需的水動力參數(shù)，以保障參數(shù)的精度。在設(shè)計控制策略時，首先采取類似于文獻[10-11]中使用的坐標變換方式，獲取水下機器人水平面運動軌跡跟蹤控制偏差系統(tǒng)。然后，為彌補文獻[10]中設(shè)計的控制策略無法跟蹤轉(zhuǎn)艏運動期望角速度為0 時的運動軌跡缺陷，本文將基于反步法和滑?？刂萍夹g(shù)提出一種新的水下機器人水平面軌跡跟蹤控制策略，并使用線性和非線性目標軌跡對控制器的有效性進行驗證。

1 模型及參數(shù)獲取

1.1 坐標系及運動方程

本文選用的坐標系如圖1 所示。在固定坐標系O?XYZ中，規(guī)定OZ的正向指向地心，OX和OY軸在水平面內(nèi)相互垂直，可選擇任意方向為軸的正方向，選取水下機器人的主航向為OX的正向。運動坐標系o?xyz建立在機器人本體上，隨機器人位置的變化而變化。為便于討論，將坐標系原點定于載體重心處。ox軸與主對稱軸一致，以艏部方向為正向；oy軸平行于基線面并垂直于ox軸，以指向右舷為正；oz軸與ox,oy軸垂直，指向底部為正。兩個坐標系滿足右手螺旋法則。

圖1 固定坐標系與運動坐標系示意圖Fig. 1 Sketch of fixed coordinate system and motion coordinate system

由給定的固定坐標系和運動坐標系，定義2 個變量 η和 υ，分別描述兩個坐標系下機器人各自由度上的運動坐標值。兩個坐標系的定義如下：

式中：x,y,z為機器人的位移；φ,θ,ψ為機器人的橫搖角、縱搖角和艏搖角；u,v,w,p,q,r分別為各自由度上的運動線速度與角速度。二者之間的關(guān)系為

式中，R為兩個坐標系的轉(zhuǎn)換矩陣，定義如下：

其中，

由上式，水下機器人六自由度運動方程為

本文主要研究的是機器人水平面運動軌跡跟蹤控制問題，僅考慮機器人在縱蕩 (x)、橫蕩 (y)和轉(zhuǎn)艏(Rz)這3 個自由度上的運動方程，其他自由度上的位移、速度及速度的一階導數(shù)均取為0。選用的仿真模型近似三平面對稱，在運動分析中部分參數(shù)矩陣僅考慮對角線部分，運動坐標系的原點為機器人的質(zhì)心。簡化后的機器人水平面運動方程為

式中，ηs={x,y,ψ}T；υs={u,v,r}T；阻尼力矩陣Ds=diag{D1,D2,D6}，其中Di=Dnoni·|υi|+Dlini，i=1,2,6；其他項也作相應(yīng)簡化?？梢姡粢獙C器人運行軌跡進行控制，需要獲取機器人的附加質(zhì)量、非線性阻尼力項和線性阻尼力項的數(shù)據(jù)。

1.2 幾何建模及其仿真數(shù)學模型水動力系數(shù)的確定

圖2 為本文水下機器人的幾何模型，機器人主體為一個艏艉近似半球形的過渡圓柱形外型結(jié)構(gòu)（總長0.4 m，直徑0.2 m），主體左右舷對稱設(shè)置有一對推進導管螺旋槳，艏艉反對稱設(shè)置有一對用于軌跡與姿態(tài)操縱的導管螺旋槳，分別用于提供進退運動的推力和產(chǎn)生轉(zhuǎn)艏運動的力矩。螺旋槳動力由臍帶纜提供，但假定臍帶纜纜徑較小，受流體干擾影響小，在運動過程中僅提供螺旋槳所需動力，對機器人運動的影響可忽略。研究中所需的機器人質(zhì)心、慣性矩等相關(guān)參數(shù)根據(jù)圖2的水下機器人幾何模型從SolidWorks 建模軟件中獲取，具體數(shù)值如表1 所示。

圖2 水下機器人模型Fig. 2 Underwater vehicle model

表1 水下機器人的相關(guān)參數(shù)Table 1 Relevant parameters of the underwater vehicle

為獲取模型的水動力參數(shù)，本文選用CFD 軟件STAR-CCM+和ANSYS AQWA 分別計算機器人模型阻尼力項參數(shù)和附加質(zhì)量參數(shù)。在阻尼力項計算中，縱（橫）蕩的運動參數(shù)采用DFBI 模型計算，轉(zhuǎn)艏運動阻尼力矩計算模型選用多參考系模型（MRF）計算。圖3 為計算這些運動參數(shù)所使用的網(wǎng)格示意圖。對于本文的縱蕩運動阻尼力參數(shù)計算，計算域大?。▁×y×z）取為4 m×3 m×2 m。為能夠捕獲流動細節(jié)，在水下機器人外部增設(shè)了2 層加密網(wǎng)格，網(wǎng)格類型為TRIM 類型。水下機器人主體網(wǎng)格類型為多面體網(wǎng)格，外部共有4 套網(wǎng)格，整體網(wǎng)格總數(shù)約為120 萬（圖3(a)）。計算域入口邊界設(shè)置流速為0～1.5 m/s 的速度入口，出口邊界設(shè)為壓力出口，其他邊界設(shè)為對稱面，水下機器人壁面設(shè)為無滑移邊界。計算模型選用Reliablek-ε湍流模型 ( RKE )。本文的橫蕩阻尼力計算方法與前述縱蕩運動的一致。

轉(zhuǎn)艏運動阻尼力參數(shù)計算所采用的計算域為高6 m、半徑3 m 的圓柱體，軸線與z軸平行（圖3(b)）。水下機器人外部共有4 層網(wǎng)格，與阻尼力項算例設(shè)計相似，另增設(shè)2 層TRIM 類型網(wǎng)格用于加密。機器人主體設(shè)定為旋轉(zhuǎn)速度為0～2.5 rad/s的旋轉(zhuǎn)域，網(wǎng)格類型為多面體網(wǎng)格，計算域邊界設(shè)置參考阻尼力項算例，計算網(wǎng)格總數(shù)約120 萬。

在計算過程中，將螺旋槳作為無轉(zhuǎn)動浮體來考慮，故在計算機器人附加質(zhì)量參數(shù)時采取了與文獻[12]相同的方法對模型進行近似簡化處理，以保證AQWA 計算效率和質(zhì)量。圖4 為所采用的模型網(wǎng)格示意圖。附加質(zhì)量參數(shù)計算采用ANSYS AQWA 軟件，水密度1 025 kg/m3，計算域（x×y×z）為10 m×10 m×40 m，網(wǎng)格總數(shù)17 362。

圖4 附加質(zhì)量計算網(wǎng)格圖Fig. 4 Grid meshes of additional mass calculation

1.3 網(wǎng)格無關(guān)性驗證計算

為降低網(wǎng)格數(shù)對計算結(jié)果精度的影響，設(shè)計了少、中、多3 種網(wǎng)格數(shù)情況。對于阻尼力（矩）計算，網(wǎng)格總數(shù)分別取為 60 萬、120 萬和300 萬，而附加質(zhì)量計算的網(wǎng)格總數(shù)則分別取為0.8 萬、1.7 萬和2.6 萬。通過對不同網(wǎng)格總數(shù)下的機器人運動阻尼力（矩）和附加質(zhì)量計算結(jié)果的對比（圖5），得到不同網(wǎng)格數(shù)下各水動力參數(shù)計算結(jié)果誤差對比結(jié)果（表2）。表2 中縱（橫）蕩阻尼力和轉(zhuǎn)艏阻尼力矩所得數(shù)據(jù)分別為流速0.9 m/s 和旋轉(zhuǎn)速度1.5 rad/s 時的計算結(jié)果，附加質(zhì)量項選取為附加質(zhì)量項m1的計算結(jié)果，取網(wǎng)格數(shù)（多）情況下的計算結(jié)果作為對比標準。

由圖5(a)和圖5(b)可見，在阻尼力（矩）計算中，3 條曲線非常接近，結(jié)合表2 中阻尼力（矩）的誤差對比，以縱蕩阻尼力和轉(zhuǎn)艏阻尼力矩為例，可見60 萬和120 萬網(wǎng)格數(shù)計算結(jié)果相較于300萬網(wǎng)格數(shù)時的結(jié)果，縱蕩阻尼力項相對誤差由2.314%減少至0.833%，轉(zhuǎn)艏阻尼力矩項由2.271%減少至0.189%，這表明在同一設(shè)計速度點下，3 種網(wǎng)格數(shù)下獲取的阻尼力（矩）值相差不大。因此，雖然增加網(wǎng)格數(shù)可獲得更精確的計算結(jié)果，但收益不大。由圖5(c)和表2 還可見，對于附加質(zhì)量計算結(jié)果，3 個參數(shù)在3 種網(wǎng)格數(shù)下的值基本保持不變，相對誤差在2%左右。因此，綜合考慮計算資源和計算精度，本文阻尼力項和附加質(zhì)量項選擇為120 萬網(wǎng)格數(shù)和1.7 萬網(wǎng)格數(shù)的計算結(jié)果。此外，本文對阻尼力項中的非線性項僅考慮了二次非線性項，對所得阻尼力（矩）參數(shù)進行二次多項式擬合后所得到如圖6 所示擬合曲線。進而可得阻尼力矩陣為：Dnons=diag{12.331·|u|,21.417·|v|,0.612·|r|}，Dlins=diag{1.012,1.541,0.137}，附加質(zhì)量矩陣為：MAs=diag{3.788,15.302,0.382}。擬合曲線表明二階非線性阻尼項遠大于線性阻尼項，只有在機器人運動速度很小的情況下，才能忽略非線性項阻尼對機器人運動的影響。

圖5 不同網(wǎng)格數(shù)下所得水動力參數(shù)對比Fig. 5 Comparison of hydrodynamic parameters under different grid numbers

表2 不同網(wǎng)格數(shù)下水動力參數(shù)計算結(jié)果對比表Table 2 Comparison of calculated hydrodynamic parameters under different grid numbers

圖6 阻尼力（矩）擬合曲線示意圖Fig. 6 Diagram of damping force (moment) fitting curve

2 水平面軌跡控制策略

通過CFD 技術(shù)獲取仿真模型的水動力參數(shù)后，為使用模型對機器人水平面軌跡跟蹤進行控制，還需對作用力矩陣Fs進行定義。本文研究的水下機器人配置有4 個螺旋槳，其中主體左右舷的兩螺旋槳僅提供x方向上的推力，艏艉端的2 個螺旋槳反對稱布置，為轉(zhuǎn)艏運動(Rz)提供一個力偶矩。因忽略臍帶纜對機器人運動的影響，且y自由度上無作用力，故本系統(tǒng)為一欠驅(qū)動系統(tǒng)。此時，作用力矩陣可表示為

式中：T為x方向上的推力；N為轉(zhuǎn)艏力矩。

結(jié)合式 (7)和式 (8)，可得如下運動系統(tǒng)：

式中，符號定義同式 (7)，并令x={x,y,ψ,u,v,r}T。參考文獻[10-11]采用的坐標變換形式，可得：

式中，z6d為轉(zhuǎn)艏運動角速度的期望值，在工程實際中該值一般變化速度不太大，可近似認為z˙6d≈0。因此，第1 個式子對時間求導，將剩余兩式代入，可得：

3 水下機器人軌跡跟蹤控制特性分析

本節(jié)在Matlab/Simulink 平臺對兩類軌跡（前者為線性軌跡，后者為非線性軌跡）下控制策略的有效性進行驗證。具體參數(shù)分別為：

3.1 線性軌跡跟蹤

已知式(27)所示線性軌跡仿真目標參數(shù)及式(29)和式(30)所示水下機器人系統(tǒng)參數(shù)，且控制參數(shù)選擇為：c1=10,η1=5,c2=15,η2=4,k=5,γ=10。但初始條件中z6d=rd=0，若依據(jù)前述控制策略分析會存在一定的跟蹤偏差，故保持其他參數(shù)不變，選取z6d=rd=0.6sin(0.1t)作為對照組進行仿真計算，結(jié)果如圖7～圖8 所示。

由圖7 可見，在設(shè)置的不同z6d條件下，e1，e3，e6響應(yīng)總是能較快收斂到0，時間約為2.2 s；e4和e5在兩種情況下也都能收斂到0，而e2在z6d=0的條件下不為0，但會收斂至一個穩(wěn)定值，而當z6d≠0時，可收斂至0，這與理論分析相符。由圖8(a)也可見，在z6d=0的條件下，仿真軌跡與目標軌跡存在一定的偏差，但兩類軌跡相互平行，表明雖然在x，y自由度上存在一定的位置偏差，但艏搖角ψ可被跟蹤上，即e3=0；而當z6d≠0時，仿真軌跡在一段時間后會與目標軌跡重合，實現(xiàn)了較好的追蹤效果。由圖8(b) 所示的螺旋槳推力示意圖可明顯看出：2 種條件對左右舷提供進退運動推力的螺旋槳影響不大，數(shù)值穩(wěn)定在7 N 左右，但對機器人艏艉端提供偏航力矩的兩個螺旋槳推力的作用較大，可見螺旋槳所需提供的推力明顯提升，而變化較為平穩(wěn)。因此，本文所提控制策略能實現(xiàn)對線性運動軌跡的跟蹤，可解決文獻[10]中軌跡跟蹤控制器無法跟蹤線性軌跡（z6d=0）的問題。

圖7 不同z6d 條件下偏差項響應(yīng)圖Fig. 7 Response of deviation under different z6d conditions

圖8 不同z6d 條件下運動軌跡和螺旋槳推力響應(yīng)圖Fig. 8 Response of trajectory and propeller thrust under different z6d conditions

3.2 非線性軌跡跟蹤

非線性軌跡仿真目標參數(shù)如式(28)所示，水下機器人系統(tǒng)參數(shù)如式(29)和式(30)所示，控制參數(shù)選擇為：c1=2,η1=2,c2=1.5,η2=2.5,k=5,γ=10，仿真結(jié)果如圖9～圖12 所示。

從圖9 可得，在非線性軌跡跟蹤過程中，誤差項ei(i=1 ～6)，都能收斂到0，在10 s 左右便收斂至0 附近，表明控制策略的控制效果較好。但同時可見部分偏差參數(shù)存在較大超調(diào)，結(jié)合圖10 軌跡跟蹤圖可見，由于在初始時刻仿真軌跡和目標軌跡存在較大偏差，使得在偏差收斂過程中有較大的起伏變動，而當偏差收斂至0 后，曲線變得平穩(wěn)。進一步觀察圖10，可見仿真軌跡能夠快速跟蹤上目標軌跡，并保持了良好的持續(xù)跟蹤效果，趨近過程中的最大超調(diào)量約10%。由圖11 可更直觀地看出，在初始時刻實際運動軌跡和目標軌跡的偏差最大，這與所設(shè)計的目標軌跡和機器人初始位置狀態(tài)相符。在運動過程中，伴隨著控制策略的控制動作，約在10 s 左右使得兩個軌跡位置偏差達到0，且后續(xù)幾乎無變化，表明此時已較好地實現(xiàn)了軌跡跟蹤，并能保持持續(xù)的跟蹤效果。由圖12 所示的螺旋槳推力示意圖可見：在0～10 s 內(nèi)，左右舷螺旋槳推力存在明顯的調(diào)節(jié)過程，而在偏差項都收斂至0 后，螺旋槳推力變化較為平緩，艏艉端螺旋槳的推力在約3 s內(nèi)便趨于穩(wěn)定，這與圖9 所示的偏差項變化趨勢是相符的。系統(tǒng)穩(wěn)定后，左右舷螺旋槳推力維持在3.2 N 左右，艏艉端螺旋槳推力維持在15 N左右。因此，對于非線性軌跡，本文設(shè)計的控制策略依舊能夠?qū)崿F(xiàn)較好的追蹤效果并保持螺旋槳推力穩(wěn)定。

圖9 偏差項響應(yīng)圖Fig. 9 Response of deviation term

圖10 軌跡跟蹤示意圖Fig. 10 Diagram of trajectory tracking

圖11 軌跡跟蹤位置偏差示意圖Fig. 11 Diagram of tracking trajectory deviation

圖12 螺旋槳推力示意圖Fig. 12 Diagram of propeller thrust

綜合上述兩類運動軌跡的仿真結(jié)果，表明本文設(shè)計的控制器可使水下機器人快速跟蹤目標軌跡，能夠保持較好的持續(xù)跟蹤效果，驗證了控制器的有效性；兩種軌跡下的螺旋槳推力在系統(tǒng)穩(wěn)定后都趨于穩(wěn)定，但線性軌跡下的螺旋槳推力存在峰值，而非線性軌跡下的推力較為平緩，這可能是因為初始時刻所設(shè)計的線性軌跡變化速度較快，使得螺旋槳必須提供一個較大推力以滿足軌跡快速趨近的要求所致。因此，在設(shè)計實際軌跡時，應(yīng)充分考慮推進器作業(yè)能力的大小。

4 結(jié) 論

本文針對水下機器人水平面軌跡跟蹤問題進行了研究。首先，考慮二階非線性阻尼力推導了水下機器人三自由度水動力數(shù)學模型，使用STAR-CCM+和ANSYS AQWA 軟件分別獲取了所需阻尼力和附加質(zhì)量參數(shù)。然后，結(jié)合模型自身的動力配置條件、計算所得的水動力參數(shù)，推導得到了機器人水平面運動軌跡跟蹤控制偏差系統(tǒng)。其后，將偏差系統(tǒng)分為3 個子系統(tǒng)，基于反步法和滑模控制技術(shù)逐個進行控制策略設(shè)計，提出了一種新的控制策略，使得能夠?qū)崿F(xiàn)對線性運動軌跡（轉(zhuǎn)艏運動期望角速度為0）和非線性運動軌跡的跟蹤。最后，使用設(shè)計的控制器，針對上述兩類運動軌跡，在Matlab/Simulink 平臺中進行了軌跡跟蹤仿真計算。得到主要結(jié)論如下：

1) 在水動力參數(shù)計算方面，本文所得阻尼項擬合曲線表明，二階非線性阻尼項遠大于線性阻尼項，只有在機器人運動速度很小的情況下才能忽略二階非線性阻尼項對運動的影響。

2) 對于線性運動軌跡（轉(zhuǎn)艏運動期望角速度為0）跟蹤問題，本文設(shè)計的控制策略可彌補所對比控制策略無法跟蹤線性運動軌跡的缺陷。對于非線性軌跡跟蹤問題，本文設(shè)計的控制策略依然可實現(xiàn)較好的跟蹤效果，趨近過程中的最大超調(diào)量約為10%，軌跡跟蹤所需時間約為10 s。

3) 在兩類運動軌跡下，螺旋槳推力在系統(tǒng)穩(wěn)定后都能趨于一個穩(wěn)定值，但線性軌跡時的螺旋槳推力曲線存在峰值，而非線性軌跡時螺旋槳推力曲線較為平緩，其原因是本文設(shè)計的軌跡參數(shù)變化較快所致。因此，在設(shè)計實際運動軌跡時，應(yīng)充分考慮推進器作業(yè)能力的大小。